Improving Alignment Using Homopolymer Collapsed Sequencing Reads

Abstract

The present disclosure provides methods, compositions, computer-implemented processes, etc. for parsing long and highly similar but not identical genomic regions to improve assembly quality, particularly for polyploid genomes. Aspects of the present disclosure relate to using exact string matching of homopolymer-folded sequence reads to determine whether two sequences overlap and thus represent the same genomic region (e.g., the same haplotype in a polyploid genome) or whether the sequences represent different genomic regions.

Description

[0001] Cross - reference to related applications

[0002] This application claims the benefit of priority of U.S. Provisional Patent Application No. 62 / 812,191, filed on February 28, 2019, the disclosure of which is hereby incorporated by reference in its entirety for all purposes. Background of the Invention

[0003] Genome sequence assembly refers to determining the nucleotide sequence of each genomic chromosome by a process of breaking each chromosome into smaller genomic fragments, "reading" the nucleotide sequence of each genomic fragment to make the fragment sequence into a read sequence, and then assembling the read sequences. Assembly requires multiple copies of genomic DNA. These multiple copies can be obtained from multiple cells from the same organism (assuming the same genomic DNA), or by replicating (e.g., PCR amplification) the genome contained in a single cell. When the same genomic locus is covered by two different fragments, these two fragments are called "overlapping". The nucleotide sequences of the overlapping fragments also overlap because they share a common subsequence. If the common subsequence shared by the overlapping fragments occurs uniquely in the genome, the overlap between these fragments can be detected from the reads of these fragments. In this case, if two reads also share a common nucleotide sequence that extends to one end of each read, it can be correctly inferred that these two reads are from a pair of overlapping genomic fragments. Two reads can be "overlapped" by superimposing the common sequence. A graph structure can be formed, where vertices (reads) are connected by edges between "overlapping" reads. Each edge represents the assertion that two reads are from genomic fragments that contain the same genomic locus. In a valid assembly, each connected component represents overlapping genomic fragments from the same chromosome. A contig can be formed from each connected component by aligning the reads and superimposing the positions corresponding to the same positions in the genome in the reads. In the absence of read errors, the nucleotide identity of each position can be correctly determined. Considering read errors, the "stacking" of many overlapping reads at each genomic position allows the draft assembly to be polished to a highly consistent accuracy using redundancy to suppress read errors.

[0004] Although the assembly process is conceptually simple, correct overlap detection across an entire genome has proven difficult, especially for genomes containing regions of long repeats. Assembly is fundamentally limited by the accuracy of detecting overlapping genomic fragments from its reads. False positive overlap errors occur when reads from two different loci are incorrectly identified as coming from the same locus. False positives can occur when two different loci have long regions of identical or nearly identical sequence. False negative overlap errors occur when reads from overlapping genomic fragments are incorrectly identified as coming from different loci. False negatives can occur when read errors obscure the common nucleotide sequence shared by overlapping genomic fragments. If not subsequently corrected, both types of overlap errors result in assembly errors. False positive overlaps can lead to chromosome fusions, or more commonly, expansion or folding of repetitive elements. False negative errors, especially systematic ones, can lead to assembly breaks where a single chromosome is represented by multiple disjoint contigs, which may be accompanied by the loss of some loci at the contig boundaries.

[0005] The present disclosure particularly addresses the challenges posed to genome assembly by the presence of highly similar but non-identical sequences in haploid and polyploid genomes. SUMMARY OF THE INVENTION

[0006] The present disclosure provides methods, compositions, computer-implemented processes, etc. for resolving long and highly similar but non-identical genomic regions to improve assembly quality, especially for polyploid genomes. At a basic level, this involves determining whether two sequences overlap, i.e., whether the sequences represent the same genomic region - the same haplotype of that region in a polyploid genome - or whether the sequences represent different genomic regions - or different haplotypes.

[0007] Aspects of the present disclosure include methods for assembling a genome or genomic region, the method comprising: obtaining a plurality of sequence reads of genomic fragments from a genome of interest; generating a homopolymer collapsed sequence (HCS) and a corresponding homopolymer encoded sequence (HES) for each of the plurality of sequence reads; generating a suffix / prefix exact string match of the HCS reads, wherein the length of the exact string match is equal to or greater than a minimum length; generating trimmed HCS reads by removing any nucleotides in each of the plurality of HCS reads that are not part of the suffix / prefix exact string match with another HCS read; generating a first directed overlap graph from the trimmed HCS reads; identifying connected components in the second directed overlap graph; generating a multiple sequence alignment for each of the connected components, wherein positions in each trimmed HCS read are labeled with consecutive integer values such that the same integer value is assigned to aligned positions in any two trimmed HCS reads; pruning merged nodes from the second directed overlap graph based on the multiple sequence alignment; generating a homopolymer collapsed consensus sequence by joining basecalls at each aligned position in the multiple sequence alignment of the trimmed HCS reads; associating a vector of homopolymer lengths with each position in the homopolymer collapsed consensus sequence, wherein: (i) the number of elements in the vector is the number of trimmed HCS reads that cover that position in the multiple sequence alignment, and (ii) each component of the vector is the length of the homopolymer in the corresponding HES at that position; assigning a consensus homopolymer length as the floor of the median of the components of the vector of homopolymer lengths associated with that position; and replacing each position in the homopolymer collapsed consensus sequence with a homopolymer string formed by N consecutive nucleotide copies at that position, wherein N is the specified consensus homopolymer length calculated for that position, to generate a homopolymer extended consensus sequence, thereby assembling the genome or genomic region of interest.

[0008] In certain embodiments, prior to generating the HCS reads, the method further comprises generating the reverse complement sequence of each of the plurality of sequence reads.

[0009] In certain embodiments, the minimum length of the overlapping region is from 0.5 kb to 10 kb. In certain embodiments, the minimum length of the overlapping region is from 5 kb to 8 kb. In certain embodiments, the minimum length of the overlapping region is from 6 kb to 7 kb. In certain embodiments, the minimum length is at least half of the average length of the HCS reads.

[0010] In certain embodiments, the plurality of sequence reads are generated in a single molecule synthesis sequencing reaction. In certain embodiments, the single molecule synthesis sequencing reaction is a single molecule real-time sequencing reaction. In certain embodiments, the plurality of sequence reads are generated in a single molecule nanopore sequencing reaction.

[0011] In some embodiments, the plurality of sequence reads are a plurality of single molecule consensus sequences (SMCS). In some embodiments, the SMCS are generated from at least 8 subreads. In some embodiments, the subreads are generated from a tandem polynucleotide substrate in a single molecule sequencing reaction. In some embodiments, the subreads are generated in a single molecule synthesis sequencing reaction. In some embodiments, the subreads are generated in a single molecule nanopore-based sequencing reaction. In some embodiments, the subreads are generated from a circular or topologically circular polynucleotide substrate in a single molecule synthesis sequencing reaction.

[0012] In some embodiments, the genome of interest is the human genome.

[0013] In some embodiments, when the genomic sample contains multiple different genomes, the method further includes generating assemblies for the multiple different genomes. In some embodiments, the sample is a metagenomic sample containing multiple microbial genomes.

[0014] In some embodiments, the HCS that is not placed into the joining component is placed into a holding bin for validating variant calls in the assembly.

[0015] In some embodiments, prior to generating the HCS, a plurality of sequence reads are preselected to map to one or more genomic regions of interest. In some embodiments, the preselection mapping is performed by a low-stringency sequence similarity search. In some embodiments, the one or more genomic regions of interest include first and second genomic loci that have high sequence similarity to each other. In some embodiments, separate consensus sequences are generated for the first and second genomic loci. In some embodiments, the one or more genomic regions of interest contain genomic loci with highly repetitive regions.

[0016] In some embodiments, the method is a method for de novo genome assembly. In some embodiments, the de novo genome assembly is a complete or partial haplotype-resolved assembly of a polyploid genome.

[0017] Aspects of the present disclosure include a system for determining a consensus sequence, comprising: a memory; an input / output; and a processor coupled to the memory, wherein the system is configured to: receive a plurality of sequence reads of genomic fragments from a genome of interest; generate a homopolymer collapse sequence (HCS) and a corresponding homopolymer encoding sequence (HES) for each of the plurality of sequence reads; generate a suffix / prefix exact string match of the HCS reads, wherein the length of the exact string match is equal to or greater than a minimum length; generate trimmed HCS reads by removing any nucleotides in each of the plurality of HCS reads that are not part of a suffix / prefix exact string match with another HCS read; generate a first directed overlap graph from the trimmed HCS reads; identify connected components in the second directed overlap graph; generate a multiple sequence alignment for each of the connected components, wherein positions in each trimmed HCS read are labeled with consecutive integer values such that the same integer value is assigned to aligned positions in any two trimmed HCS reads; prune merged nodes from the second directed overlap graph based on the multiple sequence alignment; generate a homopolymer collapse consensus sequence by joining base calls at each aligned position in the multiple sequence alignment of the trimmed HCS reads; associate a vector of homopolymer lengths with each position in the homopolymer collapse consensus sequence, wherein: (i) the number of elements in the vector is the number of trimmed HCS reads that cover that position in the multiple sequence alignment, and (ii) each component of the vector is the length of the homopolymer in the corresponding HES at that position; assign a consensus homopolymer length as the floor of the median of the components of the vector of homopolymer lengths associated with that position; and replace each position in the homopolymer collapse consensus sequence with a homopolymer string formed by N consecutive nucleotide copies, where N is the specified consensus homopolymer length calculated for that position, to generate a homopolymer extended consensus sequence; and provide the homopolymer extended consensus sequence to a user to assemble a genome of interest or a genomic region of a genome.

[0018] In certain embodiments, the system is further configured to perform the method according to any of the above embodiments and output the result of the method to a user. BRIEF DESCRIPTION OF THE DRAWINGS

[0019] Figure 1 Shows the process of generating SMCS reads from a polynucleotide substrate (a double-stranded polynucleotide with hairpin adapters at both ends).

[0020] Figure 2 Shows an example of two overlapping genomic fragments and two reads derived from these genomic fragments that share a common subsequence.

[0021] Figure 3Shown is an example of the alignment of two genomic fragments from different loci and two reads from these fragments that share a common subsequence.

[0022] Figure 4 Two reads originating from a genomic fragment containing tandem repeats and two alignments of these reads are shown.

[0023] Figure 5 Shown are the diploid genome, two genomic fragments from the maternal copy of chromosome 2, and the alignment of two reads from these fragments.

[0024] Figure 6 Two genomic fragments derived from the paternal and maternal copies of chromosome 2 and the alignment of two reads derived from these fragments are shown.

[0025] Figure 7 Two overlapping genomic fragments and two pairs of reads from these fragments are shown. The first pair has no errors, but the second read in the second pair contains a homopolymer deletion.

