Methods for haplotype phasing

Abstract

Provided herein are methods for gene allele analysis. In some cases, methods comprise identifying a plurality of unphased gene alleles in proximity ligation nucleic acid sequencing data; computing a read graph of said plurality of unphased gene alleles, wherein: a first allele of said plurality of unphased gene alleles is a first node of said read graph, a second allele of said plurality of unphased gene alleles is a second node of said read graph, said first node and said second node are connected in pairwise combination, and a weight of each connection between said first node and said second node is based on an amount of sequence overlap between said first allele and said second allele in said sequencing data; and evaluating said weight of each connection between said first node and said second node to determine which of said plurality of unphased gene alleles have the highest likelihood to be in the same phase.

Description

METHODS FOR HAPLOTYPE PHASINGCROSS REFERENCE

[0001] This application claims the benefit of U.S. Provisional Application No. 63 / 624,222, filed January 23, 2024, and U.S. Provisional Application No. 63 / 558,054, filed February 26, 2024, each of which is incorporated herein by reference in its entirety.BACKGROUND

[0002] Genotyping is the process of identifying the genetic variants of an individual organism, which is crucial for a variety of applications, including understanding genetic predispositions to diseases, personalized medicine, and population genetics studies. Determining an individual’s haplotype, or which genetic variants come from the maternal or paternal chromosome, can also be important to understanding human health and disease.SUMMARY

[0003] In an aspect, provided herein are methods for gene allele analysis. In some cases, the method comprises identifying a plurality of unphased gene alleles in proximity ligation nucleic acid sequencing data. In some cases, the method comprises computing a read graph of the plurality of unphased gene alleles, wherein: a first allele of the plurality of unphased gene alleles is a first node of the read graph, a second allele of the plurality of unphased gene alleles is a second node of the read graph, the first node and the second node are connected in pairwise combination, and a weight of each connection between the first node and the second node is based on an amount of sequence overlap between the first allele and the second allele in the sequencing data. In some cases, the method comprises evaluating the weight of each connection between the first node and the second node to determine which of the plurality of unphased gene alleles have the highest likelihood to be in the same phase. In some cases, the plurality of unphased gene alleles comprises a human leukocyte antigen (HLA) gene. In some cases, the HLA gene is HLA- A, HLA-B, HLA-C, DRB1, DQA1, DQB1, DPA1, DPB1, or a combination thereof. In some cases, the amount of sequence overlap between the first allele and the second allele is about 50 bases, about 100 bases, about 150 bases, about 200 bases, about 250 bases, about 300 bases, or about 350 bases. In some cases, the weight of each connection between the first node and the second node is determined using a Markov hitting time. In some cases, the weight of each connection between the first node and thesecond node is calculated as a symmetric hitting probability between the first node and the second node. In some cases, calculating an optimal read graph comprises solving a system of linear equations. In some cases, the proximity ligation nucleic acid sequencing data is obtained by crosslinking a sample, fragmenting nucleic acids in the sample to produce nucleic acid fragments, ligating the nucleic acid fragments to produce ligated nucleic acid fragments, reversing crosslinks, and sequencing the ligated nucleic acid fragments. In some cases, the sample is crosslinked by contacting the sample to a crosslinking agent selected from formaldehyde, psoralen, disuccinimidyl glutarate (DSG), ethylene glycol bis(succinimidyl succinate) (EGS), ultraviolet light, or a combination thereof. In some cases, the fragmenting comprises contacting the sample to an enzyme. In some cases, the enzyme is a nuclease, a restriction endonuclease, a transposase, or a combination thereof. In some cases, the nuclease is a micrococcal nuclease. In some cases, the transposase is Tn5. In some cases, fragmenting comprises non-enzymatic cleavage. In some cases, the method comprises subsequent to the ligating, adding a label to the nucleic acid fragments. In some cases, the label comprises biotin. In some cases, the label comprises an oligonucleotide. In some cases, the oligonucleotide comprises a barcode. In some cases, the cross-linking links nucleic acids to nucleic acid binding proteins in the sample. In some cases, an optimal read graph is calculated by evaluating about 20%, about 30%, about 40%, about 50%, about 60%, about 70%, about 80%, about 90%, about 99%, or greater than 99% of all possible pairwise connections between the plurality of unphased gene alleles. In some cases, an optimal read graph is calculated by evaluating a symmetric hitting probability between two alleles of the plurality of unphased gene alleles over about 20%, about 30%, about 40%, about 50%, about 60%, about 70%, about 80%, about 90%, about 99%, or greater than 99% of all possible pairwise connections between the plurality of unphased gene alleles. In some cases, the amount of sequence overlap is calculated by comparing at least one k-mer of a selected size from the first node with at least one k-mer of a same selected size from the second node. In some cases, an allele in the plurality of unphased gene alleles is a node in the read graph, and an edge connects two nodes if the at least one k-mer is shared between the two nodes. In some cases, a k-mer is about 50, about 100, about 150, about 200, about 250, about 300, or about 350 bases. In some cases, latent information about a one or more locations of heterozygous single nucleotide polymorphisms (SNPs) is encoded into the read graph. In some cases, an allele node connectivity measure between a first of two alleles and a second of the two alleles is determined using a Markov hitting time. In some cases, the allele node connectivity measure is computed bi-directionally between a first of two alleles and a secondof two alleles and between the second of two alleles and the first of two alleles. In some cases, the allele connectivity measure is computed for about 20%, about 30%, about 40%, about 50%, about 60%, about 70%, about 80%, about 90%, about 99%, or greater than 99% of all possible pairwise connections between the plurality of unphased gene alleles. In some cases, a reduced allele graph that connects the plurality of unphased gene alleles is formed, wherein each one of the plurality of unphased gene alleles is connected to each other one of the plurality of unphased gene alleles. In some cases, a weight of an edge is calculated as a symmetric hitting probability between each pair of two alleles in the reduced allele graph. In some cases, a calculation of an optimal partition of the plurality of unphased gene alleles into one or more phased haplotypes can be calculated from one or more possible partitions of disjoint sets of the plurality of unphased gene alleles. In some cases, the optimal partition is calculated using a sum of edge weights within each of the one or more possible partitions. In some cases, the optimal partition of the plurality of unphased gene alleles into the one or more phased haplotypes can be calculated from the one or more possible partitions in a bruteforce manner by calculating all possible partitions of the plurality of unphased gene alleles.INCORPORATION BY REFERENCE

[0004] All publications, patents, and patent applications mentioned in this specification are herein incorporated by reference to the same extent as if each individual publication, patent, or patent application was specifically and individually indicated to be incorporated by reference.BRIEF DESCRIPTION OF THE DRAWINGS

[0005] The novel features of the invention are set forth with particularity in the appended claims. A better understanding of the features and advantages of the present invention will be obtained by reference to the following detailed description that sets forth illustrative embodiments, in which the principles of the invention are utilized, and the accompanying drawings of which:

[0006] FIG. 1 A shows an example of sequence variants in the gene HLA-DPB 1.

[0007] FIG. IB shows an example of sequence variants in the gene HLA-DPB 1,

[0008] FIG. 2A shows read pairs that come from same physical molecule and that can establish that a pair of SNPs come from the same parental chromosome. With standardIllumina-type paired end sequencing the maximum gap that can be spanned in a single step is approximately Ikbp.

[0009] FIG. 2B shows additional SNPs, which may or may not be relevant to determining the type (e.g., they may occur in intronic or intergenic regions of no consequence to the phenotype) can be used to link paired reads together into chains to span gaps larger than the maximum insert size. In this illustration, for example, there is a G / T polymorphism at Position 2 between the type determining A / A and G / C pairs of SNPs at Positionl and Position 3. A read pair which contains the A and the G SNPs and another that contains the G and subsequent A SNPs can establish that the two A SNPs (Positionl and Positions) are on the same chromosome and correspondingly that the G and C SNPs (Positionl and Positions) are on the same chromosome.

[0010] FIG. 3 displays an example alignment for the gene HLA-DRB1, with HG002 aligned to it.

[0011] FIG. 4 displays the locations of six HLA genes on chromosome 6p 21.31.

[0012] FIG. 5 displays a general schematic of the methods described herein used for phasing HLA alleles.

[0013] FIG. 6 shows a schematic of how HLA genes are linked using overlapping k-mers.

[0014] FIG. 7 shows a schematic of allele nodes being placed in a read graph.

[0015] FIG. 8 shows a schematic of how pairwise allele node connectivity can be calculated for each possible pair of HLA genes.

[0016] FIG. 9 shows an example reduced allele graph.

[0017] FIG. 10 shows an example optimal partition for a reduced allele graph.

[0018] FIG. 11 shows a graphical representation of an optimal partition for the Int-237hlaA region.

[0019] FIG. 12 shows a graphical representation of an optimal partition for the Int-237hlaB region.

[0020] FIG. 13 shows a graphical representation of an optimal partition for the Int-223hlaD- V4 region.

[0021] FIG. 14 shows a computer system that is programmed or otherwise configured to implement methods provided herein.

[0022] FIG. 15 shows a Toy Example where k=3 with 6 basepair reads.

[0023] FIG. 16 shows a Toy Example where k=3 with 6 basepair reads.

[0024] FIG. 17 shows a Toy Example where k=3 with 6 basepair reads.

[0025] FIG. 18 shows a Toy Example where k=3 with 6 basepair reads with added long range linkage information (e.g., proximity ligation).

[0026] FIG. 19 shows a Toy Example where k=3 with 6 basepair reads with use of long range links from proximity ligation libraries.

[0027] FIG. 20 shows a Toy Example where k=3 with 6 basepair reads with use of long range links from proximity ligation libraries.

[0028] FIG. 21 shows a diagram of information stored by a node.

[0029] FIG. 22 shows another example graph.DETAILED DESCRIPTION

[0030] Traditional genotyping techniques, such as Polymerase Chain Reaction (PCR)-based approaches, have been widely used but are limited by their dependence on specific primer design. Primers are designed to target known genetic variations, which means they may miss novel or rare variations that lie outside the specific range covered by the primers. This restricts the ability of PCR-based methods to capture the full spectrum of genetic diversity and may result in incomplete or inaccurate genotype determination, particularly in cases where there are uncharacterized or unexpected genetic variations. The reliance on predefined primer sets makes PCR-based genotyping less suitable for identifying rare or novel variants, highlighting the need for alternative approaches that can provide a broader coverage of the genetic landscape.

[0031] The introduction of Next-Generation Sequencing (NGS) methods into genotyping offers the potential for more comprehensive genotyping by providing sequence data that covers every base of the entire locus or gene. NGS sequencing allows the comprehensive determination of all single nucleotide polymorphisms (SNPs) at a given genetic loci, allowing the precise characterization of both common and rare or novel alleles. However, even with these technologies, challenges persist in accurately resolving genotype ambiguities. One source of NGS genotyping ambiguities is the difficulty of assigning SNPs to a particular chromosome when they are at locations separated by SNP-free areas exceeding the read pair insert size.

[0032] The accurate determination of genotypes is crucial for various genetic studies and applications, including disease risk assessment, drug response prediction, and population genetics research. However, resolving genotype ambiguities presents a significant challenge in these endeavors.

[0033] One common issue arises when attempting to assign gene fragments, such as exons, to their respective parental chromosomes. The most significant variations are those which affect the protein constructed by the gene, although there may be other effects due to non-exonic mutations. While NGS sequencing allows the accurate determination of SNPs at every exon location, in order to determine the specific protein constructed and the type to assign to the allele, it’s necessary to establish which parental chromosome the given SNPs are located on.

[0034] This becomes particularly challenging when the fragments are separated by large SNP-free gaps (SNP deserts), which cannot be bridged using conventional sequencing methods like Illumina paired-end sequencing. These regions can occur naturally due to the characteristics of a locus and its neighboring regions in a population or may be a result of technical limitations, such as the absence of capture probes in a given region. Resolving genotype ambiguities in these regions becomes challenging due to the limited availability of informative markers. The typical mean insert size of Illumina-type shotgun sequencing is around 300 base pairs (bp), with a maximum insert size of approximately 1,000 bp. As a result, SNPs from loci separated by SNP gaps larger than these pair insert sizes cannot reliably be linked to one another to determine if they come from the same or different parental chromosome (e.g., physical DNA molecule).

[0035] Methods provided herein utilize a combination of proximity ligation assay and paired- end sequencing to overcome this limitation of normal NGS sequencing. Hi-C type proximity ligation data produces read pairs from a single parental chromosome but separated by a distance that can range from 100s of bp to chromosome length. These long-range read-pairs can be used to determine the chromosomal phasing of alleles over large ranges unavailable to typical NGS sequencing, allowing the assignment of exons or other genomic segments to the correct parental chromosome resolving genotype ambiguities that can otherwise occur.

[0036] Further, methods provided herein allow for the accurate assignment of genomic segments, such as exons, to parental chromosomes and enables the resolution of genotype ambiguities in regions with large gaps or low genetic variation, thereby facilitating a deeper understanding of an individual's genome and its implications in various genetic studies and applications.

[0037] Genotype ambiguities, where the specific genetic variants present in an individual's genome cannot be accurately determined, pose a significant challenge in various genetic applications. Resolving these ambiguities is essential for a comprehensive understanding of an individual's genetic makeup and its implications in clinical and research settings.

[0038] An illustrative example is the typing of the human leukocyte antigen (HLA)-DPBl locus, which has clinical relevance in hematopoietic stem cell transplantation and solid organ transplantation. The typing of HLA alleles involves mapping variations in each HL A gene to previously observed variations in the population. These previously identified variations are specified using a nomenclature such as DPB 1*02:01 :02, where the parts of this notation are:[gene]* [allele group]: [HLA protein]: [synonymous coding mutation]: [non-coding mutation]

[0039] The full nomenclature can have up to four comma separated fields specifying variations in the sequence at a loci. The process of typing is an attempt to unambiguously assign a unique type code for each of the two alleles present in a sample (patient). For example, DPB 1*04:02:01 :01, DPB 1*02:01 :02:01 can indicate that variations that match DPB 1*04:02:01 :01 are present on one parental chromosome and that variations that match DPB 1*02:01 :02:01 are present on the other parental chromosome.

[0040] FIG. 1 A and FIG. IB illustrate examples where an ambiguity can arise even when the precise variants present at each given genomic coordinate are known. In FIG. 1 A the sequence variants observed in codon 98 in exon 2 of the gene HLA-DPB1 are (A,G). The sequence variants observed at codon 207 in exon 3 are (C,A). Knowing only the SNPs at these loci and not their phase (e.g., which SNPs appear on the same parental chromosome) it is not possible to distinguish between Possible Type Pair 1 in the figure (DPBl*04:02:01 :01 / DPBl*02:01 :02:01) and Possible Type Pair 2 in the figure(DPB 1 * 105 :01 :01 :01 / DPB 1*416:01 :01 :01). Which pair is the correct typing for DPB1 depends on whether one parent is A / A for exon2 / exon3 respectively and the other parent is G / C for exon2 / exon3 respectively, or whether one parent is A / C for exon2 / exon3 respectively and the other parent is G / C for exon2 / exon3. The observed unphased SNPs will be the same for these two cases.

