The invention generally relates to nucleic acid sequence analysis and more specifically to the assembling of nucleic acid sequence information from a collection of short read nucleic acid subsequences.
Recent advances in sequencing technology have made possible the rapid, high-throughput and cost-effective sequencing of genomic samples. In particular, next-generation sequencing technologies have resulted in increased accuracy and a significant increase in information content. See, e.g., U.S. Pat. No. 7,282,337; U.S. Pat. No. 7,279,563; U.S. Pat. No. 7,226,720; U.S. Pat. No. 7,220,549; U.S. Pat. No. 7,169,560; U.S. Pat. No. 6,818,395; U.S. Pat. No. 6,911,345; US Pub. Nos. 2006/0252077; 2007/0070349; and 2007-0070349. These automated methods and apparatus provide for high speed and high throughput analysis of long polynucleotide sequences with simplicity, flexibility and lower cost. See, e.g., www.helicosbio.com/, particularly information on HeliScope™ Sequencer.
The most promising next-generation sequencing technologies are based upon either sequencing-by-synthesis, which utilizes the natural ability of a polymerase enzyme to incorporate a nucleotide into a primer strand in a template-dependent manner, or sequencing-by-ligation, which utilizes the natural ability of a ligase enzyme to join two fragments when correctly aligned in a template-dependent manner. Single molecule sequencing technologies provide the additional benefit of allowing detection of single nucleotide incorporation in an individual surface-bound duplex. The output of these technologies is millions of short reads, generally 15 to 100 bases in length.
One of the challenges for all next-generation sequencing technologies is to find data processing methods that allow improved sequence detection and reduced error rate.
The invention is based, in part, on the unexpected discovery that multiple short subsequences can be efficiently assembled to obtain the sequence information of a longer target nucleic acid sequence from which the short sequences (or short reads) are segments. The present invention provides methods for improving the processing of sequencing data to infer the sequence of a nucleic acid molecule that is much longer than the effective read length.
A major advantage afforded by next generation sequencing technologies is high throughput production of sequence data. However, next-generation technologies generally produce shorter sequence read lengths, e.g., less than about 100 bases in length, compared to conventional sequencing methodologies, e.g., greater than about 500 to about 1000 bases in length. The present invention provides methods for leveraging the large number of short sequence reads generated by high-throughput next generation sequencing technologies to produce longer and more accurate consensus contig reads.
Assembling short DNA or RNA sequences into longer, more accurate consensus sequences is a major challenge facing current sequencing technologies. Methods for constructing these longer consensus sequences are provided herein. These methods rely on the construction of a multi-length sequence index, statistical probability value, that can be conceptualized as a de Bruijn graph in which sequence subsequences of length (n) to (m) and the nodes in the graph are connected to each other through subsequences that are of length (n+1) to (m+1). See
In one aspect, the invention is generally related to a method for constructing a target nucleic acid sequence. The method includes: a) obtaining the sequence information of a plurality of subsequences of the target nucleic acid sequence, wherein the plurality of subsequences are segments of and together form at least substantially the complete sequence of the target nucleic acid; b) selecting an initial subsequence from the plurality of subsequences and an end base thereof and analyzing the sequence information of the plurality of subsequences to obtain a statistical probability value for the base position next to the selected end base of the initial subsequence; c) analyzing the sequence information of the plurality of subsequences to obtain a statistical probability value for the base position next to the analyzed base position in b); and d) continue to repeat c) for the next positions to construct substantially the full sequence of the target nucleic acid.
In some preferred embodiments, in steps b)-d) analyzing the sequence information of the plurality of subsequences to obtain a statistical probability value for the next base position is through a multi-scale de Bruijn graph construct. The de Bruijn graph process may utilize a single weighted matrix and/or a multiple weighted matrix. In some preferred embodiments, the subsequences are sequences having 150 or fewer base pairs, 100 or fewer base pairs, or 50 or fewer base pairs. In some preferred embodiments, the subsequences are sequences having 10 or more base pairs. In some embodiments, the target nucleic acid sequence(s) are 1,000 base pairs or longer.
In another aspect, the invention generally relates to a method for assembling the sequence of a target nucleic acid having known subsequences. The method includes: a) selecting an initial subsequence from the known subsequences and an end base thereof and analyzing the sequence information of the known subsequences to obtain a statistical probability value for the base position next to the selected end base of the initial subsequence; b) analyzing the sequence information of the known subsequences to obtain a statistical probability value for the base position next to the analyzed base position in a); and c) continue to repeat b) for the next base positions to construct the full sequence of the target nucleic acid.
In some preferred embodiments, in steps b)-c) analyzing the sequence information of the plurality of subsequences to obtain a statistical probability value for the next base position is conducted through a multi-scale de Bruijn graph construct. The de Bruijn graph process may utilize a single weighted matrix and/or a multiple weighted matrix.
In yet another aspect, the invention generally relates to a method for sequencing a target nucleic acid. The method includes: a) sequencing a plurality of subsequences of the target nucleic acid, wherein the plurality of subsequences are segments of and together form the complete sequence of the target nucleic acid sequence; and b) assembling the subsequences via a de Bruijn graph process.
The foregoing aspects and embodiments of the invention may be more fully understood by reference to the following figures, detailed description and claims.
The invention may be further understood from the following figures in which:
In general, the invention relates to methods for obtaining sequence information from a plurality of short subsequences (short reads obtained from sequencing runs). Many high-throughput sequencing technologies produce sequence read lengths that are much smaller than the genomic region of interest. For example, read lengths in many of these technologies are between about 15 base pairs and about 100 base pairs on average.
Methods described herein allow the assembly of short reads into a longer assembled sequence. In one embodiment, these methods may employ the de Bruijn graph approach to assemble short read sequence data into longer sequences. See, e.g., de Bruijn, N. G. (1946) “A Combinatorial Problem” Koninklijke Nederlandse Akademie v. Wetenschappen 49: 758-764; Flye Sainte-Marie, C. (1894) “Question 48” L'Intermédiaire Math. 1: 107-110; Good, I. J. (1946) “Normal Recurring Decimals” Journal of the London Mathematical Society 21 (3): 167-169; Zhang, et al., (1987) “On the de Bruijn-Good Graphs” Acta Math. Sinica 30 (2): 195-205.
In combinatorial mathematics, a k-ary De Bruijn sequence B(k, n) of order n is a cyclic sequence of a given alphabet A with size k for which every possible subsequence of length n in A appears as a sequence of consecutive characters exactly once. Such a sequence has the following properties:
Each B(k, n) has length kn
There are k!k(n−1)/kn distinct de Bruijn sequences B(k, n).
For example, taking A={0, 1}, there are two distinct B(2, 3): 00010111 and 11101000, one being the reverse of the other. Two of the 2048 possible B(2, 5) in the same alphabetare0000010001100101001110101101111l and 0000101001000111110111001101011.
The de Bruijn sequences can be constructed by taking a Hamiltonian path of an n-dimensional de Bruijn graph over k symbols (or equivalently, a Eulerian cycle of a (n−1)-dimensional de Bruijn graph), or via finite fields. Every four-digit sequence occurs exactly once if one traverses every edge exactly once and returns to one's starting point.
Each edge in this 3-dimensional de Bruijn graph corresponds to a sequence of four digits: the three digits that label the vertex that the edge is leaving followed by the one that labels the edge. If one traverses the edge labeled 1 from 000, one arrives at 001, thereby indicating the presence of the subsequence 0001 in the de Bruijn sequence. To traverse each edge exactly once is to use each of the 16 four-digit sequences exactly once.
For example, following the Eulerian path:
The present invention may employ nodes and edges of varying lengths.
Often samples that are subjected to DNA or RNA sequencing are comprised of many different samples. The multi-scale de Bruijn graph approach can also be used to identify sequence contexts that are multiplicatively present. For example, when constructing a consensus sequence, all possible paths from each node are followed and the resulting sequences are saved. The point at which these paths all converge represents sequence variations present within the sequenced sample. Further, the abundance of each of these variants may be correlated with the cumulative weight of each of the paths.
In one aspect, the invention is generally related to a method for constructing a target nucleic acid sequence. The method includes: a) obtaining the sequence information of a plurality of subsequences of the target nucleic acid sequence, wherein the plurality of subsequences are segments of and together form at least substantially the complete sequence of the target nucleic acid; b) selecting an initial subsequence from the plurality of subsequences and an end base thereof and analyzing the sequence information of the plurality of subsequences to obtain a statistical probability value for the base position next to the selected end base of the initial subsequence; c) analyzing the sequence information of the plurality of subsequences to obtain a statistical probability value for the base position next to the analyzed base position in b); and d) continue to repeat c) for the next positions to construct the full sequence of the target nucleic acid.
In some preferred embodiments, in steps b)-d) analyzing the sequence information of the plurality of subsequences to obtain a statistical probability value for the next base position is through a multi-scale de Bruijn graph construct. The de Bruijn graph process may utilize a single weighted matrix and/or a multiple weighted matrix.
In some preferred embodiments, the sequence information for the plurality of subsequence is obtained using a sequencing-by-synthesis process, such as a single molecule sequencing-by-synthesis process.
In some preferred embodiments, the sequence information for the plurality of subsequence is obtained using a sequencing-by-ligation process.
In some embodiments, the method of this invention is used to construct a second target nucleic acid, e.g., simultaneously or sequentially, or a third or more target nucleic acid.
The target nucleic sequence(s) may originate from a sample obtained from a single subject or from more than one subject.
In some preferred embodiments, the subsequences are sequences having 50 or fewer base pairs, 35 or fewer base pairs, 25 or fewer base pairs, or 20 or fewer base pairs. In some preferred embodiments, the subsequences are sequences having 10 or more base pairs, 15 or more base pairs, 20 or more base pairs, or 25 or more base pairs.
In some embodiments, the target nucleic acid sequence(s) are 250 base pairs or longer, 500 base pairs or longer, 1,000 base pairs or longer, 5,000 base pairs or longer, 10,000 base pairs or longer, or 50,000 base pairs or longer.
In another aspect, the invention generally relates to a method for assembling the sequence of a target nucleic acid having known subsequences. The method includes: a) selecting an initial subsequence from the known subsequences and an end base thereof and analyzing the sequence information of the known subsequences to obtain a statistical probability value for the base position next to the selected end base of the initial subsequence; b) analyzing the sequence information of the known subsequences to obtain a statistical probability value for the base position next to the base position in a); and c) continue to repeat b) for the next base positions to construct the full sequence of the target nucleic acid.
In some preferred embodiments, in steps b)-c) analyzing the sequence information of the plurality of subsequences to obtain a statistical probability value for the next base position is through a multi-scale de Bruijn graph construct. The de Bruijn graph process may utilize a single weighted matrix and/or a multiple weighted matrix.
In yet another aspect, A method for sequencing a target nucleic acid. The method includes: a) sequencing a plurality of subsequences of the target nucleic acid, wherein the plurality of subsequences are segments of and together form the complete sequence of the target nucleic acid sequence; and b) assembling the subsequences via a de Bruijn graph process.
References and citations to other documents, such as patents, patent applications, patent publications, journals, books, papers, web contents, have been made throughout this disclosure. All such documents are hereby incorporated herein by reference in their entirety for all purposes.
The representative examples which follow are intended to help illustrate the invention, and are not intended to, nor should they be construed to, limit the scope of the invention. Indeed, various modifications of the invention and many further embodiments thereof, in addition to those shown and described herein, will become apparent to those skilled in the art from the full contents of this document, including the examples which follow and the references to the scientific and patent literature cited herein. The following examples contain important additional information, exemplification and guidance which can be adapted to the practice of this invention in its various embodiments and equivalents thereof.
Examples of certain embodiments of the invention may be found in
Epoxide-coated glass slides are prepared for oligo attachment. Epoxide-functionalized 40 mm diameter #1.5 glass cover slips (slides) are obtained from Erie Scientific (Salem, N.H.). The slides are preconditioned by soaking in 3×SSC for 15 minutes at 37° C. Next, a 500-pM aliquot of 5′ aminated capture oligonucleotide (oligo dT(50)) is incubated with each slide for 30 minutes at room temperature in a volume of 80 ml. The slides are then treated with phosphate (1 M) for 4 hours at room temperature in order to passivate the surface. Slides are then stored in 20 mM Tris, 100 mM NaCl, 0.001% Triton® X-100, pH 8.0 at 4° C. until they are used for sequencing.
For the illustration of the sequencing process, see, e.g., U.S. patent application Ser. Nos. 12/043,033 (Xie et al. filed Mar. 5, 2008) and 12/113,501 (Xie et al. filed May 1, 2008) (e.g.,
Further, Virtual Terminator™ nucleotide analogs of cytosine triphosphate, guanidine triphosphate, adenine triphosphate, and uracil triphosphate, each having a cleavable cyanine-5 label (at the 7-deaza position for ATP and GTP and at the C5 position for CTP and UTP, see, e.g., U.S. patent application Ser. Nos. 11/803,339 (Siddiqi et al. filed May 14, 2007) and 11/603,945 (Siddiqi et al. filed Nov. 22, 2006), are stored separately in the buffer containing 20 mM Tris-HCl, pH 8.8, 50 μM MnSO4, 10 mM (NH4)2SO4, 10 mM HCl, and 0.1% Triton X-100, and 50 U Kienow exo-polymerase (NEB).
Sequencing proceeds as follows. First, initial imaging is used to determine the positions of duplex on the epoxide surface. The Cy3 label attached to the nucleic acid template fragments is imaged by excitation using a laser tuned to 532 nm radiation (Verdi V-2 Laser, Coherent, Santa Clara, Calif.) in order to establish duplex position. For each slide only single fluorescent molecules that are imaged in this step are counted. Imaging of incorporated nucleotides as described below is accomplished by excitation of a cyanine-5 dye using a 635-nm radiation laser (Coherent). 100 nM Cy5-dCTP is placed into the flow cell and exposed to the slide for 2 minutes. After incubation, the slide is rinsed in 1×SSC/15 mM HEPES/0.1% SDS/pH 7.0 (“SSC/HEPES/SDS”) (15 times in 60 μl volumes each, followed by 150 mM HEPES/150 mM NaCl/pH 7.0 (“HEPES/NaCl”) (10 times at 60 μl volumes). An oxygen scavenger containing 30% acetonitrile and scavenger buffer (134 μl 150 mM HEPES/100 mMNaCl, 24 μl 100 mM Trolox in 150 mM MES, pH 6.1, 10 μl 100 mM DABCO in 150 mM MES, pH 6.1, 8 μl 2M glucose, 20 μl 50 mM Nal, and 4 μl glucose oxidase (USB) is next added. The slide is then imaged (100 frames) for 2 seconds using an Inova 301K laser (Coherent) at 647 nm, followed by green imaging with a Verdi V-2 laser (Coherent) at 532 nm for 2 seconds to confirm duplex position. The positions having detectable fluorescence are recorded. After imaging, the flow cell is rinsed 5 times each with SSC/HEPES/SDS (60 μl) and HEPES/NaCl (60 μl). Next, the cyanine-5 label is cleaved off incorporated dCTP by introduction into the flow cell of 50 mM TCEP/250mM Tris, pH 7.6/100 mM NaCl for 5 minutes, after which the flow cell is rinsed 5 times each with SSC/HEPES/SDS (60 μl) and HEPES/NaCl (60 μl). The remaining nucleotide is capped with 50 mM iodoacetamide/100 mM Tris, pH 9.0/100 mM NaCl for 5 minutes followed by rinsing 5 times each with SSC/HEPES/SDS (60 μl) and HEPES/NaCl (60 μl). The scavenger is applied again in the manner described above, and the slide is again imaged to determine the effectiveness of the cleave/cap steps and to identify non-incorporated fluorescent objects.
The procedure described above is then conducted 100 nM Cy5-dATP, followed by 100 nM Cy5-dGTP, and finally 100 nM Cy5-dUTP. Uridine may be used instead of Thymidine due to the fact that the Cy5 label is incorporated at the position normally occupied by the methyl group in Thymidine triphosphate, thus turning the dTTP into dUTP. The procedure (expose to nucleotide, polymerase, rinse, scavenger, image, rinse, cleave, rinse, cap, rinse, scavenger, final image) is repeated for a total of 80-120 cycles.
Once the desired number of cycles is completed, the image stack data (i.e., the single-molecule sequences obtained from the various surface-bound duplexes) are aligned to produce the individual sequence reads. The individual single molecule sequence read lengths obtained range from 2 to 50+ consecutive nucleotides. Only the individual single molecule sequence read lengths above some predetermined cut-off depending upon the nature of the sample, e.g. greater than 20 and above, are analyzed using the method of the invention.