CSE182-L7

1. CSE182-L7 Protein Sequence Analysis using HMMs, Gene Finding

2. Domain analysis via profiles • Given a database of profiles of known domains/families, we can query our sequence against each of them, and choose the high scoring ones to functionally characterize our sequences. • What if the sequence matches some other sequences weakly (using BLAST), but does not match any known profile?

3. Psi-BLAST idea • Iterate: • Find homologs using Blast on query • Discard very similar homologs • Align, make a profile, search with profile. • Why is this more sensitive? Seq Db --In the next iteration, the red sequence will be thrown out. --It matches the query in non-essential residues

4. Pigeonhole principle again: • If profile of length m must score >= T • Then, a sub-profile of length l must score >= lT|/m • Generate all l-mers that score at least lT|/M • Search using an automaton • Multiple alignment: • Use ungapped multiple alignments only Psi-BLAST speed • Two time consuming steps. • Multiple alignment of homologs • Searching with Profiles. • Does the keyword search idea work?

5. Representation 3: HMMs • Building good profiles relies upon good alignments. • Difficult if there are gaps in the alignment. • Psi-BLAST/BLOCKS etc. work with gapless alignments. • An HMM representation of Profiles helps put the alignment construction/membership query in a uniform framework. • Also allows for position specific gap scoring. V

6. QUIZ! • Question: • your ‘friend’ likes to gamble. • He tosses a coin: HEADS, he gives you a dollar. TAILS, you give him a dollar. • Usually, he uses a fair coin, but ‘once in a while’, he uses a loaded coin. • Can you say what fraction of the times he loads the coin?

7. The generative model • Think of each column in the alignment as generating a distribution. • For each column, build a node that outputs a residue with the appropriate distribution 0.71 Pr[F]=0.71 Pr[Y]=0.14 0.14

8. A simple Profile HMM • Connect nodes for each column into a chain. Thie chain generates random sequences. • What is the probability of generating FKVVGQVILD? • In this representation • Prob [New sequence S belongs to a family]= Prob[HMM generates sequence S] • What is the difference with Profiles?

9. Profile HMMs can handle gaps • The match states are the same as on the previous page. • Insertion and deletion states help introduce gaps. • A sequence may be generated using different paths.

10. Example A L - L A I V L A I - L • Probability [ALIL] is part of the family? • Note that multiple paths can generate this sequence. • M1I1M2M3 • M1M2I2M3 • In order to compute the probabilities, we must assign probabilities of transition between states

11. Profile HMMs • Directed Automaton M with nodes and edges. • Nodes emit symbols according to ‘emission probabilities’ • Transition from node to node is guided by ‘transition probabilities’ • Joint probability of seeing a sequence S, and path P • Pr[S,P|M] = Pr[S|P,M] Pr[P|M] • Pr[ALIL AND M1I1M2M3|M] = Pr[ALIL| M1I1M2M3,M] Pr[M1I1M2M3|M] • Pr[ALIL | M] = ?

12. Formally • The emitted sequence is S=S1S2…Sm • The path traversed id P1P2P3.. • ej(s) = emission probability of symbol s in state Pj • Transition probability T[j,k] : Probability of transitioning from state j to state k. • Pr(P,S|M) = eP1(S1) T[P1,P2] eP2(S2) …… • What is Pr(S|M)?

13. Two solutions • An unknown (hidden) path is traversed to produce (emit) the sequence S. • The probability that M emits S can be either • The sum over the joint probabilities over all paths. • Pr(S|M) = ∑P Pr(S,P|M) • OR, it is the probability of the most likely path • Pr(S|M) = maxP Pr(S,P|M) • Both are appropriate ways to model, and have similar algorithms to solve them.

14. Viterbi Algorithm for HMM A L - L A I V L A I - L • Let Pmax(i,j|M) be the probability of the most likely solution that emits S1…Si, and ends in state j (is it sufficient to compute this?) • Pmax(i,j|M) = max k Pmax(i-1,k) T[k,j] ej(Si) (Viterbi) • Psum(i,j|M) = ∑ k(Psum(i-1,k) T[k,j]) ej(Si)

15. Profile HMM membership A L - L A I V L A I - L • We can use the Viterbi/Sum algorithm to compute the probability that the sequence belongs to the family. • Backtracking can be used to get the path, which allows us to give an alignment A L I L Path:M1 M2 I2 M3

16. Summary • HMMs allow us to model position specific gap penalties, and allow for automated training to get a good alignment. • Patterns/Profiles/HMMs allow us to represent families and foucs on key residues • Each has its advantages and disadvantages, and needs special algorithms to query efficiently.

17. Protein Domain databases HMM • A number of databases capture proteins (domains) using various representations • Each domain is also associated with structure/function information, parsed from the literature. • Each database has specific query mecahnisms that allow us to compare our seqeunces against them, and assign function 3D

18. Gene Finding What is a Gene?

19. Gene • We define a gene as a location on the genome that codes for proteins. • The genic information is used to manufacture proteins through transcription, and translation. • There is a unique mapping from triplets to amino-acids

20. Eukaryotic gene structure

21. Translation • The ribosomal machinery reads mRNA. • Each triplet is translated into a unique amino-acid until the STOP codon is encountered. • There is also a special signal where translation starts, usually at the ATG (M) codon.