[0041] FIG. IB illustrates another example, here the sequence variants observed at exon2 and exon3 are (TAA,CCC) and (C,A) respectively. Again, whether the correct typing is Possible Type Pair 1 (DPB 1 * 04 : 01 : 01 : 01 / DPB 1 * 105 : 01 : 01 : 01 ) or Possible Type Pair 2 (DPBl*04:02:01 :01 / DPBl*126:01) depends on whether the CCC variant in exon2 is on the same chromosome as the C variant in exon 3 or whether the CCC variant in exon2 is on the same physical chromosome as the A variant in exon3 or on the same chromosome as the C variant in exon3, and similarly on whether the TAA variant in exon2 is on the same parental chromosome as the A variant in exon3 or the C variant in exon3.

[0042] Studies have indicated that up to 26% of samples exhibit ambiguous HLA-DPB1 alleles and that the primary source of these uncertainties lies in phasing exonic heterozygouspositions between exon 2 and subsequent downstream exons such as those illustrated in FIG. 1A and FIG. IB.

[0043] Using Illumina-style paired-end sequencing it is possible to phase variants over short distances. Paired end sequencing produces two reads, one from the forward DNA (Crick) strand and one from the reverse (Watson) strand separated by a short distance (the insert size). A given sequencing run will produce a range of insert sizes with an average insert size of about 300bp with a max insert size up to 1000 bp.

[0044] Read pairs come from same physical molecule and so can establish that a pair of SNPs come from the same parental chromosome. With standard Illumina-type paired end sequencing the maximum gap that can be spanned in a single step is approximately Ikbp. This is illustrated in FIG. 2A, where the type determining SNPs A / A can be established to be on the same chromosome and the G / C SNPs can be established to be on the same chromosome implying that the proper type disambiguation is Possible Type Pair 1, that is, Typel + Type2.

[0045] In some cases, the limitation of Illumina-style NGS paired-end sequencing to establishing the phase of SNPs no more than about 1000 bp apart can be overcome by chaining together multiple such pairs across several SNPs. This is illustrated in FIG. 2B. Additional SNPs, which may or may not be relevant to determining the type (e.g., they may occur in intronic or intergenic regions of no consequence to the phenotype) can be used to link paired reads together into chains to span gaps larger than the maximum insert size. In FIG. 2B, for example, there is a G / T polymorphism at Position 2 between the type determining A / A and G / C pairs of SNPs at Positionl and Position 3. A read pair which contains the A and the G SNPs and another read pair that contains the G and subsequent A SNPs can establish that the two A SNPs (Positionl and Position3) are on the same chromosome and correspondingly that the G and C SNPs (Positionl and Position3) are on the same chromosome.

[0046] The problem that still arises with Illumina-style paired-end sequencing is that there may not always be enough intermediate SNPs to allow building chains to span large gaps. The ambiguous genotypes described in FIG. 1 A and FIG. IB are examples of this phenomena. Typing HLA-DPB1 unambiguously is challenging due to the presence of a nearly SNP-free intron approximately 4,000 base pairs (bp) in length between exon 2 and exon 3. This intron acts as a substantial barrier, preventing the determination of the parental phase for exons 3-5 relative to exons 1-2 since it is both larger than Illumina-style paired end insert sizes and does not provide any SNPs to allow chaining pairs together. As a result, itbecomes impossible to accurately assign these exons to their respective parental chromosomes using traditional sequencing approaches like Illumina sequencing. The majority of the ambiguous alleles described in Duke et al. (Duke JL, Mosbruger TL, Ferriola D, Chitnis N, Hu T, Tairis N, Margolis DJ, Monos DS. Resolving MiSeq-Generated Ambiguities in HLA-DPB1 Typing by Using the Oxford Nanopore Technology. J Mol Diagn. 2019 Sep;21(5):852-861. doi: 10.1016 / j.jmoldx.2019.04.009. Epub 2019 Jun 4. PMID: 31173929; PMCID: PMC6734860) are of this type.

[0047] Considering these factors, there is a need for alternative approaches that can bridge the gap between the limitations of short-read sequencing, such as Illumina, and the challenges associated with long-read single molecule technologies in clinical settings. The methods provided herein which utilize a proximity ligation assay combined with normal paired-end sequencing address these challenges by providing a cost-effective and practical solution for resolving genotype ambiguities in scenarios where traditional methods and long-read sequencing have limitations.

[0048] One method developed for analyzing HLA alleles may comprise aligning a plurality of sequencing reads to a reference genome, wherein at least a portion of the plurality of sequencing reads correspond to an HLA gene locus. Next, the method can comprise identifying a plurality of single nucleotide polymorphisms (SNPs) and / or indels in the plurality of sequencing reads. Alternatively, or in combination, the method can comprise identifying the plurality of SNPs by comparing the plurality of sequencing reads to a variation graph. Then, the method can comprise sorting the plurality of SNPs and indels into phase blocks. The method can then comprise aligning the plurality of sequencing reads to a variation graph to identify a plurality of HLA types. Next the method can comprise comparing the plurality of HLA types to a plurality of known HLA alleles to obtain a SNP signature for each HLA allele. Then, the method can comprise comparing the SNP signature to the phase blocks to obtain the phased HLA type.

[0049] These methods for determining phasing by building together chains of SNPs or by sorting SNPs into phase blocks can, however, still be insufficient for some applications. For example, when determining phasing for HLA genes (location of HLA genes on chromosome 6p 21.31 shown in FIG. 4), a SNP -based approach can be insufficient because the exceptionally high variability in the HLA region can cause short-read alignment to be unreliable, particularly for heterozygous SNPs. Any single reference genome will be a poor match for most individuals in the HLA region because of this variability. It is thus difficult to choose an appropriate reference genome. FIG. 3 displays an example alignment for drbl,with HG002 aligned to it. This alignment displays many misalignments. Additionally, many potential variants are not consistently called across reads. There are also coverage anomalies, such as peaks and valleys and inconsistent coverage. These issues with alignment cause significant problems for the rest of the phasing method, because the other steps in the phasing pipeline are downstream of alignment. Solving problems with alignment, then, can provide more reliable outcomes in the phasing of alleles.

[0050] Provided herein are methods for gene allele analysis. In some cases, the method comprises identifying a plurality of unphased gene alleles in proximity ligation nucleic acid sequencing data. In some cases, the method comprises computing a read graph of the plurality of unphased gene alleles. In some cases, a first allele of the plurality of unphased gene alleles is a first node of the read graph, a second allele of the plurality of unphased gene alleles is a second node of the read graph, the first node and the second node are connected in pairwise combination, and a weight of each connection between the first node and the second node is based on an amount of sequence overlap between the first allele and the second allele in the sequencing data. In some cases, the method comprises evaluating the weight of each connection between the first node and the second node to determine which of the plurality of unphased gene alleles have the highest likelihood to be in the same phase.

[0051] In various aspects of methods provided herein the plurality of sequencing reads are obtained by cross-linking the sample, fragmenting nucleic acids in the sample to produce nucleic acid fragments, ligating the nucleic acid fragments to produce ligated nucleic acid fragments, reversing crosslinks, and sequencing the ligated nucleic acid fragments. In some cases, the sample is crosslinked by contacting the sample to a crosslinking agent selected from formaldehyde, psoralen, disuccinimidyl glutarate (DSG), ethylene glycol bis(succinimidyl succinate) (EGS), ultraviolet light, or a combination thereof. In some cases, the fragmenting comprising contacting the sample to an enzyme. In some cases, the enzyme is a nuclease, a restriction endonuclease, a transposase, or a combination thereof. In some cases, the nuclease is a micrococcal nuclease. In some cases, the transposase is Tn5. In some cases, fragmenting comprises non-enzymatic cleavage. In some cases, the method further comprises subsequent to the ligating, adding a label to the nucleic acid fragments. In some cases, the label comprises biotin. In some cases, the label comprises an oligonucleotide. In some cases, the oligonucleotide comprises a barcode. In some cases, the cross-linking links nucleic acids to nucleic acid binding proteins in the sample.

[0052] In various aspects of methods provided herein, the plurality of unphased gene alleles comprises an HLA gene. In some cases, the HLA gene is HLA-A, HLA-B, HLA-C, DRB1, DQA1, DQB1, DPA1, DPB1, or a combination thereof.

[0053] In various aspects of methods provided herein, the amount of sequence overlap between the first allele and the second allele is about 50 bases, about 100 bases, about 150 bases, about 200 bases, about 250 bases, about 300 bases, or about 350 bases. In some cases, the weight of each connection between the first node and the second node is determined using a Markov hitting time. In some cases, the weight of each connection between the first node and the second node is calculated as a symmetric hitting probability between the first node and the second node. In some cases, calculating an optimal read graph comprises solving a system of linear equations. In some cases, an optimal read graph is calculated by evaluating about 20%, about 30%, about 40%, about 50%, about 60%, about 70%, about 80%, about 90%, about 99%, or greater than 99% of all possible pairwise connections between the plurality of unphased gene alleles. In some cases, an optimal read graph is calculated by evaluating a symmetric hitting probability between two alleles of the plurality of unphased gene alleles over about 20%, about 30%, about 40%, about 50%, about 60%, about 70%, about 80%, about 90%, about 99%, or greater than 99% of all possible pairwise connections between the plurality of unphased gene alleles.

[0054] In various aspects of methods provided herein, the amount of sequence overlap is calculated by comparing at least one k-mer of a selected size from the first node with at least one k-mer of a same selected size from the second node. In some cases, an allele in the plurality of unphased gene alleles is a node in the read graph, and an edge connects two nodes if at least one k-mer is shared between the two nodes. In some cases, a k-mer is about 50, about 100, about 150, about 200, about 250, about 300, or about 350 bases.

[0055] In various aspects of methods provided herein, latent information about one or more locations of heterozygous SNPs is encoded into the read graph.

[0056] In various aspects of methods provided herein, an allele node connectivity measure between a first of two alleles and a second of the two alleles is determined using a Markov hitting time. In some cases, the allele node connectivity measure is computed bi-directionally between a first of two alleles and a second of two alleles and between the second of two alleles and the first of two alleles. In some cases, the allele connectivity measure is computed for about 20%, about 30%, about 40%, about 50%, about 60%, about 70%, about 80%, about 90%, about 99%, or greater than 99% of all possible pairwise connections between the plurality of unphased gene alleles.

[0057] In various aspects of methods provided herein, a reduced allele graph that connects the plurality of unphased gene alleles is formed, wherein each one of the plurality of unphased gene alleles is connected to each other one of the plurality of unphased gene alleles. In some cases, a weight of an edge is calculated as a symmetric hitting probability between each pair of two alleles in the reduced allele graph. In some cases, a calculation of an optimal partition of the plurality of unphased gene alleles into one or more phased haplotypes can be calculated from one or more possible partitions of disjoint sets of the plurality of unphased gene alleles. In some cases, the optimal partition is calculated using a sum of edge weights within each of the one or more possible partitions. In some cases, the optimal partition of the plurality of unphased gene alleles into the one or more phased haplotypes can be calculated from the one or more possible partitions in a brute-force manner by calculating all possible partitions of the plurality of unphased gene alleles.Allele Read Graphs

[0058] In various aspects of methods provided herein a read graph as described herein is a graph as described in the field of graph theory, not a graph as in a plot or data visualization. In the field of graph theory, a graph is a mathematical structure amounting to a collection of objects in which the objects are somehow related. Each object is called a vertex, point, or node. Two nodes may be connected with an “edge” to represent a relationship between them. Visually, this may be represented as a set of dots or circles (the nodes) connected by lines (the edges).

[0059] Edges may be directed or undirected. Undirected edges means that the relationship between the two nodes is bi-directional. Directed edges mean that the relationship is unidirectional. For example, if the nodes represent people, and the edges represent whether two people know each other, then those edges would be undirected, because either both people know each other or they do not. Alternatively, if the nodes represent people, and the edges represent people who owe money to each other, then the edges may be directed or undirected, since for example person A may owe money to person B but not the other way around (directed), or both may owe money to each other (undirected).

[0060] In the methods described herein, the nodes are unphased gene alleles. The term “phasing” refers to the process of assigning alleles to the paternal and maternal chromosomes. For example, unphased gene alleles for HLA genes may be identified using a graph-guided assembler such as Kourami (Lee and Kingsford, 2018). Kourami uses an assembly technique wherein full sequences are constructed for the peptide binding domain of an HLA gene using a modified partial-order graph (POG) as a guide. The POG representationcaptures variant regions in related sequences, taking advantage of known alleles. It also can provide a framework that allows incorporation of information from sequencing reads to assemble novel alleles. Kourami can assemble novel alleles not found in a database. However, Kourami does not perform allele phasing, thus creating a need for methods to accurately assign phased haplotypes to HLA alleles.

[0061] In some cases, HLA alleles identified by Kourami or other graph-guided assemblers, along with FASTQ proximity ligation sequencing data, can be used in the method described herein. In some cases, a read graph may be built, wherein each HLA allele is represented by a node. Edges are created to connect these nodes based on overlap between reads corresponding to these alleles. A schematic for this analysis pipeline is provided in FIG. 5.

[0062] An edge connects two nodes if their corresponding reads share a A-mer . The term “k- mer” is frequently used in bioinformatics to refer to a substring of length k contained within a nucleic acid sequence. For example, a A-mer of length 5 would be a substring of 5 nucleotide bases. Aimers may be used in sequence assembly for the construction of De Bruijn graphs. In some embodiments of the present disclosure, Aimers of overlapping sequences can be used to connect the unphased alleles. A schematic showing k-mer overlap between alleles is shown in FIG. 6. An edge can connect two nodes if their corresponding reads share a A-mer of a selected size. In some cases, the selected A-mer size may be about 50, about 100, about 150, about 200, about 250, about 300, or about 350 bases.

[0063] Using Aimers, allele nodes can be placed into a read graph, with the allele nodes corresponding to the same haplotype being better connected than those on opposite haplotypes (FIG. 7). Additionally, despite no SNP typing being explicitly performed, latent information about the location of heterozygous SNPs may be encoded into the read graph. A metric called pairwise allele node connectivity may then be determined using a Markov hitting time. This metric may be calculated in both directions. For example, the pairwise allele node connectivity can be calculated for each possible pair of HLA genes (example schematic shown in FIG. 8). Once all the pairwise allele node connectivity values are calculated, a reduced allele graph may be formed which connects the HLA alleles to each other (example schematic shown in FIG. 9). The weight of each edge may be calculated as the symmetric hitting probability between the two alleles. The possible partitions of HLA alleles into haplotypes may be calculated. The optimal partition can be calculated by “brute force”, e.g., testing all partitions. The optimal partition may be found by calculating the sum of the edge weights within each partition. Calculating the sum of the edge weights may indicate the strength of connectivity between the alleles in that partition and thus mayindicate which alleles are most likely to belong to the same parent haplotype. An example schematic for an optimal partition is shown in FIG. 10.Phased de Bruijn Graphs