[0026] Figure 8 It illustrates the near orthogonality between signal—biological variation between two highly similar sequences, typically single nucleotide variation—and noise, read errors that confound the identification of overlapping genomic fragments, which are typically homopolymer insertions and deletions (indels).

[0027] Figure 9 Shown are two overlapping genomic fragments, two reads derived from these fragments, the second of which contains a homopolymer deletion, and an alignment of the homopolymer folded sequences derived from the reads.

[0028] Figure 10 An example showing how a read corrupted by homopolymers can be "perfected" by homopolymer folding. The homopolymer folded sequence of the read matches the homopolymer folded sequence of the genomic fragment from which the read was obtained, thereby masking indel errors in the read.

[0029] Figure 11 An example of filtering out homopolymer indel errors to identify a pair of overlapping reads and avoid false overlaps with reads from highly similar genomic fragments from different alleles is shown.

[0030] Figure 12 Schematic diagram showing multiple sequence alignment between exact string matches and “perfect” reads.

[0031] Figure 13 , Figure 14 and Figure 15Shows the algorithm workflow of using HCS to divide SMCS into haplotypes, calling consensus sequences for haplotypes, calling consensus sequence lengths for homopolymer regions in the consensus sequences to generate extended consensus sequences for homopolymers, and calling homozygous and heterozygous variants by comparison with a reference genome, where in some cases previously excluded HCS can be used for variant calling verification.

[0032] Figure 16 Shows how a homozygous region can induce the undesired merging of two different haplotypes into a single connected component, i.e., the haplotypes can be separated, but in the process, the haplotypes are decomposed into smaller haplotigs whose connectivity cannot be resolved in the absence of SMCS reads that fully span the homozygous region. The process of removing the merged node (i.e., node C) is sometimes referred to as "pruning" in this article.

[0033] Figure 17 Shows how SMCS reads across a homozygous region resolve haplotypes. This is also a pruning process. The process of removing the merged node (i.e., node C) is sometimes referred to as "pruning" in this article.

[0034] Figure 18 Shows a histogram of SMCS read lengths, the lengths of the HCSs derived from these reads, and the ratio of the length of each HCS to the SMCS read from which it was derived.

[0035] Figure 19 Shows the multiple sequence alignment of 11 homopolymer-folded SMCS reads from a single haplotype of SMN2.

[0036] Figure 20 Shows the multiple sequence alignment of 51 homopolymer-folded SMCS reads, which combines two haplotypes of SMN1.

[0037] Figure 21 Shows the diploid assembly of 100 SMCS reads mapped to the SMN1 and SMN2 sequences in the human genome reference GrCh38. Detailed implementation

[0038] The present disclosure particularly provides an improved method for parsing long and highly similar but distinct genomic sequences to improve the quality of genomic assembly, especially for polyploid genomes. Generally, the process includes filtering out the main forms of sequencing errors that confound genomic assembly and performing exact string matching of the filtered reads to prevent overlap of reads from highly similar genomic fragments from different loci or different haplotypes.

[0039] Definitions

[0040] The term "genomic fragment" is used herein to refer to single-stranded or double-stranded DNA molecules that are extracted from a cell and broken from the chromosome in which they reside, or alternatively, copies of such molecules formed by replication (e.g., PCR or linear amplification). Genomic fragments are identified by genomic loci - their original position in the chromosome, their nucleotide sequence, and their haplotype in a polyploid genome. When two genomic fragments share a common genomic locus and belong to the same haplotype in a polyploid genome, the two genomic fragments are "overlapping". The nucleotide sequences of overlapping genomic fragments are also overlapping; that is, the two nucleotide sequences share a common subsequence corresponding to the genomic locus shared by the overlapping genomic fragments. However, the converse is not true. Two genomic fragments whose sequences share a common subsequence are not necessarily "overlapping" because the common subsequence may occur at two different genomic loci, or in a polyploid genome, at the same locus but in different haplotypes. Genomic fragments can be derived from any source desired by the user (e.g., any animal, plant, fungus, single-celled organism, etc.). In some cases, a library of polynucleotide substrates can be derived from multiple different organisms, such as multiple different human samples or a metagenomic sample containing a mixture of different organisms. Genomic fragments can be the product of an amplification process (e.g., by PCR or linear amplification), native / non-amplified polynucleotides, or a combination of both (e.g., the polynucleotide substrate has a region of interest with amplified genomic fragments and non-amplified genomic fragments or has native strands and complementary strands produced by amplification). No limitation is intended in this regard.

[0041] The term "polynucleotide substrate" is used herein to refer to a polynucleotide that includes genomic fragments (or copies thereof) in a form that can be sequenced by a sequencing platform, regardless of the sequencing platform used. In certain embodiments, the polynucleotide substrate includes functional domains in addition to genomic fragments (e.g., synthetic or otherwise engineered sequences and / or functional portions) that facilitate obtaining and / or analyzing the sequence of the genomic fragments. Examples of such functional domains include, but are not limited to, one or more of the following: primer binding sites, binding sites for motor proteins (e.g., as employed in certain nanopore sequencing technologies), capture primer binding sites, capture moieties (e.g., cholesterol, biotin, avidin / streptavidin, etc.), sequencing primer binding sites, barcodes, registration sequences, unique molecular identifiers, detectable labels, or any other convenient sequence or portion. Such additional sequences and portions can be provided by ligating adapters to the genomic fragments, e.g., by ligation, amplification, etc., as is commonly done in the art. Libraries of polynucleotide substrates for genomic fragments of interest, e.g., whole genomes, are routinely generated and analyzed in the art.

[0042] The present disclosure uses the term "region of interest" to refer to a subset of the entire genome to which the disclosed methods can also be applied. For example, a "region of interest" can include one or more genes as a contiguous block or multiple blocks. No limitation is intended in this regard.

[0043] The present disclosure uses the term "single molecule consensus sequence" (SMCS) to refer to a consensus sequence obtained by analyzing multiple sequence reads of a genomic fragment. Each complete sequence read of a genomic fragment, which does not include any sequence of the flanking adapter polynucleotides, is referred to herein as a "sub-read". Due to differences in the construction of the polynucleotide substrate and / or the sequencing technology employed, a set of sub-reads of a region of interest can include (i) only a single strand of the polynucleotide or (ii) sub-reads of both complementary strands of the polynucleotide. For example, a polynucleotide substrate that requires sequence data may include multiple linear head-to-tail copies of a genomic fragment, which, when sequenced, provides a set of sub-reads, one for each copy, representing the same original genomic fragment (e.g., a tandem polynucleotide substrate generated by rolling circle amplification of a circular polynucleotide containing the genomic fragment). In contrast, when sequencing a double-stranded genomic fragment with hairpin adapters at both ends using a long-read synthetic sequencing method (e.g., used in sequencing the polynucleotide substrate is linear in structure but circular in topology), a set of sub-reads is generated that includes sub-reads of the forward strand of the double-stranded genomic fragment and its complementary reverse strand. The forward and reverse strand sub-reads can be analyzed to generate a consensus sequence of the genomic fragment. It should be noted that the potential sequencing method does not necessarily determine whether only single-stranded or complementary strand sub-reads are obtained. For example, rolling circle amplification of a polynucleotide can produce a linear polynucleotide substrate, which, when sequenced using nanopore sequencing technology, will return sub-reads of both complementary strands. In addition, a structurally circular double-stranded polynucleotide substrate containing a genomic fragment (similar in topology to a bacterial plasmid) sequenced using a synthetic sequencing method will return sub-reads of only one strand of the genomic fragment.

[0044] Figure 1 A schematic is provided of how to generate SMCS reads from a polynucleotide substrate in a sequencing reaction. At the Figure 1 top, a polynucleotide substrate with a double-stranded DNA genomic fragment and two terminal hairpin adapters is shown. Although only one polynucleotide substrate is shown, it should be clear that the library contains a population of Polynucleotide substrate. This polynucleotide substrate binds to a sequencing primer and a polymerase under certain conditions to form a ternary complex capable of synthesizing nucleic acids. The ternary complex undergoes sequencing in a synthesis sequencing reaction (Pacific Biosciences of California, Inc.), where the addition of each base is recorded in a single long sequencing read. Since the polynucleotide substrate is topologically circular, once the polymerase first traverses the entire polynucleotide substrate, it enters rolling circle amplification (RCA). The entire length of a single long sequencing read is referred to as a "polymerase read" and includes all sequence data from multiple passes of both the genomic fragment and the adapter. Each sub-read of the two strands of the genomic fragment in the polymerase read is identified by removing the adapter sequence. Figure 1 Each sub-read in is labeled in the order in which it is generated. (Note that sub-read 11 is still being generated). Given the Figure 1 topology of the polynucleotide substrate, odd-numbered sub-reads (i.e., sub-reads 1, 3, 5, 7, 9, and 11) represent the sequence from one strand of the double-stranded genomic fragment in the polynucleotide substrate, while even-numbered sub-reads (i.e., sub-reads 2, 4, 6, 8, and 10) represent the sequence from the other complementary strand of the double-stranded genomic fragment in the polynucleotide substrate. Sub-reads 1 to 8 are aligned in to emphasize this point (where the start of sub-read 9 is aligned because the polymerase displaces the synthesized strand from the polynucleotide substrate). After obtaining the data for the sub-reads, the SMCS read of the genomic fragment in the polynucleotide substrate is generated. The quality value (QV) of the SMCS read depends on the accuracy of the polymerase read and the number of sub-reads used to generate the SMCS. Currently, an SMCS generated from 10 sub-reads on the Figure 1 sequencing platform achieves a QV of 30 (see

[0045] Wenger, A. et al., "Highly-accurate long-read sequencing improves variant detection and assembly of a human genome," BioRxiv, doi.org / 10.1101 / 519025, January 13, 2019; hereby incorporated by reference in its entirety for all purposes).

