VIEWS: 13 PAGES: 82 CATEGORY: Business POSTED ON: 9/24/2010 Public Domain
Modern Homology Search Ming Li Canada Research Chair in Bioinformatics Computer Science. U. Waterloo I want to show you how one simple theoretical idea can make big impact to a practical field. Outline Motivation, market, biology, homology search Review dynamic programming, BLAST heuristics Optimized spaced seeds The idea How to compute them Data models, coding regions, HMM, vector seeds Multiple optimized spaced seeds PatternHunter program Mathematical theory of spaced seeds Why they are better Complexity of finding optimal spaced seeds Applications beyond bioinformatics & Open problems 1. Introduction Biology Motivation and market Definition of homology search problem Biology DNA Protein (Genotype) Phenotype Human: 3 billion bases, 30k genes. T A A T E. coli: 5 million bases, 4k genes C G T A cDNA reverse transcription A T transcription translation G C C G mRNA Protein G (A,C,G,U) (20 amino acids) C Codon: three nucleotides encode A T an amino acid. 64 codons C G 20 amino acids, some w/more codes T A A T A gigantic gold mine The trend of genetic data growth 8 Nucleotides(billion) 7 6 5 30 billion 4 3 in year 2005 2 1 0 1980 1985 1990 1995 2000 Years 400 Eukaryote genome projects underway GenBank doubles every 18 months Comparative genomics all-against-all search Software must be scalable to large datasets. Homology search Given two DNA sequences, find all local similar regions, using “edit distance” (match=1, mismatch=-1, gapopen=-5, gapext=-1). Example. Input: E. coli genome: 5 million base pairs H. influenza genome: 1.8 million base pairs Output: all local alignments (PH demo) java –jar phn.jar –i ecoli.fna –j hinf.fna –o out.txt –b Homology search vs Google search Internet search Size limit: 5 billion people x homepage size Supercomputing power used: ½ million CPU-hours/day Query frequency: Google --- 112 million/day Query type: exact keyword search --- easy to do Homology search Size limit: 5 billion people x 3 billion basepairs + millions of species x billion bases 10% (?) of world’s supercomputing power Query frequency: NCBI BLAST -- 150,000/day, 15% increase/month Query type: approximate search Bioinformatics Companies living on Blast Paracel (Celera) TimeLogic (recently liquidated) TurboGenomics (TurboWorx) NSF, NIH, pharmaceuticals proudly support many supercomputing centers mainly for homology search Hardware high cost & become obsolete in 2 years. History Dynamic programming (1970-1980) Full sensitivity, but slow. Human vs mouse genomes: 104 CPU-years BLAST, FASTA heuristics (1980-1990) Low sensitivity, still not fast enough. Human vs mouse genomes: 19 CPU-years BLAST paper was referenced 100000 times Modern technology: full sensitivity & greater speed! 2. Old Technology Dynamic programming – full sensitivity homology search BLAST heuristics --- trading sensitivity with time. Dynamic Programming: Longest Common Sequence Alignment Subsequence (LCS). (Needleman-Wunsch, 1970) s(i,j)=max V=v1v2 … vn s(i-1,j) + d(vi,-) W=w1w2 … wm s(i,j-1)+d(-,wj) s(i,j) = length of LCS s(i-1,j-1)+d(vi,wj) of V[1..i] and W[1..j] 0 (Smith-Waterman, 1981) Dynamic Programming: * No affine gap penalties, s(i-1,j) where d, for proteins, is either s(i,j)=max s(i,j-1) PAM or BLOSUM matrix. For s(i-1,j-1)+1,vi=wj DNA, it can be: match 1, mismatch -1, gap -3 (In fact: open gap -5, extension -1.) d(-,x)=d(x,-)= -a, d(x,y)=-u. When a=0, u=infinity, it is LCS. Misc. Concerns Local sequence alignment, add s[i,j]=0. Gap penalties. For good reasons, we charge first gap cost a, and then subsequent continuous insertions b<a. Space efficient sequence alignment. Hirschberg alg. in O(n2) time, O(n) space. Multiple alignment of k sequences: O(nk) Blast Popular software, using heuristics. By Altschul, Gish, Miller, Myers, Lipman, 1990. E(xpected)-value: e= dmne -lS, here S is score, m is database length and n is query length. Meaning: e is number of hits one can “expect” to see just by chance when searching a database of size m. Basic Strategy: For all 3 a.a. (and closely related) triples, remember their locations in database. Given a query sequence S. For all triples in S, find their location in database, then extend as long as e-value significant. Similar strategy for DNA (7-11 nucleotides). Too slow in genome level. Blast Algorithm Find seeded (11 bases) matches Extent to HSP’s (High Scoring Pairs) Gapped Extension, dynamic programming Report all local alignments BLAST Algorithm Example: Find seeded matches of 11 base pairs Extend each match to right and left, until the scores drop too much, to form an alignment Report all local alignments Example: 0001 1 1 011 1 11 1 1 1 11 1001101 11 10 AGCGATGTCACGCGCCCGTATTTCCGTA G | | | | | | |x| | | | | | | | | | | || TCGGATCTCACGCGCCCGGCTTACCGTG BLAST Dilemma: If you want to speed up, have to use a longer seed. However, we now face a dilemma: increasing seed size speeds up, but loses sensitivity; decreasing seed size gains sensitivity, but loses speed. How do we increase sensitivity & speed simultaneously? For 20 years, many tried: suffix tree, better programming .. 3. Optimized Spaced Seeds Why are they better How to compute them Data models, coding regions, HMM PatternHunter Multiple spaced seeds Vector seeds New thinking Three lines of existing approaches: Smith-Waterman exhaustive dynamic programming. Blast family: find seed matches (11 for Blastn, 28 for MegaBlast), then extend. Dilemma: increasing seed size speeds up but loses sensitivity; descreasing seed size gains sensitivity but loses speed. Suffix tree: MUMmer, Quasar, REPuter. Only good for precise matches, highly similar sequences. Blast approach is the only way to deal with real and large problems, but we must to solve Blast dilemma: We need a way to improve sensitivity and speed simultaneously. Goals: Improve (a) Blast, (b) Smith-Waterman. Optimal Spaced Seed (Ma, Tromp, Li: Bioinformatics, 18:3, 2002, 440-445) Spaced Seed: nonconsecutive matches and optimized match positions. Represent BLAST seed by 11111111111 Spaced seed: 111*1**1*1**11*111 1 means a required match * means “don’t care” position This seemingly simple change makes a huge difference: significantly increases hit prob. to homologous region while reducing bad hits. Formalization Given i.i.d. sequence (homology region) with Pr(1)=p and Pr(0)=1-p for each bit: 1100111011101101011101101011111011101 111*1**1*1**11*111 Which seed is more likely to hit this region: BLAST seed: 11111111111 Spaced seed: 111*1**1*1**11*111 Sensitivity: PH weight 11 seed vs Blast 11 & 10 PH 2-hit sensitivity vs Blastn 11, 12 1-hit Expect Less, Get More Lemma: The expected number of hits of a weight W length M seed model within a length L region with similarity p is (L-M+1)pW Proof: The expected number of hits is the sum, over the L-M+1 possible positions of fitting the seed within the region, of the probability of W specific matches, the latter being pW. ■ Example: In a region of length 64 with 0.7 similarity, PH has probability of 0.466 to hit vs Blast 0.3, 50% increase. On the other hand, by above lemma, Blast expects 1.07 hits, while PH 0.93, 14% less. Why Is Spaced Seed Better? A wrong, but intuitive, proof: seed s, interval I, similarity p E(#hits) = Pr(s hits) E(#hits | s hits) Thus: Pr(s hits) = Lpw / E(#hits | s hits) For optimized spaced seed, E(#hits | s hits) 111*1**1*1**11*111 Non overlap Prob 111*1**1*1**11*111 6 p6 111*1**1*1**11*111 6 p6 111*1**1*1**11*111 6 p6 111*1**1*1**11*111 7 p7 ….. For spaced seed: the divisor is 1+p6+p6+p6+p7+ … For BLAST seed: the divisor is bigger: 1+ p + p2 + p3 + … Finding Optimal Spaced Seeds (Keich, Li, Ma, Tromp, Discrete Appl. Math. 138(2004), 253-263 ) Let f(i,b) be the probability that seed s hits the length i prefix of R that ends with b. Thus, if s matches b, then f(i,b) = 1, otherwise we have the recursive relationship: f(i,b)= (1-p)f(i-1,0b') + pf(i-1,1b') where b' is b deleting the last bit. Then the probability of s hitting R is, Σ|b|=M Prob(b) f(|R|,b). Choose s with the highest probability. Improvements Brejova-Brown-Vinar (HMM) and Buhler-Keich-Sun (Markov): The input sequence can be modeled by a (hidden) Markov process, instead of iid. Multiple seeds Brejova-Brown-Vinar: Vector seeds Csuros: Variable length seeds – e.g. shorter seeds for rare query words. PatternHunter (Ma, Tromp, Li: Bioinformatics, 18:3, 2002, 440-445) PH used optimal spaced seeds, novel usage of data structures: red-black tree, queues, stacks, hashtables, new gapped alignment algorithm. Written in Java. Yet many times faster than BLAST. Used in Mouse/RAT Genome Consortia (Nature, Dec. 5, 2002), as well as in hundreds of institutions and industry. Comparison with BLAST On Pentium III 700MH, 1GB BLAST PatternHunter E.coli vs H.inf 716s 14s/68M Arabidopsis 2 vs 4 -- 498s/280M Human 21 vs 22 -- 5250s/417M Human(3G) vs Mouse(x3=9G)* 19 years 20 days All with filter off and identical parameters 16M reads of Mouse genome against Human genome for MIT Whitehead. Best BLAST program takes 19 years at the same sensitivity Quality Comparison: x-axis: alignment rank y-axis: alignment score both axes in logarithmic scale A. thaliana chr 2 vs 4 E. Coli vs H. influenza Genome Alignment by PatternHunter (4 seconds) Prior Literature Over the years, it turns out, that the mathematicians know about certain patterns appear more likely than others: for example, in a random English text, ABC is more likely to appear than AAA, in their study of renewal theory. Random or multiple spaced q-grams were used in the following work: FLASH by Califano & Rigoutsos Multiple filtration by Pevzner & Waterman LSH of Buhler Praparata et al Coding Region & HMM The Hidden Markov Model is a finite set of states, each of which is associated with a probability distribution. Transitions among the states are governed by a set of probabilities called transition probabilities. In a particular state an outcome or observation can be generated, according to the associated probability distribution. It is only the outcome, not the state, which is visible to an external observer and therefore states are ``hidden'' from the outside; hence the name Hidden Markov Model. An HMM has the following components satisfying usual probability laws: The number of states of the model, N. The number of observation symbols in the alphabet, M. A set of state transition probabilities, depending on current state. A probability distribution of the observable symbol at each of the states. Modeling coding region with HMM PatternHunter original assumption: I is a sequence of N independent Bernoulli variables X0, … XN-1 with P(Xi)=p. Coding region third position usually is less conserved. Skipping it should be a good idea (BLAT 110110110). With spaced seeds, we can model it using HMM. Brejova, Brown, Vinar: M(3) HMM: I is a sequence of N independent Bernoulli random variables X0,X1, … XN-1 where Pr(Xi=1)=pi mod 3. In particular: p0 p1 p2 human/mouse 0.82 0.87 0.61 human/fruit fly 0.67 0.77 0.40 DP algorithm naturally extends to M(3), Opt seed (weight 8) for coding region: 11011011000011 Picture of M(3) Extending M(3), Brejova-Brown-Vinar proposed M(8), above, to model dependencies among positions within codons: Σpijk = 1, each codon has conservation pattern ijk with probability pijk. Figures from Brejova-Brown-Vinar Sensitivity of all seeds Under 4 models From Brejova-Brown-Vinar Running PH Download at: www.BioinformaticsSolutions.com java –jar ph.jar –i query.fna –j subject.fna –o out.txt –b –multi 4 -Xmx512m --- for large files -j missing: query.fna self-comparison -db: multiple sequence input, 0,1,2,3 (no, query, subject, both) -W: seed weight -G: open gap penalty (default 5) -E: gap extension (default 1) -q: mismatch penalty (default 1) -r: reward for match (default 1) -model: specify model in binary -H: hits before extension -P: show progress -b: Blast output format -multi: number of seeds to be used (default 1; max 16) Multiple Spaced Seeds: Approaching Smith- Waterman Sensitivity (Li, Ma, Kisman, Tromp, PatternHunter II, J. Bioinfo Comput. Biol. 2004) The biggest problem for Blast was low sensitivity (and low speed). Massive parallel machines are built to do Smith- Waterman exhaustive dynamic programming. Spaced seeds give PH a unique opportunity of using several optimal seeds to achieve optimal sensitivity, this was not possible for Blast technology. Economy --- 2 seeds (½ time slow down) achieve sensitivity of 1 seed of shorter length (4 times speed). With multiple optimal seeds. PH II approaches Smith-Waterman sensitivity, and 3000 times faster. Finding optimal seeds, even one, is NP-hard. Experiment: 29715 mouse EST, 4407 human EST. Sensitivity and speed comparison next two slides. Finding Multiple Seeds (PH II) Computing the hit probability of given k seeds on a uniformly distributed random region is NP-hard. Finding a set of k optimal seeds to maximize the hit probability cannot be approximated within 1-1/e unless NP=P Given k seeds, algorithm DP can be adapted to compute their sensitivity (probability one of the seeds has a hit). Greedy Strategy Compute one optimal seed a. S={a} Given S, find b maximizing hit probability of S U {b}. Coding region, use M(3), 0.8, 0.8, 0.5 (for mod 3 positions) We took 12 CPU-days on 3GHz PC to compute 16 weight 11 seeds. General seeds (first 4): 111010010100110111, 111100110010100001011, 11010000110001010111, 1110111010001111. Sensitivity Comparison with Smith-Waterman (at 100%) The thick dashed curve is the sensitivity of Blastn, seed weight 11. From low to high, the solid curves are the sensitivity of PH II using 1, 2, 4, 8 weight 11 coding region seeds, and the thin dashed curves are the sensitivity 1, 2, 4, 8 weight 11 general purpose seeds, respectively Speed Comparison with Smith-Waterman Smith-Waterman (SSearch): 20 CPU- days. PatternHunter II with 4 seeds: 475 CPU-seconds. 3638 times faster than Smith-Waterman dynamic programming at the same sensitivity. Vector Seeds Definition. A vector seed is an ordered pair Q=(w,T), where w is a weight vector (w1, … wM) and T is a threshold value. An alignment sequence x1,…,xn contains a hit of the seed Q at position k if i=1..M(wi∙xk+i-1) ≥ T The number of nonzero positions in w is called the support of the seed. Comments. Other seeds are special cases of vector seeds. Experiment (Brejova PhD. thesis, 2005) Predicted sensitivity and false positive rate (prob. hit 0.25- homology region of seed length) comparison of various types of seeds in a Bernoulli region of 0.7 match probability, length 64. ----------------------------------------------------------------------------------------- Name Weight vector T Support Sensitivity False positive rate ================================================== BLAST 1111111111 10 10 41% 9.5x10 -7 PH-10 1110101100011011 10 10 60% 9.5x10 -7 BLAST-13-15 111111111111111 13 15 73% 9.2x10-7 VS-13-15 111111011011101111 13 15 85% 9.2x10-7 VS-12-13 111101100110101111 12 13 74% 6.0x10-7 VS-11-12 111011001011010111 11 12 84% 2.2x10-6 Variable weight spaced seeds M. Csuros proposal: use variable weighted spaced seeds depending on the composition of the strings to be hashed. Rarer strings are expected to produce fewer false positives, use seeds with smaller weight --- this increases the sensitivity. Open Question Can we use different patterns (of same weight) at different positions. So this can be called variable position-spaced seeds. At the same weight, we do not increase expected number of false positive hits. Can we increase sensitivity? Find these patterns? How much increase will there be? Practical enough? If the above is not the case, then prove: single optimal seed is always better than variable positions. Old field, new trend Research trend Dozens of papers on spaced seeds have appeared since the original PH paper, in 3 years. Many more have used PH in their work. Most modern alignment programs (including BLAST) have now adopted spaced seeds Spaced seeds are serving thousands of users/day PatternHunter direct users Pharmaceutical/biotech firms. Mouse Genome Consortium, Nature, Dec. 5, 2002. Hundreds of academic institutions. 4. Mathematical Theory of spaced seeds Why they are better: Spaced seeds hits first Asymptotic sensitivity: is there a “best seed” in absolute sense? NP hardness of computing sensitivity Approximation algorithm for sensitivity Why are the spaced seeds better Theorem 1. Let I be a homologous region, homology level p. For any sequence 0≤i0 < i1 … < in-1 ≤ |I|-|s|, let Aj be the event a spaced seed s hits I at ij, Bj be event the consecutive seed hits I at j. (Keich, Li, Ma, Tromp, Discrete Appl. Math. 138(2004), 253-263 ) (1) P(Uj<n Aj) ≥ P(Uj<n Bj); When ij=j, “>” holds. Proof. Induction on n. For n=1, P(A0)=pW=P(B0). Assume the theorem holds for n≤N, we prove it holds for N+1. Let Ek denote the event that, when putting a seed at I[0], 0th, 1st … k-1st 1s in the seed all match & the k-th 1 mismatches. Clearly, Ek is a partition of sample space for all k=0, .. , W, and P(Ek)=pk(1-p) for any seed. (Note, Ek for Ai and Ek for Bi have different indices – they are really different events.) It is sufficient to show, for each k≤W: (2) P(Uj=0..N Aj|Ek)≥P(Uj=0..NBj|Ek) When k=W, both sides of (2) equal to 1. For k<W since (Uj≤kBj)∩Ek=Φ and {Bk+1,Bk+2, … BN} are mutually independent of Ek, we have (3) P(Uj=0..NBj|Ek)=P(Uj=k+1..NBj) Now consider the first term in (2). For each k in {0, … W-1}, at most k+1 of the events Aj satisfy Aj∩Ek=Φ because Aj∩Ek=Φ iff when aligned at ij seed s has a 1 bit at overlapping k-th bit when the seed was at 0 – there are at most k+1 choices. Thus, there exist indices 0<mk+1<mk+2 … <mN ≤N such that Amj∩Ek≠Φ. Since Ek means all previous bits matched (except for k-th), it is clear that Ek is non-negatively correlated with Uj=k+1..N Amj, thus (4) P(Uj=0..N Aj|Ek)≥P(Uj=k+1..N Amj|Ek)≥P(Uj=k+1..NAmj) By the inductive hypothesis, we have P(Uj=k+1..N Amj) ≥ P(Uj=k+1..NBj) Combined with (3),(4), this proves (2), hence the “≥ part” of (1) of the theorem. To prove when ij=j, P(Uj=0..n-1Aj)>P(Uj=0..n-1Bj) --- (5) We prove (5) by induction on n. For n=2, we have P(Uj=0,1Aj)=2pw – p2w-Shift1>2pw-pw+1=P(Uj=0,1Bj) (*) where shift1=number of overlapped bits of the spaced seed s when shifted to the right by 1 bit. For inductive step, note that the proof of (1) shows that for all k=0,1, … ,W, P(Uj=0..n-1Aj|Ek)≥P(Uj=0..n-1Bj|Ek). Thus to prove (5) we only need to prove that P(Uj=0..n-1Aj|E0)>P(Uj=0..n-1Bj|E0). The above follows from the inductive hypothesis as follows: P(Uj=0..n-1Aj|E0)=P(Uj=1..n-1Aj)>P(Uj=1..n-1Bj)=P(Uj=0..n-1Bj|E0). ■ Corollary Corollary: Assume I is infinite. Let ps and pc be the first positions a spaced seed and the consecutive seed hit I, respectively. Then E[ps] < E[pc]. Proof. E[ps]=Σk=0..∞ kP(ps=k) =Σk=0..∞ k[P(ps>k-1)-P(ps>k)] =Σk=0..∞P(ps>k) =Σk=0..∞(1-P(Uj=0..kAj)) <Σk=0..∞(1-P(Uj=0..kBj)) =E[pc]. ■ Asymptotic sensitivity (Buhler-Keich-Sun, JCSS 2004, Li-Ma-Zhang, SODA’2006) Theorem. Given a spaced seed s, the asymptotic hit probability on a homologous region R can be computed, up to a constant factor, in time exponentially proportional to |s|, but independent of |R|. Proof. Let d be one plus the length of the maximum run of 0's in s. Let M'=|s|+d. Consider an infinite iid sequence R=R[1]R[2] ... , Prob(R[i]=1)=p. Let T be set of all length |s| strings that do not dominate s (not matched by s). Let Z[i,j] denote the event that s does not hit at any of the positions R[i], R[i+1], … , R[j]. For any ti T, let xi(n) = P[Z[1,n] | R[n..n+|s|-1] = ti ] . That is, xi(n) is the no-hit probability in the first n positions when they end with ti. We now establish a relationship between xi(n) and xi(n-M’) . Let Ci,j = P[Z[n-M'+1,n] | ti at n, tj at n-M'] x P[tj at n-M' | ti at n]. Thus Ci,j is the probability of generating tj at position n-M' (from ti at n) times the probability not having a hit in a region of length |s|+M' beginning with tj and ending with ti. Especially Ci,j ≠ 0 for any nontrivial seed of length greater than 0 (because of M’ with d spacing). Proof Continued Then, xi(n) = j=1K P[Z[1,n-M'] | tj at n-M'] x P[tj at n-M' | ti at n] x P[Z[n-M'+1,n] | ti at n, tj at n-M'] = j=1K P[Z[1,n-M'] | tj at n-M'] x Ci,j = j=1K Ci,j xj(n-M’) Define the transition matrix C=(Ci,j). It is a K x K matrix. Let x(n) = (x1(n), x2(n), … , xK(n)). Then, assuming that M' divides n, x(n) = x(1) ◦ (CT)n/M’. Because Ci,j > 0 for all i,j, C is a positive matrix. The row sum of the i-th row is the probability that a length-(|s|+M') region ending with ti does not have a hit. The row sum is lower bounded by (1-p)d (1-pW ), where (1-p)d is the probability of generating d 0's and 1-pW is the probability of generating a string in T. It is known that the largest eigenvalue of C is positive and unique, and is lower bounded by the smallest row sum and upper bounded the largest row sum of C. Let λ1>0 be the largest eigenvalue, and λ2, 0 < ||λ2|| < λ1 be the second largest. There is a unique eigenvector corresponding to λ1. x(n) / ||x(n)|| = x(1) ◦ (CT )n/M’ / ||x(1) ◦ (CT)n/M’|| converges to the eigenvector corresponding to λ1. As λ1n tends to zero, we can use standard techniques to normalize xn as y(n) = x(n) / ||x(n)||2 and then x(n+M’) = C y(n) . Proof Continued … Then (by power method) the Rayleigh quotient of λ = (y(n))T C y(n) / (y(n))T y(n) = (y(n))T x(n+M’) converges to λ1. The convergence speed depends on the ratio λ1/ |λ2|, it is known that we can upper bound the second largest eigenvalue λ2 in a positive matrix by: λ2 ≤ λ1 (K-1)/(K+1) where K= maxi,j,k,l (Ci,j Ck,l / Ci,l Ck,j). For any i,j, we have pa (1-p)b (1-p)d ≤ Ci,j ≤ pa (1-p)b where a is the number of 1's in tj, b the number of 0's in tj. K = O(1/(1-p)d ). As (1 - 1/K)K goes to 1/e, the convergence can be achieved in O(K) = O(1/(1- p)d ) steps. The time complexity is therefore upper bounded by an exponential function in |s|. QED Theorem 2. It is NP hard to find the optimal seed (Li, Ma, Kisman, Tromp, J. Bioinfo Comput. Biol. 2004) Proof. Reduction from 3-SAT. Given a 3-SAT instance C1 … Cn, on Boolean variables x1, … , xm, where Ci= l1 + l2 + l3, each l is some xi or its complement. We reduce this to seed selection problem. The required seed has weight W=m+2 and length L=2m+2, and is intended to be of the form 1(01+10)m1, a straightforward variable assignment encoding. The i-th pair of bits after the initial 1 is called the code for variable xi. A 01 code corresponds to setting xi true, while a 10 means false. To prohibit code 11, we introduce, for every 1 ≤ i ≤ m, a region Ai=1 (11)i-101(11)m-i 10m+1 1(11)i-110(11)m-i 1. Since the seed has m 0s, it cannot bridge the middle 0m+1 part, and is thus forced to hit at the start or at the end. This obviously rules out code 11 for variable xi. Code 00 automatically becomes ruled out as well, since the seed must distribute m 1s among all m codes. Finally, for each clause Ci, we introduce Ri =1(11)a-1 c1(11)m-a 1 000 1(11)b-1c2(11)m-b 1 000 1(11)c-1c3(11)m-c 1. The total length of Ri is (2m+2)+3+(2m+2)+3+(2m+2)=6m+12; it consists of three parts each corresponding to a literal in the constraint, with ci being the code for making that literal true. Thus a seed hits Ri iff the assignment encoded by the seed satisfies Ci . ■ The complexity of computing spaced seed sensitivity Theorem. It is also NP hard to compute sensitivity of a given seed. It remains to be NPC even in a uniform region, at any homology level p. Proof. Very complicated. Omitted. PTAS for computing sensitivity of a seed Trivial alg. Sample R poly times, compute the frequency s hits. Use that to approximate the real probability. When p is large, then using Chernoff bound, we have a PTAS. But when p is very small, or seed weight is very large, hitting probability is actually very small, then we do not have a PTAS. Efficient approximation of the hit probability of a seed Theorem. Let the hit probability of s be x. For any ε>0, let N= 6|R|2log |R| / ε2. Then with high probability, in polynomial time and N samples, we can output a value y such that |y-x| ≤ εx. Proof. We give a transformation to transform the low probability event ``s hits R'' dependent on p and weight of s to several events with probability sum ≥ 1, independent of p and seed weight W. Let sr be the reverse string of s. P[s hits R] = P[sr hits R] = i=0|R|-|s| P[sr hits at i, but misses at 0, … , i-1] = i=0|R|-|s| P[sr hits at i] x P[sr misses 0, … , i-1 | sr hits at i] = i=0|R|-|s| pW x P[s misses 1… i | s hits 0] (1) Because P[s hits R] ≥ pW, (1) induces that i=0|R|-|s| P[s misses 1 … i | s hits 0] ≥ 1 Now, we can do efficient sampling, and using Chernoff bounds, the rest of the proof is standard. QED Expected hit distance Corollary. E(second hit position | s hits at 0) = p-W. Proof. i=0|R|-|s| P[s misses 1… i | s hits 0] = i=0|R|-|s| j=i+1|R|-|s|+1 P[s second hits j | s hits 0] = j=0|R|-|s|+1 i=0j-1 P[s second hits j | s hits 0] = i=0|R|-|s|+1 (j-1)P[s second hits j | s hits 0] Combining with (1), we have: P[s hits R]xpW=i=0|R|-|s|+1 (j-1)P[s second hits j|s hits 0] Letting |R| ∞, the left hand is p-W, the right is E(second hit position | s hits at 0). QED Chernoff bound Consider a sequence of n independent tosses of a coin with heads prob. p. The expected number of heads is m=np. What is the probability that the number of heads deviates a lot from this? Let X be the number of successes, then: Pr[|X-m|≥εm] ≤ 2e-εεm/3 Complexity of finding the optimal spaced seeds Theorem 1 [Li-Ma]. Given a seed and it is NP-hard to find its sensitivity, even in a uniform region. Theorem 2 [Li-Ma]. The sensitivity of a given seed can be efficiently approximated by a PTAS. Theorem 3. Given a set of seeds, choose k best can be approximated with ratio 1- 1/e. Theorem 4 [Buhler-Keich-Sun, Li-Ma] The asymptotic hit probability is computable in exponential time in seed length, independent of homologous region length. Theorem 5 [L. Zhang] If the length of a spaced seed is not too long, then it strictly outperforms consecutive seed, in asymptotic hit probability. 5. Applications beyond bioinformatics Stock Market Prediction Zou-Deng-Li: Detecting Market Trends by Ignoring It, Some Days, March, 2005 A real gold mine: 4.6 billion dollars are traded at NYSE daily. Buy low, sell high. Best thing: God tells us market trend Next best thing: A good market indicator. Essentially, a “buy” indicator must be: Sensitive when the market rises Insensitive otherwise. My goal Provide a sensitive, but not aggressive market trend indicator. Learning/inference methods or complete and comprehensive systems are beyond this study. They can be used in conjunction with our proposal. Background Hundreds of market indicators are used (esp. in automated systems). Typical: Common sense: if the past k days are going up, then the market is moving up. Moving average over the last k days. When the average curve and the price curve intersect, buy/sell. Special patterns: a wedge, triangle, etc. Volume Hundreds used in automated trading systems. Note: any method will make money in the right market, lose money in the wrong market Problem Formalization The market movement is modeled as a 0-1 sequence, one bit per day, with 0 meaning market going down, and 1 up. S(n,p) is an n day iid sequence where each bit has probability p being 1 and 1-p being 0. If p>0.5, it is an up market Ik=1k is an indicator that the past k days are 1’s. I8 has sensitivity 0.397 in S(30,0.7), too conservative I8 has false positive rate 0.0043 in S(100, 0.3). Good Iij is an indicator that there are i 1’s in last j days. I811 has high sensitivity 0.96 in S(30,0.7) But it is too aggressive at 0.139 false positive rate in S(100, 0.3). Spaced seeds 1111*1*1111 and 11*11111*11 combine to have sensitivity 0.49 in S(30,0.7) False positive rate 0.0032 in S(100, 0.3). Consider a betting game: A player bets a number k. He wins k dollars for a correct prediction and o.w. loses k dollars. We say an indicator A is better than B, A>B, if A always wins more and loses less. Sleeping on Tuesdays and Fridays Spaced seeds are beautiful indicators: they are sensitive when we need them to be sensitive and not sensitive when we do not want them to be. 1.0 0.9 0.8 11111111 11*11*1*111 0.7 11*11*1*111 8-out-of-11 always beats I811 0.6 0.5 if it bets 4 dollars 0.4 for each dollar 0.3 I811 bets. It is >I8 0.2 too. 0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Two spaced seeds 1.E+00 1.E-01 Observe two spaced 1.E-02 Seeds curve vs I8, the 1.E-03 spaced seeds are 1.E-04 always more sensitive 1.E-05 in p>0.5 region, and 1111*1*1111; 11*11111*11 1.E-06 11111111 less sensitive when p<0.5 1.E-07 9-out-of-11 1.E-08 0.1 0.2 0.3 0.4 0.5 1 0.9 0.8 0.7 0.6 0.5 0.4 1111*1*1111; 11*11111*11 0.3 11111111 0.2 9-out-of-11 0.1 0 0.5 0.6 0.7 0.8 0.9 1 Two experiments We performed two trading experiments One artificial One on real data (S&P 500, Nasdaq indices) Experiment 1: Artificial data This simple HMM 4/5 generates a very artificial simple model 5000 days (bits) EVEN: P(1)=0.5 Indicators: I7, I711, 5 1/50 spaced seeds. 1/60 1/10 1/10 Trading strategy: if there is a hit, buy, and sell 5 days later. UP 1/50 DOWN P(1)=0.7 P(1)=0.3 Reward is: #(1)-#(0) in 1/60 that 5 days times the betting ratio 29/30 24/25 Results of Experiment 1. R #Hits Final MTM #Bankrupcies I7=1111111 $30 12 $679 16 I711 $15 47 $916 14 5 Spaced seeds $25 26 $984 13 Experiment 2. Historical data of S&P 500, from Oct 20, 1982 to Feb. 14, 2005 and NASDAQ, from Jan 2, ’85 to Jan 3, 2005 were downloaded from Yahoo.com. Each strategy starts with $10,000 USD. If an indicator matches, use all the money to buy/sell. While in no ways this experiment says anything affirmatively, it does suggest that spaced patterns can serve a promising basis for a useful indicator, together with other parameters such as trade volume, natural events, politics, psychology. Open Questions I have presented a simple idea of optimized spaced seeds and its applications in homology search and time series prediction. Open questions & current research: Complexity of finding (near) optimal seed, in a uniform region. Note that this is not an NP-hard question. Tighter bounds on why spaced seeds are better. Applications to other areas. Apparently, the same idea works for any time series. Model financial data Can we use varying patterns instead of one seed? When patterns vary at the same weight, we still have same number of expected hits. However, can this increase sensitivity? One seed is just a special case of this. Acknowledgement PH is joint work with B. Ma and J. Tromp PH II is joint work with Ma, Kisman, and Tromp Some joint theoretical work with Ma, Keich, Tromp, Xu, Brown, Zhang. Financial market prediction: J. Zou, X. Deng Financial support: Bioinformatics Solutions Inc, NSERC, Killam Fellowship, CRC chair program, City University of Hong Kong. Here is an example of a BLASTP match (E-value 0) between gene 0189 in C. pneumoniae and gene 131 in C. trachomatis. Query: CPn0189 Score (bits) E-value Aligned with CT131 hypothetical protein 1240 0.0 Query: 1 MKRRSWLKILGICLGSSIVLGFLIFLPQLLSTESRKYLVFSL I HKESGLSCSAEELKISW 60 MKR W KI G L +L L LP+ S+ES KYL S+++KE+GL E+L +SW Sbjct: 1 MKRSPWYKIFGYYLLVGVPLALLALLPKFFSSESGKYLFLSVLNKETGLQF EIEQLHLSW 60 Query: 61 FGRQTARKIKLTG-EAKDEVFSAEKFELDGSLLRLL I YKKPKGITLSGWSLKINEPASID 119 FG QTA+KI++ G ++ E+F+AEK + GSL RLL+Y+ PK + TL+GWSL+I+E S++ Sbjct: 61 FGSQTAKKIRIRGIDSDSEIFAAEKI IVKGSLPRLLL YRFPKALTLTGWSLQIDESLSMN 120 Etc. Note: Because of powerpoint character alignment problems, I inserted some white space to make it look more aligned. Summary and Open Questions Simple ideas: 0 created SW, 1 created Blast. 1 & 0 created PatternHunter. Good spaced seeds help to improve sensitivity and reduce irrelevant hit. Multiple spaced seeds allow us to improve sensitivity to approach Smith-Waterman, not possible with BLAST. Computing spaced seeds by DP, while it is NP hard. Proper data models help to improve design of spaced seeds (HMM) Open Problems: (a) A mathematical theory of spaced seeds. I.e. quantitative statements for Theorem 1. (b) Applications to other fields. Smaller (more solvable) questions: (perhaps good for course projects) Multiple protein seeds for Smith- Waterman sensitivity? Applications in other fields, like internet document search? Prove Theorem 1 for other models – for example, M(3) model 0.8,0.8,0.5 together with specialized seeds. More extensive study of Problem 3.