Graph-based systems and methods for mapping nucleic acid fragments

Abstract

A technical solution for mapping long nucleic acid sequence reads to a target sequence is provided. For an organism with a heterozygous genome, a directed graph representing all or some of the genome and comprising one or more nonlinear topological components is obtained. Each nonlinear topological component has a starting node and an ending node connected by at least a first branch and a second branch. One of these branches corresponds to the target sequence. The directed graph uses multiple sequence reads from a biological sample of the organism. The sequence reads overlap an unrestricted amount of overhang, provided that there is a minimum common area between every two sequence reads. A query sequence is obtained that contains at least a starting node or an ending node of a first nonlinear topological component. The directed graph is used to form a mapping of the query sequence to the directed graph.

Description

[0001] CROSS-REFERENCE TO RELATED APPLICATIONS

[0002] This application claims priority to U.S. Provisional Patent Application No. 62 / 797,022, filed on January 25, 2019, the disclosure of which is incorporated herein by reference in its entirety for all purposes. Technical Field

[0003] This specification describes the use of graph-based mapping to assemble nucleic acid fragments into coherent sequences. Background Art

[0004] Advances in biomolecular sequencing, particularly for nucleic acid and protein samples, have revolutionized the fields of cell and molecular biology. Due to the development of automated sequencing systems, it is now possible to sequence whole genomes. However, there has been a corresponding increase in the number and complexity of genomes sequenced using these new sequencing technologies. Although large quantities of high-fidelity nucleic acid sequences are now available, there are still many problems in assembling and organizing these sequences into complete genomes.

[0005] Assembling these complex genomes from sequence reads, especially diploid and higher ploidy genomes, is a difficult problem that is necessarily rooted in computational techniques. This is due to a number of factors, including the sheer size of the genomes, read size limitations of traditional sequencing technologies, and sequence variation in homologous regions of polyploid genomes. For example, while the relatively 'simple' haploid genome of Escherichia coli contains only 4.6 million base pairs in a single circular genome, the human genome contains two copies of twenty-three chromosomes, each containing over 3.2 billion base pairs.

[0006] Assuming an average sequencing read length of 250 base pairs and an average overlap between sequence reads of only 15%, 1x coverage of the human genome would require nearly 30 million sequence reads. Furthermore, the human genome includes two copies of each locus, with a sequence variation of approximately 0.1%, meaning that more than 22% of the 250 base pair sequence reads will differ from any particular reference genome used to assist the assembly process. Furthermore, the different 22% of sequence reads, representing approximately 6.5 million sequence reads, need to be phased relative to each other in order to assemble haplogroups for a diploid genome. This is only further complicated in higher-order polyploid genomes such as ferns (average 2n = 121.0; 2n range = 18–1440; Clark J et al., New Phytologist, 210: 1072–82 (2016)). It is therefore clear that genome assembly is virtually impossible without the use of computational techniques. See Kelly MJ. Computers: the best friends a human genome ever had, Genome, 31(2):1027-33(1989).

[0007] However, the computational power of general-purpose computers alone is not all that is required for efficient genome assembly, as evidenced by the massive effort required for the Human Genome Project, even with the help of computational technology (which took 13 years to complete with a global collaboration; Chial, H., DNA sequencing technologies key to the Human Genome Project, Nature Education, 1(1):219 (2008)). Instead, complex alignment algorithms are needed to handle the very large amounts of sequencing data required for genome assembly. The development of more greedy alignment algorithms has greatly accelerated the genome assembly process relative to the first assembled human genome.

[0008] Mapping sequencing reads is embodied in tools such as BLAST to BLASR, BWA-MEM, DAMAPPER, NGMLR, GraphMap, Minimap, etc. Each of these methods has its own specific requirements and uses. However, they all have a common design for mapping to linear target sequences, such as reference genomes. On the other hand, mappers such as BGREAT (as described in Limasset et al. 2016, BMC Bioinformatics 17:237) and deBGA (e.g., as described by Liu et al. 2016, Bioinformatics 32 (21): 3224-3232) are designed to process second-generation sequencing data and de Bruijn graphs, rather than linear target sequences. BlastGraph uses BLAST mapping results to cluster alignments and perform comparative genomic analysis (as described in Ye et al., 2013, Bioinformatics 29 (24): 3222-3224). However, according to experiments performed by Limasset et al., BlastGraph performs poorly even on simple nonlinear datasets, such as the Escherichia coli genome (BMC Bioinformatics 17:237, 2016). GramTools maps short reads to a population reference graph (e.g., as Maciuca et al., 2016, on the internet at dx.doi.org / 10.1101 / 059170). In general, all of these methods are designed specifically for short read sequences.

[0009] Although useful for assembling relatively simple linear genomic regions, short read sequences have many disadvantages. For example, short read sequences provide limited information in genomic regions with high GC content, repetitive regions, and regions showing structural variation (e.g., repeats) (e.g., Pollard et al., 2018, Hum Molecul Genetics 27 (R2): R234–R241). As more complex genomes are sequenced and assembled, it will become increasingly necessary to accurately and effectively use the mapper of long read sequences. Currently, a method for mapping long noise reads to arbitrary sequence graphs is a variation graph (e.g., "vg") toolkit, which provides tools for building, viewing, and manipulating variation graphs, including features such as full-length POA alignment on the graph, mapping reads, and calling variants from read comparisons (e.g., as described in Garrison et al., 2018, Nat. Biotech. 36, 875-879). It was originally designed and optimized for short reads, and was later modified to support longer sequences by dividing them into short overlapping blocks (essentially simulating short reads). However, to implement graph-based functionality, vg uses a BWT-based directed graph indexing method (GCSA2) that only indexes paths of length at most 256 in any graph (e.g., as described by Sirén et al. in Proc. ALENEX 2017, SIAM, pp. 13-27, Barcelona, ​​Spain, January 17-18, 2017. DOI: 10.1137 / 1.9781611974768.2). The upper limit on path length is intended as a heuristic to prevent "combinatorial explosion" in regions with a large number of branches. Vg also expands and unrolls cyclic directed graphs to make directed acyclic graphs (a process of "DAGing") to which partial order alignment can be applied (e.g., in Garrison et al., 2018, Nat. Biotech. 36, 875-879). Unrolling loops in a graph is also an expensive combinatorial operation that can significantly increase the size of the graph (e.g., as described in Kavya et al., 2019, J Comp Biol 26:1-15 and Naga et al., 2019, J Comp Biol 26(1)). Furthermore, while the GCSA2 indexing approach provides linear-time exact match queries independent of graph size to find super-maximal exact match seeds, algorithms based on BWT and FM indexes can often carry large hidden constant factors in the lookup process as well as significant time required to pre-build the index, which needs to be accounted for when considering maximum mapping rates for large sequences and graphs. Thus, while the vg paradigm is well suited for many applications, there are many other types of genomes that require different approaches.

[0010] An example of a use case based on graph mapping is in the context of circular genomes. When faced with sequences of circular genomes, most mappers either demonstrate a significant decrease in coverage near the end of the reference circular sequence (e.g., BLASR, as described in Chaisson and Tesler BMC Bioinformatics 13:238,2012), or maintain high sequencing coverage (they can also be placed in repeat regions) by outputting auxiliary alignments that are not directly dependent on the main alignment (e.g., Minimap2, as described in Li, Bioinformatics 15; 34(18):3094-3100,2018). One way to overcome this defect is to use a graph-based mapping method. Limiting the reference to a graph can provide seemingly transparent support for a wide range of important applications. In fact, Minimap2 and GraphMap do seek to solve the problem of circular genome mapping (e.g., see Id. and Sovi et al., 2016, Nat Comm7:11307). However, an important limitation of GraphMap in particular is that the method expects all sequences in the reference set to be circular or linear. It cannot handle a mixture of linear and circular reference genomes.

[0011] In addition to circular genomes, there are many other sequencing situations that could benefit from improved mapping methods. For example, many mappers have difficulty identifying genome topologies, such as inversions. Most mapping methods also fail to identify haplotypes from short read sequences. In addition, there is increasing interest in mapping to multi-genome maps (e.g., maps consisting of multiple human genomes for variant detection and / or diagnostic purposes, bacterial databases for species identification, or metagenomic assemblies for further assembly improvements and species separation) (e.g., discussed in Vollmers et al., 2017, PLoS ONE 12(1):e0169662).

[0012] Given the growing demand for genome sequencing and assembly of complex genomes, there is a need in the art for improved graph-based long-read mapping methods. Summary of the invention

[0013] The present disclosure addresses the shortcomings of the graph-based mapping methods presented in the background by providing an improved method for determining nucleic acid sequences based on long sequence reads. The method begins with obtaining a directed graph of a single diploid or polyploid organism of a species with a heterozygous genome, wherein the directed graph represents all or part of the heterozygous genome and includes one or more nonlinear topological components. In some embodiments, each corresponding nonlinear topological component in the one or more nonlinear topological components represents a change relative to a set of target sequences.

[0014] In some embodiments, each corresponding nonlinear topological component of the one or more nonlinear topological components has a corresponding starting node and a corresponding ending node, and the starting node and the ending node are connected by at least a corresponding first branch and a corresponding second branch, wherein one of the corresponding first branch and the corresponding second branch corresponds to the target sequence. In some embodiments, a directed graph is formed using multiple sequence reads from a biological sample of a single organism by overlapping corresponding sequence reads in the multiple sequence reads, and the sequence reads have a random error greater than ten percent relative to the target sequence, wherein the amount of protrusion between the overlapping first sequence read and the second sequence is not limited, provided that there is a minimum common region between the first and second sequence reads.

[0015] The method is performed by obtaining a query sequence from a biological sample in an electronic format. The query sequence is (i) a sequence read in a plurality of sequence reads, or (ii) a scaffold formed by a plurality of sequence reads. In addition, the query sequence at least comprises a corresponding start node or a corresponding end node of a first nonlinear topological component in one or more nonlinear topological components.

[0016] The method continues by using the directed graph to form a mapping of the query sequence to the directed graph. The mapping of the query sequence includes (i) a first alignment component to a first portion of the directed graph, (ii) a second alignment component to a second portion of the directed graph, and (iii) a path in the directed graph describing a relationship between the first portion and the second portion of the directed graph. In this case, one of the first portion and the second portion is in a first or second branch of the first nonlinear topology component, and the other of the first portion and the second portion is not in a branch of the first nonlinear topology component.

[0017] Other embodiments relate to systems, portable consumer devices, and computer-readable media associated with the methods described herein.

[0018] Where applicable, any embodiment disclosed herein may be applied to any other aspect.

[0019] Other aspects and advantages of the present disclosure will become apparent to those skilled in the art from the following detailed description, wherein only illustrative embodiments of the present disclosure are shown and described. It should be understood that the present disclosure is susceptible of other and different embodiments, and that its several details may be modified in various obvious respects, all without departing from the present disclosure. Therefore, the drawings and description are to be regarded as illustrative in nature, and not restrictive.

[0020] Incorporated by Reference