[0045] As described above, any method of generating an SMCS of a genomic fragment using a single-molecule sequencing platform can be used in the assembly methods disclosed herein. Thus, the term SMCS can be used for data obtained using any single-molecule sequencing platform, e.g., in single-molecule real-time sequencing by Pacific Biosciences Sequencing of polynucleotide substrates, such as genomic fragments used in nanopore sequencing platforms such as those from Oxford Nanopore Technologies, Genia, or any other convenient single-molecule sequencing platform. For example, SMCS reads can be generated using subreads from nanopore-based single-molecule sequencing data for concatemers formed from multiple copies of genomic fragments (e.g., as described in Volden et al., PNAS 2018, v115(39), p. 9726-9731 “Improving nanopore read accuracy with the R2C2 method enables the sequencing of highly multiplexed full-length single-cell cDNA”, which is incorporated herein by reference in its entirety) or polynucleotide substrates with unique molecular identifiers (UMIs). There is no intention to limit in this regard. Thus, any consensus sequence generated from single-molecule sequencing data of multiple subreads from a single genomic fragment or its copy / copies is included in this term. For sequencing, SMCS refers to a consensus sequence determined using subreads obtained from a single polynucleotide substrate that is sequenced in a single zero-mode waveguide (ZMW) in a sequencing chip (as described above Figure 1 ). For nanopore sequencing platforms, SMCS represents a consensus sequence determined using subreads from a single original genomic fragment that is sequenced in a single nanopore, e.g., a single polynucleotide substrate that contains joined complementary strands and / or repeats derived from a single original genomic fragment (the “concateners” as described above), or from multiple nanopores, e.g., separate copies of the same original genomic fragment that are sequenced in multiple different nanopores, where for example each copy is labeled with a UMI. For examples of single-molecule sequencing platforms and methods, see the following U.S. patents and U.S. patent application publications, each of which is incorporated herein by reference: US8324914, US2013 / 0244340, US2015 / 0119259, US2010 / 0196203, US2011 / 0229877, US2016 / 0162634, US7315019, US2009 / 0087850, and US2018 / 0023134.

[0046] As used herein, the term "homopolymer collapsed sequence" or "HCS" refers to a sequence derived from a parental sequence in which each instance of multiple consecutive identical nucleotides in the parental sequence is replaced by a single nucleotide of the same type. For example, the HCS of the polynucleotide sequence AATGGGCCG is ATGCG. Thus, "homopolymer collapse", "collapsed homopolymer", etc. are used to describe the process of creating an HCS from a parental sequence (which is not an HCS).

[0047] A "homopolymer indel error" is a sequencing error in which a nucleotide identical to an adjacent and correct nucleotide in the sequence read is inserted or deleted in the sequence read. For example, inserting an incorrect G next to a correct G in a sequence read where the correct read is a single G, resulting in a GG read, is a homopolymer indel error. As another example, deleting one G from a 4-G long sequence, resulting in a GGG read instead of the correct GGGG read, is also a homopolymer indel error. A homopolymer indel error can insert or delete more than one nucleotide identical to an adjacent and correct nucleotide in the sequence read, such as a homopolymer indel of 2, 3, or 4 nucleotides. As described herein, homopolymer indel errors in the original sequence read are filtered out by the process of forming the corresponding HCS (i.e., homopolymer collapse). Thus, homopolymer collapse converts a sequencing read containing a homopolymer indel error (i.e., a sequencing read different from the genomic fragment from which it originated) into a sequence (an HCS) identical to the HCS of the genomic fragment from which the sequence originated.

[0048] A "perfect" sequence read is a sequence read whose homopolymer collapsed sequence (HCS) is identical to the HCS of the genomic fragment from which it originated. Indel errors in homopolymers in a sequence read are masked by homopolymer collapse. If the only error in a sequence read is a homopolymer indel, the read is made perfect by homopolymer collapse.

[0049] Genome assembly problems

[0050] As described above, genome assembly relies on the correct overlap of sequence reads derived from different genomic fragments. When sequence reads from two independent genomic fragments share a common nucleotide sequence that extends to one end of each read (performing a "dovetail" alignment), it can be correctly inferred that these two reads are from a pair of overlapping genomic fragments. Thus, two sequencing reads can be overlapped by superimposing this common sequence. Figure 2A simple diagram is provided that shows how two genomic fragments (A and B in the second inset) of a chromosome from a haploid genome containing the same locus (shown in the top inset) overlap. In this diagram, genomic fragment A includes nucleotides 123,000 to 133,000 from chromosome 2 (Chr2: 123000 - 133000), while genomic fragment B includes nucleotides 127,000 to 137,000 from chromosome 2 (Chr2: 127000 - 137000). Both of these genomic fragments contain nucleotides 127,000 to 133,000 (locus Chr2: 127000 - 133000). Thus, when these genomic fragments are sequenced (sequences a and b in the bottom inset), their respective sequence reads will contain a common overlapping subsequence, namely the sequence of Chr2: 127000 - 133000, which allows them to be superimposed during the genome assembly process.

[0051] When the common subsequence of two genomic fragments occurs only once in the genome, the overlap between the genomic fragments can be correctly inferred from the sequence reads of these fragments (as Figure 2 shown). However, since genomes typically contain many repetitive elements where the same or highly similar sequences occur at multiple different loci in the genome, the genome assembly process can be confounded. For example, many repetitive elements (even those that are not identical) share such a high sequence similarity that their differences are not easily detectable from their sequencing reads. In addition, long contiguous regions of repetitive sequences, such as long stretches of a 5 - base sequence, can lead to assembly errors. Thus, the basic assumption that sequence reads sharing a common sequence must necessarily originate from genomic fragments of the same locus can be invalidated by repetitive elements. Therefore, detecting the identical (or nearly identical) sequence regions shared between a pair of reads is necessary but not sufficient for the two reads to represent overlapping genomic fragments.

[0052] Tandem repeats and interspersed repeats are particularly troublesome regions that can lead to errors or breaks in assembly. Tandem repeats consist of multiple consecutive copies of a repetitive sequence motif, while interspersed repeats include sequences that occur at two or more non - adjacent positions in the genome. Figure 3 An example showing how interspersed repeats can have a negative impact on genome assembly is presented. Figure 3The upper small figure in [ ] shows genomic fragments that contain the same nucleotide subsequence but are derived from different loci in the genome. Specifically, genomic fragment A ends with the subsequence 127000 - 133000 (starting somewhere upstream), and genomic fragment C starts with the same subsequence 257000 - 263000 (ending somewhere downstream). The sequence reads of these genomic fragments (a and c in the lower small figure) can overlap in this same subsequence region. However, this overlap can lead to incorrect inferences about the underlying genome. Specifically, this overlap results in the deletion of nucleotides 1330001 to 262999 in the genome assembly. Figure 4 Shows an example of how tandem repeats can have a negative impact on genome assembly. In this figure, genomic fragments D and E include a common subsequence within a tandem repeat sequence, which has a total of 4 copies of the same nucleotide sequence, spanning nucleotides 124000 - 136000. The sequence reads of these genomic fragments (d and e in the lower small figure) can be aligned either to delete one repeat, thus folding the repeat region (middle small figure) or to add a repeat, thus expanding the repeat region (lower small figure).

[0053] The biological origin of tandem repeats and dispersed repeats is typically one or more duplication events, followed by evolutionary divergence - the independent accumulation of mutations in descendants. When two different loci share the same sequence, correct assembly is possible only if at least one read completely spans each occurrence of the common sequence. In such cases, the sequences flanking the common sequence can be used to distinguish the loci. More commonly, the genome contains repetitive elements that were originally created by the repeated duplication or insertion of the same element in an ancestral organism, but these elements have undergone mutations over multiple generations, resulting in sequence differences between them. In the absence of sequencing errors when reading genomic fragments, the differences between genomic fragments will be detected, allowing different loci to be distinguished.

[0054] In dispersed repeats, the region flanking one repeat has low similarity to the corresponding region flanking the second repeat. Therefore, it is possible to construct a contiguous assembly that bridges the dispersed repeat with two overlapping reads within the repeat, where one of the overlapping reads starts upstream of the dispersed repeat and the second read extends downstream from the dispersed repeat. In the case where the tandem repeat region consists of identical copies of a repeat motif, contiguous assembly requires reads to completely span the entire tandem repeat block because the correct registration between two reads anchored on opposite sides of the tandem repeat block cannot be determined. Specifically, bridging the tandem repeat block with two reads from opposite sides, rather than completely spanning the region with a single read, will result in the expansion or folding of the number of repeat units in the tandem repeat region (as shown in Figure 4 ).

[0055] In addition to the problem of repetitive elements occurring at different loci, there is also the problem of homozygosity in polyploid genomes, which contain multiple homologous copies of each chromosome. This is shown in the top inset of Figure 5 , where paternal chromosomes are denoted by ♂ and maternal chromosomes are denoted by ♀. The human genome is an example of a highly homozygous diploid genome, with less than 0.1% difference between homologous chromosomes. The desired assembly of a polyploid genome is a set of contigs, where each contig represents a complete chromosome and each homologous chromosome is represented by a different contig. As shown in Figure 5 , genomic segments A and B include the common locus 127,000 - 133,300 from the maternal chromosome 2. Their respective sequence reads a and b thus include the common subsequence of this shared maternal genomic locus, i.e., the sequence of locus 127,000 - 133,000. The overlap of these sequence reads (shown in the bottom inset) accurately reflects the underlying genomic structure.

[0056] At a given locus, two alleles, i.e., homologous loci on two different homologous chromosomes, are said to be homozygous if they have the same sequence at that locus. As shown in Figure 6 , genomic segments A and C include a homozygous locus in chromosome 2: nucleotides 127,000 - 133,000 of the maternal chromosome 2 and nucleotides 127,000 - 133,000 of the paternal chromosome. Their respective sequence reads a and c thus include the common subsequence of this homozygous genomic locus, i.e., the sequence of locus 127,000 - 133,000 of the maternal and paternal chromosomes. The overlap of these sequence reads (shown in the bottom inset) does not accurately reflect the underlying genomic structure. This incorrect overlap can lead to assembly errors that merge or break the maternal and paternal contigs at that locus. A correct diploid assembly requires not only determining whether two reads are from the same genomic locus but also whether they are from the same haplotype of that locus. The assumption of genomic segments caused by reads sharing a common sequence from the same allele can be invalidated by homozygous regions. When the common sequences are the same, assembly is limited by the length of the repetitive elements and the homozygous locus relative to the read length. Long reads are needed to span long regions where the two haplotypes are the same and extend to regions where there is enough variation between the haplotypes to easily distinguish them.

[0057] The ability to distinguish highly similar, different sequences depends not only on the length of these sequences but also on the degree of similarity. Noisy reads may need to be very long to fully span moderately similar regions that extend over long distances in the genome. However, if the accuracy is sufficient to distinguish intermediate regions with only moderate similarity, then even moderately long, highly accurate reads can also assemble the same region by spanning many shorter regions of the same sequence, thus anchoring the ends of the reads.

