Method for sequencing polynucleotides

This invention relates to computational methods in molecular biology, and more specifically to methods for determining the sequence of a polynucleotide.

Baines, W., and Smith, G. C., J Theor. Biology, 135:303-307 (1988).

Ben-Dor, A., Pe'er, I., Shamir, R., and Sharan, R., In Proceedings of the Tenth International Conference on Combinatorial Pattern Matching (CPM '99), 88-100 (1999), New York: ACM Press.

Cargill, M., Altshuler, D., Ireland, J., Sklar, P., Ardlie, K., Patil, N., Lane, C. R., Lim, E. P., Kalyanaraman, N., Nemesh, J., Ziaugra, L, Friedland, L., Rolfe, A., Warrington, J., Lipshutz, R., Daley, G. Q., and Lander, E. S., Nature Genetics, 22:231-238 (1999).

Drmanac, R., and Crkvenjakov, R., Yugoslav Patent Application 570 (1987).

Durbin, R., Eddy, S., Krogh, A., and Mitchison, G., Biological Sequence Analysis: Probabilistic Models of proteins and Nucleic Acids, Cambridge University Press, (1998).

Eddy, S. R., Current Opinions in Structural Biology, 6(3):361-365 (1996).

Hirschberg, D. S., Communications of the ACM, 18,6:341-343 (1975).

Jukes, T. H., and Cantor, C. R., Mammalian Protein Metabolism, New York: Academic Press 21-123 (1969).

Khrapko, K. R., Lysov, Y. P., Khorlyn, A. A., Shick, V. V., Florentiev, V. L., and Mirzabekov, A. D., FEBS Letters, 256:118-122 (1989).

Kimura, M., Journal of Molecular Evolution, 16:111-120 (1980).

Krogh, A., Brown, M. Mian, S., Sjolander, M., and Haussler, D., Applications to protein modeling. Technical Report UCSC-CRL-93-32, Department of Computer and Information Sciences, University of California at Santa Cruz (1993).

Krogh, A., Brown, M., Mian S., Sjolander, M. and Haussler, D., Appliations to protein modeling 235(5):1501-1531 (1994).

Lysov, Y, Floretiev, V., Khorlyn, A., Khrapko, K., Shick, V, and Mirzabekov, A., Dokl, Acad Sci., USSR, 303:1508-1511 (1988).

Macevices, S.C., International Patent Application PS US89 04741 (1989). National Center for Biotechnology Information, 2000, A database of single nucleotide polymorphisms, World Wide Web address: ncbi.nlm.nih.gov/SNP./

Pevzner, P. A., and Lipshutz, R. J., Mathematical Foundations of Computer Science, LNCS 841:143-158 (1994).

Pevzner, P. A., Lysov, Y. P., Khrapko, K. R., Belyavsky, A. V., Florentiev, V. L., and Mirzabekov, A. D., J. Biomol. Struct. Dyn. 7:63-73 (1989).

Preparata, F., Frieze, A., and Upfal, E., Journal of Computational Biology 6(3-4):361-368 (1999).

Skiena, S. S., and Sundaram, G., J. Comput. Biol. 2:333-353 (1995).

Smith, T. F., and Waterman, M. S., Journal of Molecular Biology 147(1):195-197 (1981).

Southern, E. M., Maskos, U., and Elder, J. K., Genomics 13:1008-1017 (1992).

Southern E., UK patent Application GB 8,810,400 (1988).

Southern, E. M., Trends in Genetics 12:110-115 (1996).

Wang, D. G., Fan, J., Siao, C., Berno, A., Young P., Sapolsky, R., Ghandour, G., Perkins, N., Winchester, E., Spencer, J., Kruglyak, L., Stein, L., Hsie, L., Topalogiou, T., Hubbell, E., Robinson, E., Mittrnann, M., Morris, M.S., Shen, N., Kilburn, D., Rioux, J., Nusbaum, C., Lipshutz, R., Chee, M., and Lander, E. S., Science 280:1077-1082.

Yang, Z., Molecular Biology and Evolution 10:1396-1401.

Sequencing by hybridization (SBH) is a method for sequencing a polynucleotide such as a DNA molecule (Bains & Smith 1988, Lysov et al. 1988, Southern 1988, Drmanac and Crkvenjakov 1987, Macevics 1989). In this method, a chip, or microarray is used consisting of a surface upon which all possible oligonucleotide probes of a particular length k (referred to herein as “k-mers”) are immobilized (Southern 1996). The DNA molecule whose sequence is to be determined, referred to as the “target molecule”, is allowed to hybridize to the k-mers on the chip. The target molecule and the k-mers on the chip may all be single stranded molecules. Alternatively, a double stranded target may first be cut into fragments having single stranded “sticky ends”, and the k-mers on the chip may be the sticky ends of double stranded molecules. Ideally, a single stranded target or the sticky end of a double stranded target hybridizes to a k-mer on the chip if and only if the sequence complementary to the k-mer occurs somewhere in the target sequence or the sticky end. Thus, in principle, it is possible to experimentally determine the “k-spectrum” of the target (the set of all k-long substrings present in the target). In practice, however, the data are ambiguous due to the ability of the target to bind to k-mers that are only partially complementary to one of its substrings. Thus, any binarization of the hybridization signal will contain errors.