[0021] All publications, patents, and patent applications cited herein are incorporated by reference in their entirety. In the event of a conflict in a terminology herein and that of an incorporated reference, the terminology herein controls. BRIEF DESCRIPTION OF THE DRAWINGS

[0022] In the figures of the accompanying drawings, there are illustrated by way of example and not limitation, embodiments disclosed herein. Like reference numerals refer to corresponding parts throughout the several views of the drawings.

[0023] Figure 1 is an example block diagram illustrating a computing device according to some embodiments of the present disclosure.

[0024] Figure 2A and 2B An example flow chart of a method of mapping a query sequence using a directed graph according to some embodiments of the present disclosure is illustrated, wherein optional steps are depicted in dashed boxes.

[0025] Figure 3 An example of phasing haplogroups across an assembled genome alignment according to some embodiments of the present disclosure is illustrated.

[0026] Figure 4A , 4B 4C together show example diagrams according to some embodiments of the present disclosure.

[0027] Figure 5 An example of mapping to a circular genome according to some embodiments of the present disclosure is illustrated.

[0028] Fig. 6A and 6B Shown together are examples of maps with genomic inversions according to some embodiments of the present disclosure.

[0029] Fig. 7A and 7B Together, examples of mapping using known insert sequences (e.g., plasmids, mitochondrial DNA, or targeted genomic regions) are shown according to some embodiments of the present disclosure.

[0030] Fig. 8A and 8B Together, examples of transcriptome mapping via edges defining splice sites are shown according to some embodiments of the present disclosure. DETAILED DESCRIPTION

[0031] Reference will now be made in detail to the embodiments, examples of which are shown in the accompanying drawings. In the following detailed description, many specific details are set forth in order to thoroughly understand the present disclosure. However, it will be apparent to one of ordinary skill in the art that the present disclosure can be implemented without these specific details. In other examples, well-known methods, procedures, components, circuits, and networks are not described in detail to avoid unnecessarily obscuring various aspects of the embodiments.

[0032] Embodiments described herein provide various technical solutions using long read sequence reads.For single diploid or polyploid organisms with heterozygous genomes, a directed graph representing all or part of a genome is obtained.A query sequence (e.g., a sequence of interest) from a biological sample of a corresponding biological species is also obtained.Then a directed graph is used to form a mapping of the query sequence to the directed graph.Therefore, the method can detect structural variants and single nucleotide variants.The method provided herein needs to use a long sequence read.

[0033] definition

[0034] As used herein, the term "about" or "approximately" may mean within an acceptable error range for a particular value determined by a person of ordinary skill in the art, which may depend in part on how the value is measured or determined (e.g., limitations of the measurement system). For example, according to the practice in the art, "about" may mean within or above 1 standard deviation. "About" may mean a range of ±20%, ±10%, ±5%, or ±1% of a given value. The term "about" or "approximately" may mean within an order of magnitude, within 5 times, or within 2 times of a value. Where specific values ​​are described in the application and claims, unless otherwise stated, it should be assumed that the meaning of the term "about" is within an acceptable error range for the specific value. The term "about" may have the meaning commonly understood by a person of ordinary skill in the art. The term "about" may refer to ±10%. The term "about" may refer to ±5%.

[0035] As disclosed herein, the term "subject" refers to any living or inanimate organism, including but not limited to humans (e.g., men, women, fetuses, pregnant women, children, etc.), non-human animals, plants, bacteria, fungi, or protists. Any human or non-human animal can be selected from, including but not limited to mammals, reptiles, birds, amphibians, fish, ungulates, ruminants, bovines (e.g., cattle), equines (e.g., horses), caprines and ovines (e.g., sheep, goats), swine (e.g., pigs), camelids (e.g., camels, llamas, alpacas), monkeys, apes (e.g., gorillas, chimpanzees), ursids (e.g., bears), poultry, dogs, cats, rodents, fish, dolphins, whales, and sharks. The subject can be male or female (e.g., male, female, or child) at any stage.

[0036] As used herein, the term "biological sample", "patient sample" or "sample" refers to any sample taken from a subject that can reflect a biological state associated with the subject. Examples of biological samples include, but are not limited to, blood, whole blood, plasma, serum, urine, cerebrospinal fluid, feces, saliva, sweat, tears, pleural effusion, pericardial fluid or peritoneal fluid of a subject. A biological sample may include any tissue or material derived from a living or dead subject. A biological sample may be a cell-free sample. A biological sample may include a nucleic acid (e.g., DNA or RNA) or a fragment thereof. A sample may be a liquid sample or a solid sample (e.g., a cell or tissue sample). A biological sample may be a body fluid, such as blood, plasma, serum, urine, vaginal fluid, fluid from a hydrocele (e.g., testicular), vaginal washings, pleural effusion, ascites, cerebrospinal fluid, saliva, sweat, tears, sputum, bronchoalveolar lavage fluid, fluid discharged from the nipple, fluid aspirated from different parts of the body (e.g., thyroid, breast), etc. A biological sample may be a stool sample. Biological samples can be processed to physically disrupt tissue or cellular structures (e.g., centrifugation and / or cell lysis), thereby releasing intracellular components into a solution that can further contain enzymes, buffers, salts, detergent solutions, etc., which can be used to prepare the sample for analysis.

[0037] As used herein, the terms "nucleic acid" and "nucleic acid molecule" are used interchangeably. The term refers to any compositional form of nucleic acid, such as deoxyribonucleic acid (DNA, such as complementary DNA (cDNA), genomic DNA (gDNA), etc.), ribonucleic acid (RNA, such as information RNA (mRNA), short inhibitory RNA (siRNA), ribosomal RNA (rRNA), transfer RNA (tRNA), microRNA, RNA highly expressed by the fetus or placenta, etc.) and / or DNA or RNA analogs (e.g., containing base analogs, sugar analogs and / or non-natural main chains, etc.), RNA / DNA hybrids and polyamide nucleic acids (PNA), all of which can be in single-stranded or double-stranded form. Unless otherwise limited, nucleic acids can contain known natural nucleotide analogs, some of which can function in a manner similar to naturally occurring nucleotides. Nucleic acids can be in any form (e.g., linear, circular, supercoiled, single-stranded, double-stranded, etc.) for performing the processes herein. Nucleic acid can be or can be from plasmid, phage, autonomous replication sequence (ARS), centromere, artificial chromosome, chromosome or can be replicated or replicated in vitro or host cell, cell, in certain embodiments nucleus or cytoplasm other nucleic acid.In some embodiments, nucleic acid can be from single chromosome or its fragment (for example, nucleic acid sample can be from a chromosome of sample obtained from diploid organism).In certain embodiments, nucleic acid includes nucleosome, fragment or part of nucleosome or nucleosome-like structure.Nucleic acid sometimes comprises protein (for example, histone, DNA binding protein etc.).Nucleic acid analyzed by methods described herein is sometimes basically separated and not associated with protein or other molecules.Nucleic acid also includes derivatives, variants and analogs ("sense" or "antisense", "positive" strand or "negative" strand, "forward" reading frame or "reverse" reading frame) and double-stranded polynucleotides of RNA or DNA synthesized, replicated or amplified from single strand.Deoxyribonucleotides include deoxyadenosine, deoxycytidine, deoxyguanosine and deoxythymidine. For RNA, the base cytosine is replaced with uracil and the sugar 2' position includes a hydroxyl moiety. The nucleic acid can be prepared using a nucleic acid obtained from a subject as a template.

[0038] As used herein, the term "termination position" or "end position" (or just "end") refers to the genomic coordinates or genomic identity or nucleotide identity of the outermost base, for example, at the end of a cell-free DNA molecule, such as a plasma DNA molecule. The end position can correspond to either end of a DNA molecule. In this way, if the start and end of a DNA molecule are referred to, both can correspond to the end position. In some cases, an end position is the genomic coordinates or nucleotide identity of the outermost base on one end of a cell-free DNA molecule, which is detected or determined by an analytical method, such as massively parallel sequencing or next-generation sequencing, single-molecule sequencing, double-stranded or single-stranded DNA sequencing library preparation protocols, polymerase chain reaction (PCR) or microarray. In some cases, such in vitro techniques can change the true in vivo physical end of the cell-free DNA molecule. Therefore, each detectable end can represent a biologically true end or the end is one or more nucleotides inward or extended from the original end of the molecule, such as the 5' blunt end and 3' filling of the overhang of a non-blunt-ended double-stranded DNA molecule by Klenow fragment. The genomic identity or genomic coordinates of the end position can be derived from the alignment of the sequence reads with a human reference genome, such as hg19. It can come from a directory of indices or codes representing the original coordinates of the human genome. It can refer to a position or nucleotide identity on a cell-free DNA molecule, through but not limited to target-specific probes, mini-sequencing, DNA amplification reads. The term "genomic position" can refer to a nucleotide position in a polynucleotide (e.g., a gene, a plasmid, a nucleic acid fragment, a viral DNA fragment). The term "genomic position" is not limited to a nucleotide position within a genome (e.g., in a haploid chromosome set of a gamete or microorganism, or in each cell in a multicellular organism).

[0039] As used herein, the term "fragment" (e.g., DNA fragment) refers to a portion of a polynucleotide or polypeptide sequence comprising at least three consecutive nucleotides. Nucleic acid fragments can retain the biological activity and / or some characteristics of the parent polynucleotide. In one example, nasopharyngeal carcinoma cells can deposit fragments of Epstein-Barr virus (EBV) DNA into a subject, such as a patient's bloodstream. These fragments can include one or more BamHI-W sequence fragments, which can be used to detect the level of tumor-derived DNA in plasma. The BamHI-W sequence fragment corresponds to a sequence that can be recognized and / or digested using the Bam-HI restriction enzyme. The BamHI-W sequence can refer to a sequence 5'-GGATCC-3'.

[0040] As used herein, "local maximum" may refer to a genomic position (e.g., nucleotide) at which the maximum of a parameter of interest is obtained when compared to an adjacent position, or to the value of a parameter of interest at such a genomic position. For example, the range of an adjacent position may be from 50bp to 2000bp. Examples of parameters of interest include, but are not limited to, the number of fragments that terminate at a genomic position, the number of fragments that overlap with the position, or the ratio of fragments that cover a genomic position greater than a threshold size. When a parameter of interest has a periodic structure, many local maxima may occur. A global maximum is a specific one of a local maximum. Similarly, a "local minimum" refers to a genomic position at which the minimum of a parameter of interest is obtained when compared to an adjacent position, or to the value of a parameter of interest at such a genomic position.