[0058] When the accuracy of two reads is so high that the number of differences between the reads is significantly higher than the expected number of read errors, reads generated from two different but highly similar sequences can be distinguished. However, it is also possible to determine that two reads come from genomic fragments of different nucleotide sequences with even higher similarity by examining the types of differences between the two reads. For example, read errors in many long-read platforms are mainly indels. For example, in Figure 7 , when analyzing the error-free reads a and b of genomic fragments A and B, they correctly overlap (the upper right small figure), while when analyzing the error-free sequence read a and the error-containing sequence read b*, the homopolymer deletion error in b* (i.e., the removal of "T" from the "TT" homopolymer) results in the failure of the reads that should overlap to overlap (based on the fact that they are derived from overlapping genomic fragments A and B in chromosome 2). Contrary to this main form of sequence read error, the true (or biological) differences between two highly similar genomic loci or heterozygous alleles are usually single nucleotide substitutions. Therefore, if the only difference between a pair of reads is a homopolymer indel, it can be inferred that the reads come from genomic fragments with the same sequence and the difference is a sequence read error. On the contrary, if the only difference between a pair of reads is a single nucleotide substitution, it can be inferred that the reads come from highly similar but non-overlapping genomic fragments, such as genomic fragments from different alleles. In extreme cases, we can separate two reads derived from genomic fragments that differ only at one position into their different haplotypes (i.e., the two reads differ by a single nucleotide variant, or SNV).

[0059] Noise filtering: True biological variation vs. sequencing read errors

[0060] An important aspect of noise filtering is to identify and utilize situations where the signal and noise are essentially in orthogonal directions in certain coordinate spaces. Regarding the genome assembly process, the signal we are considering is the true biological variation between repetitive sequence elements or haplotypes (e.g., SNVs), while the noise is the sequencing read errors (e.g., homopolymer indels).

[0061] The relationship between these signal and noise vectors is as Figure 8As shown. In this figure, two approximately orthogonal vectors representing signal and noise are shown, where the signal vector represents the biological difference, which can be used to identify when two genomic segments do not overlap and thus belong to different genomic loci and / or haplotypes (in this case SNVs), and the noise vector represents sequence read errors, which prevent the identification of two overlapping genomic segments and thus belong to the same genomic locus and / or haplotype (in this case homopolymer indels). In the genome, most of the biological differences between highly similar sequences belonging to different haplotypes and / or different genomic loci are single nucleotide variants (SNVs). In many sequencing platforms, read errors are mainly homopolymer indels (see Table 1 in Wenger, A. et al., "Highly-accurate long-read sequencing improves variant detection and assembly of a human genome", BioRxiv, doi.org / 10.1101 / 519025 on January 13, 2019; hereby incorporated by reference in its entirety for all purposes). In contrast, nucleotide substitution errors that may be mistaken for biological SNVs are relatively few. The difference between biological variation and read errors provides an opportunity for filtering. The approximate orthogonality between the signal and noise depicted in this figure means that noise can be suppressed without significantly reducing the signal strength.

[0062] The assembly process involves finding pairs of reads (R1, R2) that form a long overlapping alignment, where the suffix of R1 aligns with the prefix of R2 and vice versa. Alignments that exceed a defined threshold in length and certain sequence similarity are assumed to be true overlaps and are used for assembly. When the reads are error-free (i.e., no noise), the alignment of the suffix and prefix is an exact string match. Gusfield et al. (Gusfield, Dan, Gad M. Landau, and Baruch Schieber. "An efficient algorithm for the all pairs suffix-prefix problem." Information Processing Letters 41.4 (1992): 181-185; hereby incorporated by reference in its entirety) describe an algorithm using a suffix tree that solves the all pairs suffix-prefix problem, and its time complexity is linear in the sum of the inputs (i.e., the sum of the read lengths) and the sum of the outputs (i.e., the square of the number of reads). Since detecting pairwise overlaps between reads is considered the rate-limiting step in genome assembly, methods that accelerate this step result in significantly faster assembly.

[0063] We use well-characterized differences between haplotypes in a typical human genome as a model of biological variation (“signal”). The human genome consists of approximately three billion positions where homologous sequences on the paternal and maternal chromosomes are aligned. For this rough analysis of the mutation rate, we ignore differences in the male sex chromosomes (X and Y). In a typical person, there are approximately three million single nucleotide variants (SNVs; one nucleotide replacing another) and approximately 300,000 insertion and deletion variants (indels). The corresponding ratios of SNVs and indels are 1 to 1,000 and 1 to 10,000 (see Chaisson, Mark JP et al. “Multi-platform discovery of haplotype-resolved structural variation in human genomes.” bioRxiv (2018): 193144; hereby incorporated by reference in its entirety).

[0064] When assembling a genome or genomic region from a set of sequencing reads, it is important to overlap two reads (“overlap with a tail”) when the prefix of one read and the suffix of a second read are from the same genomic fragment. To prevent spurious overlap with a tail of reads with identical sequences (i.e., sequences from different positions in the genome that are identical to each other) in two different genomic fragments, the overlap length can be set to exceed the length of all (or most) such identical genomic fragments, e.g., from about 1,000 to about 7,000 nucleotides. It should be noted that the adjustment of the overlap length parameter can be done by the user to address specific issues known to be associated with the genome being sequenced and / or the sequencing platform being used, and thus, there is no strict threshold for the expected overlap length. Generally, increasing the minimum overlap length parameter increases the specificity of overlap detection while decreasing the sensitivity. Assemblies formed with higher sensitivity (i.e., with a lower minimum overlap length) have higher continuity but may result in joining two reads that originate from non-overlapping genomic fragments. In addition, even when the overlap of sequence reads is correctly determined (i.e., it is not the result of a sequence read error), two reads from different haplotypes that do not themselves overlap may still be joined to a third read that overlaps a homozygous region shared by the two haplotypes. For example, two reads with homozygous suffix regions can both overlap the same third read, the prefix of which includes all or part of that homozygous region. In this case, two different haplotypes may be undesirably merged into one joined component. Fortunately, these merges can generally be resolved in subsequent steps of the assembly process, e.g., by trimming the joined component of the third read to break such haplotype merges.

[0065] In the genome assembly methods described herein, we wish to avoid overlapping two sequencing reads that do not share the same contiguous subsequence for a sufficient length. Any differences between sequencing reads indicate that the two reads are derived from polynucleotide substrates containing different, non-overlapping genomic fragments, which may be either different regions of the genome or different haplotypes from the same region. In either case, incorrectly overlapping such sequencing reads will introduce errors in the assembly, such as in diploid assembly.

[0066] In some cases, two independent genomic fragments, i.e., genomic fragments that occur at different positions in the genome or are different haplotypes at the same locus, may be identical in length beyond the length threshold used to score overlapping sequencing reads. When such genomic fragments occur at different genomic positions, incorrect overlap of the sequencing reads derived from these genomic fragments will result in assembly errors. When such genomic fragments occur in different haplotypes at the same genomic position, incorrect overlap of the sequencing reads derived from these haplotypes will result in the merging of the two haplotypes, thus leading to the end of a contiguous phased block (a phased block is a region in genome assembly where haplotype sequences are separable, e.g., maternal and paternal sequences are resolved). The relative phase of two different phased blocks interrupted by a homozygous block cannot be determined. Without additional information at a scale longer than the read lengths provided, it is not possible to avoid incorrect overlaps caused by identical sequences.

[0067] Our current goal is to detect the smallest possible sequence differences between two genomic fragments with high sensitivity and specificity, i.e., single substitutions or indels within two sequence reads (e.g., two SMCS reads).

[0068] Filtering out noise (i.e., sequencing read errors) can successfully detect potential biological variations and prevent the types of assembly and consensus errors described above. The resulting assembly is more accurate, more contiguous, and has improved haplotype resolution in terms of both the length and consistency accuracy of contiguous phased blocks.

[0069] In many sequencing platforms, homopolymer indels pose a significant challenge. Consider a genomic sequence containing five consecutive A's (i.e., AAAAA). If the positions of the five A's cannot be distinguished in the read, there are five ways to generate the read sequence AAAA, i.e., by deleting any one of the five A's. Similarly, there are six ways to generate the read sequence AAAAAA, i.e., by inserting an A before the first A, after the last A, or between any two A's. Since the degeneracy of indels increases linearly with the length of the homopolymer, the single-pass error rate also increases with the length of the homopolymer.

[0070] Consensus sequences of homopolymers (e.g., SMCS reads) are particularly error-prone because the single-pass error rate in these regions is high compared to non-homopolymer indel errors (e.g., substitutions). Thus, the error distribution in consensus reads is significantly biased towards homopolymer indels and away from other types of errors. The enrichment of homopolymer indel errors as the major error type in consensus sequence reads increases with the length of the homopolymer region and the number of reads used to generate the consensus. For SMCS reads, the higher the number of subreads, the higher the proportion of homopolymer indel errors in the total sequence errors. For example, in an SMCS read formed by 10 subreads by a nucleic acid sequencing instrument of Pacific Biosciences, approximately 99% of the errors are homopolymer indels.

[0071] The prevalence of homopolymer indel errors means that a high read coverage (a combination of single-molecule and multi-molecule reads) is required to reliably determine the length of long homopolymers. However, the concentration of SMCS read errors in a single channel (i.e., homopolymer indels) provides an opportunity for noise filtering in the genome assembly process.

[0072] Recall that haplotype variations in the human genome are 90% SNVs and 10% indels. Approximately one-quarter of these occur in homopolymers. Thus, only a few true human haplotype variations (signals) are homopolymer indels. Therefore, when we observe that two aligned reads (e.g., SMCS reads) differ only in indels in the homopolymer region, the difference is likely a read error (noise) and the reads are from the same genomic fragment.