The goal of SBH is to determine the target sequence from the target spectrum. However, even if the target spectrum were error free, the target sequence is not uniquely determined by the spectrum. If the number of sequences consistent with the spectrum is large, there is no satisfactory method to select the true sequence. Theoretical analysis and simulations (Southern et al., 1992, Pevzner and Lipshutz 1994) have shown that even when the spectrum is errorless and the correct multiplicity of each k-mer in the target sequence is known, the average length of a uniquely reconstructible target sequence using a chip of 8-mers is only about two hundred nucleotides, far below the length of a DNA molecule that may be sequenced by electrophoresis.

Let Σ=(A,C,G,T) designate the set of nucleotides composing a DNA molecule. M=4 is the “alphabet size”. A DNA sequence is a string over Σ which is denoted herein between braces (< >). The k-spectrum of a target sequence T of length L, T=<t 1 , t 2 , . t L >, is the set of all k-long substrings (k-mers) of T. For each k-mer {overscore (x)}=<x 1 , x 2 , . x k >in εΣ k , we define T({overscore (x)})to be 1 if {overscore (x)} is a substring of T, and 0 otherwise. We denote K=M k , the number of k-mers. A hybridization experiment measures, for each k-mer {overscore (x)} in εΣ k , an intensity of its hybridization with the target.

The result of an SBH experiment may be described by a graph in which each candidate target sequence is represented as a path in a graph (Pevzner et al., 1989). The graph is a directed de-Bruijn graph G(V,E) whose vertices are labeled by all the (k-1)-mers (the set of vertices V=Σ k−1 ), and its edges are labeled by k-mers, (the set of edges E=Σ k ). The edge labeled <x 1 , x 2 . . . X k >connects the vertex <x 1 , x 2 . . . x k−1 >to the vertex <x 2 . . . x k >. There is a 1—1 correspondence between L-long candidate target sequences and (L−k+1)- long paths in G, whose edge labels comprise the target spectrum. Hereafter, we interchangeably refer to edges and their labels, and also to sequences and their corresponding paths.

Since k-mers may reoccur in the target sequence, the paths do not have to be simple. When the spectrum is perfect and the multiplicities of the k-mers in the spectrum are known, every solution is an Eulerian path (Pevzner et al. 1989). In practice, however, the spectrum is not perfect and the multplicities are not known.

Alternative chip designs (Bains and Smith 1988, Khrapko et al 1989, Pevzner et al 1991, Preparata et al. 1999, Ben-Dor et al. 1999), as well as interactive protocols (Skiena and Sundararn 1995) havebeen suggested, often assuming additional information, in order to reduce the ambiguity of the hybridization-based reconstruction.

Nucleotide sequences from different sources may resemble each other, due to a common ancestral gene. This phenomenon is encountered within a species, between duplicated regions within a genome, and between individuals within a population. Small differences in sequences, referred to as “Single Nucleotide Polymorphisms” or SNPs, efficiently serve as genetic markers that are useful in medicine. Thus the detection and genotyping of SNPs has become an important task of human geneticists. The evolution of homologous sequences from a common ancestral gene is mainly due to nucleotide substitution. Insertions and deletions of nucleotides are also known to have occurred during evolution of homologous sequences, though at lower rates.

A DNA molecule having a known sequence and known to be homologous to a target molecule has not yet been used to reduce the ambiguity of SBH data in order to determine the target sequence.

In the following description and set of claims, two parameters are considered to be equivalent to each other if they are proportional to each other.

The present invention provides a method for sequencing a target sequence. In accordance with the invention, experimental spectrum data obtained from a DNA chip is combined with sequence information of a reference DNA molecule. The reference molecule is preferably a molecule believed to be homologous with the target. For example, the target sequence may be a mutant gene and the reference sequence the previously sequenced normal gene. As another example, the target sequence may be a human gene and the reference sequence the homologous gene in another organism. A score is defined for each sequence in a set of candidate target sequences based upon a simultaneous comparison of the candidate sequence with the spectrum and with the reference sequence. A candidate target sequence is then selected having a essentially maximal score. Calculating the score does not require knowledge of the multiplicities of the k-mers in the k-spectrum. Moreover, unlike all prior art algorithms, the score does not assume that the spectrum is perfect.

The invention therefore provides a novel probabilistic method that handles imperfect hybridization data with unknown multiplicities. Thus, in accordance with the invention the hybridization of the target T with the k-mer on the DNA chip complementary to {overscore (x)} is described by probabilities P 0 ({overscore (x)}) and P 1 ({overscore (x)}) of the observed hybridization signal when T({overscore (x)})=0, and T({overscore (x)})=1, respectively The results of a hybridization experiment are described by the “probabilistic spectrum” (PS ) defined as the pair (P 0 ,P 1 ) of functions P i : Σ k →[0, 1]. If the experiment were perfect, i.e., if P 0 ({overscore (x)}) and P 1 ({overscore (x)}) are either 0 or 1 with P 0 ({overscore (x)})+P 1 ({overscore (x)}) 1, then the PS would represent the k-spectrum. In practice, however, P 0 ({overscore (x)}) and P 1 ({overscore (x)})are both positive There is thus a chance 1−P 0 ({overscore (x)}) for a false positive (a k-mer ({overscore (x)}) not occurring in T, whose complementary sequence produces a hybridization signal indicative of hybridization) and a chance 1−P 1 ({overscore (x)}) for a false negative (a k-mer ({overscore (x)}) occurring in T, whose complementary sequence produces a signal indicative of no hybridization). (When handling probabilities, some of which are perfect, problems of division by zero might occur. This is avoided by implicitly perturbing probabilities 0 and 1 to ε and 1− ε .)