[0041] As used herein, the term "ratio" of a nucleic acid molecule (e.g., DNA or RNA) terminated at a position may relate to the frequency at which a nucleic acid molecule terminates at that position. The ratio may be based on the number of nucleic acid molecules terminated at that position normalized relative to the number of nucleic acid molecules analyzed. The ratio may be based on the number of nucleic acid molecules terminated at that position normalized relative to the number of nucleic acid molecules terminated at different positions. The ratio may be based on the number of nucleic acid molecules from a first sample terminated at that position normalized relative to the number of nucleic acid molecules from a second sample (e.g., a reference sample). The ratio may be based on the number of nucleic acid molecules from a first sample terminated at a first group of positions (e.g., genomic positions) normalized relative to the number of nucleic acid molecules from a second sample (e.g., a reference sample) terminated at a second group of positions. Therefore, the ratio may correspond to the frequency at which how many nucleic acid molecules terminate at a certain position, and in some cases has nothing to do with the periodicity of the position with a local maximum in the number of nucleic acid molecules terminated at that position.

[0042] As used herein, the term "relative abundance" may refer to the ratio of the first amount of a nucleic acid fragment with a specific characteristic (e.g., terminating at one or more specified coordinates / termination positions or a specified length aligned with a specific region in a genome) to the second amount of a nucleic acid fragment with a specific characteristic (e.g., terminating at one or more specified coordinates / termination positions or a specified length aligned with a specific region in a genome). In an example, relative abundance may refer to the ratio of the number of DNA fragments terminating at a first group of genomic positions to the number of DNA fragments terminating at a second group of genomic positions. In some aspects, "relative abundance" may be a separation value that associates the amount (a value) of a DNA molecule terminated in a window of a genomic position with the amount (other values) of a DNA molecule terminated in another window of a genomic position. The two windows may overlap, but the sizes may be different. In other implementations, the two windows may not overlap. In addition, a window may be the width of a nucleotide, which is therefore equivalent to a genomic position.

[0043] As used herein, the term "reference genome" refers to any specific known, sequenced or characterized genome of any organism or virus, whether partial or complete, which can be used to reference the identified sequence from the subject. Exemplary reference genomes for human subjects and many other organisms are provided in online genome browsers hosted by the National Center for Biotechnology Information ("NCBI") or the University of California, Santa Cruz (UCSC). "Genome" refers to the complete genetic information of an organism or virus expressed in a nucleic acid sequence. As used herein, a reference sequence or reference genome is typically an assembled or partially assembled genome sequence from an individual or multiple individuals. In some embodiments, a reference genome is an assembled or partially assembled genome sequence from one or more human individuals. A reference genome can be considered as a representative example of a species gene set. In some embodiments, a reference genome comprises a sequence assigned to a chromosome. Exemplary human reference genomes include, but are not limited to, NCBI build 34 (UCSC equivalent: hg16), NCBI build 35 (UCSC equivalent: hg17), NCBI build 36.1 (UCSC equivalent: hg18), GRCh37 (UCSC equivalent: hg19), and GRCh38 (UCSC equivalent: hg38).

[0044] As used herein, the term "sequence read" or "read" refers to a nucleotide sequence generated by any sequencing process described herein or known in the art. Reads can be generated from one end of a nucleic acid fragment (a "single-end read"), and sometimes reads are generated from both ends of a nucleic acid fragment (e.g., a double-ended read, a double-ended read, or a paired read). The length of a sequence read is generally associated with a particular sequencing technology. In some embodiments, the mean, median, or average length of the sequence reads is between about 5,000 bp and 50,000 bp long (e.g., about 5,000 bp, about 7,500 bp, about 10,000 bp, about 12,500 bp, about 15,000 bp, about 20,000 bp, about 25,000 bp, about 30,000 bp, about 35,000 bp, about 40,000 bp, about 45,000 bp, about 50,000 bp, about 55,000, about 60,000, about 65,000, about 70,000, about 75,000, or about 80,000). In some embodiments, the sequence reads have a mean, median, or average length of about 1000 bp, 2000 bp, 5000 bp, 10,000 bp, 50,000 bp, or more. A sequence read (or sequencing read) can refer to sequence information (e.g., a string of nucleotides) corresponding to a nucleic acid molecule. For example, a sequence read can correspond to a string of nucleotides (e.g., up to about 50,000 bp) of a portion of a nucleic acid fragment, can correspond to a string of nucleotides at one or both ends of a nucleic acid fragment, or can correspond to nucleotides of the entire nucleic acid fragment. Sequence reads can be obtained in a variety of ways, for example, using sequencing techniques or using probes, for example, in hybridization arrays or capture probes, or amplification techniques, such as polymerase chain reaction (PCR) or linear amplification using a single primer or isothermal amplification.

[0045] As used herein, the terms "sequencing", "sequence determination", etc. used herein generally refer to any and all biochemical methods that can be used to determine the order of biological macromolecules such as nucleic acids or proteins. For example, sequencing data can include all or part of the nucleotide bases in a nucleic acid molecule such as a DNA fragment.

[0046] As used herein, the term "single nucleotide variant", "single nucleotide polymorphism" or "SNP" refers to a nucleotide replacement with a different nucleotide at the position (e.g., site) of a nucleotide sequence, such as an individual's sequence read. The substitution from the first nucleobase X to the second nucleobase Y can be expressed as "X>Y". For example, cytosine to thymine SNP can be expressed as "C>T". The term "het-SNP" refers to a heterozygous SNP, in which the genome is at least diploid and at least one of two or more homologous sequences, but not all, exhibits a specific SNP. Similarly, "hom-SNP" is a homologous SNP, in which each homologous sequence of a polyploid genome has the same variant compared to a reference genome. As used herein, the term "structural variant" or "SV" refers to a large (e.g., greater than 1kb) region of a genome that has undergone physical transformation, such as inversion, insertion, deletion or duplication (e.g., see Spielmann et al., 2018, Nat Rev Genetics 19: A review of human genome SVs in 453–467).

[0047] As used herein, the terms "folded reads", "folded sequence reads" or "folded regions" refer to multiple sequence reads derived from similar regions of a polynucleotide molecule and combined into a resulting sequence. For example, as part of sequence determination, multiple reads from a polynucleotide can be compared to each other (e.g., to provide a consensus sequence of the original polynucleotide). In addition, the resulting folded reads can include information about nucleotide mismatches, insertions, or deletions. The use of folded reads in nucleotide sequencing provides a higher degree of certainty in sequencing results (e.g., rather than relying on a single sequence read for each region of the polynucleotide being sequenced).

[0048] As used herein, the term "node" refers to a sequence read, and the term "branch" refers to the overlap between two separate sequence reads (e.g., nodes). A directed graph includes a plurality of nodes and a plurality of branches. Folded reads are aligned to the directed graph such that any folded read can be represented in sequence by a subset of edges and corresponding vertices (e.g., nodes). As described in Compeau et al., branches

[0049] (e.g., directed edges) can indicate the directionality of overlap between nodes, thereby helping to determine the overall topological organization of multiple nodes (Nat Biotechnol 29(11):987-991,2011). This is the basis for forming a directed graph. De Bruijn graphs are usually derived from sequence reads with predefined overlaps, and this method is very suitable for short read sequencing. This can be achieved by dividing the reads into specified lengths (e.g., k-polymers) or by limiting the length of the edge (Id.). Typically, for reads that are cut into k-polymers, the overlap will be k-1. The term "k-polymer" used herein refers to a set length of a DNA sequence. The unidirectional edges of a directed graph represent the sequence of k nucleotide bases (e.g., "k-polymers") in the target region, and these edges are connected by vertices (or nodes).

[0050] As used herein, the term "non-linear topological component" or "non-linear component" refers to a complex nucleic acid sequence. In particular, such components can be sequence reads comprising one or more variants, wherein different options for mapping sequences at one or more positions are indicated. Alternatively, a non-linear component can include circular nucleic acids. In some cases, these variants can be structural variations (SVs), single nucleotide variations (SNPs), small nucleotide variations, or regions with elevated numbers of het-SNPs that are simulated as larger substitution events.

[0051] As used herein, the term "haplogroup" refers to a haplotype contig. A contig contains a set of continuous overlapping nucleic acid segments representing a common region of a nucleic acid (e.g., a genome). And a haplotype includes several polymorphisms that are usually inherited together (e.g., those that are close to each other in a set of regions of a genome, such as chromosomes). Formally, a haplogroup is a "contig of clones with the same haplotype," as described in Makoff and Flomen, 2007, Genome Bio 8(6): R114.

[0052] As used herein, the term "phasing" refers to identifying alleles in a diploid or polyploid genome. In particular, phasing is used to determine which alleles belong together to a single chromosome or a specific region of the genome. Choi et al., 2018, PLOS Genetics 14(4): e1007308 describes a method for phasing the human genome, which is incorporated herein by reference in its entirety.

[0053] The terms used herein are for the purpose of describing a particular situation only, and the terms are not intended to be limiting. As used herein, the singular forms "a", "an", and "the" are intended to include the plural forms as well, unless the context clearly indicates otherwise. In addition, if the terms "include", "comprising", "having", "having", "with", or variations thereof are used in the detailed description and / or claims, such terms are intended to be inclusive in a manner similar to the term "comprising".

[0054] Several aspects are described below with reference to example applications for illustration. It should be understood that many specific details, relationships and methods are set forth to provide a comprehensive understanding of the features described herein. However, those of ordinary skill in the relevant art will readily recognize that the features described herein can be practiced without one or more specific details or by other methods. The features described herein are not limited by the order of the actions or events described, because some actions can occur in different orders and / or simultaneously with other actions or events. In addition, not all actions or events described need to be used to implement the method according to the features described herein.

[0055] Exemplary System Implementation

[0056] Now combine Figure 1 Details of an exemplary system are described. Figure 1 1 is a block diagram illustrating a system 100 according to some embodiments. In some embodiments, the device 100 includes one or more processing units (CPUs) 102 (also referred to as processors or processing cores), one or more network interfaces 104, a user interface 106, non-persistent storage 111, permanent storage 112, and one or more communication buses 114 for interconnecting these components. One or more communication buses 114 optionally include circuit systems (sometimes referred to as chipsets) that interconnect and control the communication between system components. Non-persistent storage 111 typically includes high-speed random access memory, such as DRAM, SRAM, DDR RAM, ROM, EEPROM, flash memory, while permanent storage 112 typically includes CD-ROM, digital versatile disk (DVD) or other optical storage device, tape cassette, magnetic tape, magnetic disk storage device or other magnetic storage device, magnetic disk storage device, optical disk storage device, flash memory device or other non-volatile solid-state storage device. Persistent storage 112 optionally includes one or more storage devices arranged away from CPU 102. Persistent storage 112 and one or more non-volatile storage devices within non-persistent storage 112 include non-temporary computer-readable storage media. In some embodiments, the non-persistent storage 111 or (alternatively) the non-transitory computer-readable storage medium (sometimes in combination with the persistent storage 112) stores the following programs, modules, and data structures, or a subset thereof:

[0057] An optional operating system 116, which contains routines for handling various basic system services and for performing hardware-related tasks;

[0058] Optional network communication module (or instructions) 118, which is used to connect the system 100 to other devices or to a communication network;

[0059] a graph module 120 for mapping long sequence reads to a directed graph, comprising a directed graph 122 of an organism having a heterozygous genome, the heterozygous genome comprising at least one non-linear topological component 124, each non-linear topological component comprising a start node 126 and an end node 128, and each pair of nodes being connected by at least a corresponding first branch 130-1 and a corresponding second branch 130-2; and