[0073] This property provides the basis for methods to suppress read errors (noise) to reveal subtle biological sequence variations (signals). Specifically, the sequence alignment method described herein eliminates the confounding effect of homopolymer indel errors by reducing homopolymer runs in sequence reads to a single base of the same type (referred to as homopolymer folding) before alignment. Reads that differ only by homopolymer indels become identical after homopolymer folding and can be paired by exact string matching. For example, in Figure 9 , the reads a and b* shown in the upper right inset (identical to Figure 7 ) can be correctly overlapped by first converting them to their homopolymer-folded forms, which masks the homopolymer indel error in b* (see Figure 8(lower right inset in). Due to the highly skewed error distribution (mainly homopolymer indels), sequence reads that align by exact string matching over most of their length (e.g., 100, 200, 300, 400, 500, 750, 1000, 2000, 3000, 4000, 5000 bases or more) after homopolymer folding are assumed to be from the same genomic fragment and to overlap. Combinations of many such exact sequence overlaps form the basis of the draft assembly.

[0074] In the current polyploid genome assembly process, the draft assembly is "polished" to resolve inconsistencies in the multiple sequence alignment of the aligned reads, resulting in a consensus sequence for each haplotype. In many cases, polishing a polyploid genome assembly involves an iterative process of partitioning the reads into haplotypes and then calling a consensus sequence for each partition.

[0075] In contrast to this iterative polishing process, the draft assembly produced by exact string matching of overlapping homopolymer-folded reads described herein is largely haplotype-resolved, except that long homozygous regions not spanned by a single sequence read may lead to haplotype merging. In the exact string matching-based method described herein, different haplotype blocks are formed by removing sequence reads that fall entirely within regions of overlap where all aligned positions are identical (i.e., for each position in the sequence read, if there is only one base represented at that position in all overlapping reads, that read is removed). Once these reads are removed, at each position within a haplotype block, all reads belonging to that haplotype have the same nucleotide at that position. Thus, for a diploid genome, there should be at most two haplotype blocks for a genomic region. At this point, the consensus for each haplotype block is trivial because, by its construction, homopolymer-folded reads mapped to the same haplotype are identical at each aligned position. Thus, the homopolymer-folded consensus sequence for each haplotype is determined by simply calling the consensus base at each aligned position. Thus, in genomic regions where multiple haplotype blocks are formed, this process results in a consensus sequence for each distinct haplotype, which, by definition, differ at one or more positions. For genomic homozygous regions that interrupt haplotype blocks, since no single sequence read spans the entire homozygous region, reads from both haplotypes produce a single consensus sequence.

[0076] After dividing the genome into heterozygous and homozygous regions and assigning a shared homopolymer fold sequence to each haplotype in each homozygous region and phased block, the remaining step is to generate a complete polyploid assembly by re-extending the homopolymer fold sequences by making a consensus call on the haplotype-resolved homopolymer lengths. As described herein, when each sequence read is folded, the length of its homopolymer is recorded. For each homopolymer, the set of lengths of that homopolymer in the aligned reads of a given haplotype is used to determine the consensus length call (an example of this process is described below).

[0077] As noted elsewhere herein, aspects of the present disclosure employ single molecule consensus (SMCS) reads, which are formed by obtaining multiple individual reads derived from a single original polynucleotide fragment (e.g., a single genomic fragment) and combining them to form a single consensus sequence of that original polynucleotide fragment. As with multimolecule consensus, where reads from different original polynucleotide fragments are aligned and analyzed, the redundancy in the multiple reads used to generate SMCS provides a mechanism for suppressing read noise (i.e., sequencing errors). Unlike multimolecule consensus, the multiple reads used to form SMCS reads are known to be from the same original polynucleotide fragment, thus eliminating the possibility of mapping errors. This allows SMCS reads to be "polished" to high accuracy before overlapping with other SMCS reads. The high accuracy of SMCS reads may be sufficient to distinguish sequences derived from mutually different but highly similar genomic fragments that cannot be distinguished by lower accuracy single-pass reads.

[0078] Errors in SMCS reads are a direct result of errors in the single-pass reads from which they are derived. In platforms where indels are the primary error type (in single-pass reads), indels will also be the primary error type in SMCS reads. Error types that occur less frequently in single-pass reads (e.g., substitutions) tend to be quickly "weeded out" from SMCS reads. In general, as the number of subreads increases, each type of single-pass error is exponentially weeded out from SMCS reads. The exponential factor that determines the incidence of a particular error type in SMCS reads is the incidence of that error type in single-pass reads. Thus, when comparing error rates in SMCS reads, variations in the error rates of various types of single-pass reads are amplified.

[0079] Computer-implemented analysis

[0080] Aspects of the methods presented herein can be embodied, in whole or in part, as software that is recorded on a fixed medium for use in a computer (or computer system). The computer can be any electronic device having at least one processor (e.g., a CPU, etc.), a memory, an input / output (I / O), and a data repository. The CPU, memory, I / O, and data repository can be connected by one or more system buses or using any type of communication connection. The computer can also include a network interface for wired and / or wireless communication. In one embodiment, the computer can include a personal computer (e.g., a desktop computer, a laptop computer, a tablet computer, etc.), a server, a client computer, or a wearable device. In another embodiment, the computer can include any type of information device for interacting with remote data applications and can include devices such as an Internet-enabled television, a mobile phone, etc.

[0081] The processor controls the operation of the computer and can read information (e.g., instructions and / or data) from the memory and / or data repository and execute the instructions accordingly to implement the exemplary embodiments. The term processor is intended to include one processor, multiple processors, or one or more processors having multiple cores.

[0082] For example, the I / O can include any type of input device, such as a keyboard, a mouse, a microphone, etc., and any type of output device, such as a monitor and a printer. In an embodiment where the computer includes a server, the output device can be coupled to a local client computer.

[0083] Generally, the present disclosure provides computer-implemented methods that employ a homopolymer collapse sequence (HCS) to improve alignment sequences, determine consensus sequences, map sequences to a reference, and / or sequence assembly processes, such as in de novo assembly of a genome. As defined above, an HCS is a sequence derived from a parental sequence in which each instance of multiple consecutive identical nucleotides in the parental sequence is replaced by a single nucleotide of the same type. For example, the HCS of the polynucleotide sequence AATGGGCCG is ATGCG. It should be noted that each HCS stores the length of each collapsed homopolymer, so this information is not lost. These stored homopolymer lengths are used for downstream analysis, e.g., making consensus homopolymer length calls for haplotype resolution to refine a draft genome assembly.

[0084] As described herein, homopolymer folding allows for a significant improvement in sequence analysis when applied to sequencing platforms where the primary type of sequencing error is homopolymer indel errors. As defined above, homopolymer indel errors are errors in which a nucleotide identical to an adjacent and correct nucleotide in a sequencing read is inserted or deleted. Applying homopolymer folding to a sequencing read containing a homopolymer indel error and a reference sequence (or the polynucleotide substrate sequence from which it is derived) to which it is compared results in a perfect match between the sequences. In other words, the homopolymer indel errors are masked and thus do not negatively impact sequence alignment algorithms. Additionally, homopolymer folding of multiple sequencing reads allows for computer-implemented contig and genome assembly using exact string matching rather than relying on similarity thresholds or error-tolerant algorithms that use short k-mer seeds (e.g., k < 30) and linked exact matches.

[0085] The differences between the homopolymer folding / exact string matching method detailed herein and the k-mer matching method are as follows. In current practice, k-mer matching is used to identify short common subsequences shared by two reads, which may be part of the overlapping region between the two reads. However, even if the aligned region contains sequence differences between the two reads, i.e., differences between the regions of the sequence between the identified perfect k-mer matches, the two reads can be judged to overlap (i.e., will have originated from overlapping genomic fragments). Thus, k-mer matching is error-tolerant. In contrast, exact string matching is not error-tolerant and thus is not simply k-mer matching as currently performed with longer k-values. Instead, exact string matching only judges two reads to overlap when the overlapping regions between the two reads are identical, i.e., when there are no differences between the reads across the entire overlapping region. Because exact string matching is not error-tolerant, exact string matching of homopolymer-folded sequences can significantly speed up the alignment, consensus, and assembly processes (as described below). Additionally, for SMCS reads and other read types where homopolymer indels are the primary error type (e.g., nanopore sequencing), exact string matching has higher sensitivity and specificity and can be used to identify true overlaps between the genomic sequences from which a pair of reads was obtained.

[0086] In some embodiments of the present disclosure, the sequence reads employed are single molecule consensus sequence (SMCS) reads, which can be derived from any sequencing platform on which SMCS reads can be generated, e.g., a sequencing platform and a nanopore sequencing platform. Generally speaking, SMCS reads are consensus sequences generated by analyzing multiple single-pass sequence reads that are derived from the same original polynucleotide substrate molecule, e.g., by repeated sequencing of the original polynucleotide substrate (such as sequencing) or by sequencing multiple copies of the original polynucleotide substrate (as in sequencing concatemers generated by rolling circle amplification or other means using nanopore sequencing). (See, e.g., Figure 1 and the description above it.) Note that in sequencing applications, sequencing of concatemers can be achieved by generating polynucleotide substrates, each of which includes a concatemer derived from a single polynucleotide substrate, and / or by generating multiple polynucleotide substrates, each of which includes a copy from the same original polynucleotide substrate. In addition, certain nanopore sequencing methods can be used to sequence topologically circular polynucleotide substrates, such as the technology from Genia, now part of Roche (see Fuller et al., 2016, PNAS 113(19):5233-8, which is hereby incorporated by reference in its entirety). Thus, no limitation is intended in this regard.

[0087] It should be noted here that while SMCS reads are described for use in the subject methods, the methods described herein are not limited to SMCS reads. In fact, the methods described herein are applicable to any sequence reads for which homopolymer indel errors are a significant or major sequence read error type and thus a confounding problem in genome assembly, including single-pass sequence reads. No limitation is intended in this regard.

[0088] Current algorithms for read mapping and alignment involve a fast screening step based on detecting one or more perfect k-mer matches between sequences, followed by a dynamic programming step to find the best sequence alignment. The fast screening step involves a trade-off between specificity and sensitivity, which is regulated by the choice of k, the length of the k-mer. The larger the k value, the less likely it is for two sequences to overlap randomly. The smaller the k value, the less likely it is for sequencing read errors to mask a match to the correct target (i.e., the locus from which the read is derived or another read from the same locus). Reducing the number of differences between a sequencing read and its target (e.g., other sequencing reads, reference sequences, etc.) means that a larger k value can be used without losing sensitivity to correct matches. However, as described above, current k-mer alignment algorithms are error-tolerant, so some form of polishing is required to arrive at a consensus in the overlapping regions of sequence reads, which may include sequence differences outside the aligned k-mer regions.

[0089] Dynamic programming is a method for exploring all alignments between two sequences in time proportional to the product of the sequence lengths. If the sequences are error-free, the alignment can be found in time proportional to the length of the longer sequence (i.e., linear time). By classifying the HCS of sequence reads as error-free, e.g., the HCS of SMCS reads, we can take advantage of this feature of dynamic programming by requiring exact string matching to align the sequences (instead of using current k-mer matching).

[0090] The signals and noise in sequencing data are not completely "orthogonal". For example, while the vast majority of read errors (noise) in sequencing platforms are homopolymer indels, there are occasionally cases where genomic fragments have biological homopolymer indel differences (signals), e.g., a genomic fragment from the first haplotype at a genomic locus will differ from a genomic fragment from the second haplotype at the same genomic locus in the length of the homopolymer sequence. Based on our current understanding of the human genome, the sequences of two 5kb genomic fragments with 99.9% similarity will on average differ by approximately five nucleotide substitutions and approximately 0.5 indels. For indels, approximately 0.4 of the 0.5 indels occur outside homopolymers and approximately 0.1 of the 0.5 indels occur within homopolymers. A 5kb error overlap occurs between two SMCS reads when there are no substitution differences in the 5kb overlapping region, no indels outside homopolymers, and no one or more homopolymer indels. Thus, even if the genomic fragments from which the SMCS reads are derived are highly similar, a 5kb error overlap between SMCS reads is extremely unlikely. The vast majority of error overlaps result in failure to identify heterozygous variants, leading to phase block collapse, which most commonly occurs outside coding regions. Error overlaps that result in incorrect genome assembly can occur within repetitive regions, where a large number of repetitive elements have very high sequence similarity, such as centromeres, but are otherwise extremely unlikely. Even so, the ability to detect single-base differences (most commonly substitutions) between genomic fragments has greatly increased the average length of phase blocks in highly homozygous genomes (such as the human genome).

[0091] In some embodiments, the present disclosure takes advantage of the unique properties of long SMCS reads (e.g., 10 - 15 kb or longer) that can be generated from long-read sequencing techniques, such as those that produce reads of 50 kb, 75 kb, 100 kb, 150 kb or longer. Specifically, the long read lengths result in a large number of subreads (e.g., 4, 5, 6, 7, 8, 9, or 10 subreads or more) being obtained from an original polynucleotide substrate of approximately 10 - 15 kb in length, which can be used to generate SMCS reads with an accuracy of 99% to 99.99% or higher. In some embodiments, the polynucleotide substrates analyzed according to the present disclosure are derived from genomic DNA samples, where in some cases the genomic DNA samples are from polyploid organisms, such as plant, fungal, animal, or human genomes. In other cases, the sample is a metagenomic sample containing a variety of different microorganisms such as bacteria, protozoa, yeast, or other single-celled organisms. These SMCS reads greatly reduce homopolymer indel errors, including substitution errors (errors where one base is changed to a different base, e.g., reading the polynucleotide substrate sequence AGCTG as AGATG) and indel errors where a nucleotide base is inserted or deleted that is different from two adjacent bases (e.g., reading the polynucleotide substrate AGCTG as ATGCTG or ACTG). For sequencing, we found that all types of errors decrease exponentially with the number of passes.

[0092] Based on the above discussion, it is clear that most errors in SMCS reads (e.g., generated from approximately 4 - 10 subreads or more) are homopolymer indels. Since most biological variations are single nucleotide variations (one base replacing another), the SMCS read error types show very low overlap with true biological variations. Thus, removing homopolymer indels in SMCS reads by homopolymer folding (thereby generating HCS reads) preferentially removes errors based on the sequencing platform while leaving true biological variations. Therefore, filtering out these errors will improve many downstream sequence analysis algorithms, from mapping and alignment to de novo genome assembly. Once any desired downstream alignment of the HCS reads is complete, the folded homopolymers of each HCS read can be extended (based on their length in the original SMCS read). The extended homopolymer regions of the SMCS reads can then be analyzed to determine the consensus length at each distinct position. These consensus homopolymer lengths can then be added back to any consensus sequences generated from the processes using the HCS reads (e.g., assembly, alignment, and / or any generated consensus sequences).

[0093] The following figures and their descriptions are intended to illustrate certain embodiments of the methods disclosed herein and are not intended to be limiting. For example, while the following description relates to HCS from SMCS reads, HCS from single-pass sequence reads can be employed where homopolymer indel errors are the primary or significant error type.

[0094] Figure 11 Shows an example of aligning SMCS read pairs after filtering out homopolymer indels, which represent the vast majority of sequencing errors. The shaded blocks represent homopolymer indel errors, which are the major error type in SMCS. The solid blocks in SMCS3 represent single nucleotide variants (SNVs), which identify SMCS3 as originating from a different haplotype than SMCS1 and SMCS2. Homopolymer indel errors are masked by homopolymer folding and ignored when determining whether two reads are from the same haplotype. Since the only difference between SMCS1 and SMCS2 is homopolymer indels, and homopolymer indels are assumed to be read errors during the assembly overlap step, it is assumed that SMCS1 and SMCS2 originate from the same haplotype (the same genomic fragment). In contrast, single nucleotide substitution differences are considered to be real biological differences between haplotypes.

[0095] Figure 12 Shows a toy example of a multiple sequence alignment formed by pairwise exact string matching of SMCS reads. Pairwise exact string matching can be simply characterized by integer offsets. Multiple sequence alignments are generally very complex and trivial for exactly string-matching reads from the same haplotype. Exact string matching is transitive, while offsets are additive.

[0096] Figures 13 to 15 Shows an embodiment of a sequence analysis pipeline that uses homopolymer folding and exact alignment mapping to separate SMCS reads into haplotypes. Although these figures depict haplotype separation for a diploid genome (e.g., the human genome), this analysis pipeline is applicable to any sequence analysis that requires separating SMCS reads into groups of sequences originating from the same original genomic / polynucleotide substrate, e.g., in metagenomic sequence analysis. The analysis pipeline also handles genomes with higher ploidy, such as tetraploid (n = 4), hexaploid (n = 6), or octaploid (n = 8). No limitation is intended in this regard.

[0097] In Figure 13In the first step of the middle pipeline, select SMCS reads that map to a specific region of the reference genome. This step is not an essential feature of the algorithm but is used here to construct a problem of limited scale, namely, the haplotype-resolved assembly of the highly similar SMN1 and SMN2 loci, allowing for an easy-to-understand demonstration of the algorithm's utility. This initial mapping can be performed with relatively low stringency to maximize the number of SMCS reads available for downstream analysis, as reads that are mis-mapped to this region can be easily filtered out during the assembly process. One or more regions can be selected by the user, for example, regions that are associated with or predicted to be associated with a phenotype (e.g., a disease phenotype). Once a subset of SMCS reads that map to the region(s) of interest (or regions of interest) (represented as a "heap of disorganized SMCS reads" in Figure 13 are obtained, they are converted to HCS reads and proceed to a "whole-to-whole" pairwise alignment with strict filtering, as described herein ( Figure 13 step 2). For example, the alignment can be filtered such that the aligned region is (1) at least 1 / 4 to 1 / 2 the length of the average sequence read length (or a threshold minimum length that is predicted to span a homozygous region in the genome under study, e.g., ~1 kb to ~5 kb), and (2) an exact match between the suffix of one read and the prefix of the other read. The alignments on the right side of step 2 meet these criteria and are processed in step 3, with the aligned regions represented by right-facing arrows. All pairwise alignments that do not meet these criteria are discarded or placed in a reservoir. SMCS reads that contain any read errors other than homopolymer indels will not form an exact string match with other reads and will also be placed in the reservoir. The alignments on the left side of step 2 are placed in the reservoir because they have multiple mismatches (represented by "*") in the aligned region. In step 3, an overlap layout algorithm is used to compare and separate the aligned regions of all pairwise alignments that meet this filtering requirement (represented by arrows), where pairwise alignments that have an exact overlap in their respective aligned regions are separated into the same group (or haplotype, as Figure 13Among them; haplotypes 1 and 2). Reads belonging to different haplotypes are determined by treating the alignments between reads as vertices and edges in a graph respectively and finding the connected components of the graph. In this case, each alignment between a pair of reads indicates that the two reads may belong to the same haplotype, but also provides the relative offset between the read start positions, which will require arranging the corresponding positions of the aligned sequences. These pairwise offsets can be used to arrange a set of connected reads along a common coordinate axis, as shown in step 3. In this case, each small graph contains a set of reads belonging to the same haplotype. Therefore, at any given position in the multiple sequence alignment, all reads covering that position have the same base call at that position. Regions forming pairwise alignments that do not overlap with any other regions forming pairwise alignments are placed into a reservoir. These orphan pairwise alignment regions may come from SMCS reads that were mis-mapped to the region of interest in step 1 and / or SMCS reads that may come from polynucleotide contaminants or sample preparation artifacts (e.g., unintentional mixing of the initial genomic DNA sample or generation of chimeric polynucleotide substrates and / or amplification products during sample preparation, etc.). The criteria for putting pairwise alignments (and / or their SMCS reads) into the reservoir can be determined by the user and can be based on known information about the genomic sample, such as the ploidy or expected number of organisms in a metagenomic sample, sample preparation details, etc. In this way, reads can be grouped by haplotype according to the differences observed in the pairwise alignments.

[0098] As Figure 14 shown, a consensus sequence is then generated for each haplotype or overlapping sequence group (step 4). The consensus sequence of a haplotype is determined by reading the base calls at each position in the sequence. The consensus sequence here represents the consensus sequence of the homopolymer folds of each haplotype / group. After generating the consensus sequence from the HCS, the homopolymer fold regions can be extended in step 5 to generate a consensus sequence with homopolymer extension. This process involves attaching the homopolymer lengths observed and recorded at each folded position of each read, converting a set of aligned homopolymer fold reads (HCS) into a set of aligned homopolymer extended reads (HES). Note that the alignments of these reads are retained because we "extend" each homopolymer, not by a string of repeated nucleotides to represent the homopolymer, but as a base call and a repeat number. For example, a homopolymer of 4 A's is represented by "A4" rather than "AAAA" (the top HES read in step 5). Figure 14The right small figure shows two positions in the multiple sequence alignment where the (extended) homopolymer lengths in the reads are inconsistent. In this example, to form homopolymer length calls at these positions, we find the floor of the median. We use the floor of the median because the common homopolymer length must be an integer value. We choose the floor instead of the ceiling because shorter homopolymers occur more frequently than longer ones. By calling the homopolymer length at each position in the homopolymer collapsed consensus sequence, we form a run length encoding representation of the homopolymer extended consensus sequence. Now, we expand each run length encoded homopolymer into a string of repeated nucleotides, e.g., converting "A4" to "AAAA", to produce the final homopolymer extended consensus sequence, as Figure 14 shown.

[0099] An example of homopolymer extension is as follows. First, a vector of homopolymer lengths is associated with each position in the homopolymer collapsed sequence, where (i) the number of elements in the vector is the number of trimmed HCSs in the multiple sequence alignment that cover that position, and (ii) each component of the vector is the homopolymer length observed in the original reads at that position in the HCS. For example, in Figure 14 , the vector for the "A" nucleotide at position 2 in the HCS is from the corresponding position in the HES and is thus: 4, 4, 4, 4, 3, 4. Next, the consensus homopolymer length at each position in the homopolymer collapsed sequence is calculated as the floor of the median of the components of the homopolymer length vector associated with that position, e.g., the floor of the median of the lengths derived from the corresponding positions in the HES. In Figure 14 , the value is 4 because the floor of the median of the series 3, 4, 4, 4, 4, 4 is 4. Finally, each position in the homopolymer collapsed sequence is replaced with a homopolymer string N of the same nucleotide, where N is the consensus homopolymer length calculated for that position.

[0100] As Figure 15 shown, once the consensus sequences of homopolymer extension are called in step 5, these consensus sequences are compared with the genomic reference sequence (e.g., the genomic region used to select the initial SMCS reads) in step 6 to call any heterozygous variants (represented as 1, 2, and 3) and / or homozygous variants (represented as 4). In some embodiments, if there are low coverage regions in the consensus sequence, reads in the reservoir can be used to confirm the variant calls. This is shown in Figure 15 as the dashed arrow for variant 3 from the HCS reads in the reservoir, which supports the call of variant 3 in the haplotype 2 consensus. Note that variant positions can occur in homopolymer regions because they have been extended. Analyzing the reads in the reservoir by extending the homopolymer regions may also help to determine the consensus homopolymer length, if beneficial.

[0101] Perfect reads, such as SMCS reads whose errors are completely masked by homopolymer collapses (as defined above), participate in diploid assembly by exact string matching with other perfect reads. When the prefix of the polynucleotide substrate HCS that gives rise to one read is the suffix of the polynucleotide substrate HCS that gives rise to another read, the two perfect reads overlap during assembly, forming a perfect dovetail alignment. This alignment is required to produce an accurate genome assembly.

[0102] Accordingly, to maintain the accuracy of genome assembly, we desire to exclude from the assembly process reads that are not yet perfect. The requirement that an overlap be formed only when the HCSs of two SMCSs exactly match has the effect of excluding many reads with errors not masked by homopolymer collapses. Except for rare coincidences, reads containing (unmasked) errors near either end will not exactly match any other read.

[0103] However, we must also consider the case of reads with a single (unmasked) error. Roughly speaking, such a read has one perfect half that will overlap with other perfect reads, but the other half with the error will not overlap with other reads. The read is retained in the analysis because it forms a perfect dovetail assembly with a perfect read. One possible outcome is that such a read will terminate a contig in the assembly because only one side of the read forms a perfect dovetail alignment. Another possibility is that such a read will result in a "spur", which resembles a unique haplotype variant that forms a branch separate from other perfect reads.

[0104] To avoid the adverse effects of including this type of imperfect SMCS read during alignment, we remove these errors from such SMCS reads by trimming any positions at the ends of the reads that cannot overlap with any other read before the layout step of the assembly. This quality control step ensures that all bases used in the assembly process are represented at that position by at least two separate SMCS reads. In embodiments where the threshold overlap length is at least half the average read length, positions not covered by at least one overlap may be at the ends of the reads.

[0105] Before the layout step (e.g., Figure 13 step 3 in ), we first generate a graph that represents the pairwise overlaps between reads. Each read is represented by a vertex in the graph. Each overlap between a pair of reads is represented by an edge between the corresponding vertices. In an ideal situation, the connected components of the graph would represent a chromosome (e.g., one haplotype of the genome). In a diploid genome, each paternal chromosome has one component and each maternal chromosome has one component. Different chromosomes will be represented by different connected components.

[0106] However, in many cases, due to fragmentation during assembly, the chromosome is represented by multiple connected components. Fragmentation may be caused by systematic and / or random coverage loss, leaving some positions not covered by any reads. In the currently disclosed algorithms, the continuity of the assembly at a position requires that the position be covered by at least two perfect SMCS reads.

[0107] In addition to fragmentation, the connected components may represent the merging of segments from multiple chromosomes. Most commonly, the merged connected components are caused by homozygous regions shared by two or more haplotypes. For example, as Figure 16 shown, reads A and B belong to different haplotypes and contain one or more positions where the haplotypes differ (represented by the "x" positions), and thus do not overlap. However, both reads A and B overlap with a third read C. The overlap between A and C contains only homozygous positions, i.e., positions where the two haplotypes have the same sequence. Similarly, the overlap between B and C contains only homozygous positions. In this case, reads A and B, which belong to different haplotypes, are merged into the same connected component through their mutual overlap with read C in the genomic homozygous region. In Figure 16 it, reads D and E that vary at position "y" overlap with the other end of read C in a similar manner. Thus, read C contains only homozygous positions at this locus in the genome; it contains neither x nor y.

[0108] In Figure 16 this alignment scheme results in the figure labeled "merged haplotypes". The subgraph of the connected components is induced by removing the edges representing the overlaps that contain only homozygous positions (e.g., by removing node C from the graph) to separate such merged haplotypes (or "connected components"). This process is called pruning. For example, the overlaps between A and C and between B and C will be removed, and the overlaps between D and C and between E and C will also be removed. If there is no SMCS read that contains both positions x and y, then removing read C divides the graph into four connected components, as shown in the "separated but unresolved haplotypes" box. For a diploid genome, there are two possible layouts for these homologous haplotypes: 1) A is connected to D through C, and B is connected to E through C (as shown in the upper right layout of Figure 16 ); or 2) A is connected to E through C, and B is connected to D through C (as shown in the lower right layout of Figure 16 ). Thus, the homozygous region between the flanking resolved haplotype regions induces a haplogroup break rather than a contig break (i.e., the haplotypes cannot be resolved, but the contigs through this region remain intact).