[0064] In various aspects of method provided herein, in some cases de Bruijn graphs are used to phase a haplotype, such as an HLA haplotype.

[0065] Traditional sequence alignment to a reference genome (e.g., hg38) can be prone to problems in regions with high complexity & diversity in the population. For example, the HLA locus features multiple genes (HLA-A, B, C ... DRB1, etc.). Each gene has a diverse variety of alleles (HLA*A:01:01, *A:02:01, ...) that are defined by the presence of a unique set of intra- and inter-genic SNVs & indels. Furthermore, the genes are located in a region of chromosome 6 with complex structural variation. There are 8 major structural haplotypes for the HLA locus, with many more minor structural variants branching off these major haplotypes. For regions with such high diversity, it can be a challenge to select an appropriate reference sequence. Every SNP / indel / structural variant in the reference genome that is not present in the sequencing data will induce a variant call at that site (and vice versa). If the sequencing data and the reference are too dissimilar (e.g., a dense cluster of small variants, the wrong structural haplotype), reads can be misaligned (poor quality mappings) or not aligned at all (read dropout). A reference-free approach sidesteps these issues by reframing the alignment problem as the presence of a set of kmers describing the sequence(s) of interest (e.g., HLA gene alleles) in the sequencing dataset where a kmer is a short string of fixed length (e.g., 23bp) comprised of DNA nucleotides (A, C, T, G). Furthermore, kmers can be used to construct a de Bruijn graph capable of storing a huge variety of genes / alleles / structural haplotypes in a single condensed graph structure. Typically built from a large reference database (e.g., FASTA sequences from IMGT / HLA database) to encode known / ob served sequence variability. One can also construct de Bruijn graphs from the raw sequencing data itself for reference-free assembly and variant detection. No information about where these sequences are located in a reference genome is required, but can be added for visualization purposes or identifying kmers belonging to a genomic neighborhood.

[0066] In an aspect, methods of haplotype phasing using phased de Bruijn graphs can include collecting all reads mapping to a locus of interest (and their mates). In some cases, the method comprises using single-ended reads, build an uncolored, de Bruijn graph (e.g., forward & reverse, k = 23). In some cases, the method comprises pruning the resulting graph. For example by removing kmers with support less than the min-threshold and morethan the max-threshold; collecting statistics and removing top / bottom 5%; using user-defined thresholds (e.g. 200x library -> require 50-100x support); removing disconnected ends (chimeric reads), which is a more targeted approach than across-the-board removal of poorly support kmers; removing all kmer branches that do not reconnect back to graph; pruning poorly supported bubbles, which is another approach that is more targeted; identifying bubbles by finding kmers along bubble path with weakest read support if path’s minimal kmer support is lower than threshold, then remove bubble from graph; require minimum support threshold for phased links; or a combination of above approaches. In some cases, the method comprises identifying all bubbles in graph, where for all kmers, search for branch points (start of the bubble) where a given kmer is connected to multiple adjacent kmers; from each branch point, find first shared kmer downstream (where branches reconnect, breadth- first search); backtrack to enumerate all unique paths back to source (tag kmers during BFS to know which paths to follow in reverse); and assign each bubble path an index (bpi), and tag corresponding kmers with their bpi (or store kmers & their bubble-indices in db). In some cases, the method comprises connecting disparate bubble paths using proximity-linked read pairs. In some cases, each bubble in graph corresponds to a variant (SNP, indel, etc.) present within the data and read pairs are used to make direct connections between (potentially distant) kmers in the de Bruijn graph where increment link’s support for each read pair supporting it. During this stage, it can be optional to only link kmers if they belong to an enumerated bubble in graph (i.e., only link variants, not regions of homology) if it reduces computational / storage needs. In some cases, bubbles are linked by choosing bubble with greatest read support; assigning each unique path (i.e. set of kmers) through bubble a different “color”, where in the case of a bubble-within-bubble (multi-allelic site), only color the kmers unique to each path; for each colored kmer, find linked kmers using bubble-to- bubble (b2b) connections and assign them the same color (repeat this step for all newly colored kmers, finding their linked kmers using b2b connections & coloring them the same color); repeat the previous step all linked kmers are visited & colored; removing all colored bubble paths from the set of bubble paths to link, and repeat steps 1-5 (using new colors) until all bubbles in graph are colored. In some cases, the method comprises resolving linking ambiguities (bubble path gets 2+ colors), then: removing all colors for path (discard links); removing colors with weak support (threshold); and retaining color with the most support, discard all others.

[0067] Simple bi-allelic bubbles in a de Bruijn can be searched for more quickly by implementing constraints on which paths to follow when searching for bubbles. Theminimum read support for any node along a given bubble’s path can be set by a user-defined fraction of the read support of the source node for said bubble. For example, if a source node’s read support is 100, and the user-defined fraction is 0.25, then only nodes with read support of 25 or more are followed to define each path of a bubble starting at that source node. This drastically helps identify simple bubbles for source nodes with high read support (e.g., greater than 1000 reads). Without such a constraint, source nodes with high read support can have numerous potential paths induced by low frequency sequencing errors (e.g., 2% sequencing error for 1000 reads could cause 20 random errors at any given site, each potentially generating a new kmer node in the de Bruijn graph). These additional paths can complicate bubble finding, causing a potentially exponential increase in the number of paths to follow. Often such high read support nodes can fail to produce any simple bubbles, instead being classified as a “super bubble” with all of the additional nested bubbles. By constraining to just nodes with relatively similar support to the source node, many fewer paths are followed, and simple bubbles are quickly identified even when node read support is high.

[0068] In various aspects of method provided herein, in some cases the phased de Bruijn graph is useful for various applications, including by not limited to phasing a set of alleles or sequences-of-interest using graph by breaking apart sequences into same-sized kmers, finding kmers linked to bubbles in graph, identifying which color(s) each sequence belongs to, where sequences with same color get phased together. The phased de Bruijn graph is also useful for determining a phased kmer copy number, where given a sequence of interest, this can provide min / max coverage across its set of kmers. If restricted to considering only variant kmers (i.e., kmer inside of bubbles), this will provide an estimate of allele-specific copy number. Over a set of such sequences-of-interest (including a “control” sequence of known / stable copy number), can compute a normalized copy number ratio vs. that control sequence. In some cases, if no control sequence is provided, then normalized coverage is computed against the average read support across all sequences-of-interest. Finally, one can kmer-ize a reference genome and store the coordinates of each kmer and use these positioned-kmers to anchor / linearize haplotypes stored in graph. This could also be used to help visualize results in a more convenient / interpretable manner.

[0069] FIG. 15 shows a toy example (k=3, 6 bp reads) by building a de Bruijn graph from input sequencing reads treated as single-ended, unpaired reads. In practice, shotgun-seq reads will be ~150bp, and kmer size will be sized to be as small as possible while achieve as much uniqueness as needed across the sequence space of the problem at hand. Generally, k is chosen to be an odd number (to avoid reverse-complement palindromes), somewhere in therange of 20-30bp. An example kmer size is 23bp. As kmer size is reduced, total number unique kmers decreases (less memory) and graph becomes more & more lossy. As kmer size is increased, total number unique kmers increases (more memory) and more information is retained in graph. Multiple kmer sizes can be tested to determine the ’’sweet spot” that balances memory demands with graph information. Lastly, lossy graphs are likely to still retain the information needed to solve the phasing problem.

[0070] FIG. 16 shows a toy example (k=3, 6 bp reads) where we learn that SNPs / MNPs / Indels create bubbles in the graph. Shared paths indicate regions of homology between input sequences.

[0071] FIG. 17 shows a toy example (k=3, 6 bp reads) where we see four unique paths through the graph or four different sequences It is noted that each variant creates a bubble in the graph and each bubble doubles the number of potential paths or sequences.

[0072] FIG. 18 shows a toy example (k=3, 6 bp reads) where long range linkage information is added from proximity ligation data. Here, the linked reads between bubble paths are informative for phasing purposes. All kmers from R1 can be linked to all kmers from R2 but to reduce memory requirements can only store links between bubble kmers. The kmer to kmer link supports tallied and stored in graphs or an adjacent database.

[0073] FIG. 19 shows a toy example (k=3, 6 bp reads) where long range links are used from proximity ligation libraries. Here paths are eliminated that are incompatible with long read phasing.

[0074] FIG. 20 shows a toy example (k=3, 6 bp reads) where long range links are used from proximity ligation libraries. Here multiple sequences of interest are broken up into kmers and linked kmers are used to determine their phasing.

[0075] FIG. 21 shows a summary of information stored in a node. Specifically, Nodes are stored in a dictionary indexed by two keys: key 1 is kmer (AGC in example) and key is = reverse-comp kmer (GCT in example). Each node stores: read support (int); bubble path index (null if not in bubble); assigned phase color(s) (int); optionally a set of phase-linked kmers; optionally pre-nucleotides (G in example); and optionally post-nucleotides (A, T in example). The database is separated to store phase information: phase-linked kmers and read support (e.g., kl, k2, n_pairs). This separate database used to reduce redundancy & memory requirements.

[0076] FIG. 22 shows a more realistic initial graph that includes branches to nowhere (e.g., chimeric reads), bubbles within bubbles (multi -allelic), weakly supported bubbles (e.g., sequencing errors), multiple ends and ambiguous phasing. Many techniques for pruninggraph prior to analysis include removing kmers with support less than the minimum threshold and / or support greater than the maximum threshold; collecting stats, remove kmers with “extreme” support (e.g. top / bottom 5%), user-defined thresholds; removing disconnected ends (chimeric reads), this is a more targeted approach than across-the-board removal of poorly support kmers, including remove all kmer branches that do not reconnect back to graph and prune poorly-supported bubbles. Another approach that is more targeted includes identifying bubbles, finding kmers along bubble path with weakest read support. Remove bubble path from graph if path’s minimal kmer support is lower than threshold. Then require phased links to have support above a minimum threshold. Alternatively, one could remove all ambiguous links where a kmer from one bubble path is linked to kmers belonging to multiple paths of a different bubble.

[0077] In various aspects of method provided herein, misphasing bubbles can be identified and removed from the graph. A misphasing bubble is a bubble where both paths of the bubble are linked (i.e., via linked reads) to the same single path of another bubble. This can arise if one or both paths of a bubble exist in more than one copy within the region being analyzed. Misphased bubbles can give rise to ambiguous phasing and unstable haplotype solutions. To avoid these problems, prior to solving the phasing matrix, either the misphasing bubble can be removed entirely, including each path and all its constitutive kmers, or the set of kmers that are contributing to the misphasing bubble can be identified and only those kmers removed.