[0060] A mapping module 139 comprising a mapping of query sequences 140 to directed graphs 122, each mapping comprising one or more query sequences 142, each query sequence comprising at least a first alignment 144, a second alignment 146, and a path 148 of a topological relationship between the first and second alignments.

[0061] In some embodiments, directed graph 122 is derived from a de novo assembly of sequence reads or contigs. In some embodiments, directed graph 122 is derived from a reference genome (e.g., one or more existing annotations). In some embodiments, the sequence reads used to form the directed graph are derived from a biological sample of an organism, and relative to the target sequence set, the sequence reads have a total random error.

[0062] In some embodiments, the query sequence 142 comprises a sequence read in a plurality of sequence reads, or a scaffold (e.g., a contig, a unitig, or a haplogroup) formed by a plurality of sequence reads. In some embodiments, the query sequence 142 describes at least a start node or an end node of a corresponding nonlinear topological component 124.

[0063] In various embodiments, one or more of the above-mentioned elements are stored in one or more previously mentioned storage devices, and correspond to the instruction group for performing the above-mentioned functions. The above-mentioned modules, data or programs (such as instruction sets) do not need to be implemented as separate software programs, processes, data sets or modules, and therefore each subset of these modules and data can be combined or otherwise rearranged in various embodiments. In some embodiments, non-permanent memory 111 optionally stores a subset of the above-mentioned modules and data structures. In addition, in some embodiments, the memory stores other modules and data structures not described above. In some embodiments, one or more of the above-mentioned elements are stored in a computer system outside the computer system of the visualization system 100, and the computer system can be addressed by the visualization system 100 so that the visualization system 100 can retrieve all or part of such data when needed.

[0064] although Figure 1 A "system 100" is depicted, but the diagram is intended more as a functional description of various features that may be present in a computer system than as a schematic diagram of the architecture of the embodiments described herein. In practice, and as recognized by one of ordinary skill in the art, items shown separately may be combined and some items may be separated. Furthermore, although Figure 1 Certain data and modules are depicted in non-persistent storage 111 , but some or all of these data and modules may reside in persistent storage 112 .

[0065] Although referenced Figure 1 A system according to the present disclosure is disclosed, but now reference is made to Figure 2A and 2B The method according to the present disclosure is described in detail.

[0066] Use graphs when mapping long sequence reads or contigs to target sequences.

[0067] Box 202. A method for sequencing nucleic acids by graph-based mapping is provided.

[0068] In some embodiments, the target sequence will be the sequence to which the query is mapped.The target sequence can include any type of sequence, such as a reference genome, or an assembled contig, unit group, or scaffold.

[0069] Box 204. Using computer system 100, obtain a directed graph of a single diploid or polyploid organism of a species having a heterozygous genome, wherein the directed graph represents all or part of the heterozygous genome and includes one or more nonlinear topological components (e.g., bubbles or rings). In some embodiments, each corresponding nonlinear topological component in the one or more nonlinear topological components represents a variation relative to a set of target sequences (e.g., a structural variation, a single nucleotide variation (SNP), a small nucleotide variation, or a region that is heavy in a het-SNP and simulated as a larger substitution event).

[0070] In some embodiments, the directed graph represents at least ten percent of the heterozygous genome, and the directed graph includes ten or more non-linear topological components.

[0071] In some embodiments, the directed graph represents at least fifteen percent of the heterozygous genome, and the directed graph includes twenty-five or more non-linear topological components.

[0072] In some embodiments, the directed graph represents at least five percent of the heterozygous genome, at least ten percent of the heterozygous genome, at least twenty percent of the heterozygous genome, at least thirty percent of the heterozygous genome, at least forty percent of the heterozygous genome, at least fifty percent of the heterozygous genome, at least seventy-five percent of the heterozygous genome, at least ninety percent of the heterozygous genome, or at least one hundred percent of the heterozygous genome.

[0073] In some embodiments, directed graph includes multiple nodes, for example, each node can represent sequence reads.Therefore, the total number of nodes in directed graph will be directly proportional to the size of genome or the genome part mapped.For example, in some embodiments where complex diploid and higher order polyploid genomes of large proportion are mapped, directed graph will have at least 10,000 nodes.In some embodiments, directed graph will have at least 100,000 nodes.In some embodiments, directed graph will have at least 500,000 nodes.In other embodiments, directed graph will have at least 1 million, 2 million, 3 million, 4 million, 5 million or more nodes.

[0074] In some embodiments, the directed graph includes two or more nonlinear topological components, five or more nonlinear topological components, ten or more nonlinear topological components, twenty-five or more nonlinear topological components, fifty or more nonlinear topological components, or one hundred or more nonlinear topological components.

[0075] In some embodiments, for example, in the case of mapping a large proportion of complex diploid or higher order polyploid genomes, the directed graph will have at least 1,000 nonlinear topological components. In some embodiments, the directed graph will have at least 10,000 nonlinear topological components. In some embodiments, the directed graph will have at least 100,000 nonlinear topological components. In other embodiments, the directed graph will have at least 500,000, 1 million, 2 million, 3 million, 4 million, 5 million or more nonlinear topological components. In some embodiments, for each heterozygous allele in a subject with at least one sequence read mapped to each possible allele, there is a nonlinear element in the graph. That is, it is observed that at least one sequence read in each SNP and each INDEL has been mapped to allele 1, and at least one read has been mapped to allele 2.

[0076] In some embodiments, the directed graph includes any combination of the above percentages of heterozygous genomes, numbers of nodes, and numbers of non-linear topological components.

[0077] In some embodiments, each corresponding nonlinear topological component of the one or more nonlinear topological components has a corresponding starting node and a corresponding ending node, and the starting node and the ending node are connected by at least a corresponding first branch and a corresponding second branch. In such an embodiment, the corresponding first branch and the corresponding second branch correspond to (e.g., map to) a target sequence. In some embodiments, the group of target sequences represents one or more of the following: (i) a reference sequence of a species or (ii) an assembled contig, unit group, haplogroup or scaffold. Assembly can be performed from scratch or by mapping multiple sequencing reads to another target sequence.

[0078] In some embodiments, the corresponding first and second branches of the nonlinear topological components in one or more nonlinear topological components represent the first allele and the second allele of the variation in the heterozygous genome, respectively. In some embodiments, there may be more branches and / or more alleles. For example, in some embodiments, the genome of the organism is polyploid rather than just diploid, and there may be as many alleles as chromosome pairs. In some embodiments, one or more nonlinear topological components include three corresponding branches and three corresponding alleles, four corresponding branches and four corresponding alleles, five corresponding branches and five corresponding alleles, or six corresponding branches and six corresponding alleles.

[0079] A directed graph is formed using a plurality of sequence reads from a biological sample of a single organism by overlapping corresponding sequence reads in a plurality of sequence reads, the sequence reads having a random error greater than ten percent relative to a target sequence, wherein the amount of overhang between the overlapping first sequence read and the second sequence is unrestricted, provided that there is a minimal common region between the first and second sequence reads (e.g., thereby distinguishing the methods described herein from a de Bruijn graph). In some embodiments, the plurality of sequence reads have an average sequence length greater than 10,000 base pairs.

[0080] In some embodiments, the plurality of sequence reads have an average length of greater than 2,500 base pairs, greater than 5,000 base pairs, greater than 10,000 base pairs, greater than 15,000 base pairs, greater than 20,000 base pairs, greater than 50,000 base pairs, or greater than 100,000 base pairs.

[0081] Block 218. Obtain a query sequence from the biological sample in electronic format using computer system 100. The query sequence is (i) a sequence read from the plurality of sequence reads, or (ii) a scaffold (e.g., a contig, a haplogroup, or a haplogroup) formed by the plurality of sequence reads. Here, the query sequence comprises at least a corresponding start node or a corresponding end node of a first nonlinear topological component of the one or more nonlinear topological components.

[0082] In some embodiments, the query sequence is a sequence read in the plurality of sequence reads. In some embodiments, the query sequence is a scaffold formed by the common content of two or more sequence reads in the plurality of sequence reads.

[0083] In some embodiments, the query sequence is obtained by single molecule real-time (SMRT) sequencing, which generally provides long sequence reads. In some embodiments, SMRT sequencing is performed as described by Rhoads and Au (Genomics Proteomics Bioinformatics 13: 278-289, 2015). This long read sequencing provides enough information for de novo whole genome sequencing alone (e.g., as described in Chin et al., 2013, Nat Methods 10 (6): 563-569).

[0084] Box 224. The method continues by forming a mapping of the query sequence to the directed graph using the directed graph. The mapping of the query sequence includes (i) a first alignment component to a first portion of the directed graph, (ii) a second alignment component to a second portion of the directed graph, and (iii) a path in the directed graph that describes a relationship between the first portion and the second portion of the directed graph. In this case, one of the first portion and the second portion is in a first or second branch of a first nonlinear topological component, and the other of the first portion and the second portion is not in a branch of the first nonlinear topological component. In some embodiments, each nonlinear topological component has more than two branches, and the same query can also span multiple nonlinear topological components (e.g., mapped to a large overlapping group of a reference graph).

[0085] In some embodiments, the nonlinear topological component includes at least 2 branches, at least 3 branches, at least 4 branches, at least 5 branches, at least 6 branches, at least 7 branches, at least 8 branches, at least 9 branches, or at least 10 branches.

[0086] In some embodiments, mapping provides a starting position and an ending position in the target sequence for the query sequence. In some embodiments, one or more programs individually or collectively further include instructions for using a sequence alignment algorithm and the starting position and the ending position to compare the query sequence with the target sequence.

[0087] In some embodiments, the sequence alignment algorithm is a dynamic programming global alignment algorithm, a dynamic programming semi-global alignment algorithm, or a dynamic programming local alignment algorithm.

[0088] Different alignment algorithms have different efficiencies or are adjusted to better detect specific genomic features. In fact, there are many different types of algorithms, as discussed below. GASSST (Global Alignment Short Sequence Search Tool) aims to improve the detection of insertion and deletion (e.g., as described in Rizk and Lavenier, Bioinformatics 26 (20): 2534–2540, 2010). The Smith-Waterman alignment (BWA-SW) method of Burrows-Wheeler Aligner aims to improve efficiency using long read sequences (e.g., as described in Li and Durbin, Bioinformatics 26 (5): 589–595, 2010). Parasail aims to provide faster local pairwise alignment within the sequence (e.g., as described in Daily, BMC Bioinformatics 17: 81, 2016). Skewer is particularly focused on using next-generation sequencing paired end reads (e.g., as described in Jiang et al., 2014, BMC Bioinformatics 15: 182). DNACLUST aims to improve the clustering efficiency of sequences in highly similar DNA regions (e.g., as described in Ghodsi et al., 2011, BMC Bioinformatics 12:271). Local alignment algorithms have a long history, including universal methods such as BLAST and Smith-Waterman (Altschul et al., 1990, J Mol Biol 215:403-410; Smith and Waterman, 1981, J Mol Biol 147 (1):195-197). However, there are also recently developed methods that attempt to improve the initial algorithm. For example, Bowtie 2 improves the performance of gap alignment (e.g., as described in Langmead and Salzberg, Nat Methods 9 (4):357-359,2012). Li and Homer provided a comparison of different sequence alignment algorithms and, in particular, concluded that long-read sequencing can obtain more detailed genomic information than short-read-based sequencing (Briefings in Bioinformatics 11(5):473-483, 2010).