[0109] Figure 17 shows the same as Figure 16Situations related to those described in , except that the set of sequence reads (shown in the upper left) includes reads F and G, each of which spans positions x and y, i.e., spans the homozygous region. These reads can be used to resolve the two haplotypes. If F overlaps with reads A and D (meaning it contains the same variants as reads A and D at positions x and y) and read G overlaps with reads B and E (meaning it includes the same variants as reads B and E at positions x and y), then removing the edges connected to the vertex associated with read C (i.e., pruning) results in a graph with two connected components, one for each contiguous haplotype (as shown on the right).

[0110] Example: SMN1 / SMN2 genomic region

[0111] In the following example, the survival of the motor neuron 1 and 2 loci (SMN1 and SMN2) is analyzed according to an embodiment of the present disclosure. SMN1 and SMN2 are part of a 500 kb inverted duplication on chromosome 5q13, with SMN1 being the telomeric copy and SMN2 being the centromeric copy. These genes encode the same protein, SMN. This duplicated region contains at least four genes and repetitive elements, making it prone to rearrangement and deletion. The repetitive and complex nature of the sequence also makes it difficult to determine the organization of this genomic region. Mutations in the telomeric copy SMN1 are associated with spinal muscular atrophy (also known as Werdnig-Hoffmann disease or Kugelberg-Welander disease); mutations in the centromeric copy SMN2 do not cause disease. The centromeric copy may be a modulator of the disease caused by mutations in the telomeric copy. Mutations in SMN1 and SMN2 result in embryonic lethality. The key sequence difference between the two genes is a single nucleotide in exon 7, which is thought to be an exon splicing enhancer. The 9 exons of the telomeric and centromeric copies have historically been designated as exons 1, 2a, 2b, and 3 - 8. Gene conversion events are thought to possibly involve the two genes, resulting in different copy numbers of each gene. The protein encoded by this gene is localized in the cytoplasm and nucleus. Within the nucleus, the protein is localized in subnuclear bodies called gems, which are found near the coiled bodies containing a high concentration of small nuclear ribonucleoproteins (snRNPs). This protein forms hetero-oligomeric complexes with proteins such as SIP1 and GEMIN4, and also interacts with several proteins known to be involved in snRNP biogenesis, such as hnRNP U protein and small nucleolar RNA-binding protein. Two transcript variants encoding different isoforms have been described.