[0078] In various aspects of method provided herein, normalization of kmer counts in the phasing matrix can improve the ability to phase entities with different numbers of associated kmers. The phasing matrix is the input matrix to phasing algorithm (e.g., MAXCUT) and the matrix contains the linked read support between any two entities being phased (e.g., alleles, bubbles). For example, Matrix(A, B) = (# Linked Reads connecting A Kmers & B Kmers) / (# A Kmers * # B Kmers), where A and B are entities being phased into haplotypes and can be either HLA alleles, bubbles, or combinations thereof. This change to the phasing matrix can enable greatly improved HLA haplotype phasing performance, verified across a full cohort of samples with known truth (e.g., via family trios).Phasing Singleton Genes and Ambiguous Types

[0079] In some cases, an individual could have one or more singleton gene alleles which lack opposing partners to phase against. For example, in the HLA region with the HLA-DRB3, - DRB4, and -DRB5 alleles, an individual could have two alleles from one of these genes (HLA-DRB3*01 + HLA-DRB3*02), or one allele each from two of these genes (e.g. HLA-DRB3*01+ HLA-DRB5*02). In the latter case, each allele didn’t have an opposing partner to phase against, which can pose a challenge for some phasing methods. One approach to address this challenge is to use a “meta-gene” (e.g., in this example labelled DRB345: HLA- DRB345:DRB3*01 + HLA-DRB345:DRB5*02). This enables the phasing method to phase singleton alleles to HLA haplotypes without any other modifications, because any combination of DRB3 / 4 / 5 alleles would have partners to phase against since they were all belonged to the same “meta-gene.” Alternatively, rather than a “meta-gene”, a “pseudoallele” can be used to act as a phasing partner for singleton alleles. In this approach, every singleton allele, defined as an allele belonging to a gene that is present only once amongst all input gene types, gets a phase mate labeled: “GENE*XXX” (e.g., HLA-DRB5*XXX). No other changes were necessary in phasing algorithm to enable these singleton alleles to be placed in correct haplotypes.

[0080] In some cases, an input can contain ambiguous genotyping information. For example, in the HLA region, for genes like MICA or MICB, an HLA typing program (e.g., Kourami) can provide information with an excessively long list of alleles. In an example, for each of the two alleles of MICA, the following information was provided: MICA*242;MICA*241;MICA*009:01 :09;MICA*009:01 :08;MICA*009:01 :07;MICA*009:0 1 :06;MICA*009:01 :05;MICA*009:01 :04;MICA*009:01 :03;MICA*009:01 :02;MICA*009:0 1 :01;MICA*057:01 :04;MICA*057:01 :03;MICA*057:01 :02;MICA*057:01 :01;MICA*237; MIC A*235 ;MIC A* 199;MIC A* 194;MIC A* 193 ;MIC A* 192;MIC A* 191 ;MIC A*224;MIC A* 223 ;MIC A*221 ;MIC A* 186;MIC A* 182;MIC A* 181 ;MIC A* 180;MIC A*219;MIC A*217;MI CA*212;MICA*210;MICA*176;MICA*175;MICA*174;MICA*216:01 :02;MICA*173;MIC A*216:01 :01;MICA*105Q;MICA*004:04;MICA*004:03;MICA*205;MICA*004:02;MICA *204;MICA*049:01 :02;MICA*203;MICA*049:01 :01;MICA*009:02:06;MICA*009:02:05;MIC A* 009 : 02 : 04;MIC A* 009 : 02 : 03 ;MIC A* 163 ;MIC A* 009 : 02 : 02;MIC A* 130: 02;MIC A* 00 9:02:01;MICA*130:01;MICA*098:01;MICA*160;MICA*159;MICA*154;MICA*166:01 :03 ;MIC A* 166:01 : 02;MIC A* 166:01 :01 ;MIC A* 149 ;MIC A* 027 : 01 : 03 ;MIC A* 148 ;MIC A* 027 : 01 :02;MICA*147;MICA*027:01 :01;MICA*146;MICA*145;MICA*144;MICA*141;MICA* 140;MICA*208:01 :03;MICA*208:01 :02;MICA*208:01 :01;MICA* 151 :01 :03;MICA* 151 :01 :02;MICA*151 :01 :01;MICA*138;MICA*016:01 :03;MICA*136;MICA*016:01 :02;MICA*0 16:01 :01;MICA*134;MICA*097;MICA*094;MICA*128;MICA*127;MICA*126;MICA*00 9:03 :02;MICA*009:03 :01 ;MICA* 123;MICA* 122;MICA* 120;MICA*087;MICA*085;MIC A*082;MIC A*080;MIC A* 117;MIC A* 116;MIC A* 115 ;MIC A*008 : 01 : 17;MIC A* 114;MIC A *019:01 :04;MICA*008:01 : 16;MICA*113;MICA*019:01 :03;MICA*008:01 : 15;MICA*019:01 :02;MICA*008:01 : 14;MICA*019:01 :01;MICA*008:01 : 13;MICA*008:01 : 12;MICA*077; MICA*008:01 : 1 l;MICA*076;MICA*008:01 : 10;MICA*074;MICA*004:01 : 13;MICA*004:0 1 : 12;MICA*004:01 : 1 l;MICA*004:01 : 10;MICA* 109;MICA*016:03;MICA* 108;MICA*016 :02;MICA*008:01 :09;MICA*008:01 :08;MICA*008:01 :07;MICA* 104;MICA*008:01 :06;MI CA*103;MICA*008:01:05;MICA*102;MICA*008:01:04;MICA*008:01:03;MICA*008:01:0 2;MICA*067;MICA*008:01 :01;MICA*004:01 :09;MICA*004:01 :08;MICA*004:01 :07;MIC A*004:01 :06;MICA*004:01 :05;MICA*004:01 :04;MICA*004:01 :03;MICA*004:01 :02;MIC A*004:01 :01;MICA*058;MICA*056;MICA*051;MICA*027:02;MICA*048;MICA*042;MI CA*064N;MICA*033;MICA*032;MICA*031;MICA*028;MICA*124:01:02;MICA*124:01: 01;MICA*024;MICA*008:23:02;MICA*008:23:01;MICA*006;MICA*005;MICA*008:31; MICA*049:02;MICA*008:29;MICA*008:28;MICA*008:27;MICA*008:26;MICA*008:24; MICA*008:22;MICA*008:21;MICA*008:20;MICA*008: 19;MICA*008: 18;MICA*008:17; MICA*008: 16;MICA*008: 15;MICA*008: 14;MICA*008: 13;MICA*008: 1 l;MICA*008:09; MICA*008:08;MICA*008:06;MICA*008:05;MICA*008:03;MICA*008:02;MICA* 119:01 :0 2;MIC A* 119:01 :01 ;MIC A*291 ;MIC A*284;MIC A*283 ;MIC A*282;MIC A*280;MIC A*088 N;MIC A*277;MIC A*274;MIC A*273 ;MIC A*271 ;MIC A*019 : 03 ;MIC A*019 : 02;MIC A*008 :04:07;MICA*008:04:06;MICA*008:04:05;MICA*008:04:04;MICA*008:04:03;MICA*008: 04:02;MICA*008:04:01;MICA*269;MICA*268;MICA*267;MICA*264;MICA*263;MICA* 262;MICA*261;MICA*009:08;MICA*009:07;MICA*009:05;MICA*009:04;MICA*255;MI CA*009:01 :l l;MICA*009:01 : 10;MICA*248. Such a list of alleles can be aggregated under a single “meta-allele” and given a unique label consisting of the first allele in the list followed by the number of additional ambiguous alleles in the list: “FIRST ALLELE+N”, where N is the number of additional alleles in the list. For the above example, the label could be “MICA*242+228”. All sequences for the ambiguous alleles belonging to each “meta-allele” can be kmerized as usual, and any shared kmers between meta-alleles for the same gene (e.g. MICA) can be removed to ensure that only kmers unique to each meta-allele are used for phasing.

[0081] In some cases, imbalanced alleles can be removed. For example, in the HL A region, the allele sequences in the IMGT / HLA database consist of either full genomic sequences or partial sequences (e.g. just the exons, no intronic / intergenic sequence). Often for common, well-characterized genes, both genomic and partial sequences are available. Rarer alleles may only have partial sequences available. This can present a problem when the input types for a given gene are a mix of common and rare alleles such that the common allele has a full genomic sequence, including introns, and / or the rare allele has just exonic sequence. In sucha case, the common allele will get many more kmers associated with it compared to the rare allele. Furthermore, in such a case it can be unclear whether the intronic kmers are actually unique to the common allele and would not be observed in the rare allele. Without full genomic sequences for the rare allele, one cannot know which kmers are unique to both alleles. Therefore, such “imbalanced” alleles for each gene can be identified. When one allele of a pair has only partial sequences available in IMGT / HLA, the partner allele can also only use partial sequences from IMGT / HLA.Long Range Haplotype Phasing

[0082] In various aspects of methods provided herein, methods for generating read sets are provided, including phased read-sets, for applications including genome assembly and haplotype phasing, using long-read or short-read sequencing technologies. Example techniques include but are not limited to proximity ligation techniques such as Hi-C, Chicago, Micro-C, and Omni-C. Nucleic acid molecules can be bound (e.g., in a chromatin structure), cleaved to expose internal ends, re-attached at junctions to other exposed ends, freed from binding, and sequenced. This technique can produce nucleic acid molecules comprising multiple sequence segments. The multiple sequence segments within a nucleic acid molecule can have phase information preserved while being rearranged relative to their natural or starting position and orientation. Sequence segments on either side of a junction can be confidently considered to come from the same phase of a sample nucleic acid molecule.

[0083] Nucleic acid molecules, including high molecular weight DNA, can be bound or immobilized on at least one nucleic acid binding moiety. For example, DNA assembled into in vitro chromatin aggregates and fixed with formaldehyde treatment are consistent with methods herein. Nucleic acid binding or immobilizing approaches include, but are not limited to, in vitro or reconstituted chromatin assembly, native chromatin, DNA-binding protein aggregates, nanoparticles, DNA-binding beads, or beads coated using a DNA-binding substance, polymers, synthetic DNA-binding molecules or other solid or substantially solid affinity molecules. In some cases, the beads are solid phase reversible immobilization (SPRI) beads (e.g., beads with negatively charged carboxyl groups such as Beckman-Coulter Agencourt AMPure XP beads).

[0084] Nucleic acids bound to a nucleic acid binding moiety such as those described herein can be held such that a nucleic acid molecule having a first segment and a second segment separated on the nucleic acid molecule by a distance greater than a read distance on a sequencing device (10 kb, 50 kb, 100 kb or greater, for example) are bound togetherindependent of their common phosphodiester bonds. Upon cleavage of such a bound nucleic acid molecule, exposed ends of the first segment and the second segment may ligate to one another. In some cases, the nucleic acid molecules are bound at a concentration such that there is little or no overlap between bound nucleic acid molecules on a solid surface, such that exposed internal ends of cleaved molecules are likely to re-ligate or become reattached only to exposed ends from other segments that were in phase on a common nucleic acid source prior to cleavage. Consequently, a DNA molecule can be cleaved, and cleaved exposed internal ends can be re-ligated, for example at random, without loss of phase information.

[0085] A bound nucleic acid molecule can be cleaved to expose internal ends through one of any number of enzymatic and non-enzymatic approaches. For example, a nucleic acid molecule can be digested using a restriction enzyme, such as a restriction endonuclease that leaves a single stranded overhang. Mbol digest, for example, is suitable for this purpose, although other restriction endonucleases are contemplated. Lists of restriction endonucleases are available, for example, in most molecular biology product catalogues. Other non-limiting techniques for nucleic acid cleavage include using a transposase, tagmentation enzyme complex, topoisomerase, nonspecific endonuclease, DNA repair enzyme, RNA-guided nuclease, fragmentase, or alternate enzyme. Transposase, for example, can be used in combination with unlinked left and right borders to create a sequence-independent break in a nucleic acid that is marked by attachment of transposase-delivered oligonucleotide sequence. Physical methods can also be used to generate cleavage, including mechanical methods (e.g., sonication, shear), thermal methods (e.g., temperature change), or electromagnetic methods (e.g., irradiation, such as UV irradiation).

[0086] Immobilization of nucleic acids at this stage can keep the cleaved nucleic acid molecule fragments in close physical proximity, such that phase information for the initial molecule is preserved. A benefit of the fixation, e.g., to chromatin aggregates, is that separate regions of a common nucleic acid molecule can be held together independent of their phosphodiester backbone, such that their phase information is not lost upon cleavage of the phosphodiester backbone. This benefit is also conveyed through alternate scaffolds to which a nucleic acid molecule is attached prior to cleavage.

[0087] Optionally, single stranded “sticky” end overhangs are modified to prevent reannealing and re-ligation. For example, sticky ends are partially filled-in, such as by adding one nucleotide and a polymerase. In this way, the entire single-stranded end cannot be filled in, but the end is modified to prevent re-ligation with a formerly complementary end. In theexample of Mbol digestion, which leaves a 5’ GATC 5-prime overhang, only the Guanosine nucleotide triphosphate is added. This results in only a “G” fill-in of the first complementary base (“C”) and result in a 5’ GAT overhang. This operation renders the free sticky ends incompatible for re-ligation to one another but preserves sticky ends for downstream applications. Alternately, blunt ends are generated through completely filling in the overhangs, restriction digest with blunt-end generating enzymes, treatment with a singlestrand DNA exonuclease, or nonspecific cleavage. In some cases, a transposase is used to attach adapter ends having blunt or sticky ends to the exposed internal ends of the DNA molecule.

[0088] Optionally, a “punctuation oligonucleotide” is introduced. This punctuation oligonucleotide marks cleavage / re-ligation sites. Some punctuation oligonucleotides have single-stranded overhangs on both ends that are compatible with the partially filled-in overhangs generated on the exposed nucleic acid sample internal ends. An example of a punctuation oligonucleotide is shown below. In some cases, the double-stranded oligonucleotide having single-stranded overhangs is modified, such as by 5’ phosphate removal at its 5’ ends, so that it cannot form concatemers during ligation. Alternately, blunt punctuation oligonucleotides are used, or cleavage sites are not marked using a distinct punctuation oligonucleotide. In some systems, such as when a transposase is used, punctuation is accomplished through addition of transpososome border sequences, followed by ligation of border sequences to one another or to a punctuation oligo. An example punctuation oligo is presented below. However, alternate punctuation oligos are consistent with the disclosure herein, varying in sequence, length, overhang presence or sequence, or modification such as 5’ de-phosphorylation.

[0089] In some cases, the double-stranded region of the punctuation oligonucleotide will vary. A relevant feature of the punctuation oligonucleotide is the sequence of its overhang, allowing ligation to the nucleic acid sample but optionally modified precluding auto-ligation or concatemer formation. It is often desired that the punctuation oligonucleotide comprise sequence that does not occur or is less likely to occur in a target nucleic acid molecule, such that it is easily identified in a downstream sequence reaction. Punctuation oligos are optionally barcoded, for example with a known barcode sequence or with a randomly generated unique identifier sequence. Unique identifier sequences can be designed to make it highly unlikely for multiple junctions in a nucleic acid molecule or in a sample to be barcoded with the same unique identifier.

[0090] Cleaved ends can be attached to one another directly or through an oligo (e.g., a punctuation oligo), for example using a ligase or similar enzyme. Ligation can proceed such that the free single-stranded ends of an immobilized high-molecular weight nucleic acid molecule are ligated directly or to the punctuation oligonucleotide. Because the punctuation oligonucleotide, if utilized, can have two ligatable ends, this ligation can effectively chain regions of the high molecular weight nucleic acid molecule together. Alternative approaches resulting in affixing a punctuating sequence or molecule between two exposed ends can also be employed, as can approaches for directly connecting two exposed ends without punctuation.

[0091] Nucleic acids can then be liberated from the nucleic acid binding moiety. In the case of in vitro chromatin aggregates, this can be accomplished by reversing the cross-links, or digesting the protein components, or both reversing the crosslinking and digesting protein components. A suitable approach is treatment of complexes with proteinase K, though many alternatives are also contemplated. For other binding techniques, suitable methods can be employed, such as the severing of linker molecules or the degradation of a substrate.

[0092] Nucleic acid molecules resulting from such techniques can have a variety of relevant features. Sequence segments within a nucleic acid molecule can be rearranged relative to their natural or starting positions and orientations, but with phase information preserved. Consequently, sequence segments on either side of a junction can be confidently assigned to a common phase of a common sample molecule. Thus, segments far removed from one another on a molecule can be, by such techniques, brought together or in proximity such that portions or the entirety of each segment is sequenced in a single run of a single molecule sequencing device, allowing definitive phase assignment. Alternately, in some cases originally adjacent segments can become separated from one in the resultant nucleic acid. In some cases, the nucleic acid molecules can be re-ligated such that at least about 50%, 55%, 60%, 65%, 70%, 75%, 80%, 85%, 90%, 91%, 92%, 93%, 94%, 95%, 96%, 97%, 98%, 99%, 99.9%, 99.99%, 99.999%, or 100% of re-ligations are between segments that were in phase on a common nucleic acid source prior to cleavage.

[0093] Another relevant feature of the resultant molecules is that, in some cases, most or all the original molecular sequence is preserved, though in some cases rearranged, in the final punctuated or rearranged molecule. For example, in some cases no more than 1%, 2%, 3%, 4%, 5%, 10%, 15%, or 20% of the original molecule is lost in producing the resultant molecule or molecules. Consequently, in addition to being useful as a phase determinant, the resultant molecule retains a substantial proportion of the original molecule sequence, suchthat the resultant molecule is optionally used to concurrently generate sequence information such as contig information useful in de novo sequencing or as independent verification of previously generated contig information.

[0094] Another feature of libraries of some resultant molecules is that cleavage junctions are not common to multiple members of a population of resultant molecules. That is, that different copies of the same starting nucleic acid molecule can end up with different patterns of junction and rearrangement. Random cleavage junctions can be generated with a nonspecific cleavage molecule, or through variation in restriction endonuclease selection or digestion parameters.

[0095] A consequence of having molecule-specific cleavage sites is that in some cases punctuation oligonucleotides are optionally excluded from the process that results in the ‘punctuation molecule’ re-shuffling and re-ligation to no ill effect. By aligning segments of three or more reshuffled molecules, one observes that cleavage sites are readily identified by their absence in the majority of other members of a library. That is, when three or more reshuffled molecules are locally aligned, a segment can be found to be common to all of the molecules, but the edges of the segment can vary among the molecules. By noting where segment local sequence similarity ends, one can map cleavage junctions in an ‘unpunctuated’ rearranged nucleic acid molecule.

[0096] The resulting nucleic acid molecules can be sequenced, for example on a long-read sequencer. The resulting sequence reads contain segments that alternate between nucleic acid sequence from the original input molecule and, if they are used, sequences of the punctuation oligo. These reads can be processed by a computer to split sequence data from each read using the punctuation oligonucleotide sequence or are otherwise processed to identify junctions. The sequence segments within each read can be segments from a single input high molecular weight DNA molecule. The original nucleic acid molecule can comprise a genome sequence or fraction thereof, such as a chromosome. The sets of segment reads can be discontinuous in the original nucleic acid molecule but reveal long-range, haplotype-phased data. These data can be used for de novo genome assembly and phasing heterozygous positions in the input genome. Sequence between junctions indicates contiguous nucleic acid sequence in the source nucleic acid sample, while sequence across a junction is indicative of a nucleic acid segment that is in phase in the nucleic acid sample but that may be far removed in the arranged scaffold from the adjacent segment.

[0097] Junctions can be identified by a variety of approaches. If punctuation oligos are used, junctions can be identified at reads containing the punctuation oligo sequence. Alternately,junctions can be identified by comparison to a second sequence source (and, often, a third sequence source) for a nucleic acid molecule, such as a previously generated contig sequence dataset or a second, independently generated DNA chain molecule having independently derived junctions. As the sequence is aligned, for example, the quality or confidence of alignment to a particular location can indicate where one segment ends and another begins. If restriction enzymes are used to generate cleavages, sequences containing the restriction enzyme recognition site can be evaluated for potentially containing a junction. Note that not every restriction enzyme recognition site may contain a junction, as some restriction enzyme recognition sites may not have been physically accessible by the enzyme while the nucleic acid was bound to the support, for example. Statistical information can also be employed in identifying junctions; for example, the length segments between junctions may be predicted to be of a certain average value or to follow a certain distribution.

[0098] A benefit of the manipulations herein is that they can preserve molecular phase information while bringing nonadj acent regions of the molecule in proximity such that they are included in a single nucleic acid molecule at a distance suitable for sequencing in a single read, such as a long read. Thus, regions that are separated in the starting sample by greater than the distance of a single long read operation (for example 10 kb, 15 kb, 20 kb, 30 kb, 50 kb, 100 kb or greater) are brought into local proximity such that they are within the distance covered by a single read of a long-range sequencing reaction. Thus, regions that are separated by more than the range of the sequencing technology for a single read in the original sample are read in a single reaction in the phase-preserved, rearranged molecule.

[0099] Resultant rearranged molecules can be sequenced, and their sequence information mapped to independently or concurrently generated sequence reads or contig information, or to a reference genome sequence (for example, the sequence of the human genome).Segments adjacent on the resultant rearranged molecule reads are presumed to be in phase. Accordingly, when these segments are mapped to disparate contigs or long range sequence reads, the reads are assigned to a common phase of a common molecule in the sequence assembly.

[0100] Alternately, if multiple independently generated resultant rearranged molecules are sequenced concurrently, phased sample data is optionally generated from these molecules alone, such that segment sequences separated by junctions are inferred to be in phase, while sequences not separated by junctions are inferred to represent stretches of nucleic acids contiguous in the sample itself and useful for, for example, de novo sequence determination as well as being useful for phase determination. However, additionally or as an alternative,multiple independently generated resultant rearranged molecules sequenced concurrently can still be compared to independently generated scaffold or contig information

[0101] Methods and compositions presented herein can preserve long-range phase information, particularly for molecule segments separated by greater than the length of a read in a sequencing technology (10 kb, 20 kb, 50 kb, 100 kb, 500 kb or greater, for example), while providing such nonadj acent segments in a rearranged or often ‘punctuated’ molecule where the segments are adjacent or close enough to be covered by a single read.

[0102] In some instances, resultant rearranged molecules are combined with native molecules for sequencing. The native molecules can be recognized and utilized informatically by the lack of punctuation sequences, if employed. Native molecules are sequenced using short or long read technology, and their assembly is guided by the phase information and segment sequence information generated through sequencing of the rearranged molecule or library.Haplotype Phasing

[0103] In various aspects of methods provided herein, it often important to know which allelic variants are linked on the same chromosome. This is known as the haplotype phasing. Short reads from high-throughput sequence data rarely allow one to directly observe which allelic variants are linked. Computational inference of haplotype phasing can be unreliable at long distances. The disclosure provides one or more methods that allow for determining which allelic variants are linked using allelic variants on read pairs. In some cases, phasing with methods of the present disclosure is conducted without imputation.

[0104] In various embodiments, the methods and compositions of the disclosure enable the haplotype phasing of diploid or polyploid genomes with regard to a plurality of allelic variants. The methods described herein can thus provide for the determination of linked allelic variants that are linked based on variant information from read pairs and / or assembled contigs using the same. Examples of allelic variants include, but are not limited to, those that are known from the lOOOgenomes, UK10K, HapMap and other projects for discovering genetic variation among humans. Disease association to a specific gene can be revealed more easily by having haplotype phasing data as demonstrated, for example, by the finding of unlinked, inactivating mutations in both copies of SH3TC2 leading to Charcot-Marie-Tooth neuropathy (Lupski JR, Reid JG, Gonzaga- Jauregui C, et al. N. Engl. J. Med. 362: 1181-91, 2010) and unlinked, inactivating mutations in both copies of ABCG5 leading to hypercholesterolemia 9 (Rios J, Stein E, Shendure J, et al. Hum. Mol. Genet. 19:4313-18, 2010).

[0105] Humans are heterozygous at an average of 1 site in 1,000. In some cases, a single lane of data using high-throughput sequencing methods can generate at least about 150,000,000 read pairs. Read pairs can be about 100 base pairs long. From these parameters, one-tenth of all reads from a human sample is estimated to cover a heterozygous site. Thus, on average one-hundredth of all read pairs from a human sample is estimated to cover a pair of heterozygous sites. Accordingly, about 1,500,000 read pairs (one-hundredth of 150,000,000) provide phasing data using a single lane. With approximately 3 billion bases in the human genome, and one in one-thousand being heterozygous, there are approximately 3 million heterozygous sites in an average human genome. With about 1,500,000 read pairs that represent a pair of heterozygous sites, the average coverage of each heterozygous site to be phased using a single lane of a high-throughput sequence method is about (IX), using a typical high-throughput sequencing machine. A diploid human genome can therefore be reliably and completely phased with one lane of a high-throughput sequence data relating sequence variants from a sample that is prepared using the methods disclosed herein. In some examples, a lane of data can be a set of DNA sequence read data. In further examples, a lane of data can be a set of DNA sequence read data from a single run of a high-throughput sequencing instrument.

[0106] As the human genome comprises two homologous sets of chromosomes, understanding the true genetic makeup of an individual requires delineation of the maternal and paternal copies or haplotypes of the genetic material. Obtaining a haplotype in an individual is useful in several ways. First, haplotypes are useful clinically in predicting outcomes for donor-host matching in organ transplantation and are increasingly used as a method to detect disease associations. Second, in genes that show compound heterozygosity, haplotypes provide information as to whether two deleterious variants are located on the same allele, greatly affecting the prediction of whether inheritance of these variants is harmful. Third, haplotypes from groups of individuals have provided information on population structure and the evolutionary history of the human race. Lastly, recently described widespread allelic imbalances in gene expression suggest that genetic or epigenetic differences between alleles may contribute to quantitative differences in expression. An understanding of haplotype structure will delineate the mechanisms of variants that contribute to allelic imbalances.Sequencing

[0107] In various embodiments of methods provided herein, suitable sequencing methods described herein or other suitable methods will be used to obtain sequence information fromnucleic acid molecules within a sample. Sequencing can be accomplished through classic Sanger sequencing methods. Sequence can also be accomplished using high-throughput systems some of which allow detection of a sequenced nucleotide immediately after or upon its incorporation into a growing strand, e.g., detection of sequence in real time or substantially real time. In some cases, high-throughput sequencing generates at least 1,000, at least 5,000, at least 10,000, at least 20,000, at least 30,000, at least 40,000, at least 50,000, at least 100,000 or at least 500,000 sequence reads per hour; where the sequencing reads can be at least about 50, about 60, about 70, about 80, about 90, about 100, about 120, about 150, about 180, about 210, about 240, about 270, about 300, about 350, about 400, about 450, about 500, about 600, about 700, about 800, about 900, or about 1000 bases per read.

[0108] Sequencing can be whole-genome, with or without enrichment of particular regions of interest. Sequencing can be targeted to particular regions of the genome. Regions of the genome that can be enriched for or targeted include but are not limited to single genes (or regions thereof), gene panels, gene fusions, human leukocyte antigen (HLA) loci (e.g., Class I HLA- A, B, and C; Class II HLA-DRB1 / 3 / 4 / 5, HLA-DQA1, HLA-DQB1, HLA-DPA1, HLA-DPB1), exonic regions, exome, and other loci. Genomic regions can be relevant to immune response, immune repertoire, immune cell diversity, transcription (e.g., exome), cancers (e.g., BRCA1, BRCA2, panels of genes or regions thereof such as hotspot regions, somatic variants, SNVs, amplifications, fusions, tumor mutational burden (TMB), microsatellite instability (MSI)), cardiac diseases, inherited diseases, and other diseases or conditions. A variety of methods can be used to enrich for or target regions of interest, including but not limited to sequence capture. In some cases, Capture Hi-C (CHi-C) or CHi- C-like protocols are employed, employing a sequence capture operation (e.g., by target enrichment array) before or after library preparation.

[0109] In some embodiments, high-throughput sequencing involves the use of technology available by Illumina’s Genome Analyzer IIX, MiSeq personal sequencer, or HiSeq systems, such as those using HiSeq 2500, HiSeq 1500, HiSeq 2000, or HiSeq 1000 machines. These machines use reversible terminator-based sequencing by synthesis chemistry. These machines can do 200 billion DNA reads or more in eight days. Smaller systems may be utilized for runs within 3, 2, 1 days or less time.

[0110] In some embodiments, high-throughput sequencing involves the use of technology available by ABI Solid System. This genetic analysis platform that enables massively parallel sequencing of clonally-amplified DNA fragments linked to beads. The sequencing methodology is based on sequential ligation with dye-labeled oligonucleotides.[OHl] The next generation sequencing can comprise ion semiconductor sequencing (e.g., using technology from Life Technologies (Ion Torrent)). Ion semiconductor sequencing can take advantage of the fact that when a nucleotide is incorporated into a strand of DNA, an ion can be released. To perform ion semiconductor sequencing, a high-density array of micromachined wells can be formed. Each well can hold a single DNA template. Beneath the well can be an ion sensitive layer, and beneath the ion sensitive layer can be an ion sensor. When a nucleotide is added to a DNA, H+ can be released, which can be measured as a change in pH. The H+ ion can be converted to voltage and recorded by the semiconductor sensor. An array chip can be sequentially flooded with one nucleotide after another. No scanning, light, or cameras can be required. In some cases, an IONPROTON™ Sequencer is used to sequence nucleic acid. In some cases, an IONPGM™ Sequencer is used. The Ion Torrent Personal Genome Machine (PGM). The PGM can do 10 million reads in two hours.

[0112] In some embodiments, high-throughput sequencing involves the use of technology available by Helicos BioSciences Corporation (Cambridge, Massachusetts) such as the Single Molecule Sequencing by Synthesis (SMSS) method. SMSS is unique because it allows for sequencing the entire human genome in up to 24 hours. Finally, SMSS is described in part in US Publication Application Nos. 20060024711; 20060024678; 20060012793; 20060012784; and 20050100932.

[0113] In some embodiments, high-throughput sequencing involves the use of technology available by 454 Lifesciences, Inc. (Branford, Connecticut) such as the PicoTiterPlate device which includes a fiber optic plate that transmits chemiluminescent signal generated by the sequencing reaction to be recorded by a CCD camera in the instrument. This use of fiber optics allows for the detection of a minimum of 20 million base pairs in 4.5 hours.

[0114] Methods for using bead amplification followed by fiber optics detection are described in Marguiles, M., et al. “Genome sequencing in microfabricated high-density picolitre reactors,” Nature, doi: 10.1038 / nature03959; and well as in US Publication Application Nos. 20020012930; 20030068629; 20030100102; 20030148344; 20040248161; 20050079510, 20050124022; and 20060078909.

[0115] In some embodiments, high-throughput sequencing is performed using Clonal Single Molecule Array (Solexa, Inc.) or sequencing-by-synthesis (SBS) utilizing reversible terminator chemistry. These technologies are described in part in US Patent Nos. 6,969,488; 6,897,023; 6,833,246; 6,787,308; and US Publication Application Nos. 20040106110; 20030064398; 20030022207; and Constans, A., The Scientist 2003, 17(I3):36.

[0116] The next generation sequencing technique can comprise real-time (SMRT™) technology by Pacific Biosciences. In SMRT, each of four DNA bases can be attached to one of four different fluorescent dyes. These dyes can be phospho linked. A single DNA polymerase can be immobilized with a single molecule of template single stranded DNA at the bottom of a zero-mode waveguide (ZMW). A ZMW can be a confinement structure which enables observation of incorporation of a single nucleotide by DNA polymerase against the background of fluorescent nucleotides that can rapidly diffuse in an out of the ZMW (in microseconds). It can take several milliseconds to incorporate a nucleotide into a growing strand. During this time, the fluorescent label can be excited and produce a fluorescent signal, and the fluorescent tag can be cleaved off. The ZMW can be illuminated from below. Attenuated light from an excitation beam can penetrate the lower 20-30 nm of each ZMW. A microscope with a detection limit of 20 zepto liters (20x 10'21liters) can be created. The tiny detection volume can provide 1000-fold improvement in the reduction of background noise. Detection of the corresponding fluorescence of the dye can indicate which base was incorporated. The process can be repeated.

[0117] In some cases, the next generation sequencing is nanopore sequencing (see, e.g., Soni GV and Mell er A. (2007) Clin Chem 53: 1996-2001). A nanopore can be a small hole, of the order of about one nanometer in diameter. Immersion of a nanopore in a conducting fluid and application of a potential across it can result in a slight electrical current due to conduction of ions through the nanopore. The amount of current which flows can be sensitive to the size of the nanopore. As a DNA molecule passes through a nanopore, each nucleotide on the DNA molecule can obstruct the nanopore to a different degree. Thus, the change in the current passing through the nanopore as the DNA molecule passes through the nanopore can represent a reading of the DNA sequence. The nanopore sequencing technology can be from Oxford Nanopore Technologies; e.g., a GridlON system. A single nanopore can be inserted in a polymer membrane across the top of a microwell. Each microwell can have an electrode for individual sensing. The microwells can be fabricated into an array chip, with 100,000 or more microwells (e.g., more than 200,000, 300,000, 400,000, 500,000, 600,000, 700,000, 800,000, 900,000, or 1,000,000) per chip. An instrument (or node) can be used to analyze the chip. Data can be analyzed in real-time. One or more instruments can be operated at a time. The nanopore can be a protein nanopore, e.g., the protein alpha-hemolysin, a heptameric protein pore. The nanopore can be a solid-state nanopore made, e.g., a nanometer sized hole formed in a synthetic membrane (e.g., SiNx, or SiO2). The nanopore can be a hybrid pore (e.g., an integration of a protein pore into a solid-state membrane). The nanopore can be a nanoporewith an integrated sensor (e.g., tunneling electrode detectors, capacitive detectors, or graphene-based nano-gap or edge state detectors (see e.g., Garaj et al. (2010) Nature vol. 67, doi: 10.1038 / nature09379)). A nanopore can be functionalized for analyzing a specific type of molecule (e.g., DNA, RNA, or protein). Nanopore sequencing can comprise “strand sequencing” in which intact DNA polymers can be passed through a protein nanopore with sequencing in real time as the DNA translocates the pore. An enzyme can separate strands of a double stranded DNA and feed a strand through a nanopore. The DNA can have a hairpin at one end, and the system can read both strands. In some cases, nanopore sequencing is “exonuclease sequencing” in which individual nucleotides can be cleaved from a DNA strand by a processive exonuclease, and the nucleotides can be passed through a protein nanopore. The nucleotides can transiently bind to a molecule in the pore (e.g., cyclodextran). A characteristic disruption in current can be used to identify bases.

[0118] Nanopore sequencing technology from GENIA can be used. An engineered protein pore can be embedded in a lipid bilayer membrane. “Active Control” technology can be used to enable efficient nanopore-membrane assembly and control of DNA movement through the channel. In some cases, the nanopore sequencing technology is from NABsys. Genomic DNA can be fragmented into strands of average length of about 100 kb. The 100 kb fragments can be made single stranded and subsequently hybridized with a 6-mer probe. The genomic fragments with probes can be driven through a nanopore, which can create a current-versus- time tracing. The current tracing can provide the positions of the probes on each genomic fragment. The genomic fragments can be lined up to create a probe map for the genome. The process can be done in parallel for a library of probes. A genome-length probe map for each probe can be generated. Errors can be fixed with a process termed “moving window Sequencing By Hybridization (mwSBH).” In some cases, the nanopore sequencing technology is from IBM / Roche. An electron beam can be used to make a nanopore sized opening in a microchip. An electrical field can be used to pull or thread DNA through the nanopore. A DNA transistor device in the nanopore can comprise alternating nanometer sized layers of metal and dielectric. Discrete charges in the DNA backbone can get trapped by electrical fields inside the DNA nanopore. Turning off and on gate voltages can allow the DNA sequence to be read.

[0119] The next generation sequencing can comprise DNA nanoball sequencing (as performed, e.g., by Complete Genomics; see e.g., Drmanac et al. (2010) Science 327: 78-81). DNA can be isolated, fragmented, and size selected. For example, DNA can be fragmented (e.g., by sonication) to a mean length of about 500 bp. Adaptors (Adi) can be attached to theends of the fragments. The adaptors can be used to hybridize to anchors for sequencing reactions. DNA with adaptors bound to each end can be PCR amplified. The adaptor sequences can be modified so that complementary single strand ends bind to each other forming circular DNA. The DNA can be methylated to protect it from cleavage by a type IIS restriction enzyme used in a subsequent step. An adaptor (e.g., the right adaptor) can have a restriction recognition site, and the restriction recognition site can remain non-methylated. The non-methylated restriction recognition site in the adaptor can be recognized by a restriction enzyme (e.g., Acul), and the DNA can be cleaved by Acul 13 bp to the right of the right adaptor to form linear double stranded DNA. A second round of right and left adaptors (Ad2) can be ligated onto either end of the linear DNA, and all DNA with both adaptors bound can be PCR amplified (e.g., by PCR). Ad2 sequences can be modified to allow them to bind each other and form circular DNA. The DNA can be methylated, but a restriction enzyme recognition site can remain non-methylated on the left Adi adaptor. A restriction enzyme (e.g., Acul) can be applied, and the DNA can be cleaved 13 bp to the left of the Adi to form a linear DNA fragment. A third round of right and left adaptor (Ad3) can be ligated to the right and left flank of the linear DNA, and the resulting fragment can be PCR amplified. The adaptors can be modified so that they can bind to each other and form circular DNA. A type III restriction enzyme (e.g., EcoP15) can be added; EcoP15 can cleave the DNA 26 bp to the left of Ad3 and 26 bp to the right of Ad2. This cleavage can remove a large segment of DNA and linearize the DNA once again. A fourth round of right and left adaptors (Ad4) can be ligated to the DNA, the DNA can be amplified (e.g., by PCR), and modified so that they bind each other and form the completed circular DNA template.

[0120] Rolling circle replication (e.g., using Phi 29 DNA polymerase) can be used to amplify small fragments of DNA. The four adaptor sequences can contain palindromic sequences that can hybridize, and a single strand can fold onto itself to form a DNA nanoball (DNB™) which can be approximately 200-300 nanometers in diameter on average. A DNA nanoball can be attached (e.g., by adsorption) to a microarray (sequencing flowcell). The flow cell can be a silicon wafer coated with silicon dioxide, titanium and hexamethyldisilazane (HMDS) and a photoresist material. Sequencing can be performed by unchained sequencing by ligating fluorescent probes to the DNA. The color of the fluorescence of an interrogated position can be visualized by a high-resolution camera. The identity of nucleotide sequences between adaptor sequences can be determined.

[0121] In some embodiments, high-throughput sequencing can take place using AnyDot.chips (Genovoxx, Germany). In particular, the AnyDot.chips allow for lOx - 50xenhancement of nucleotide fluorescence signal detection. AnyDot. chips and methods for using them are described in part in International Publication Application Nos. WO 02088382, WO 03020968, WO 03031947, WO 2005044836, PCT / EP 05 / 05657, PCT / EP 05 / 05655; and German Patent Application Nos. DE 101 49 786, DE 102 14 395, DE 103 56 837, DE 10 2004 009 704, DE 10 2004 025 696, DE 10 2004 025 746, DE 10 2004 025 694, DE 10 2004 025 695, DE 10 2004 025 744, DE 10 2004 025 745, and DE 10 2005 012 301.

[0122] Other high-throughput sequencing systems include those disclosed in Venter, J., et al. Science 16 February 2001; Adams, M. et al. Science 24 March 2000; and M. J. Levene, et al. Science 299:682-686, January 2003; as well as US Publication Application No. 20030044781 and 2006 / 0078937. Overall, such systems involve sequencing a target nucleic acid molecule having a plurality of bases by the temporal addition of bases via a polymerization reaction that is measured on a molecule of nucleic acid, e.g., the activity of a nucleic acid polymerizing enzyme on the template nucleic acid molecule to be sequenced is followed in real time. Sequence can then be deduced by identifying which base is being incorporated into the growing complementary strand of the target nucleic acid by the catalytic activity of the nucleic acid polymerizing enzyme at each step in the sequence of base additions. A polymerase on the target nucleic acid molecule complex is provided in a position suitable to move along the target nucleic acid molecule and extend the oligonucleotide primer at an active site. A plurality of labeled types of nucleotide analogs are provided proximate to the active site, with each distinguishable type of nucleotide analog being complementary to a different nucleotide in the target nucleic acid sequence. The growing nucleic acid strand is extended by using the polymerase to add a nucleotide analog to the nucleic acid strand at the active site, where the nucleotide analog being added is complementary to the nucleotide of the target nucleic acid at the active site. The nucleotide analog added to the oligonucleotide primer as a result of the polymerizing operation is identified. The steps of providing labeled nucleotide analogs, polymerizing the growing nucleic acid strand, and identifying the added nucleotide analog are repeated so that the nucleic acid strand is further extended, and the sequence of the target nucleic acid is determined.Hi-C Methods Using Micrococcal Nuclease (MNase)

[0123] Additionally, provided herein are methods of obtaining phased genotype information, such as HLA type information that may comprise obtaining a stabilized biological sample comprising a nucleic acid molecule complexed to at least one nucleic acid binding protein; contacting the stabilized biological sample to a micrococcal nuclease (MNase) to cleave the nucleic acid molecule into a plurality of segments; and attaching a first segment and a secondsegment of the plurality of segments at a junction. Use of MNase in methods herein may provide specific information about where DNA binding proteins are bound to the chromatin with up to single base pair resolution because, for example, MNase can cleave all base pairs not bound to a DNA binding protein. In addition, use of MNase digestion may allow for creation of contact maps and topologically associated domains to decipher three-dimensional chromatin structural information. In some cases, the MNase may be coupled or fused to an immunoglobulin binding protein or fragment thereof, such as a Protein A, a Protein G, a Protein A / G, or a Protein L.

[0124] For example, MNase Hi-C methods can provide locations of protein binding or genome contact interactions at a resolution of less than or equal to about 1 bp, 2 bp, 3 bp, 4 bp, 5 bp, 6 bp, 7 bp, 8 bp, 9 bp, 10 bp, 20 bp, 30 bp, 40 bp, 50 bp, 60 bp, 70 bp, 80 bp, 90 bp, 100 bp, 200 bp, 300 bp, 400 bp, 500 bp, 600 bp, 700 bp, 800 bp, 900 bp, 1000 bp, 2000 bp, 3000 bp, 4000 bp, 5000 bp, 6000 bp, 7000 bp, 8000 bp, 9000 bp, 10 kb, 20 kb, 30 kb, 40 kb, 50 kb, 60 kb, 70 kb, 80 kb, 90 kb, or 100 kb. In some cases, protein binding sites, protein footprints, contact interactions, or other features can be mapped to within 1000 bp, within 900 bp, within 800 bp, within 700 bp, within 600 bp, within 500 bp, within 400 bp, within300 bp, within 200 bp, within 190 bp, within 180 bp, within 170 bp, within 160 bp, within150 bp, within 140 bp, within 130 bp, within 120 bp, within 110 bp, within 100 bp, within 90 bp, within 80 bp, within 70 bp, within 60 bp, within 50 bp, within 40 bp, within 30 bp, within20 bp, within 10 bp, within 9 bp, within 8 bp, within 7 bp, within 6 bp, within 5 bp, within 4 bp, within 3 bp, within 2 bp, or within 1 bp.

[0125] In certain aspects, methods involving a MNase digestion operation may further comprise subjecting a plurality of segments to size selection to obtain a plurality of selected segments. In some cases, the plurality of selected segments can be from about 145 to about 600 bp. In some cases, the plurality of selected segments can be from about 100 to about 2500 bp. In some cases, the plurality of selected segments can be from about 100 to about 600 bp. In some cases, the plurality of selected segments can be from about 600 to about 2500 bp. In some cases, the plurality of selected segments can be from about 100 bp to about 600 bp, from about 100 bp to about 700 bp, from about 100 bp to about 800 bp, from about 100 bp to about 900 bp, from about 100 bp to about 1000 bp, from about 100 bp to about 1100 bp, from about 100 bp to about 1200 bp, from about 100 bp to about 1300 bp, from about 100 bp to about 1400 bp, from about 100 bp to about 1500 bp, from about 100 bp to about 1600 bp, from about 100 bp to about 1700 bp, from about 100 bp to about 1800 bp, from about 100 bp to about 1900 bp, from about 100 bp to about 2000 bp, from about 100 bpto about 2100 bp, from about 100 bp to about 2200 bp, from about 100 bp to about 2300 bp, from about 100 bp to about 2400 bp, or from about 100 bp to about 2500 bp.

[0126] In another aspect of methods involving a MNase digestion operation as provided herein, the methods may further comprise preparing a sequencing library from the plurality of segments. In some embodiments, the method may further comprise subjecting the sequencing library to a size selection to obtain a size-selected library. In some cases, the size-selected library may be from about 350 bp to about 1000 bp in size. In some cases, the size-selected library may be from about 100 bp to about 2500 bp in size, for example, from about 100 bp to about 350 bp, from about 350 bp to about 500 bp, from about 500 bp to about 1000 bp, from about 1000 to about 1500 bp, from about 2000 bp to about 2500 bp, from about 350 bp to about 1000 bp, from about 350 bp to about 1500 bp, from about 350 bp to about 2000 bp, from about 350 bp to about 2500 bp, from about 500 bp to about 1500 bp, from about 500 bp to about 2000 bp, from about 500 bp to about 3500 bp, from about 1000 bp to about 1500 bp, from about 1000 bp to about 2000 bp, from about 1000 bp to about 2500 bp, from about 1500 bp to about 2000 bp, from about 1500 bp to about 2500 bp, or from about 2000 bp to about 2500 bp.

[0127] In another aspect, methods involving a MNase digestion operation as provided herein can further comprise analyzing the plurality of segments to obtain a QC value. In some cases, a QC value may be selected from a chromatin digest efficiency (CDE) and a chromatin digest index (CDI). A CDE can be calculated as the proportion of segments having a desired length. For example, in some cases, the CDE can be calculated as the proportion of segments from 100 bp to 2500 bp in size prior to size selection. In some cases, a sample may be selected for further analysis when the CDE value is at least 65%. In some cases, a sample may be selected for further analysis when the CDE value is at least about 50%, at least about 55%, at least about 60%, at least about 65%, at least about 70%, at least about 75%, at least about 80%, at least about 85%, at least about 90%, or at least about 95%.

[0128] A CDI can be calculated as a ratio of a number of mononucleosome-sized segments to a number of dinucleosome-sized segments prior to size selection. For example, a CDI may be calculated as a logarithm of the ratio of fragments having a size of 600-2500 bp versus fragments having a size of 100-600 bp. In some cases, a sample may be selected for further analysis when the CDI value is greater than -1.5 and less than 1. In some cases, a sample may be selected for further analysis when the CDI value is greater than about -2 and less than about 1.5, greater than about -1.9 and less than about 1.5, greater than about -1.8 and less than about 1.5, greater than about -1.7 and less than about 1.5, greater than about -1.6 andless than about 1.5, greater than about -1.5 and less than about 1.5, greater than about -1.4 and less than about 1.5, greater than about -1.3 and less than about 1.5, greater than about - 1.2 and less than about 1.5, greater than about -1.1 and less than about 1.5, greater than about -2 and less than about 1.5, greater than about -1 and less than about 1.5, greater than about - 0.9 and less than about 1.5, greater than about -0.8 and less than about 1.5, greater than about -0.7 and less than about 1.5, greater than about -0.6 and less than about 1.5, greater than about -0.5 and less than about 1.5, greater than about -2 and less than about 1.4, greater than about -2 and less than about 1.3, greater than about -2 and less than about 1.2, greater than about -2 and less than about 1.1, greater than about -2 and less than about 1, greater than about -2 and less than about 0.9, greater than about -2 and less than about 0.8, greater than about -2 and less than about 0.7, greater than about -2 and less than about 0.6, or greater than about -2 and less than about 0.5.

[0129] In another aspect, stabilized biological samples used in methods involving a MNase digestion operation as provided herein may comprise biological material that has been treated with a stabilizing agent. In some cases, the stabilized biological sample may comprise a stabilized cell lysate. Alternatively, the stabilized biological sample may comprise a stabilized intact cell. Alternatively, the stabilized biological sample may comprise a stabilized intact nucleus. In some cases, contacting the stabilized intact cell or intact nucleus sample to a MNase may be conducted prior to lysis of the intact cell or the intact nucleus. In some cases, cells and / or nuclei may be lysed prior to attaching a first segment and a second segment of a plurality of segments at a junction.

[0130] In another aspect, methods involving a MNase digestion operation as provided herein may be conducted on small samples containing few cells or small amounts of nucleic acid. For example, in some cases, the stabilized biological sample may comprise fewer than 3,000,000 cells. In some cases, the stabilized biological sample may comprise fewer than2,000,000 cells. In some cases, the stabilized biological sample may comprise fewer than1,000,000 cells. In some cases, the stabilized biological sample may comprise fewer than500,000 cells. In some cases, the stabilized biological sample may comprise fewer than400,000 cells. In some cases, the stabilized biological sample may comprise fewer than300,000 cells. In some cases, the stabilized biological sample may comprise fewer than200,000 cells. In some cases, the stabilized biological sample may comprise fewer than100,000 cells. In some cases, the stabilized biological sample comprises fewer than 50,000 cells. In some cases, the stabilized biological sample comprises fewer than 40,000 cells. In some cases, the stabilized biological sample comprises fewer than 30,000 cells. In somecases, the stabilized biological sample comprises fewer than 20,000 cells. In some cases, the stabilized biological sample comprises fewer than 10,000 cells. In some cases, the stabilized biological sample comprises about 10,000 cells. In some cases, the stabilized biological sample may comprise less than 10 pg DNA. In some cases, the stabilized biological sample may comprise less than 9 pg DNA. In some cases, the stabilized biological sample may comprise less than 8 pg DNA. In some cases, the stabilized biological sample may comprise less than 7 pg DNA. In some cases, the stabilized biological sample may comprise less than 6 pg DNA. In some cases, the stabilized biological sample may comprise less than 5 pg DNA. In some cases, the stabilized biological sample may comprise less than 4 pg DNA. In some cases, the stabilized biological sample may comprise less than 3 pg DNA. In some cases, the stabilized biological sample may comprise less than 2 pg DNA. In some cases, the stabilized biological sample comprises less than 1 pg DNA. In some cases, the stabilized biological sample comprises less than 0.5 pg DNA.

[0131] In another aspect, methods involving a MNase digestion operation herein may be conducted on individual or single cells. For example, methods herein may be conducted on cells distributed into individual partitions. Example partitions include, but are not limited to, wells, droplets in an emulsion, or surface positions (e.g., array spots, beads, etc.) comprising distinct patches of differentially sequenced linker molecules as described elsewhere herein. Additional partitions are also contemplated and consistent with the methods, compositions, and systems disclosed herein.

[0132] In additional aspects, stabilized biological samples used in methods involving a MNase digestion operation herein may be further treated with an additional nuclease, such as a DNase to create fragments of DNA. In some cases, the DNase may be non-sequence specific. In some cases, the DNase may be active for both single-stranded DNA and doublestranded DNA. In some cases, the DNase may be specific for double-stranded DNA. In some cases, the DNase may preferentially cleave double-stranded DNA. In some cases, the DNase may be specific for single-stranded DNA. In some cases, the DNase may preferentially cleave single-stranded DNA. In some cases, the DNase can be DNase I. In some cases, the DNase can be DNase II. In some cases, the DNase may be selected from one or more of DNase I and DNase II. In some cases, the DNase may be coupled or fused to an immunoglobulin binding protein or fragment thereof, such as a Protein A, a Protein G, a Protein A / G, or a Protein L. Other suitable nucleases are also within the scope of this disclosure.

[0133] In additional aspects, stabilized biological samples as provided herein for use in methods involving a MNase digestion operation can be treated with a crosslinking agent. In some cases, the crosslinking agent may be a chemical fixative. In some cases, the chemical fixative comprises formaldehyde, which has a spacer arm length of about 2.3-2.7 angstrom (A). In some cases, the chemical fixative comprises a crosslinking agent with a long spacer arm length, For example, the crosslinking agent can have a spacer length of at least about 3 A, 4 A, 5 A, 6 A, 7 A, 8 A, 9 A, 10 A, 11 A, 12 A, 13 A, 14 A, 15 A, 16 A, 17 A, 18 A, 19 A, or 20 A. The chemical fixative can comprise ethylene glycol bis(succinimidyl succinate) (EGS), which has a spacer arm with length about 16.1 A. The chemical fixative can comprise disuccinimidyl glutarate (DSG), which has a spacer arm with length about 7.7 A. In some cases, the chemical fixative comprises formaldehyde and EGS, formaldehyde and DSG, or formaldehyde, EGS, and DSG. In some cases where multiple chemical fixatives are employed, each chemical fixative is used sequentially; in other cases, some or all of the multiple chemical fixatives are applied to the sample at the same time. The use of crosslinkers with long spacer arms can increase the fraction of read pairs with large (e.g., > 1 kb) read pair separation distances. DSG is membrane-permeable, allowing for intracellular crosslinking. DSG can increase crosslinking efficiency compared to disuccinimidyl suberate (DSS) in some applications. EGS has NHS ester reactive groups at both ends and can be reactive towards amino groups (e.g., primary amines). EGS is membrane-permeable, allowing for intracellular crosslinking. EGS crosslinks can be reversed, for example, by treatment with hydroxylamine for 3 to 6 hours at pH 8.5; in an example, lactose dehydrogenase retained 60% of its activity after reversible crosslinking with EGS. In some cases, the chemical fixative may comprise psoralen. In some cases, the crosslinking agent may be ultraviolet light, chlormethine, cyclophosphamide, chlorambucil, uramustine, melphalan, bendamustine, bis(2-chloroethyl)ethylamine, bis(2-chloroethyl)methylamine, tris(2-chloroethyl)amine, isofamide, carmustine, lomustine, streptozocin, busulfan, cisplatin, carboplatin, cicycloplatin, eptaplatin, lobaplatin, miriplatin, nedaplatin, oxaliplatin, picoplatin, satraplatin, triplatin tetranitrate, procarbazine, altretamine, dacarbazine, mitozolomide, temozolomide, mitomycin C, nitrous acid, formaldehyde, acetylaldehyde, doxorubicin, daunorubicin, epirubicin, or idarubicin. In some cases, the crosslinking agent comprises an intercalating agent, an antibiotic, or a minor groove binding agent. In some cases, the stabilized biological sample may be a crosslinked paraffin-embedded tissue sample.

[0134] In further aspects, methods involving a MNase digestion operation provided herein may comprise contacting the plurality of selected segments to an antibody. In some cases, an immunoglobulin binding protein or fragment thereof tethered to an oligonucleotide adaptor may be targeted to the antibody bound to a plurality of selected segments.

[0135] In additional aspects, methods involving a MNase digestion operation provided herein may comprise attaching a first segment and a second segment of a plurality of segments at a junction. In some cases, attaching may comprise filling in sticky ends using biotin tagged nucleotides and ligating the blunt ends. In some cases, attaching may comprise contacting at least the first segment and the second segment to a bridge oligonucleotide. In some cases, attaching may comprise contacting at least the first segment and the second segment to a barcode. In some embodiments, bridge oligonucleotides herein may be from at least about 5 nucleotides in length to about 50 nucleotides in length. In some embodiments, bridge oligonucleotides herein may be from about 15 to about 18 nucleotides in length. In some embodiments, bridge oligonucleotides may be about 5, about 6, about 7, about 8, about 9, about 10, about 11, about 12, about 13, about 14, about 15, about 16, about 17, about 18, about 19, about 20, about 25, about 30, about 35, about 40, about 45, or about 50 nucleotides in length. In some embodiments, bridge oligonucleotides herein may comprise a barcode.

[0136] In further aspects of methods involving a MNase digestion operation herein, methods can comprise obtaining at least some sequence on each side of the junction to generate a first read pair. For example, the methods may comprise obtaining at least about 50 bp, at least about 100 bp, at least about 150 bp, at least about 200 bp, at least about 250 bp, or at least about 300 bp of sequence on each side of the junction to generate a first read pair.

[0137] In additional aspects of methods involving a MNase digestion operation herein, methods can comprise mapping the first read pair to a set of contigs and determining a path through the set of contigs that represents an order and / or orientation to a genome.

[0138] In further aspects of methods involving a MNase digestion operation herein, methods can comprise mapping the first read pair to a set of contigs; and determining, from the set of contigs, a presence of a structural variant or loss of heterozygosity in the stabilized biological sample.

[0139] In additional aspects of methods involving a MNase digestion operation herein, methods can comprise mapping the first read pair to a set of contigs and assigning a variant in the set of contigs to a phase.

[0140] In further aspects of methods involving a MNase digestion operation herein, methods can comprise mapping the first read pair to a set of contigs; determining, from the set ofcontigs, a presence of a variant in the set of contigs, and conducting an operation selected from one or more of: (1) identifying a disease stage, a prognosis, or a course of treatment for the stabilized biological sample; (2) selecting a drug based on the presence of the variant; or (3) identifying a drug efficacy for the stabilized biological sample.Computer systems

[0141] The present disclosure provides computer systems that are programmed to implement methods of the disclosure. FIG. 14 shows a computer system 1401 that is programmed or otherwise configured to obtain a phased HLA type. The computer system 1401 can regulate various aspects of HLA type phasing of the present disclosure, such as, for example, synthesize SNPs and indels obtained paired end read information to produce phase information for an HLA genomic locus. The computer system 1401 can be an electronic device of a user or a computer system that is remotely located with respect to the electronic device. The electronic device can be a mobile electronic device.

[0142] The computer system 1401 includes a central processing unit (CPU, also “processor” and “computer processor” herein) 1405, which can be a single core or multi core processor, or a plurality of processors for parallel processing. The computer system 1401 also includes memory or memory location 1410 (e.g., random-access memory, read-only memory, flash memory), electronic storage unit 1415 (e.g., hard disk), communication interface 1420 (e.g., network adapter) for communicating with one or more other systems, and peripheral devices 1425, such as cache, other memory, data storage and / or electronic display adapters. The memory 1410, storage unit 1415, interface 1420 and peripheral devices 1425 are in communication with the CPU 1405 through a communication bus (solid lines), such as a motherboard. The storage unit 1415 can be a data storage unit (or data repository) for storing data. The computer system 1401 can be operatively coupled to a computer network (“network”) 1430 with the aid of the communication interface 1420. The network 1430 can be the Internet, an internet and / or extranet, or an intranet and / or extranet that is in communication with the Internet. The network 1430 in some cases is a telecommunication and / or data network. The network 1430 can include one or more computer servers, which can enable distributed computing, such as cloud computing. The network 1430, in some cases with the aid of the computer system 1401, can implement a peer-to-peer network, which may enable devices coupled to the computer system 1401 to behave as a client or a server.

[0143] The CPU 1405 can execute a sequence of machine-readable instructions, which can be embodied in a program or software. The instructions may be stored in a memory location,such as the memory 1410. The instructions can be directed to the CPU 1405, which can subsequently program or otherwise configure the CPU 1405 to implement methods of the present disclosure. Examples of operations performed by the CPU 1405 can include fetch, decode, execute, and writeback.

[0144] The CPU 1405 can be part of a circuit, such as an integrated circuit. One or more other components of the system 1401 can be included in the circuit. In some cases, the circuit is an application specific integrated circuit (ASIC).

[0145] The storage unit 1415 can store files, such as drivers, libraries, and saved programs. The storage unit 1415 can store user data, e.g., user preferences and user programs. The computer system 1401 in some cases can include one or more additional data storage units that are external to the computer system 1401, such as located on a remote server that is in communication with the computer system 1401 through an intranet or the Internet.

[0146] The computer system 1401 can communicate with one or more remote computer systems through the network 1430. For instance, the computer system 1401 can communicate with a remote computer system of a user (e.g., a researcher desiring phased HLA haplotypes for an individual). Examples of remote computer systems include personal computers (e.g., portable PC), slate or tablet PC’s (e.g., Apple® iPad, Samsung® Galaxy Tab), telephones, Smart phones (e.g., Apple® iPhone, Android-enabled device, Blackberry®), or personal digital assistants. The user can access the computer system 1401 via the network 1430.

[0147] Methods as described herein can be implemented by way of machine (e.g., computer processor) executable code stored on an electronic storage location of the computer system 1401, such as, for example, on the memory 1410 or electronic storage unit 1415. The machine executable or machine readable code can be provided in the form of software. During use, the code can be executed by the processor 1405. In some cases, the code can be retrieved from the storage unit 1415 and stored on the memory 1410 for ready access by the processor 1405. In some situations, the electronic storage unit 1415 can be precluded, and machine-executable instructions are stored on memory 1410.

[0148] The code can be pre-compiled and configured for use with a machine having a processer adapted to execute the code or can be compiled during runtime. The code can be supplied in a programming language that can be selected to enable the code to execute in a pre-compiled or as-compiled fashion.

[0149] Aspects of the systems and methods provided herein, such as the computer system 1401, can be embodied in programming. Various aspects of the technology may be thoughtof as “products” or “articles of manufacture” typically in the form of machine (or processor) executable code and / or associated data that is carried on or embodied in a type of machine readable medium. Machine-executable code can be stored on an electronic storage unit, such as memory (e.g., read-only memory, random-access memory, flash memory) or a hard disk. “Storage” type media can include any or all of the tangible memory of the computers, processors or the like, or associated modules thereof, such as various semiconductor memories, tape drives, disk drives and the like, which may provide non-transitory storage at any time for the software programming. All or portions of the software may at times be communicated through the Internet or various other telecommunication networks. Such communications, for example, may enable loading of the software from one computer or processor into another, for example, from a management server or host computer into the computer platform of an application server. Thus, another type of media that may bear the software elements includes optical, electrical, and electromagnetic waves, such as used across physical interfaces between local devices, through wired and optical landline networks and over various air-links. The physical elements that carry such waves, such as wired or wireless links, optical links, or the like, also may be considered as media bearing the software. As used herein, unless restricted to non-transitory, tangible “storage” media, terms such as computer or machine “readable medium” refer to any medium that participates in providing instructions to a processor for execution.

[0150] Hence, a machine readable medium, such as computer-executable code, may take many forms, including but not limited to, a tangible storage medium, a carrier wave medium or physical transmission medium. Non-volatile storage media include, for example, optical or magnetic disks, such as any of the storage devices in any computer(s) or the like, such as may be used to implement the databases, etc. shown in the drawings. Volatile storage media include dynamic memory, such as main memory of such a computer platform. Tangible transmission media include coaxial cables; copper wire and fiber optics, including the wires that comprise a bus within a computer system. Carrier-wave transmission media may take the form of electric or electromagnetic signals, or acoustic or light waves such as those generated during radio frequency (RF) and infrared (IR) data communications. Common forms of computer-readable media therefore include for example: a floppy disk, a flexible disk, hard disk, magnetic tape, any other magnetic medium, a CD-ROM, DVD or DVD- ROM, any other optical medium, punch cards paper tape, any other physical storage medium with patterns of holes, a RAM, a ROM, a PROM and EPROM, a FLASH-EPROM, any other memory chip or cartridge, a carrier wave transporting data or instructions, cables or linkstransporting such a carrier wave, or any other medium from which a computer may read programming code and / or data. Many of these forms of computer readable media may be involved in carrying one or more sequences of one or more instructions to a processor for execution.

[0151] The computer system 1401 can include or be in communication with an electronic display 1435 that comprises a user interface (UI) 1440 for providing, for example, synthesizing SNP and indel information obtained from paired end reads of an HLA genomic region. Examples of UI’s include, without limitation, a graphical user interface (GUI) and web-based user interface.

[0152] Methods and systems of the present disclosure can be implemented by way of one or more algorithms. One or more algorithms can be implemented by way of software upon execution by the central processing unit 1405.

[0153] In some cases, the graph alignment operation takes less than 30 hours, 25 hours, 20 hours, 15 hours, 10 hours, 5 hours, 4 hours, 3 hours, 2 hours, 1 hour, 50 minutes, 40 minutes, 30 minutes, 20 minutes, 10 minutes, or 5 minutes of CPU time. In some cases, the entire computational method takes less than 30 hours, 25 hours, 20 hours, 15 hours, 10 hours, 5 hours, 4 hours, 3 hours, 2 hours, 1 hour, 50 minutes, 40 minutes, 30 minutes, 20 minutes, 10 minutes, or 5 minutes of CPU time.Definitions

[0154] Unless defined otherwise, all technical and scientific terms used herein have the same meaning as commonly understood to one of ordinary skill in the art to which this disclosure belongs. Although any methods and reagents similar or equivalent to those described herein can be used in the practice of the disclosed methods and compositions, the example methods and materials are now described.

[0155] As used herein and in the appended claims, the singular forms “a,” “and,” and “the” include plural referents unless the context clearly dictates otherwise. Thus, for example, reference to “contig” includes a plurality of such contigs and reference to “probing the physical layout of chromosomes” includes reference to one or more methods for probing the physical layout of chromosomes and equivalents thereof known to those skilled in the art, and so forth.

[0156] Also, the use of “and” means “and / or” unless stated otherwise. Similarly, “comprise,” “comprises,” “comprising,” “include,” “includes,” and “including” are interchangeable and not intended to be limiting.

[0157] It is to be further understood that where descriptions of various embodiments use the term “comprising,” those skilled in the art may understand that in some specific instances, an embodiment can be alternatively described using language “consisting essentially of’ or “consisting of.”

[0158] The term “sequencing read” as used herein, refers to a fragment of DNA in which the sequence has been determined.

[0159] The term “subject” as used herein can refer to any eukaryotic or prokaryotic organism.

[0160] The term “read pair” or “read-pair” as used herein can refer to two or more elements that are linked to provide sequence information. In some cases, the number of read-pairs can refer to the number of mappable read-pairs. In other cases, the number of read-pairs can refer to the total number of generated read-pairs.

[0161] The term “stabilized” as used herein can describe a sample that has been preserved or otherwise protected from degradation. In some cases, a stabilized sample is crosslinked or treated with a fixative or crosslinking agent. In some cases, a stabilized sample is treated with formaldehyde, formalin, paraformaldehyde, glutaraldehyde, osmium tetroxide, or the like.

[0162] The term “about” as used herein can describe a number, unless otherwise specified, as a range of values including that number plus or minus 10% of that number.EXAMPLES

[0163] The following examples are included for illustrative purposes only and are not intended to limit the scope of the invention.Example 1: Optimal Partition for Int-237hlaA

[0164] Haplotype phasing of the Int-237hlaA region was performed. First, proximity ligation sequencing of the HLA genes was performed. The FASTQ files from this sequencing experiment were then input into the graph-guided assembler Kourami, which identified HLA gene alleles of interest. A read graph was then built, connecting the alleles using overlapping k-mers as previously described. The pairwise allele connectivity using a Markov hitting time metric was then calculated between each possible pair of alleles. A reduced allele graph was then created, with the alleles as nodes and the weight of each edge as the symmetric hitting probability between each pair of alleles. The optimal partition of alleles into haplotypes was then calculated by “brute-force”, e.g., testing all possible partitions of alleles into haplotypes and selecting the best one based on the sum of edge weights in that partition. The weightsindicated the strength of connectivity between the alleles in that partition, which was a proxy for the more strongly-connected alleles belonging to the same haplotype.

[0165] FIG. 11 shows a graphical representation of the optimal partition for the Int-237hlaA region calculated using the methods described herein. This figure shows aggregated connectivity values between each of the eight HLA genes. The color of an edge indicates whether the pairwise phasing was concordant (gray) or non-concordant (red) with an HLA pipeline using an approach of sorting SNPs into phase blocks as described herein. The width of an edge indicates the strength of the phasing. The methods described herein produced a phasing of the HLA region for Int-237hlaA that was identical to that produced with an HLA pipeline using an approach of sorting SNPs into phase blocks as described herein.Example 2: Optimal Partition for Int-237hlaB

[0166] Haplotype phasing of the Int-237hlaB region was performed. First, proximity ligation sequencing of the HLA genes was performed. The FASTQ files from this sequencing experiment were then input into the graph-guided assembler Kourami, which identified HLA gene alleles of interest. A read graph was then built, connecting the alleles using overlapping k-mers as previously described. The pairwise allele connectivity using a Markov hitting time metric was then calculated between each possible pair of alleles. A reduced allele graph was then created, with the alleles as nodes and the weight of each edge as the symmetric hitting probability between each pair of alleles. The optimal partition of alleles into haplotypes was then calculated by “brute-force”, e.g., testing all possible partitions of alleles into haplotypes and selecting the best one based on the sum of edge weights in that partition. The weights indicated the strength of connectivity between the alleles in that partition, which was a proxy for the more strongly-connected alleles belonging to the same haplotype.

[0167] FIG. 12 shows a graphical representation of the optimal partition for the Int-237hlaB region calculated using the methods described herein. This figure shows aggregated connectivity values between each of the eight HLA genes. The color of an edge indicates whether the pairwise phasing was concordant (gray) or non-concordant (red) with an HLA pipeline using an approach of sorting SNPs into phase blocks as described herein. The width of an edge indicates the strength of the phasing. The methods described herein produced a phasing of the HLA region for Int-237hlaB that was identical to that produced with an HLA pipeline using an approach of sorting SNPs into phase blocks as described herein. It was found that the genes HLA-B, HLA-C, and DPA1 are homozygous.Example 3: Optimal Partition for Int-223hlaD-V4

[0168] Haplotype phasing of the Int-223hlaD-V4 region was performed. First, proximity ligation sequencing of the HLA genes was performed. The FASTQ files from this sequencing experiment were then input into the graph-guided assembler Kourami, which identified HLA gene alleles of interest. A read graph was then built, connecting the alleles using overlapping k-mers as previously described. The pairwise allele connectivity using a Markov hitting time metric was then calculated between each possible pair of alleles. A reduced allele graph was then created, with the alleles as nodes and the weight of each edge as the symmetric hitting probability between each pair of alleles. The optimal partition of alleles into haplotypes was then calculated by “brute-force”, e.g. testing all possible partitions of alleles into haplotypes and selecting the best one based on the sum of edge weights in that partition. The weights indicated the strength of connectivity between the alleles in that partition, which was a proxy for the more strongly-connected alleles belonging to the same haplotype.

[0169] FIG. 13 shows a graphical representation of the optimal partition for the Int-223hlaD- V4 region calculated using the methods described herein. This figure shows aggregated connectivity values between each of the eight HLA genes. The color of an edge indicates whether the pairwise phasing was concordant (gray) or non-concordant (red) with an HLA pipeline using an approach of sorting SNPs into phase blocks as described herein. The width of an edge indicates the strength of the phasing. It was found that the gene DRB 1 (green) had no pipeline call. For DRB1, the edge colors correspond to the optimal partition phasing.

[0170] While preferred embodiments of the present invention have been shown and described herein, it will be obvious to those skilled in the art that such embodiments are provided by way of example only. Numerous variations, changes, and substitutions will now occur to those skilled in the art without departing from the invention. It should be understood that various alternatives to the embodiments of the invention described herein may be employed in practicing the invention. It is intended that the following claims define the scope of the invention and that methods and structures within the scope of these claims and their equivalents be covered thereby.

Claims

CLAIMSWHAT IS CLAIMED IS:

1. A method for gene allele analysis, comprising:(a) identifying a plurality of unphased gene alleles in proximity ligation nucleic acid sequencing data;(b) computing a read graph of said plurality of unphased gene alleles, wherein: (i) a first allele of said plurality of unphased gene alleles is a first node of said read graph, (ii) a second allele of said plurality of unphased gene alleles is a second node of said read graph, (iii) said first node and said second node are connected in pairwise combination, and (iv) a weight of each connection between said first node and said second node is based on an amount of sequence overlap between said first allele and said second allele in said sequencing data; and(c) evaluating said weight of each connection between said first node and said second node to determine which of said plurality of unphased gene alleles have the highest likelihood to be in the same phase.

2. The method of claim 1, wherein said plurality of unphased gene alleles comprises a human leukocyte antigen (HLA) gene.

3. The method of claim 2, wherein said HLA gene is HLA- A, HLA-B, HLA-C, DRB1, DQA1, DQB1, DPA1, DPB1, or a combination thereof.

4. The method of any one of claims 1 to 3, wherein said amount of sequence overlap between said first allele and said second allele is about 50 bases, about 100 bases, about 150 bases, about 200 bases, about 250 bases, about 300 bases, or about 350 bases.

5. The method any one of claims 1 to 4, wherein said weight of each connection between said first node and said second node is determined using a Markov hitting time.

6. The method of any one of claims 1 to 5, wherein said weight of each connection between said first node and said second node is calculated as a symmetric hitting probability between said first node and said second node.

7. The method of any one of claims 1 to 6, wherein calculating an optimal read graph comprises solving a system of linear equations.

8. The method of any one of claims 1 to 7, wherein said proximity ligation nucleic acid sequencing data is obtained by crosslinking a sample, fragmenting nucleic acids in said sample to produce nucleic acid fragments, ligating said nucleic acid fragments toproduce ligated nucleic acid fragments, reversing crosslinks, and sequencing said ligated nucleic acid fragments.

9. The method of claim 8, wherein said sample is crosslinked by contacting said sample to a crosslinking agent selected from formaldehyde, psoralen, disuccinimidyl glutarate (DSG), ethylene glycol bis(succinimidyl succinate) (EGS), ultraviolet light, or a combination thereof.

10. The method of claim 8 or claim 9, wherein said fragmenting comprises contacting said sample to an enzyme.

11. The method of claim 10, wherein said enzyme is a nuclease, a restriction endonuclease, a transposase, or a combination thereof.

12. The method of claim 11, wherein said nuclease is a micrococcal nuclease.

13. The method of claim 11, wherein said transposase is Tn5.

14. The method of claim 8 or claim 9, wherein fragmenting comprises non-enzymatic cleavage.

15. The method of any one of claims 8 to 14, further comprising, subsequent to said ligating, adding a label to said nucleic acid fragments.

16. The method of claim 15, wherein said label comprises biotin.

17. The method of claim 15 or claim 16, wherein said label comprises an oligonucleotide.

18. The method of claim 17, wherein said oligonucleotide comprises a barcode.

19. The method of any one of claims 8 to 18, wherein said cross-linking links nucleic acids to nucleic acid binding proteins in said sample.

20. The method of any one of claims 1 to 19, wherein an optimal read graph is calculated by evaluating about 20%, about 30%, about 40%, about 50%, about 60%, about 70%, about 80%, about 90%, about 99%, or greater than 99% of all possible pairwise connections between said plurality of unphased gene alleles.

21. The method of claim 20, wherein an optimal read graph is calculated by evaluating a symmetric hitting probability between two alleles of said plurality of unphased gene alleles over about 20%, about 30%, about 40%, about 50%, about 60%, about 70%, about 80%, about 90%, about 99%, or greater than 99% of all possible pairwise connections between said plurality of unphased gene alleles.

22. The method of any one or claims 1 to 21, wherein said amount of sequence overlap is calculated by comparing at least one k-mer of a selected size from said first node with at least one k-mer of a same selected size from said second node.

23. The method of claim 22, wherein an allele in said plurality of unphased gene alleles is a node in said read graph, and an edge connects two nodes if said at least one k-mer is shared between said two nodes.

24. The method of claim 23, wherein a k-mer is about 50, about 100, about 150, about 200, about 250, about 300, or about 350 bases.

25. The method of any one of claims 1 to 24, wherein latent information about a one or more locations of heterozygous single nucleotide polymorphisms (SNPs) is encoded into said read graph.

26. The method of any one of claims 1 to 25, wherein an allele node connectivity measure between a first of two alleles and a second of said two alleles is determined using a Markov hitting time.

27. The method of claim 26, wherein said allele node connectivity measure is computed bi-directionally between a first of two alleles and a second of two alleles and between said second of two alleles and said first of two alleles.

28. The method of claim 26 or claim 27, wherein said allele connectivity measure is computed for about 20%, about 30%, about 40%, about 50%, about 60%, about 70%, about 80%, about 90%, about 99%, or greater than 99% of all possible pairwise connections between said plurality of unphased gene alleles.

29. The method of claim 28, wherein a reduced allele graph that connects said plurality of unphased gene alleles is formed, wherein each one of said plurality of unphased gene alleles is connected to each other one of said plurality of unphased gene alleles.

30. The method of claim 29, wherein a weight of an edge is calculated as a symmetric hitting probability between each pair of two alleles in said reduced allele graph.

31. The method of claim 30, wherein a calculation of an optimal partition of said plurality of unphased gene alleles into one or more phased haplotypes can be calculated from one or more possible partitions of disjoint sets of said plurality of unphased gene alleles.

32. The method of claim 31, wherein said optimal partition is calculated using a sum of edge weights within each of said one or more possible partitions.

33. The method of claim 32, wherein said optimal partition of said plurality of unphased gene alleles into said one or more phased haplotypes can be calculated from said one or more possible partitions in a brute-force manner by calculating all possible partitions of said plurality of unphased gene alleles.-SO-

Images