[0089] A directed graph is formed when long sequence reads or overlapping groups are mapped to a target sequence. In some embodiments of the present disclosure, a nucleic acid sequencing method for constructing a directed graph of a single diploid or polyploid organism of a species with a heterozygous genome is provided. The method is performed by mapping a plurality of sequence reads obtained from a biological sample of an organism to the target sequence of the species. A plurality of sequence reads have a random error greater than ten percent relative to the target sequence. A directed graph represents all or part of a heterozygous genome and includes one or more nonlinear topological components (e.g., bubbles or rings). Each corresponding nonlinear topological component in one or more nonlinear topological components represents a variation relative to a reference sequence (wherein the variation is, for example, a structural variation, a single nucleotide variation (SNP), a small nucleotide variation, or one of the concentrated regions that are heavy and simulated as larger substitution events in het-SNPs). Each corresponding nonlinear topological component in one or more nonlinear topological components has a corresponding starting node and a corresponding terminating node, and the starting node and the terminating node are at least connected by a corresponding first branch and a corresponding second branch. One of the corresponding first branch and the corresponding second branch corresponds to the target sequence, and the other of the corresponding first branch and the corresponding second branch corresponds to a variation identified by mapping in one or more sequence reads in the plurality of sequence reads.

[0090] In some embodiments, mapping provides a starting position and an ending position in a target sequence for a sequence read in a plurality of sequence reads. In such embodiments, the sequence reads are aligned to the target sequence using a sequence alignment algorithm and the starting position and the ending position of the sequence reads.

[0091] In some embodiments, the directed graph represents at least ten percent of the species genome, at least twenty percent of the species genome, at least thirty percent of the species genome, at least forty percent of the species genome, at least fifty percent of the species genome, at least sixty percent of the species genome, or at least seventy percent of the species genome.

[0092] In some embodiments, directed graph includes multiple nodes, for example, each node can represent sequence reads.Therefore, the total number of nodes in directed graph will be directly proportional to the size of genome or the genome part mapped.For example, in some embodiments where complex diploid and higher order polyploid genomes of large proportion are mapped, directed graph will have at least 10,000 nodes.In some embodiments, directed graph will have at least 100,000 nodes.In some embodiments, directed graph will have at least 500,000 nodes.In other embodiments, directed graph will have at least 1 million, 2 million, 3 million, 4 million, 5 million or more nodes.

[0093] In some embodiments, for example, in the case of mapping a large proportion of complex diploid or higher order polyploid genomes, the directed graph will have at least 1,000 nonlinear topological components. In some embodiments, the directed graph will have at least 10,000 nonlinear topological components. In some embodiments, the directed graph will have at least 100,000 nonlinear topological components. In other embodiments, the directed graph will have at least 500,000, 1 million, 2 million, 3 million, 4 million, 5 million or more nonlinear topological components. In some embodiments, for each heterozygous allele in a subject with at least one sequence read mapped to each possible allele, there is a nonlinear element in the graph. That is, it is observed that at least one sequence read in each SNP and each INDEL has been mapped to allele 1, and at least one read has been mapped to allele 2.

[0094] Use case 1: Graph formation from phased de novo assembly of long sequencing reads.

[0095] Typical mapping methods are often unable to handle single haplogroups. Using conventional mappers and only the primary assembly as reference, haplogroups are broken down into primary and one or more supplementary alignments in regions of structural variation - map bubbles. For example, some mappers (e.g., BLASR with the "--bestn 1" option selected) will report only the primary alignment (see Chaisson and Tesler BMC Bioinformatics 13:238, 2012). However, this obscures some of the underlying genetic variation.

[0096] A method for nucleic acid sequencing is provided, comprising constructing a haplogroup assembly map for a single diploid or polyploid organism of a species having a heterozygous genome. In some embodiments, the haplogroup assembly map represents the heterozygous genome and comprises one or more nonlinear topological components (e.g., bubbles or rings). Each corresponding nonlinear topological component of the one or more nonlinear topological components has a corresponding starting node and a corresponding ending node, and the starting node and the ending node are connected by at least a corresponding first branch and a corresponding second branch. The corresponding first and second branches represent first and second phased haplogroups corresponding to the first and second homologous chromosomes in the heterozygous genome, respectively.

[0097] In some embodiments, there may be more branches and / or more alleles. For example, in some embodiments, the genome of the organism is polyploid rather than just diploid, and there may be as many alleles as chromosome pairs. In some embodiments, one or more nonlinear topological components include three corresponding branches and three corresponding alleles, four corresponding branches and four corresponding alleles, five corresponding branches and five corresponding alleles, or six corresponding branches and six corresponding alleles.

[0098] The sequencing method continues by obtaining haplogroups derived from a plurality of sequence reads from a biological sample obtained from a single diploid or polyploid organism in an electronic format. The haplogroup comprises at least a corresponding start node or a corresponding end node of a first nonlinear topological component in one or more nonlinear topological components. The method continues by forming a map of the haplogroup using a haplogroup assembly graph. The map of the haplogroup comprises (i) a first alignment component to a first portion of the haplogroup assembly graph, (ii) a second alignment component to a second portion of the haplogroup assembly graph, and (iii) a path in the haplogroup assembly graph describing a relationship between the first portion and the second portion of the haplogroup assembly graph. One of the first portion and the second portion is in a first or second branch of the first nonlinear topological component, and the other of the first portion and the second portion is not in a branch of the first nonlinear topological component.

[0099] In some embodiments, the mapping provides a starting position and an ending position in the target sequence for the haplogroup. The one or more programs, individually or collectively, further include instructions for aligning the haplogroup with the target sequence of the species using a sequence alignment algorithm and the starting position and the ending position.

[0100] In some embodiments, the haplogroup assembly graph represents at least ten percent of the heterozygous genome, and the haplogroup assembly graph comprises ten or more non-linear topological components. In some embodiments, the haplogroup assembly graph represents at least fifteen percent of the heterozygous genome, and the haplogroup assembly graph comprises twenty-five or more non-linear topological components.

[0101] In some embodiments, the haplogroup assembly map represents at least five percent of the heterozygous genome, at least ten percent of the heterozygous genome, at least twenty percent of the heterozygous genome, at least thirty percent of the heterozygous genome, at least forty percent of the heterozygous genome, at least fifty percent of the heterozygous genome, at least seventy-five percent of the heterozygous genome, at least ninety percent of the heterozygous genome, or at least one hundred percent of the heterozygous genome.

[0102] In some embodiments, the haplogroup assembly graph includes two or more nonlinear topological components, five or more nonlinear topological components, ten or more nonlinear topological components, twenty-five or more nonlinear topological components, fifty or more nonlinear topological components, or one hundred or more nonlinear topological components.

[0103] In some embodiments, the haplogroup assembly graph includes any combination of the above percentages of heterozygous genomes and numbers of non-linear topological components.

[0104] In some embodiments, a haplogroup assembly graph is formed using the plurality of sequence reads by overlapping corresponding sequence reads in the plurality of sequence reads, wherein the amount of overhang between the overlapping first sequence read and the second sequence is unrestricted, provided that there is a minimum common region between the first and second sequence reads (e.g., to distinguish a de Bruijn graph with a predefined overlap length).

[0105] In some embodiments, the minimum consensus region comprises at least 1 base, at least 2 bases, at least 3 bases, at least 4 bases, at least 5 bases, at least 6 bases, at least 7 bases, at least 8 bases, at least 9 bases, at least 10 bases, at least 15 bases, at least 20 bases, at least 30 bases, at least 40 bases, at least 50 bases, at least 60 bases, at least 70 bases, at least 80 bases, at least 90 bases, at least 100 bases, at least 500 bases or at least 1000 bases. In some embodiments, each consensus region has the same overlapping length. In some embodiments, each consensus region has different overlapping lengths. In some embodiments, the first group of consensus regions has the first overlapping length and the second group of consensus regions has the second overlapping length.

[0106] Figure 3 An example of aligning the output of phased haplogroups in an assembled genome containing structural variation between haplotypes is illustrated. Haplogroup phasing is suitable for the FALCON-unzip assembly method, which is intended as a diploid assembler (e.g., as described in Chin et al., 2016, Nat Methods 13(12): 1050-1054). In some embodiments, the disclosed methods are applied herein to enhance FALCON-unzip assembly by enabling accurate haplogroup placement and construction of less fragmented haplogroup graphs (as well as longer phased haplogroups). In some embodiments, the resulting assembly is output in a graphical format to retain as much information as possible about divergent haplotypes (e.g., as opposed to the FASTA output format that only provides linear contigs).

[0107] exist Figure 3 , the main contig 302 is shown as having two related contigs 304 and 306, generating respective bubbles 310 and 312. In some embodiments, the sequence reads that are best mapped within a single contig produce a continuous alignment. In some embodiments, reads that map across one or more branch positions are reported as split alignments with accompanying alignment paths. In some embodiments, a score is assigned to each alignment. In some embodiments, each score indicates the quality of the corresponding alignment. In some embodiments, only alignments with scores above a threshold are reported. In some embodiments, if more than one alignment has a related score above a threshold, each such alignment is reported.