[0112] Figures 18 to 20 Preliminary results of the diploid assembly of the SMN1 and SMN2 regions from a set of SMCS reads are shown. Figure 21 The final results are shown. The data and assembly process are described in more detail below.

[0113] We first obtained human genome (HG002) DNASMCS reads from a narrow band (+ / - 1 kb) fragment centered at 13.5 kb (these reads are described in Wenger, A. et al. "Highly-accurate long-read sequencing improves variant detection and assembly of a human genome" BioRxiv, doi.org / 10.1101 / 519025; for all purposes, it is hereby incorporated by reference in its entirety). We used a subset of these reads that were mapped to SMN1 or SMN2 with relatively low stringency by minimap2 (some may map to both because of their very high sequence similarity). A histogram of the lengths of the SMN-mapped SMCS reads selected for this analysis is shown in the upper left inset of Figure 18 . This led to the selection of 154 SMCS reads.

[0114] Next, we made reverse-complementary copies of each SMCS read, forming a set of 308 SMCS reads. Note that the initial set of SMCS reads represents reads of genomic fragments from both strands of the genome. We consider two genomic fragments to be "overlapping" if one fragment overlaps the reverse complement of the other. By making two "mirror" copies for each read, we will form two mirror assemblies from the set of reads. The genomic reference (arbitrary) represents one of the two strands, so we retain the assembly corresponding to the reference strand. Then, we generated homopolymer-folded sequences (HCS) for each SMCS read. A histogram of the lengths of the HCS is shown in the lower left inset of Figure 18 . The average length of the HCS was 9.5 kb. A histogram of the ratio between the HCS and SMCS lengths is shown in Figure 18Homopolymer folding reduces most SMCS reads to 69-70% of their original length. For comparison, folding of strings generated by randomly and independently drawing four letters with equal probability reduces the strings to 75% of their original length. Next, we performed an all-to-all pairwise alignment of these 308 SMCS reads. A total of 494 alignments were formed between pairs of reads. A pair of reads has an alignment between them if the suffix of one read is the same as the prefix of the other read and the length of this common subsequence is longer than the minimum overlap length. Here, we chose a minimum overlap length of 6kb, a value just over half of the longest HCS in the set. These candidate alignments were then divided into groups based on their connectivity. For this step, the aligned reads are represented as a graph where reads are vertices and alignments are directed edges. Reads whose suffix matches the prefix of another read point to this directed edge.

[0115] The graph resulting from 494 alignments between 308 reads has twelve connected assemblies between 200 reads—six pairs in which the members of the pair are mirror images of each other. The other 108 reads are singleton reads that do not overlap with any other read. Most likely, these singleton reads failed to overlap with the other reads because they were corrupted by one or more read errors. Because we chose the minimum overlap to be greater than half the length of any HCS, a single read error at the midpoint of a read will result in a read that fails to overlap with any other read—that is, except in the extremely unlikely event that another read is identical at 6000 or more positions, none of which contains an error except one of exactly the same type at exactly the same position. More commonly, read errors at either end of a read exclude it from assemblies built from overlapping reads based on exact string matches.