22. Translation • The ribosomal machinery reads mRNA. • Each triplet is translated into a unique amino-acid until the STOP codon is encountered. • There is also a special signal where translation starts, usually at the ATG (M) codon. • Given a DNA sequence, how many ways can you translate it?

23. ATG 5’ UTR 3’ UTR exon intron Translation start Acceptor Donor splice site Transcription start Gene Features

24. Gene identification • Eukaryotic gene definitions: • Location that codes for a protein • The transcript sequence(s) that encodes the protein • The protein sequence(s) • Suppose you want to know all of the genes in an organism. • This was a major problem in the 70s. PhDs, and careers were spent isolating a single gene sequence. • All of that changed with better reagents and the development of high throughput methods like EST sequencing

25. Expressed Sequence Tags • It is possible to extract all of the mRNA from a cell. • However, mRNA is unstable • An enzyme called reverse transcriptase is used to make a DNA copy of the RNA. • Use DNA polymerase to get a complementary DNA strand. • Sequence the (stable) cDNA from both ends. • This leads to a collection of transcripts/expressed sequences (ESTs). • Many might be from the same gene AAAA TTTT AAAA TTTT

26. EST sequencing • The expressed transcript (mRNA) has a poly-A tail at the end, which can be used as a template for Reverse Transcriptase. • This collection of DNA has only the spliced message! • It is sampled at random and sequenced from one (3’/5’) or both ends. • Each message is sampled many times. • The resulting collection of sequences is called an EST database AAAA TTTT AAAA TTTT

27. EST Sequencing • Often, reverse transcriptase breaks off early. Why is this a good thing? • The 3’ end may not have a much coding sequence. • We can assemble the 5’ end to get more of the coding sequence

28. Project Input • EST clustering and assembly • Given a collection of EST (3’/5’) sequences, your goal is to cluster all ESTs from the same gene, and produce a consensus. • Note that all the 3’ ESTs should line up at the 3’ end. • 5’ and 3’ ESTs from the same clone should have the same clone ID, which should allow us to recruit them Output

29. Project Extra credit • Some genes may be alternatively spliced and may have multiple transcripts • Can you deconvolute the information back from ESTs? ATG

30. ATG 5’ UTR 3’ UTR exon intron Translation start Acceptor Donor splice site Transcription start Computational Gene Finding • Given Genomic DNA, identify all the coordinates of the gene

31. Gene Finding: The 1st generation • Given genomic DNA, does it contain a gene (or not)? • Key idea: The distributions of nucleotides is different in coding (translated exons) and non-coding regions. • Therefore, a statistical test can be used to discriminate between coding and non-coding regions.

32. Coding versus Non-coding • You are given a collection of exons, and a collection of intergenic sequence. • Count the number of occurrences of ATGATG in Introns and Exons. • Suppose 1% of the hexamers in Exons are ATGATG • Only 0.01% of the hexamers in Intons are ATGATG • How can you use this idea to find genes?

33. Generalizing I E AAAAAA AAAAAC AAAAAG AAAAAT Compute a frequency count for all hexamers. Use this to decide whether a sequence is an exon/intron

34. Coding versus non-coding • Fickett and Tung (1992) compared various measures • Measures that preserve the triplet frame are the most successful. • Genscan: 5th order Markov Model • Conservation across species

35. Coding vs. non-coding regions Compute average coding score (per base) of exons and introns, and take the difference. If the measure is good, the difference must be biased away from 0.

36. Coding differential for 380 genes

37. Other Signals ATG AG GT Coding

38. Coding region can be detected • Plot the coding score using a sliding window of fixed length. • The (large) exons will show up reliably. • Not enough to predict gene boundaries reliably Coding

39. Other Signals • Signals at exon boundaries are precise but not specific. Coding signals are specific but not precise. • When combined they can be effective ATG AG GT Coding

40. The second generation of Gene finding • Ex: Grail II. Used statistical techniques to combine various signals into a coherent gene structure. • It was not easy to train on many parameters. Guigo & Bursett test revealed that accuracy was still very low. • Problem with multiple genes in a genomic region

42. HMMs and gene finding • HMMs allow for a systematic approach to merging many signals. • They can model multiple genes, partial genes in a genomic region, as also genes on both strands.

43. The Viterbi Algorithm

44. HMMs and gene finding • The Viterbi algorithm (and backtracking) allows us to parse a string through the states of an HMM • Can we describe Eukaryotic gene structure by the states of an HMM? • This could be a solution to the GF problem.

45. An HMM for Gene structure

46. Generalized HMMs, and other refinements • A probabilistic model for each of the states (ex: Exon, Splice site) needs to be described • In standard HMMs, there is an exponential distribution on the duration of time spent in a state. • This is violated by many states of the gene structure HMM. Solution is to model these using generalized HMMs.

47. Length distributions of Introns & Exons

48. Generalized HMM for gene finding • Each state also emits a ‘duration’ for which it will cycle in the same state. The time is generated according to a random process that depends on the state.

49. Forward algorithm for gene finding qk j i Duration Prob.: Probability that you stayed in state qk for j-i+1 steps Emission Prob.: Probability that you emitted Xi..Xj in state qk (given by the 5th order markov model) Forward Prob: Probability that you emitted I symbols and ended up in state qk

50. HMMs and Gene finding • Generalized HMMs are an attractive model for computational gene finding • Allow incorporation of various signals • Quality of gene finding depends upon quality of signals.