[0108] In some embodiments, quality scores are reported in the SAM format (eg, as described in the 'Sequence Alignment / Map Format Specification', available on the Internet at samtools.github.io / hts-specs / SAMv1.pdf.

[0109] In some embodiments, the threshold score is at least 1, at least 5, at least 10, at least 15, at least 20, at least 25, at least 30, at least 35, at least 40, at least 45, at least 50, or at least 55.

[0110] In some embodiments, the quality score is calculated using the sum of the quality scores of each base across all alignments in a given sequence read and the sum of the quality scores of each mismatched base (e.g., as described in Li et al., 2008, Genome Res 18(11):1851-1858, which is incorporated herein by reference in its entirety). In some embodiments, the quality score is a Phred score, which is a log-scaled probability that a base call is incorrect (e.g., as described in Ruffalo et al., 2012, Bioinformatics 28(5):i349-i355).

[0111] This approach provides increased coverage near the ends of contigs (enabling better refinement), improved read phasing, and increased ability to resolve complex repetitive structures in the assembly.

[0112] Use case 2: Phasing reads by mapping to haplogroup maps.

[0113] A: Phasing of long reads using haplogroup maps.

[0114] Resolving haplotype scaffolds by combining the FALCON-unzip haplogroup assembly and long-range sequencing data (e.g., 10x, Dovetail, Bionano) is a common procedure used by groups such as the Vertebrate Genome Project (VGP) consortium. Mapping long-range reads should be more accurate when the context of the assembly graph is available. Unzip assemblies are highly contiguous, but often discard connectivity information when aligning long-range data. Often, short paired-end reads map only to contigs / haplogroups. In the case of slightly repetitive regions (which are abundant on larger genomes), mapping can become ambiguous. Graph context can restrict mapping results to a region where there is a path (of acceptable length) between candidate target sequences. A related approach uses Hi-C data in combination with an assembly graph in GFA format to accurately build assemblies (see Ghurye et al. bioRxiv 261149; on the internet at doi.org / 10.1101 / 261149).

[0115] Another nucleic acid sequencing method is provided, which includes constructing a haplogroup assembly map for a single diploid or polyploid organism of a species having a heterozygous genome. The haplogroup assembly map represents the heterozygous genome and includes one or more nonlinear topological components (e.g., bubbles or rings). Each corresponding nonlinear topological component of the one or more nonlinear topological components has a corresponding starting node and a corresponding ending node, and the starting node and the ending node are connected by at least a corresponding first branch and a corresponding second branch. The corresponding first and second branches represent first and second phased haplogroups corresponding to the first and second homologous chromosomes in the heterozygous genome, respectively.

[0116] In some embodiments, for example in a polyploid genome, there may be 2 or more branches. In some embodiments, there may be 3 or more branches, 4 or more branches, 5 or more branches, 6 or more branches, 7 or more branches, 8 or more branches, 9 or more branches, or 10 or more branches.

[0117] The method continues by obtaining a plurality of sequence reads from a biological sample obtained from a single diploid or polyploid organism in an electronic format. A first sequence read in the plurality of sequence reads comprises at least a corresponding start node or a corresponding end node of a first nonlinear topological component in one or more nonlinear topological components. In a subsequent embodiment, the method may also include the start node and the end node and the entire bubble (or ring or even multiple nearby bubbles). The method continues using the haplogroup assembly graph to form a mapping of each corresponding sequence read in the plurality of sequence reads. The mapping of the first sequence read is reported as a split alignment, which includes (i) a first alignment component to a first portion of the haplogroup assembly graph, (ii) a second alignment component to a second portion of the haplogroup assembly graph, and (iii) a path describing a relationship between the first portion and the second portion of the haplogroup assembly graph in the haplogroup assembly graph, wherein one of the first portion and the second portion is in a first or second branch of the first nonlinear topological component, and the other of the first portion and the second portion is not in a branch of the first nonlinear topological component.

[0118] In some embodiments, the mapping provides a starting position and an ending position of a corresponding sequence read in a target sequence, and wherein one or more programs individually or collectively further include instructions for aligning the corresponding sequence read to the target sequence using a sequence alignment algorithm and the starting position and the ending position.

[0119] In some embodiments, a haplogroup assembly graph is formed using a plurality of sequence reads by overlapping corresponding sequence reads in the plurality of sequence reads. The amount of overhang between the overlapping first sequence read and the second sequence is not restricted, provided that there is a minimum common region between the first and second sequence reads. This is different from a de Bruijn graph, which has a specific k-mer length overlap for each sequence read. The details of the overlap length are discussed above and will not be repeated here.

[0120] B: Phasing of paired-end short sequence reads using haplogroup maps.

[0121] One aspect of the present disclosure provides a method for nucleic acid sequencing, which includes constructing a haploid assembly map for a single diploid or polyploid organism having a heterozygous genome. The haploid assembly map represents the heterozygous genome and includes one or more nonlinear topological components (bubbles or rings). Each corresponding nonlinear topological component in the one or more nonlinear topological components has a corresponding starting node and a corresponding ending node, and the starting node and the ending node are at least connected by a corresponding first branch and a corresponding second branch. The corresponding first and second branches represent the first and second phased haplogroups corresponding to the first and second homologous chromosomes in the heterozygous genome, respectively.

[0122] The method continues by obtaining a plurality of paired-end sequence reads from a biological sample obtained from a single diploid or polyploid organism in an electronic format. A first group of paired-end sequence reads in the plurality of paired-end sequence reads at least comprises a corresponding start node or a corresponding end node of a first nonlinear topological component in one or more nonlinear topological components (e.g., bubbles or rings). The method continues by using a haplogroup assembly graph to form a mapping of each corresponding paired-end sequence read in the plurality of paired-end sequence reads. The mapping of the first group of paired-end sequence reads is reported as a split alignment, which includes (i) a first alignment component to a first portion of a haplogroup assembly graph representing a first sequence read in the first group of paired-end sequence reads, (ii) a second alignment component to a second portion of a haplogroup assembly graph representing a second sequence read in the first group of paired-end sequence reads, and (iii) a path describing a relationship between the first portion and the second portion of the haplogroup assembly graph in the haplogroup assembly graph. One of the first portion and the second portion is in a first or second branch of the first nonlinear topological component, and the other of the first portion and the second portion is not in a branch of the first nonlinear topological component.

[0123] In some embodiments, mapping provides a starting position and an ending position for each corresponding group of double-ended sequence reads in the target sequence. One or more programs, individually or collectively, further include instructions for aligning each corresponding group of double-ended sequence reads with the target sequence using a sequence alignment algorithm and the starting position and the ending position.

[0124] In some embodiments, a haplogroup assembly graph is formed using the plurality of paired-end sequence reads by overlapping corresponding paired-end sequence read groups in the plurality of paired-end sequence reads. The amount of overhang between the overlapping first set of paired-end sequence read groups and the second set of paired-end sequence read groups is not limited, provided that there is a minimum common region between the first set of paired-end sequence read groups and the second set of paired-end sequence read groups.

[0125] In some embodiments, the minimum consensus region comprises at least 1 base, at least 2 bases, at least 3 bases, at least 4 bases, at least 5 bases, at least 6 bases, at least 7 bases, at least 8 bases, at least 9 bases, at least 10 bases, at least 15 bases, at least 20 bases, at least 30 bases, at least 40 bases, at least 50 bases, at least 60 bases, at least 70 bases, at least 80 bases, at least 90 bases, at least 100 bases, at least 500 bases or at least 1000 bases. In some embodiments, each consensus region has the same overlapping length. In some embodiments, each consensus region has different overlapping lengths. In some embodiments, the first group of consensus regions has the first overlapping length and the second group of consensus regions has the second overlapping length.

[0126] like Figure 4A ,4B As shown in 4C, the haplogroup assembly map may include one or more bubbles representing phased haplogroups (e.g., 420, 422, and 432). The regions between each phased haplogroup (e.g., 406, 412, and 418) are generally folded sequences that do not contain enough information to be successfully phased.

[0127] After unzipping the assembly (e.g., using FALCON-unzip), a haplogroup map can be used as an intermediate output. In some embodiments, the haplogroup map is a directed acyclic graph (DAG) consisting of bubbles with at least two branches (one for each haplotype) and collapsed regions (e.g., regions of the genome where there is insufficient information for phasing). Each bubble branch in the map (e.g., 420, 422) represents a phased portion of a diploid genome. Reads can be mapped to the FALCON-unzip haplogroup map and phased by selecting the branch with the highest alignment score in the map as the most likely destination (e.g., phase).

[0128] Figure 4B and 4C Two possible graphs of the same assembly are illustrated, along with potential mappings of two long-range paired-end reads 402 and 404. In some embodiments, the orientations of the two paired-end reads 402 and 404 are known. In such embodiments, it can be determined Figure 4B , where the location of each bubble 420 and 422 is known with respect to each other. In some embodiments, the first mapping location may be preferred due to the connectivity information provided by the known locations of 402 and 404. In some embodiments, the connectivity between 402 and 404 is unknown, so only each individual bubble 420 and 422 is determined. In some embodiments, read segment 404 is repeated and is also mapped to a second map location (e.g., in Figure 4C 430 in the area).

[0129] In some embodiments, this method is used to improve results for resolving haplotype scaffolds, provide higher quality assemblies, and improve FALCON decompression results.

[0130] Use case 3: Mapping to a circular genome.

[0131] Circular genomes can be simulated in a graph context by simply adding an edge connecting the end of a chromosome to the beginning of the same chromosome. In this case, a graph-aware mapper will know that the alignment can safely extend beyond the end of the sequence and continue from its own beginning. Specifying a circular context via a graph definition allows for mixed references of linear and circular chromosomes in input FASTA sequence sets (e.g., eukaryotic genomes).

[0132] Therefore, one aspect of the present disclosure provides a method for nucleic acid sequencing, in which a graph of a species with a circular genome is obtained. The graph represents a circular genome and includes a segment line of a sequence specifying the circular genome. The beginning of the sequence defines the beginning of the segment line, the end of the sequence defines the end of the segment line, and a sideline representing a circular boundary in the graph, the sideline specifying the end of the segment line is connected to the beginning of the segment line. The method continues to obtain multiple sequence reads from a biological sample obtained from an organism of a species in an electronic format. The first sequence read in the multiple sequencing reads spans the circular boundary. The method continues to use the graph to form a mapping of each sequencing read in the multiple sequencing reads. The mapping of the first sequencing read is reported as a split alignment, including (i) a first alignment component to the end portion of the sequence, (ii) a second alignment component to the initial portion of the sequence, and (iii) a path describing the relationship between the end portion of the sequence in the graph and the initial portion of the sequence in the graph.

[0133] Another aspect of the present disclosure provides a method for nucleic acid sequencing, in which a graph of a species having a circular genome is obtained. The graph represents a circular genome and includes a segment line specifying a sequence of the circular genome. The beginning of the sequence defines the beginning of the segment line, and the end of the sequence defines the end of the segment line. The edge represents a circular boundary in the graph, and the end of the edge specifies that the segment line is linked to the beginning of the segment line. The method continues to obtain multiple double-end sequencing reads from a biological sample obtained from an organism of the species in an electronic format. The first group of double-end sequencing reads in the multiple double-end sequencing reads includes a circular boundary. The method continues to use the graph to form a mapping of each double-end sequencing read in the multiple double-end sequence reads. The mapping report of the first group of double-end sequence reads is a split alignment, including (i) an alignment of the first sequence read in the first group of double-end sequence reads with the end portion of the sequence, (ii) an alignment of the second sequence read in the first group of double-end sequence reads with the initial portion of the sequence, and (iii) a path describing the relationship between the end portion and the initial portion in the graph.

[0134] Figure 5 A circular reference genome is illustrated. A circular reference is the simplest form of a graph, with only one edge 504 connecting the end of a sequence 508 to its beginning 506. A graph-based mapper can detect reads that should be aligned across gaps and build the correct alignment. This ensures that high coverage of mapped reads is retained at the front and back of such contigs, enabling better refinement of circular contigs. In some embodiments, the FALCON assembler includes a feature to naturally output circular edges for circular contigs. This means that if a graph-based mapper is used with such contigs, the circular edges of the contigs will be automatically output. Figure 1 If used together, Arrow or other refinement tools will have nearly uniform available read coverage at their disposal (this is currently a problem when using BLASR as the default aligner).

[0135] In some embodiments, the reads across the circular boundary will report two mapping positions. In some embodiments, the first mapping position occurs near the beginning 506 of chromosome 502, until reaching the end 508. In some embodiments, the second mapping position occurs near the end 508 of chromosome, and ideally arrives at the beginning 506. In some additional embodiments, a mapping path is provided to describe the relationship between the first and second mappings. This method ensures good coverage of the mapping reads in the entire circular genome, and allows better refinement of the sequence.

[0136] Use case 4: Mapping using genomic inversions.

[0137] Inversions occur as natural structural variations (SVs) between different genomes. Many current mappers perform poorly for sequences containing inversions. For example, Chen et al. describe structural variation predictions for various mapping tools, and each method analyzed has a relatively low success rate for inversions (Nat Methods 6(9):677-681, 2009). The detection of inversions emerges seamlessly from the definition of the graph-based mapping described herein. In some embodiments of the present disclosure, a portion of the sequence mapped to the reverse complementary strand of the reference can be connected to the remainder of the query fragment by an implicit edge.

[0138] Therefore, one aspect of the present disclosure provides a method for nucleic acid sequencing, wherein a query sequence representing a portion of a genomic sequence of an organism of a species containing a genomic inversion is obtained in an electronic format. The query sequence includes a first query portion and a second query portion. In the method, the query sequence is mapped to a target sequence of the species. Mapping includes mapping the first query portion to a first target portion of a first chain of the target sequence and mapping the second query portion to a second target portion of a second chain of the target sequence, wherein the first chain is an inversion complement of the second chain, thereby forming a split alignment, which includes (i) a first alignment component, which includes an alignment of the first query portion with the first target portion of the first chain, (ii) a second alignment component, which includes an alignment of the second query portion with the second target portion of the second chain, and (iii) an alignment path describing the relationship between the first target portion and the second target portion.

[0139] In some embodiments, the query sequence includes a third query portion. In some embodiments, the query sequence includes a fourth query portion, a fifth query portion, a sixth query portion, a seventh query portion, an eighth query portion, a ninth query portion, or a tenth query portion. In some embodiments, there can be multiple inversions in the same query. In some embodiments, a query can include 1 or more inversions, 2 or more inversions, 3 or more inversions, 4 or more inversions, 5 or more inversions, 10 or more inversions, 15 or more inversions, or 20 or more inversions.

[0140] The mapping includes mapping the third query portion to the third reference portion of the first strand, and the split alignment further includes a third alignment component comprising an alignment of the third query portion to the third reference portion of the first strand, and the alignment path further describes the relationship of the first reference portion and the second reference portion to the third reference portion.

[0141] In some embodiments, the query sequence is a single sequence read obtained from a biological sample of an organism. In some embodiments, the query sequence is a contig formed from multiple sequence reads obtained from a biological sample of an organism.

[0142] Fig. 6A and 6B An example of inversion mapping is shown. The query sequence 600 has one or more anchors of a reference genome. Fig. 6A As shown, in some embodiments, query sequence 600 includes three anchor sequences 602 , 604 , and 606 .

[0143] In some embodiments, the query sequence comprises at least 2 anchor sequences, at least 3 anchor sequences, at least 4 anchor sequences, at least 5 anchor sequences, at least 6 anchor sequences, at least 7 anchor sequences, at least 8 anchor sequences, at least 9 anchor sequences, or at least 10 anchor sequences.

[0144] like Figure 6B As shown, two of the anchor sequences 602 and 606 are mapped to the forward strand 610 of the reference sequence, and the third anchor sequence 604 is mapped to the corresponding reverse strand 608 of the reference sequence. Unlike described herein for other use cases, the map of this embodiment cannot be predetermined by the reference sequence or other input. Instead, in this embodiment, the map is implicitly determined. In some embodiments, the map determined by this method is subsequently analyzed for structural variants.

[0145] Use Case 5: SMRTbell TM The reference simulation is shown in Fig.

[0146] In some embodiments, a SMRTbell template is used. As described by Travers et al., in some embodiments, a SMRTbell template is composed of a sequence of interest, an adapter sequence, and an optional barcode sequence (Nucleic Acids Res 38 (15): e159, 2010). In some embodiments, a known region of the genome is targeted for sequencing. In such an embodiment, a closely related reference sequence is used to model the insert and construct a circular graph. In some embodiments, the SMRTbell template is used to improve sequencing efficiency because the adapter and / or barcode sequence contribute to the alignment of the expanded long read sequence (e.g., a long read sequence that passes through the circular template multiple times). In some embodiments, the long read sequence passes through the circular template at least 1 time, at least 2 times, at least 5 times, at least 10 times, at least 20 times, or at least 50 times.

[0147] In one aspect according to use case 5, a nucleic acid sequencing method is provided, comprising obtaining a circular graph. The circular graph comprises a first segment specifying an insert sequence in a first direction, a second segment specifying an adapter sequence in the first direction, a third segment specifying an insert sequence in a second direction, a fourth segment specifying an adapter sequence in the second direction, a first edge connecting an end of the first segment to a beginning of the second segment, a second edge connecting an end of the second segment to a beginning of the third segment, a third edge connecting an end of the third segment to a beginning of the fourth segment, and a fourth edge connecting an end of the fourth segment to a beginning of the first segment. The second direction is the reverse complement of the first direction.

[0148] The method proceeds by obtaining a query sequence read from a biological sample in an electronic format. The sequence read spans a sequence represented by a circular graph. The method then proceeds to form a map of the sequence read using the graph. The map of the sequence read is reported as an alignment of each instance of the inserted sequence in the sequence read.

[0149] In some embodiments, the mapping also includes an alignment of each instance of the adaptor sequence.

[0150] In some embodiments, the sequence reads span the ring graph multiple times. In some embodiments, multiple times is three or more times. In some embodiments, multiple times is five or more times. In some embodiments, multiple times is twenty or more times.

[0151] SMRTbell TM Constructs, such as Fig. 7AAs shown, it is composed of the required insertion sequence (e.g., 702 and 710), adapters (e.g., 706 and 714) and optional barcode sequences (e.g., 704, 708, 712, and 716). In some embodiments, the known region of the genome is targeted for sequencing, and a closely related reference sequence is used to model the insert and build a circular graph from the circular construct. In some embodiments, the barcode sequences are the reverse complements of each other (e.g., 712 / 716 and 704 / 708 pairs). In some embodiments, the barcode sequence is symmetrical (e.g., 716 will be the reverse complement of 708 or 712). In some embodiments, the barcode sequence is not symmetrical (e.g., 716 is not the reverse complement of 708 or 712). In some embodiments, the required insert will exist in forward (e.g., 702) and reverse (e.g., 710) sequences. In some embodiments, one or more regions of the circular construct are double-stranded (e.g., 719). In some embodiments, the adapter sequences (e.g., 714 and 706) are the reverse complements of each other. In some embodiments, the adapter sequences (e.g., 714 and 706) are not the reverse complements of each other. In some embodiments, each adapter sequence forms a hairpin region (e.g., 718 and 717) at each end of the circular construct.

[0152] In some embodiments of the present disclosure, the expansion and alignment of ZMW reads along a known amplicon reference are modeled in a graph context. In the simplest case (e.g., without the optional barcode sequences 704, 708, 712, and 716), a graph based on a circular construct will include a first segment (e.g., 714), a second segment (e.g., 702), a third segment (e.g., 706), a fourth segment (e.g., 710), and an edge between each two of these segments. In some embodiments, two of the segments will be the reverse complement of each other (e.g., 702 and 710).

[0153] In some embodiments, from SMRTbell TM The sequencing reads of the construct include ZMW reads 720, such as Figure 7B As shown. Figure 7B As can be seen in FIG. 1 , in some embodiments—depending on how long the ZMW read is—multiple copies of the insert sequence (e.g., 722-A, 722-B, and 722-C) will exist. In some embodiments, each copy of the insert sequence is referred to as a sub-read, and multiple sub-reads are compared to each other, thereby improving the accuracy of the overall sequencing results.

[0154] Use case 6: Transcriptome mapping by edges defining splice sites.

[0155] Transcriptome mapping is a natural application of graph-based mappers. Therefore, in some embodiments, the disclosed methods are directed to circRNA mapping. circRNA is a circular molecule. In some embodiments, sequencing such circular molecules produces exons that appear to be disordered but are actually circular arrangements of spliced ​​exons. In some embodiments, adding circular edges to the transcriptome enables mapping of such RNA-seq data. Li et al. describe an example of circRNA mapping, but the method relies on k-mer analysis and is not suitable for long sequence reads used in the present disclosure (BMC Genomics 19(Suppl 6):572,2018).

[0156] Therefore, according to the method of use case 6, a nucleic acid sequencing method is provided, the method comprising obtaining a transcriptome annotation that defines the mapping of each corresponding exon in a plurality of exons in a reference sequence of a species, and constructing a graph corresponding to the transcriptome annotation. The graph comprises the reference sequence of the species, the first fragment, the second fragment, and a link between the 3' end of the first exon and the 5' end of the second exon.

[0157] In some embodiments, the first fragment represents the first exon in a plurality of exons. In some embodiments, an exon is a sequence represented as a node in the figure. In some embodiments, an edge represents a splice site (e.g., those positions where two fragments are virtually connected together). In some embodiments, the first fragment also includes the corresponding coordinates of the 5' end and 3' end of the first exon in the reference sequence. In some embodiments, the second fragment represents the second exon in a plurality of exons and includes the corresponding coordinates of the 5' end and 3' end of the second exon in the reference sequence.

[0158] Fig. 8A and 8B An example of transcriptome mapping is illustrated. In some embodiments, the transcriptome is described by a graph marking a splice region along a reference genome (e.g., 810). A first fragment (e.g., a single node or a group of nodes) represents a first exon 802 in a plurality of exons (e.g., 802, 804, 806, and 808). In some embodiments, the exons as sequences are represented as nodes or node sets in the graph. Edges represent splice sites where two fragments should be virtually connected together. The first fragment also includes the corresponding coordinates of the 5' end and 3' end of the first exon in the reference sequence. The second fragment represents a second exon in a plurality of exons and includes the corresponding coordinates of the 5' end and 3' end of the second exon in the reference sequence.

[0159] The method proceeds by obtaining sequence reads of mRNA from a biological sample of an organism of the species in an electronic format. The sequence reads include a first exon and a second exon (e.g., 802 and 804). In some embodiments, the sequence reads include at least 2 exons, at least 3 exons, at least 4 exons, at least 5 exons, at least 6 exons, at least 7 exons, at least 8 exons, at least 9 exons, or at least 10 exons. The method then proceeds to use the graph to form a mapping of the sequence reads to a reference sequence 810.

[0160] In some embodiments, the linking allows mapping even when the order of the exons in the sequence reads differs from the linear order in which they appear in the reference genome. In some embodiments, a third exon that is not in the sequence reads is between the first exon and the second exon in the reference sequence.

[0161] In some embodiments, the query sequence is obtained from a circular RNA molecule (circRNA), and the order of exons in the query sequence is not linear with their corresponding order in the reference sequence.

[0162] Transcriptome annotations can be represented in the form of a graph, connecting adjacent splice sites. By applying a graph-based mapper, there is no need to explicitly construct transcriptome sequences before mapping Isoseq / RNA-seq reads.

[0163] In some embodiments, the expected result of transcriptome mapping is a splice alignment with a predefined number of CIGAR (Concise Idiosyncratic Gapped Alignment Report) operations in the same alignment, or an alignment path along a single exon region of a genome. In some embodiments, the set of possible CIGAR operations includes at least an alignment match (e.g., a sequence match or mismatch), an insertion to a reference, a deletion from a reference, a region skipped from a reference, a sequence match, or a sequence mismatch.

[0164] Conclusion

[0165] Multiple instances can be provided for the components, operations or structures described herein as single instances. Finally, the boundaries between the various components, operations and data storages are arbitrary to some extent, and specific operations are described in the context of specific illustrative configurations. Other functional allocations are contemplated, and can fall within the scope of the (multiple) embodiments. In general, the structure and function presented as a separate component in the example configuration can be implemented as a combined structure or assembly. Similarly, the structure and function presented as a single component can be implemented as a separate component. These other variations, modifications, additions and improvements fall within the scope of the (one or more) embodiments.

[0166] It will also be understood that although the terms first, second, etc. may be used herein to describe various elements, these elements should not be limited by these terms. These terms are only used to distinguish one element from another element. For example, without departing from the scope of the present disclosure, the first subject may be referred to as the second subject, and similarly, the second subject may be referred to as the first subject. Although both the first subject and the second subject are subjects, these subjects are not the same subject.

[0167] The terms used in this disclosure are only used to describe specific embodiments and are not intended to limit the present invention. As used in the specification of the present invention and the appended claims, the singular forms "one", "an" and "said" are also intended to include plural forms, unless the context clearly indicates otherwise. It should also be understood that the term "and / or" used herein refers to and encompasses any and all possible combinations of one or more of the associated listed items. It should be further understood that when used in this specification, the term "comprises and / or comprising" specifies the presence of stated features, integers, steps, operations, elements and / or components, but does not exclude the presence or addition of one or more other features, integers, steps, operations, elements, components and / or groups thereof.

[0168] As used herein, the term "if" may be interpreted to mean "when" or "upon" or "in response to determining" or "in response to detecting," depending on the context. Similarly, the phrases "if it is determined" or "if [a stated condition or event] is detected" may be interpreted to mean "upon determining" or "in response to determining" or "upon detecting (a stated condition or event)" or "in response to detecting (a stated condition or event)," depending on the context.

[0169] The foregoing description includes example systems, methods, techniques, instruction sequences, and computing machine program products that embody illustrative embodiments. For purposes of explanation, many specific details are set forth in order to provide an understanding of the various embodiments of the subject matter of the present invention. However, it will be apparent to those skilled in the art that the subject matter of the present invention can be practiced without these specific details. In general, well-known illustrative examples, protocols, structures, and techniques are not shown in detail.

[0170] For the purpose of explanation, the foregoing description has been described with reference to specific embodiments. However, the above illustrative discussion is not intended to be exhaustive or to limit the embodiments to the precise forms disclosed. In view of the above teachings, many modifications and variations are possible. These embodiments are selected and described in order to best explain these principles and their practical applications, thereby enabling other persons skilled in the art to best use these embodiments and various embodiments with a variety of modifications suitable for the intended specific use.

Claims

1. A method for assembling long sequence reads to determine a nucleic acid sequence, comprising: In a computer system comprising one or more processors and memory, the memory stores one or more programs that individually or collectively include instructions for: obtaining a directed graph of a single diploid or polyploid organism of a species having a heterozygous genome, wherein the directed graph represents all or a portion of the heterozygous genome and comprises one or more non-linear topological components, wherein each corresponding non-linear topological component in the one or more non-linear topological components represents a variation relative to a set of target sequences, Each corresponding nonlinear topological component in the one or more nonlinear topological components has a corresponding starting node and a corresponding ending node, wherein the starting node and the ending node are connected by at least a corresponding first branch and a corresponding second branch, wherein one of the corresponding first branch and the corresponding second branch corresponds to a target sequence in the target sequence group, and forming the directed graph using a plurality of sequence reads from a biological sample of a single organism by overlapping corresponding sequence reads of a plurality of sequence reads having a random error greater than ten percent relative to a target sequence in a target sequence set, wherein an amount of overhang between an overlapping first sequence read and a second sequence is not restricted, provided that a minimum common region exists between the first sequence read and the second sequence read; obtaining a query sequence from the biological sample in an electronic format, wherein the query sequence is (i) a sequence read from the plurality of sequence reads, or (ii) a scaffold formed by the plurality of sequence reads, wherein the query sequence comprises at least the corresponding start node or the corresponding end node of a first nonlinear topological component from among the one or more nonlinear topological components; forming a mapping of the query sequence to the directed graph using the directed graph, wherein the mapping of the query sequence comprises (i) a first alignment component to a first portion of the directed graph, (ii) a second alignment component to a second portion of the directed graph, and (iii) a path describing a relationship between a first portion and a second portion of the directed graph in the directed graph, wherein one of the first portion and the second portion is in a first or second branch of a first nonlinear topological component and the other of the first portion and the second portion is not in a branch of the first nonlinear topological component; as well as The mapping is used to determine the nucleic acid sequence of the query sequence.

2. The nucleic acid sequencing method according to claim 1, wherein the corresponding first branch and second branch of the nonlinear topological component in the one or more nonlinear topological components respectively represent the first allele and the second allele of the variation in the heterozygous genome.

3. The nucleic acid sequencing method according to claim 1 or 2, wherein the query sequence is a sequence read among the plurality of sequence reads. 4 . The nucleic acid sequencing method according to claim 1 , wherein the query sequence is a scaffold formed by a common content of two or more sequence reads in the plurality of sequence reads. The nucleic acid sequencing method according to claim 1 , wherein the plurality of sequence reads have an average sequence length greater than 10,000 base pairs.

6. A nucleic acid sequencing method according to claim 1, wherein the mapping provides a starting position and an ending position in the target sequence for the query sequence, and wherein the one or more programs individually or collectively further include instructions for aligning the query sequence with the target sequence using a sequence alignment algorithm and the starting position and the ending position.

7. The nucleic acid sequencing method according to claim 1, wherein the sequence alignment algorithm is a dynamic programming global alignment algorithm, a dynamic programming semi-global alignment algorithm or a dynamic programming local alignment algorithm.

8. The nucleic acid sequencing method according to claim 1, wherein The directed graph represents at least ten percent of the heterozygous genome, and The directed graph includes ten or more non-linear topological components.

9. The nucleic acid sequencing method according to claim 1, wherein The directed graph represents at least fifty percent of the heterozygous genome, and The directed graph includes twenty-five or more non-linear topological components.

10. The nucleic acid sequencing method according to claim 1, wherein the directed graph comprises at least 500,000 nodes.

11. The nucleic acid sequencing method according to claim 1, wherein the directed graph comprises at least 100,000 non-linear topological components.

12. An executable software product stored on a computer readable medium comprising program instructions for assembling long sequence reads to determine a nucleic acid sequence, the program instructions being executed on at least one processor, comprising: obtaining a directed graph of a single diploid or polyploid organism of a species having a heterozygous genome, wherein the directed graph represents all or a portion of the heterozygous genome and comprises one or more non-linear topological components, wherein each corresponding non-linear topological component in the one or more non-linear topological components represents a variation relative to a set of target sequences, Each corresponding nonlinear topological component in the one or more nonlinear topological components has a corresponding starting node and a corresponding ending node, wherein the starting node and the ending node are connected by at least a corresponding first branch and a corresponding second branch, wherein one of the corresponding first branch and the corresponding second branch corresponds to a target sequence in the target sequence group, and forming the directed graph using a plurality of sequence reads from a biological sample of a single organism by overlapping corresponding sequence reads of a plurality of sequence reads having a random error greater than ten percent relative to a target sequence in a target sequence set, wherein an amount of overhang between an overlapping first sequence read and a second sequence is not restricted, provided that a minimum common region exists between the first sequence read and the second sequence read; obtaining a query sequence from the biological sample in an electronic format, wherein the query sequence is (i) a sequence read from the plurality of sequence reads, or (ii) a scaffold formed by the plurality of sequence reads, wherein the query sequence comprises at least the corresponding start node or the corresponding end node of a first nonlinear topological component from among the one or more nonlinear topological components; forming a mapping of the query sequence to the directed graph using the directed graph, wherein the mapping of the query sequence includes (i) a first alignment component to a first portion of the directed graph, (ii) a second alignment component to a second portion of the directed graph, and (iii) a path in the directed graph describing a relationship between the first portion and a second portion of the directed graph, wherein one of the first portion and the second portion is in a first or second branch of a first non-linear topological component, and the other of the first portion and the second portion is not in a branch of the first non-linear topological component; as well as The mapping is used to determine the nucleic acid sequence of the query sequence.

13. A system for assembling long sequence reads to determine a nucleic acid sequence, comprising: Memory; as well as a processor coupled to the memory, configured to: obtaining a directed graph of a single diploid or polyploid organism of a species having a heterozygous genome, wherein the directed graph represents all or a portion of the heterozygous genome and comprises one or more non-linear topological components, wherein each corresponding non-linear topological component in the one or more non-linear topological components represents a variation relative to a set of target sequences, Each corresponding nonlinear topological component in the one or more nonlinear topological components has a corresponding starting node and a corresponding ending node, wherein the starting node and the ending node are connected by at least a corresponding first branch and a corresponding second branch, wherein one of the corresponding first branch and the corresponding second branch corresponds to a target sequence in the target sequence group, and forming the directed graph using a plurality of sequence reads from a biological sample of a single organism by overlapping corresponding sequence reads of a plurality of sequence reads having a random error greater than ten percent relative to a target sequence in a target sequence set, wherein an amount of overhang between an overlapping first sequence read and a second sequence is not restricted, provided that a minimum common region exists between the first sequence read and the second sequence read; obtaining a query sequence from the biological sample in an electronic format, wherein the query sequence is (i) a sequence read from the plurality of sequence reads, or (ii) a scaffold formed by the plurality of sequence reads, wherein the query sequence comprises at least the corresponding start node or the corresponding end node of a first nonlinear topological component from among the one or more nonlinear topological components; forming a mapping of the query sequence to the directed graph using the directed graph, and wherein the mapping of the query sequence comprises (i) a first alignment component to a first portion of the directed graph, (ii) a second alignment component to a second portion of the directed graph, and (iii) a path describing a relationship between the first portion and a second portion of the directed graph in the directed graph, wherein one of the first portion and the second portion is in a first or second branch of a first non-linear topological component, and the other of the first portion and the second portion is not in a branch of the first non-linear topological component; as well as The mapping is used to determine the nucleic acid sequence of the query sequence.

Images