The probability of obtaining a specific spectrum PS when T is used as the target is referred to as the “experimental likelihood”. The experimental likelihood is calculated assuming that the hybridization results of the target to different k-mer probes are mutually independent. In one embodiment of the invention, an experimental likelihood L e ({circumflex over (T)}) is used that does not assume knowledge of the multiplicities of each k-mer in the sequence. L e ({circumflex over (T)}) is given by:

L e ⁡ ( T ^ ) = Prob ⁡ ( PS ❘ T ^ ) = ∏ x → ∈   ⁢ ∑ k ⁢   ⁢ P T ^ ⁡ ( x → ) ⁡ ( x → ) ( 1 )

Taking logarithms and defining

ω ⁡ ( x → ) = log ⁢ P 1 ⁡ ( x → ) P 0 ⁡ ( x → )

we can write:

log ⁢   ⁢ P T ^ ⁡ ( x → ) ⁡ ( x → ) = { log ⁢   ⁢ P 0 ⁡ ( x → ) if ⁢   ⁢ T ^ ⁡ ( x → ) = 0 log ⁢   ⁢ P 0 ⁡ ( x → ) + ω ⁡ ( x → ) if ⁢   ⁢ T ^ ⁡ ( x → ) = 1. ( 2a )

Hence,

log ⁢   ⁢ L e ⁡ ( T ^ ) = ∑ x → ∈ ∑     ⁢   ⁢ k   ⁢   ⁢ log ⁢   ⁢ P 0 ⁢   ⁢ ( x → ) + ∑ T ^ ⁡ ( x → ) = 1 ⁢   ⁢ ω ⁡ ( x → ) . ( 2b )

The first term is a constant (independent of TP), and is omitted hereafter.

In another embodiment, an approximate likelihood {tilde over (L)}({circumflex over (T)}) is used, that is defined as follows: Let p=e 0 . . . e L−k be the path in G corresponding to {circumflex over (T)} and define

log ⁢ L ~ e ⁡ ( T ^ ) = ∑ i = 0 L - k ⁢   ⁢ ω ⁡ ( e i ) . ( 3 )

{tilde over (L)} e =L e ({circumflex over (T)}) for a path in which all edges have a multiplicity of 1, and is otherwise an approximation to L e ({circumflex over (T)}). {tilde over (L)}({circumflex over (T)}) has the advantage of being easily computable in a recursive manner:

log{tilde over (L)} e ( e 0 , . . . e l )=log{tilde over (L)} e ( e 0 , . . . e l−1 )+ω( e l )  (4)

In yet another einbodinent, an experimental likelihood L e ({circumflex over (T)}) is used that takes into account the multiplicities of edges. In this case, the probabilistic spectrum consists of probabilities P i ({overscore (x)}), denoting the probability of the observed hybridization signal when the multiplicity of x in the target is i. L e ({circumflex over (T)}) is defined by:

L e _ ⁡ ( T ^ ) = Prob ⁡ ( PS ❘ T ^ ) = ∏ x → ∈   ⁢ ∑ k ⁢ P T ^ _ ⁡ ( x → ) ⁡ ( x → ) ( 4b )

where {circumflex over (T)} ({overscore (x)}) is the multiplicity of {overscore (x)} in {circumflex over (T)}.

Thus in its first aspect, the invention provides a method for obtaining a candidate sequence, the candidate nucleotide sequence being indicative of a sequence of a target polynucleotide molecule T, T producing a hybridization signal I({overscore (x)}) upon incubating T with a polynucleotide {overscore (x)} for each polynucleotide {overscore (x)} in a set E of polynucleotides, the method comprising the steps of:

(a) for each polynucleotide {overscore (x)} in the set E of polynucleotides, obtaining a probability P 0 ({overscore (x)}) of the hybridization signal I({overscore (x)}) when the sequence {overscore (x)} is not complementary to a subsequence of T and a probability P 1 ({overscore (x)}) of the hybridization signal when the sequence x is complementary to a subsequence of T; so as to obtain a probabilistic spectrum (PS) of T; (b) assigning a score to each of a plurality of candidate nucleotide sequences, the score being based upon the probabilistic spectrum and upon at least one reference nucleotide sequence H; and (c) selecting one or more candidate nucleotide sequences having an essentially maximal score.

In its second aspect, the invention provides a program storage device readable by machine, tangibly embodying a program of instructions executable by the machine to perform method steps for obtaining a candidate nucleotide sequence, the candidate nucleotide sequence being indicative of a sequence of a target polynucleotide molecule T, T producing a hybridization signal I({overscore (x)}) upon incubating T with a polynucleotide x for each polynucleotide {overscore (x)} in a set E of polynucleotides, the method comprising the steps of:

(a) for each polynucleotide x in the set E of polynucleotides, obtaining a probability P 0 ({overscore (x)}) of I({overscore (x)}) when the sequence {overscore (x)} is not complementary to a subsequence of T and a probability P 1 ({overscore (x)}) of I({overscore (x)}) when the sequence {overscore (x)} is complementary to a subsequence of T; so as to obtain a probabilistic spectrum (PS) of T; (b) assigning a score to each of a plurality of candidate nucleotide sequences, the score being based upon the probabilistic spectrum and upon at least one reference nucleotide sequence H; and (c) selecting a candidate nucleotide sequence having an essentially maximal score.

In its third aspect the invention provides a computer program product comprising a computer useable medium having computer readable program code embodied therein for obtaining a candidate nucleotide sequence, the candidate nucleotide sequence being indicative of a sequence of a target polynucleotide molecule T, T producing a hybridization signal l({overscore (x)}) upon incubating T with a polynucleotide {overscore (x)} for each polynucleotide {overscore (x)} in a set E of polynucleotides, the computer program product comprising:

(a) for each polynucleotide {overscore (x)} in the set E of polynucleotides, computer readable program code for causing the computer to obtain a probability P 0 ({overscore (x)}) of I({overscore (x)}) the sequence x is not complementary to a subsequence of T and a probability P 1 ({overscore (x)}) of I({overscore (x)}) when the sequence {overscore (x)} is complementary to a subsequence of T; (b) computer readable program code for causing the computer to assign a score to each of a plurality of candidate nucleotide sequences, the score being based upon the probabilistic spectrum and upon at least one reference nucleotide sequence H; and (c) computer readable program code for causing the computer to select a candidate nucleotide sequence having an essentially maximal score.

First Embodiment

In this embodiment, the unknown target sequence T=<t 1 . . . t 1 >has a known, homologous reference sequence H=<h 1 . . . h 1 >. H and T are known to differ from each other by nucleotide substitutions without insertions or deletions (indels). This would be the case, for instance, when the target T is a mutant sequence whose wild type sequence H has already been sequenced, and one expects that nucleotide substitutions are the only cause of variability between H and T (statistically, substitutions are much more prevalent than indels (Wang et al. 1998)). A set of M×M position specific substitution matrices M (1) , . . . , M (1) are used, where for each position j along the sequence:

M (j) [i, i ′]=Prob( t j =i|h j =i ′) (5) for nucleotides i and i′εΣ.

for nucleotides i and i′εΣ.

The matrices M (j) may be the same for all j, or may different for different positions j. The matrices M (j) are used to calculate a distribution on the space of possible target sequences. This “prior distribution for ungapped homology”, D U , is given, for each candidate target sequence T by:

D u ⁡ ( T ^ ) = Prob ⁡ ( T ^ ❘ H ) = ∏ j = 1 l ⁢   ⁢ M ( j ) ⁡ [ t j , h j ] ( 6 )

One may recursively compute:

D u ( t 1 . . . t j )=( t 1 . . . t j−1 ). M (j) [t j , h j ]  (7)

We denote L (j) [x,y]≡log M (j) [x,y] .

The probability of a candidate target sequence {circumflex over (T)}. given the probability spectrum PS and the reference sequence H is:

Prob ⁡ ( T ^ ❘ H , PS ) = Prob ⁡ ( H ) · Prob ⁡ ( T ^ ❘ H ) · Prob ⁡ ( PS ❘ H , T ^ ) Prob ⁡ ( H , PS ) ( 8 )

Given {circumflex over (T)}, the hybridization signal is independent of H:

P rob ( PS|H,{circumflex over (T)} )= P rob ( PS|{circumflex over (T)} )

Thus, omitting the constant

Prob ⁡ ( H ) Prob ⁡ ( H , PS )

we can write:

P rob ( {circumflex over (T)}|H, PS )≅ D u ( {circumflex over (T)} )· L e ( {circumflex over (T)} )  (9 a )

P rob ( {circumflex over (T)}|H,PS )≅ D u ( {circumflex over (T)} )·{tilde over (L)} e ( {circumflex over (T)} )  (9 b )

or

P rob ( {circumflex over (T)}|H,PS )≅ D u ( {circumflex over (T)} )· L e ( {circumflex over (T)} )  (9 c )

Taking logarithms, the following “ungapped scores” of a candidate target are obtained:

S core 1 u ( {circumflex over (T)} )=log L e ( {circumflex over (T)} )+log D u ( {circumflex over (T)} )  (10 a )

S core 2 u ( {circumflex over (T)} )=log {tilde over (L)} e ( {circumflex over (T)} )+log D u ({circumflex over (T)})  (10 b )

S core 3 u ( {circumflex over (T)} )=log L e ( {circumflex over (T)} )+log D u ( {circumflex over (T)} )  (10 c )

With Score u l , Score u 2 or Score u 3 , the higher the score of a sequence {circumflex over (T)}, the more likely it is to be the target sequence. The highest scoring candidate sequence may be determined by any method known in the art. In the search for the highest scoring candidate sequence, complexity is preferably reduced by deleting from the graph edges for which {tilde over (L)} e ({circumflex over (T)}) L e ({circumflex over (T)}) or L e ({circumflex over (T)}) is less than a predetermined constant. Isolated vertices corresponding to highly improbable (k−1)-mers, are also preferably deleted from the graph.

For example, using {tilde over (L)} e ({circumflex over (T)}), the search for a high scoring candidate sequence may be performed by the following algorithm referred to herein as “Algorithm A ”. In accordance with Algorithm A, for each vertec {right arrow over (y)}=(y 1 . . . y k−1 εΣ k−1 , and integer j=k−1, k, k+1, . . . , l , let S u [{right arrow over (y)},j] be the maximum score of a j-long sequence ending with {right arrow over (y)} aligned to h 1 . . . h j .

Initialize, for each {right arrow over (y)}:

S u ⁡ [ y → , k - 1 ] = ∑ j = 1 k - 1 ⁢ L ( j ) ⁡ [ y j , h j ] ( 11 )

Loop over j=k, . . . , l, and for each vertex {right arrow over (y)}=y 1 . . . y k−1 recursively update

S u ⁡ [ y → , j ] = L ( j ) ⁡ [ y k - 1 , h j ] + max e = ( z → , y → ) ∈ E ⁢ { S u ⁡ [ z , j - 1 ] + ω ⁡ ( e ) } ( 12 ⁢ a )

Finally, return:

MAX ⁢   ⁢ Score u = max y → ∈ V ⁢ S u ⁡ [ y → , l ] ( 12 ⁢ b )

A sequence T* attaining the maximal score is found from the matrix S u as is known in the art, for example, by saving trace-back pointers:

P ⁡ [ y → , j ] = arg ⁢   ⁢ max Z → = < Z 0 ⁢ Z 1 ⁢ … ⁢   ⁢ Z K - 1 > , E = ( Z → , y → ) ∈ c ⁢ { S u ⁡ [ z , j - 1 ] + ω ⁡ ( e ) } ( 13a ) MAXPtr = argmax y → ∈ V ⁢   ⁢ S u ⁡ [ y _ , l ] ( 13b )

The maximum-scoring path in the graph is then followed, by setting:

Z 1 =MAXPtr, and for all j =k, . . . ,l: Z j−1 =P[Z j ,j]. Denote Z j =<Z j 1 z j 2 . . . z j k−1 >.

The final result is the sequence of nucleotides <z k−1 1, z k−1 2, . . . z k−1 k−1 , z k k−1 , z k+1 k−1 , . . . z 1 k−1 >

The time complexity is O(lK), since the maximization in ( 12 a ),( 13 a ) is a maximum of only a constant number (four) of terms. Although the complexity is exponential in k, it is constant for a given microarray (currently feasible values re k=8 or 9). Moreover, the complexity scales linearly with the size of the hybridization experimental results, which are part of the input.

Space complexity requires a more elaborate analysis. When naively using this algorithm, it requires O(lK) memory space, which is quite high for current technology microarrays. We now detail how we can modify the algorithm to reduce space complexity.

Observe, that this algorithm consists of two computations: Computing the optimal score (equations (11),(12 a ) and (12 b )), and reconstructing the optimal sequence (equations (13 a ) and (13 b )). The first task, of computing the optimal score alone, is space-efficient: it can be accomplished using space which is linear in the (effective) size of the hybridization experimental data, that is, O([K]) space.

By following the paradigm of Hirschberg (Hirschberg 1975), for example, for linear-space pair-wise alignment, a version of the algorithm is obtained which requires only linear space. The reduced space complexity is traded for time complexity, which increases by an O(log 1) factor.

For each position j =l, l '1, . . . k, k−1, the score of the entire sequence is decomposed. The total score is represented as a sum of two expressions: the contribution of its (j−k+l)-prefix, which equals the score of this prefix computed by S u , plus the contribution of the corresponding suffix. Formally, for each vertex {overscore (y)} y 1 . . . y k−1 εV , let R u [{overscore (y)}j ] be the maximum contribution to the score of a (l−j+k−l )-long sequence beginning with {overscore (y)} aligned to h j−k+2 . . . h l .

Initialize, for each {overscore (y)}:

R u [{right arrow over (y)},l ]=0  (14)

Loop over j=l−1, 1−2, . . . k−1, and for each vertex {overscore (y)}=y 1 . . . y k−1 recursively update:

R u ⁡ [ y -> , j ] = max e = ( y -> , z -> ) ∈ E ⁢ { R u ⁡ [ z -> , j + 1 ] + ω ⁡ ( e ) + L ( j + 1 ) ⁡ [ z k - 1 , h j + 1 ] } ( 15 )

Observe that, for all k−l ≦j ≦l

MAX ⁢   ⁢ Score u = max y -> ∈ V ⁢ { S u ⁡ [ y -> , j ] + R u ⁡ [ y -> , j ] } ( 16 )

Equation (16) can be used to decompose the problem into two similar problems, of half its size. Recursively solving these. sub-problems gives a divide-and-conquer approach for finding the optimal sequence. The linear space algorithm is therefore as follows:

1. If the length l of the target is smaller than some constant C, for example, 25 nucleotides:

Solve the problem directly, according to the dynamic program of Equations (11), (12a), (12b), (13a) and (13b).

Otherwise,

2. Set

m = l + k - 1 2 .

3. For each j =k−l, . . . ,m:

Compute S u [{overscore (y)},j] (following equations (11) and (12a)) for all {right arrow over (y)}, re-using space.

4. For each j =l, l−1, . . . , m:

Compute R u [{overscore (y)},j] (following equations (14) and (15)) for all {overscore (y)}, re-using space.

5. Find

y -> m = arg ⁢   ⁢ max y -> ∈ V ⁢ { S u ⁡ [ y -> , m ] + R u ⁡ [ y -> , m ] ,

thereby computing: MAXScore u , by (16).

6. Recursively compute:  (a) The optimal sequence aligned to h l . . . h m ending with {overscore (Y)} m .  (b) The optimal sequence aligned to (h m . . . h l ) beginning with {overscore (Y)} m

Observe, that for each {overscore (y)}, j, the values of S u [{overscore (y)},j] and R u {overscore (y)}j, are computed a total of log l times. Thus the algorithm takes O([K]l log 1) time and O([K]) space,. using the effective spectrum.

Second Embodiment: substitutions and deletions

In this embodiment, the unknown target sequence T=t 1 . . . t 1 · differs from the reference H=(h, . h,), by substitutions and deletions only, without insertions.

Denote the probability of initiating a gap right before h j (aligning h j to space) is 2 α j . Similarly, β j is the logarithm of the probability for gap extension at h j . Also define β j =log(1−2 β j ) αj=log(1−2 α j ). To overcome boundary problems at the ends of the sequence, we extend the alphabet by including left and right space characters: {overscore (Σ)}=Σ∪, . We augment the reference sequence by the string k on its left and k on the right. We extend the substitution matrix by using probabilities that force alignment of each of and to itself. Formally, we define:

Σ k - 1 _ = Σ k - 1 ⁢ ⋃ { x -> ⁢   ⁢ z ->  ⁢ x -> = ⊳ j , z -> ∈ Σ k - 1 - j } ⋃ { z -> ⁢   ⁢ x ->  ⁢ z -> ∈ Σ j , x -> = ⊲ k - 1 - j } ( 17 )

We arbitrarily set ω({overscore (y)}) to 0 for each {overscore (y)}ε{overscore (Σ k−1 )} |Σ k−1 . Thus, the weighted de-Bruijn graph is naturally extended over {overscore (Σ k−1 )}, and so is [G]=([V], [E]), its effective subgraph. Hereafter, we use the notation [G] for the extended graph. As with the previous embodiment, in order to reduce complexity, edges for which {tilde over (L)} e ({overscore (T)}) or L e ({overscore (T)}) is less than E are preferably deleted from the graph. Isolated vertices corresponding to highly improbable (k−1)-mers, are also preferably deleted from the graph.

The search for a high scoring candidate sequence may be performed by the following algorithm referred to herein as “Algorithm B”. In accordance with Algorithm B, for each {overscore (y)}=y 1 . . . y k−1 ε[V], j=k=1, k, k+1, . . . , S d [{overscore (y)},j] is defined as the maximum score of aligning a sequence ending with {overscore (y)} to h l . . . h j where h j is aligned to a gap (and y −1 is aligned to some h 1 . . . h j ). Further T d [{overscore (x)},j] is defined as the maximum score of aligning a sequence ending with y 1 . . . y k−1 , to h l . . . h j where he aligned to y k−1 . Initialize, for each {overscore (y)}:

S d [{right arrow over (y)}, k− 1]=−∞;  (19)

T d ⁡ [ y -> , k - 1 ] = { 0 y -> = ⊳ k - 1 - ∞ otherwise ( 20 )

Loop over j =k, . . . , 1, and for each {right arrow over (y)}=y 1 , . . . y k−l ε, [V], recursively update:

S d [{right arrow over (y)}, j]= max { T d [{right arrow over (y)}, j− 1]

+αa j , S d [{right arrow over (y)}, j− 1]+β j   (21)

T d ⁡ [ y -> , j ] = L ( j ) ⁡ [ y k - 1 , h j ] + max e = ( z , -> ⁢   ⁢ y -> ) ∈ E ⁢ { ω ⁡ ( e ) + max ⁢ { T d ⁡ [ z -> , j - 1 ] + α ^ j S d ⁡ [ z -> , j - 1 ] + β ^ j } } ( 22 )

Finally, return:

MAX Score d =T d [ k−1 , l]   (23)

The complexity of this algorithm is still O(l[K]) and a linear space variant can be obtained, as described in the previous embodiment. A sequence T* attaining the maximal score is then formed from the matrix T d as is known in the art, for example, by saving trace-back pointers to follow the maximally scoring path in analogous manner to that described in the previous embodiment.

Third Embodiment: Substitutions, Deletions and Insertions.

In this embodiment, a target sequence is determined when the target is known to be obtained from the reference by substitutions, insertions and deletions. The algorithm is an extension of the dynamic programs of the previous embodiments.

Denote by T j the target prefix whose last nucleotide is aligned to h j in the reference sequence. Further denote by a j (respectively b j ) the log-probability of initiating (extending) an insertion in the target after T j , and define â j =1−a j , {circumflex over (b)} j =1=b j .

Consider the weighted graph (G,ω). Define the K×K matrix W as follows:

W ⁡ [ x -> , y -> ] = { 2 w ⁡ ( y -> ) The ⁢   ⁢ ( k - 1 ) - suffix ⁢   ⁢ of ⁢   ⁢ x ->   is ⁢   ⁢ the ⁢   ⁢ ( k - 1 ) - prefix ⁢   ⁢ of ⁢   ⁢ y -> 0 Otherwise ( 24 )

W ′ [{right arrow over (x)}, {right arrow over (y)}]is thus the probability of moving from {right arrow over (x )} to {right arrow over (y)} along i edges. The probability of an insertion of length i after T j is a j b i j {circumflex over (b)} j . Suppose that the prefix T j ends with {overscore (x)}. Then a j b 1 i=1 {circumflex over (b)} j W ′ [{right arrow over (x)},{right arrow over (y)}] is the probability of T j+1 ending with {right arrow over (y)} and being i nucleotides longer than T j . The matrix W′ governing the stochastic progression from T j to T j+1 is calculated as follows:

W ′ = ⁢ a ^ j ⁢ b j ⁢ W 2 ⁢ b ^ j + a j ⁢ b j 2 ⁢ W 3 ⁢ b ^ j ⁢   ⁢ … ⁢ ( 25 ) = ⁢ a ^ j ⁢ W + a j ⁢ b j ⁢ b ^ j ⁢ W 2 ⁢ ∑ i ≥ 2 ⁢ b j i - 2 ⁢ W i - 2 ⁢ ( 26 ) = ⁢ a ^ j ⁢ W + a j ⁢ b j ⁢ b ^ j ⁢ W 2 ⁡ ( I - b j ⁢ W ) - 1 ⁢ ( 27 )

A new weighted graph (G′,ω′) is now defined as follows. The vertex set of G is also the vertex set of G′. The edge set E′ of G′ is the set of all pairs {right arrow over (x)}, {right arrow over (y)} with W′[{right arrow over (x)},{right arrow over (y)}]>0. Each such edge e=[{overscore (x)}{overscore (,)}y] is associated with a weight w′(e)=log W′[{right arrow over (x)}, {right arrow over (y)}].

The search for a high scoring candidate sequence may be performed by the following algorithm referred to herein as “Algorithm C”. In accordance with Algorithm C, Algorithm B of the second embodiment is applied to (G′, ω′) instead of (G, ω).

In contrast to G, degrees in G′ are not bounded by 4. Therefore, computing each dynamic program cell has complexity O(K) in the worst case, with the total complexity of the algorithm being O(l|E′|). Again, considering only the effective size of the graph allows more efficient computation, taking O(l|[E′]|).

Fourth Embodiment: Substitutions, Deletions and Insertions.

In this embodiment, homology between nucleotide sequences is described by Hidden Markov Models (HMMs) using a set Q of Markov chain states with a predefined set of allowed transitions between them. For each level (position along the sequence) j=1, . . . , l Q , Q includes three states: M j (match), I j (insert), and D j (delete). M j and D j can be reached from the three (j−1) (th) level states. I j can be reached from the three (j)-(th) level states (including a self-loop). Transition probabilities are as described in previous sections, e.g., a j =ProbS(M j , I j ). Additionally, each insert or match state, q, induces a vector of emission probabilities M q , where M q [i ] is the probability that the target nucleotide is i. We denote L q [i]≡0 for q=D j , L q [i] ≡log M q [i] otherwise. We write lpb(X)≡log Prob(X) for short.

The search for a high scoring candidate sequence may be performed by the following algorithm referred to herein as “Algorithm D”. In accordance with Algorithm D, a three dimensional array S is defined, where for each qεQ,{right arrow over (y)}=y 1 . . . y k−1 ε[V]=r=k, . . . , L, S[q,{right arrow over (y)},r] is defined as the maximum score of an r-long sequence ending with (y 1 . . . y k−1 >, whose alignment to the profile ends in q. Thus, initialize:

S[q start. k−1 , k −1]=0  (28)

S[q,{right arrow over (y)},k− 1]=−∞for other values of {right arrow over (y)},q   (29)

Loop over r=k, . . . 1, and for each {right arrow over (y)}=y 1 . . . y k−1 ε[V] p r≦l Q , recursively update:

S ⁡ [ q , y -> , r ] = L q ⁡ [ y k - 1 ] + max ⁢   e = ( E , q ) ∈ E q ′  ⁢ q ′ ∣ ⁢   → q ⁢ { S ⁡ [ q ′ , z -> , r - 1 ] + lpb ( q ′ ∣ → q ) + ω ⁡ ( e ) } ( 30 )

Finally, return:

MAX ⁢   ⁢ Score = max l ⁢ { S ⁡ [ q end ⁢   ⁢ 1 ⁢ ⊲ k - 1 , l ] } ( 31 )

A sequence T * is maximal score is then found in a manner similar to that described in the previous embodiments.

This algorithm requires O(l Q ·[K]·L) time and space, where L is an upper bound on the size of the target sequence. As with the previous embodiments, the complexity of this algorithm can be reduced to O(l Q ·[K] ·L log L) time and O(l Q ·[K]) memory. Furthermore, one can consider the dynamic program as filling a l Q ×L matrix, with a [K]-long vector in each matrix cell. Since all values far from the main diagonal of this matrix should be negligible, preferably only values within a distance less than a predetermined constant R to the main diagonal are calculated, reducing the complexity to O(R(l Q +L) ·[K] ·log L) time and O(R(l Q +L)·[K]) space.

Fifth Embodiment: Summation over all paths

In this embodiment the graph nodes (HMM states and k-mers) that are most likely to be visited at a certain position along the target sequence are obtained. The “Forward-Backward” algorithm is used (see, e.g., Durbin et al., 1998) providing the likelihood summed over all paths entering a node, instead of the likelihood of the maximum path. The only change to the equation presented thus far is that max operators are changed into log-sum-of-exponents. More specifically, equations (12a), (12b), (15), (16), (20), (21), (29), and (30) are re-written, respectively, as follows:

S u ⁡ [ y ⇀ , j ] = L ( j ) ⁡ [ y k - 1 , h j ] + log ⁢ ∑ e = ( z ⇀ , y ⇀ ∈ E )   ⁢ exp ⁢   ⁢ ( S u ⁡ [ z , j - 1 ] + ω ⁡ ( e ) ) ( 12 ⁢ a ′ ) MAXScore u = log ⁢ ∑ y ⇀ ∈ V   ⁢ exp ⁡ ( S u ⁡ [ y ⇀ , l ] ) ( 12 ⁢ b ′ ) R u ⁡ [ y ⇀ , j ] = log ⁢ ∑ e = ( y ⇀ , z ⇀ ) ∈ E   ⁢ exp ⁡ ( R u ⁡ [ z ⇀ , j + 1 ] + ω ⁡ ( e ) + L ( j + 1 ) ⁡ [ z k - 1 , h j + 1 ] ) ( 15 ′ ) MAXScore u = log ⁢ ∑ y ⇀ ∈ V   ⁢ exp ⁡ ( S u ⁡ [ y ⇀ , j ] + R u ⁡ [ y ⇀ , j ] ) ( 16 ′ ) S d ⁡ [ y ⇀ , j ] = log ⁡ ( exp ⁡ ( T d ⁡ [ y ⇀ , j - 1 ] + α j ) + exp ⁡ ( S d ⁡ [ y ⇀ , j - 1 ] + β j ) ) ( 20 ′ ) T d ⁡ [ y ⇀ , j ] = L ( j ) ⁡ [ y k - 1 , h j ] + log ⁢ ∑ e = ( z ⇀ , y ⇀ ) ∈ E   ⁢ exp ⁡ ( ω ⁡ ( e ) ) + log ⁡ ( exp ⁡ ( T d ⁡ [ z ⇀ , j - 1 ] + α ^ j ) + exp ⁡ ( S d ⁡ [ z ⇀ , j - 1 ] + β ^ j ) ) ( 21 ′ ) S ⁡ [ q , y ⇀ , r ] = L q ⁡ [ y k - 1 ] + ∑ e = ( E , q ) ∈ E q ′ ❘ q ′ ↦ q   ⁢ exp ⁡ ( S ⁡ [ q ′ , z ⇀ , r - 1 ] + l ⁢   ⁢ p ⁢   ⁢ b ⁡ ( q ′ ↦ q ) + ω ⁡ ( e ) ) ( 29 ′ ) MAXScore = log ⁢ ∑ l   ⁢ exp ⁡ ( S ⁡ [ q end ⁢   ⁢ 1 ⁢ ⊲ k - 1 , l ] ) ( 30 ′ )

Sixth Embodiment: Enhancements

In this embodiment the exact likelihood calculated according to Equation 10a of several top-scoring candidates found using the approximated likelihood (Equation 10b) is calculated. These sequences are then re-ranked. This 2-phase filtering is more discriminative than approximated scoring, while still tractable using the formulae presented.

If the score of a dynamic programming cell is very low, that cell probably does not participate in the maximum solution. This allows discarding such cells, without risking loss of the true optimum. Computing time and space may thus be saved.

The invention may be used for simultaneously re-sequencing several short targets, instead of a single long sequence. This scenario arises when considering many exons of a single gene. The invention may also be generalized to deal with DNA chips that do not contain the set of all k-mers.

When the set of oligonucleotides on the microarray is not the set of all k-mers, a graph is constructed consisting, as vertices, instead of all the (k−1)-mers, all the prefixes and suffixes of oligonucleotides on the microarray. Edges in this graph connect two vertices if there is one base pair suffix (suffix) addition to one of them, that makes the other its proper suffix (prefix). The scoring mechanism remains the same. This also applies for oligonucleotides containing “gaps” or “universal bases” (Preparata et al., 1999).

The invention may be used also for sequencing polypeptides. Given a polypeptide chain homologous to a target, and given a collection of probabilities of occurrence of sub-chains along the target, our algorithms will find the optimal candidate target sequence.

The invention was implemented and tested on simulated data. Nucleotide substitutions were equiprobable and insertions and deletions were not allowed. As a reference sequence, prefixes of the gene-rich human mitochondrial sequence, (Accession Number J01415) were used. For each reference sequence, the following were generated:

1. A target sequence generated according to a prescribed probability q of substitution, defining the matrix M as 1-q on the diagonal and q/3 elsewhere. 2. An 8-spectrum for the target was generated using the probabilistic spectrum defined by P i ({right arrow over (x)})=1 −p if T({right arrow over (x)})=i, where p is a fixed probability.

All probabilistic parameters were constant, i.e., position/k-mer independent. For each 8-spectrum and target sequence, candidate sequences were scored using Eq. (10), and a candidate sequence of maximal score was found.

The algorithm was implemented in C++ and executed on Linux and SGI machines. Running times, on a Pentium 3, 600 MHz machine, were roughly 0.12l log l seconds for an l-long reference sequence (ranging from roughly 7 minutes for a 500 bp-long sequence to 2.5 hours for 6 Kb). Only the main memory was used, with the application consuming at most 40 Mb. The graph was not reduced to its effective size. This would have reduced both space and time dramatically, at the expense of possibly missing the truly maximal scoring sequence.

The performance of the algorithm was quantified by the following FIGS. of merit:

1. Full success rate —The fraction of runs for which the target sequence was perfectly reconstructed. 2. ε-success rate —The fraction of runs for which the target sequence length 1 was reconstructed with fewer than ε•1 nucleotide errors. 3. Average sequencing error —The fraction of nucleotide errors.

Table 1 presents results for a scenario of distinct, but closely related sequences, e.g., orthologous genes in a pair of primates. We assume perfect hybridization data with 97% sequence similarity (that is q=0.03). The results show that sequences of length up to 2000 can be reconstructed almost perfectly. The non-monotonicity of the figures of merit with respect to the target length is due to sequence contents.

Table 2 presents results for a scenario of SNP-genotyping. The rate of SNPs is assumed to be 1:700 (Wang et al. 1998), and p=2% was used. The results show that a high success rate is achievable even in the presence of spectrum errors.

Table 1: Results on simulated date, for different sequence lengths, assuming 97% sequence similarity between the target and the reference, and perfect hybridization data.

It will also be understood that the system according to the invention may be a suitably programmed computer. Likewise, the invention contemplates a computer program being readable by a computer for executing the method of the invention. The invention further contemplates a machine-readable memory tangibly embodying a program of instructions executable by the machine for executing the method of the invention.