[0116] The process of determining the connected components of the graph also generates a layout of the reads within each component. Components are formed by performing a breadth-first traversal from an arbitrary read and assigning that read an arbitrary coordinate value of zero. The prefix of each newly arriving read in the traversal matches the suffix of a read that has already arrived, so that the coordinates of each new read are at least as large as a read that was already part of the traversal. When a traversal from a new read touches a read that has already been assigned to a component, the two components are merged. The coordinates of all reads in the newly touched component are increased by a fixed offset, making the coordinates in the merged component self-consistent. Figure 19The upper small figure shows the layout of component 3, which consists of 11 HCSs from the SMCS reads. This layout covers approximately 20 kb, but after trimming, there are only 17,577 bases. The slightly thicker lines extending from the ends of four HCSs (arrows) show the regions of the reads that were trimmed because these regions do not overlap with any other HCSs in the set, most likely due to read errors. These trimmed base pairs do not contribute to the assembly. Positions at the left and right ends of the layout that are not represented by at least 2 HCS reads (covered by only one of the HCS reads) are trimmed and not used to form the consensus sequence. Figure 19 The bottom small figure shows the number of variant base calls in the multiple sequence alignment caused by the HCS layout. In this case, at each alignment position covered by at least two reads (i.e., excluding the trimmed regions), each read covering that position has the same base call. The reads provide a consistent consensus sequence for each base call in the consensus sequence. Another way to describe this multiple sequence alignment is that each constituent HCS is a correct (exact) substring of the consensus sequence of the homopolymer fold. Figure 19 The zero values in the variant curve in the bottom small figure correspond to positions covered only by the trimmed regions of the reads.

[0117] In contrast to Figure 19 the case shown, Figure 20 shows a connected component containing two merged haplotypes. In Figure 20 , 54 HCSs form a connected component that spans nearly 40 kb before trimming. Figure 20 The variant curve of this component shown in the bottom figure of Figure 20 indicates that while most positions among all reads at that position are consistent, some positions contain inconsistent reads. At each inconsistent position, the reads can be divided into two groups, which are defined by the base calls they contain at that position. These two groups represent two different haplotypes. Three reads (

[0118] Figure 21Shows the final diploid assembly, where the consensus sequences representing each contig are mapped to the sequences of SMN1 and SMN2 that appear in the human genome reference GrCh38. Prior to the assembly process, most individual SMCS reads could not be reliably mapped to SMN1 or SMN2 because the similarity between the reference sequences was higher than the similarity between the SMCS reads and either reference sequence. Although the SMCS reads have high accuracy, it remains ambiguous whether the genomic fragments from which the SMCS reads were obtained originated from the SMN1 locus or the SMN2 locus. Even the homopolymer collapses of the SMCS reads, which eliminate most read errors, do not resolve this ambiguity for most reads. Figure 21 The mapping shown in Figure 21 is possible only because several nucleotides in exons 7 and 8 distinguish SMN1 from SMN2. This allows us to map a limited number of reads to the correct locus, but only in that region. However, because we have a diploid assembly, the joining of these "mappable" reads with other reads belonging to the same haplotype anchors the entire haplogroup to the correct locus. The variant positions marked between the SMCS reads and the reference allow us to call variants across the entire length of both loci. The consistency of these variants among multiple aligned reads provides strong evidence for the correctness of these variant calls. At many positions, heterozygous variants are evident, and the two haplotypes can be clearly identified.

[0119] The SMN1 and SMN2 loci are difficult to assemble because of their high similarity and high homozygosity. In the high-quality assemblies recently obtained from these reads, current assemblers were unable to map the reads to any exons from either SMN1 or SMN2. [This is noted in Wenger et al. Figure 2 c, where the SMN1 and SMN2 exons are listed as 0% mappable (Wenger, A. et al. "Highly-accurate long-read sequencing improves variant detection and assembly of a human genome" BioRxiv, doi.org / 10.1101 / 519025; January 13, 2019; which is hereby incorporated by reference in its entirety for all purposes).]

[0120] It will be apparent to those of ordinary skill in the relevant art that other suitable modifications and adaptations can be made to the methods and compositions described herein without departing from the scope of the invention or any of its embodiments. Having now described the invention in detail, the invention will be more clearly understood by reference to the following examples, which are included herein for illustrative purposes only and are not intended to limit the invention.

[0121] Although the above invention has been described in some detail for purposes of clarity and understanding, it will be apparent to those skilled in the art that various changes in form and detail may be made without departing from the true scope of the invention. For example, all of the above techniques and devices may be used in various combinations. All publications, patents, patent applications, and / or other documents cited in this application are hereby incorporated by reference in their entirety for all purposes to the extent as if each individual publication, patent, patent application, and / or other document were specifically and individually indicated to be incorporated by reference for all purposes.

Claims

1. A method of assembling a genome or a genomic region, the method comprising: obtaining a plurality of sequence reads of genomic fragments from a genome of interest, wherein each sequence read comprises a base call sequence; generating a homopolymer collapse sequence (HCS) and a corresponding homopolymer encoding sequence (HES) for each sequence read of the plurality of sequence reads, thereby generating a plurality of HCS reads and a plurality of HES reads; generating a plurality of suffix / prefix exact string matches of the plurality of HCS reads, wherein the length of the exact string match is equal to or greater than a minimum length; generating a plurality of trimmed HCS reads by removing any nucleotides from each of the plurality of HCS reads that are not part of a suffix / prefix exact string match with another HCS read; generating a first directed overlap graph from the trimmed HCS reads; identifying connected components in the first directed overlap graph; generating a multiple sequence alignment for each connected component, thereby generating one or more multiple sequence alignments, wherein positions in each trimmed HCS read are labeled with consecutive integer values such that the same integer value is assigned to aligned positions in any two trimmed HCS reads; pruning merged nodes from the first directed overlap graph based on one or more multiple sequence alignments of the trimmed HCS reads; generating a homopolymer collapsed consensus sequence by joining base calls at each aligned position in a first multiple sequence alignment of one or more multiple sequence alignments of the trimmed HCS reads; associating a vector of homopolymer lengths with each position in the homopolymer collapsed consensus sequence, wherein: (i) the number of elements in the vector is the number of trimmed HCS reads covering that position in the multiple sequence alignment, and (ii) each component of the vector is the length of the homopolymer in the corresponding homopolymer encoding sequence at that position; assigning a consensus homopolymer length as the floor of the median of the components of the vector of homopolymer lengths associated with that position to each position in the homopolymer collapsed consensus sequence; and replacing each position in the homopolymer collapsed consensus sequence with a homopolymer string formed by N consecutive nucleotide copies, where N is the assigned consensus homopolymer length calculated for that position, to generate a homopolymer extended consensus sequence, thereby assembling a genome or a genomic region.

2. The method according to claim 1, wherein the minimum length of the overlapping region is from 0.5 kb to 10 kb.

3. The method according to claim 2, wherein the minimum length of the overlapping region is from 5 kb to 8 kb.

4. The method according to claim 3, wherein the minimum length of the overlapping region is from 6 kb to 7 kb.

5. The method according to claim 1, wherein the minimum length is at least half the length of the average length of the HCS reads.

6. The method according to any one of claims 1 to 5, wherein the plurality of sequence reads are generated in a single molecule synthesis sequencing reaction.

7. The method according to claim 6, wherein the single molecule synthesis sequencing reaction is a single molecule real-time (SMRT) sequencing reaction.

8. The method according to any one of claims 1 to 5, wherein the plurality of sequence reads are generated in a single molecule nanopore sequencing reaction.

9. The method according to any one of claims 1 to 5, wherein the plurality of sequence reads are a plurality of single molecule consensus sequences (SMCS).

10. The method according to claim 9, wherein the SMCS are generated from at least 8 subreads.

11. The method according to claim 10, wherein the subreads are generated from a tandem polynucleotide substrate in a single molecule sequencing reaction.

12. The method according to claim 11, wherein the subreads are generated in a single molecule sequencing by synthesis reaction.

13. The method according to claim 11, wherein the subreads are generated in a single molecule nanopore-based sequencing reaction.

14. The method according to claim 10, wherein the subreads are generated from a circular or topologically circular polynucleotide substrate in a single molecule sequencing by synthesis reaction.

15. The method according to any one of claims 1 to 5, wherein the genome of interest is the human genome.

16. The method according to any one of claims 1 to 5, wherein the genomic sample comprises a plurality of different genomes, and the method further comprises generating assemblies for the plurality of different genomes.

17. The method according to claim 16, wherein the sample is a metagenomic sample comprising a plurality of microbial genomes.

18. The method according to any one of claims 1 to 5, wherein HCS not placed into the connection assembly is placed into a holding bin for validating variant calls in the assembly.

19. The method according to any one of claims 1 to 5, wherein prior to generating the HCS, a plurality of sequence reads are preselected to map to one or more genomic regions of interest.

20. The method according to claim 19, wherein the preselection mapping is performed by a low-stringency sequence similarity search.

21. The method according to claim 19, wherein the one or more genomic regions of interest comprise first and second genomic loci having high sequence similarity to each other.

22. The method according to claim 21, wherein separate consensus sequences are generated for the first and second genomic loci.

23. The method according to claim 19, wherein the one or more genomic regions of interest comprise a genomic locus having a highly repetitive region.

24. The method according to any one of claims 1 to 5, wherein the method is a method for de novo genome assembly.

25. The method according to claim 24, wherein the de novo genome assembly is a complete or partial haplotype-resolved assembly of a polyploid genome.

26. A system for assembling a genome or genomic region, comprising: a memory; an input / output; and a processor coupled to the memory, wherein the system is configured to: receive a plurality of sequence reads of genomic fragments from a genome of interest, wherein each sequence read comprises a base call sequence; Generate a homopolymer collapse sequence (HCS) and a corresponding homopolymer encoding sequence (HES) for each of the plurality of sequence reads, thereby generating a plurality of HCS reads and a plurality of HES reads; Generate a plurality of suffix / prefix exact string matches of the plurality of HCS reads, wherein the length of the exact string match is equal to or greater than a minimum length; Generate a plurality of trimmed HCS reads by removing any nucleotides in each of the plurality of HCS reads that are not part of a suffix / prefix exact string match with another HCS read; Generate a first directed overlap graph from the plurality of trimmed HCS reads; Identify the connected components in the first directed overlap graph; Generate a multiple sequence alignment for each connected component, thereby generating one or more multiple sequence alignments, wherein the positions in each trimmed HCS read are labeled with consecutive integer values such that the same integer value is assigned to aligned positions in any two trimmed HCS reads; Prune merged nodes from the first directed overlap graph based on the multiple sequence alignment; Generate a homopolymer collapse consensus sequence by concatenating base calls at each alignment position in a first multiple sequence alignment among one or more multiple sequence alignments of the trimmed HCS reads; Associate a vector of homopolymer lengths with each position in the homopolymer collapse consensus sequence, wherein: (i) the number of elements in the vector is the number of trimmed HCS reads that cover that position in the first multiple sequence alignment, and (ii) each component of the vector is the length of the homopolymer in the corresponding homopolymer encoding sequence at that position; Assign a consensus homopolymer length as the floor of the median of the components of the vector of homopolymer lengths associated with that position; Replace each position in the homopolymer collapse consensus sequence with a homopolymer string formed by N consecutive nucleotide copies, where N is the assigned consensus homopolymer length calculated for that position, to generate a homopolymer extended consensus sequence; and Provide the homopolymer extended consensus sequence to a user to assemble a genome or genomic region.

Images