Gene Prediction: Statistical Methods

1. Gene Prediction: Statistical Methods (Lecture for CS498-CXZ Algorithms in Bioinformatics) Sept. 20, 2005 ChengXiang Zhai Department of Computer Science University of Illinois, Urbana-Champaign Many slides are taken/adapted from http://www.bioalgorithms.info/slides.htm and Ying Xu’s lecture

2. Approaches to Gene Prediction • Similarity-based approaches: • Exploit the fact that many genes are conserved across species • Can be highly reliable • Only good for finding unknown genes • Statistical approaches • Exploit statistical characteristics of coding regions and non-coding regions and other knowledge about genes • Can potentially detect new genes • May not be reliable • They can/should be combined • Currently no principled approaches for doing this Given a new genome, identify “known genes” first Learn from “known genes” to identify new gene

3. Gene Prediction Analogy • Newspaper written in unknown language • Certain pages contain encoded message, say 99 letters on page 7, 30 on page 12 and 63 on page 15. • How do you recognize the message? You could probably distinguish between the ads and the story (ads contain the “$” sign often) • Statistics-based approach to Gene Prediction tries to make similar distinctions between exons and introns.

4. Statistical Approach: Metaphor in Unknown Language Noting the differing frequencies of symbols (e.g. ‘%’, ‘.’, ‘-’) and numerical symbols could you distinguish between a story and the stock report in a foreign newspaper?

5. A Few Basic Questions • What is exactly a gene for the purpose of prediction? • In Prokaryotes, gene = mRNAProtein • In Eukaryotes, gene = Exon (coding region) • What does a gene look like? • Where does it start? • Where does it end? • What is the codon usage inside a gene (exon)? • What is the codon usage outside a gene (intron)? • … • How do we exploit such knowledge to identify genes?

6. Statistical Characteristics of a Gene • Gene starts with a start codon • Gene ends at a stop codon • Splicing signals • Codon usage distributions • …

7. Gene Structure

8. Genetic Code and Stop Codons UAA, UAG and UGA correspond to 3 Stop codons that (together with Start codon ATG) delineate Open Reading Frames

9. Six Frames in a DNA Sequence • stop codons • start codons CTGCAGACGAAACCTCTTGATGTAGTTGGCCTGACACCGACAATAATGAAGACTACCGTCTTACTAACAC CTGCAGACGAAACCTCTTGATGTAGTTGGCCTGACACCGACAATAATGAAGACTACCGTCTTACTAACAC CTGCAGACGAAACCTCTTGATGTAGTTGGCCTGACACCGACAATAATGAAGACTACCGTCTTACTAACAC CTGCAGACGAAACCTCTTGATGTAGTTGGCCTGACACCGACAATAATGAAGACTACCGTCTTACTAACAC GACGTCTGCTTTGGAGAACTACATCAACCGGACTGTGGCTGTTATTACTTCTGATGGCAGAATGATTGTG GACGTCTGCTTTGGAGAACTACATCAACCGGACTGTGGCTGTTATTACTTCTGATGGCAGAATGATTGTG GACGTCTGCTTTGGAGAACTACATCAACCGGACTGTGGCTGTTATTACTTCTGATGGCAGAATGATTGTG GACGTCTGCTTTGGAGAACTACATCAACCGGACTGTGGCTGTTATTACTTCTGATGGCAGAATGATTGTG

10. The “Sly Fox” & Effect of Base Deletion • In the following string THE SLY FOX AND THE SHY DOG • Delete 1, 2, and 3 nucleotifes after the first ‘S’: THE SYF OXA NDT HES HYD OG THE SFO XAN DTH ESH YDO G THE SOX AND THE SHY DOG • Which of the above makes the most sense?

11. Splicing Signals and Exon Boudnaries Exons are interspersed with introns and typically flanked by GT and AG

12. Splicing mechanism (http://genes.mit.edu/chris/)

13. Donor and Acceptor Sites: GT and AG dinucleotides • The beginning and end of exons are signaled by donor and acceptor sites that usually have GT and AC dinucleotides • Detecting these sites is difficult, because GT and AC appear very often Donor Site Acceptor Site GT AC exon 1 exon 2

14. Donor and Acceptor Sites: Motif Logos Donor: 7.9 bits Acceptor: 9.4 bits (Stephens & Schneider, 1996) (http://www-lmmb.ncifcrf.gov/~toms/sequencelogo.html)

15. Codon Usage in Human Genome Biased codon usage in exons allows us to distinguish exons from introns

16. Codon Frequencies • Coding sequences are translated into protein sequences • We found the following – the dimer frequency in protein sequences is NOT evenly distributed The average frequency is 5% Some amino acids prefer to be next to each other Some other amino acids prefer to be not next to each other shewanella

17. Dicodon Frequencies • The biased (uneven) dimer frequencies are the foundation of many gene finding programs! • Basic idea of gene finding – if a dimer has lower than average dimer frequency; this means that proteins prefer not to have such dimers in its sequence; otherwise proteins prefer to have such dimers Hence if we see a dicodon encoding this dimer, we may want to bet against this dicodon being in a coding region!

18. General Steps for Gene Prediction • Identify candidate exons in Open Reading Frames (ORFs) • Determine ORFs: An ORF starts with a start codon and ends at a stop codon • Determine sites for receptors/donors • Evaluate the potential of a candidate exon for coding (Exploit codon usages)

19. Step 1: Identify ORFs

20. Long vs.Short ORFs • Long open reading frames may be a gene • At random, we should expect one stop codon every (64/3) ~= 21 codons • However, genes are usually much longer than this • A basic approach is to scan for ORFs whose length exceeds certain threshold • This is naïve because some genes (e.g. some neural and immune system genes) are relatively short

21. Prediction of Translation Starts • Translation start: ATG • How to predict a translation start • Collect a set of experimentally validated translation starts with flanking regions and align them up …. ATG …… GCCATGGCGA ….. ACGATGCTGT …. GACATGGTAC … AGGATGGGCT … GCGATGTGGC …

22. Prediction of Translation Starts • Certain nucleotides prefer to be in certain position around start “ATG” and other nucleotides prefer not to be there • The “biased” nucleotide distribution is information! It is a basis for translation start prediction • Question: which one is more probable to be a translation start? ATG A C G T -4 -3 -2 -1 +3 +4 +5 +6 CACC ATG GC TCGA ATG TT

23. A C G T Prediction of Translation Starts • Mathematical model: Fi (X): frequency of X (A, C, G, T) in position i • Score a string by i log (Fi (X)/0.25) CACC ATG GC TCGA ATG TT log (6/25) + log (6/25) + log (15/25) + log (7/25) + log (13/25) + log (14/25) = -(0.62 + 0.62 + 0.22 + 0.55 + 0.28 + 0.25) = -2.54 log (58/25) + log (49/25) + log (40/25) + log (50/25) + log (43/25) + log (39/25) = 0.37 + 0.29 + 0.20 + 0.30 + 0.24 + 0.20 = 1.60 The model captures our intuition!

24. Prediction of Translation Starts • Build a mathematical model, based on collected translation start sequence • For each candidate translation start sequence, apply the model and get a score • If the score if larger than zero, predict it is a “translation start”; the higher score, the higher the probability the prediction is true ATG ……

25. Step 2: Identify Exon Boundaries

26. { translation start, acceptor site } { translation stop, donor site } Prediction of Splice Junction Sites • A start exon starts with a translation start and ends with a donor site • An internal exon starts with an acceptor site and ends with a donor site • A terminal exon starts with an acceptor site and ends with a stop codon exon Accurate prediction of exons/genes requires accurate prediction of splice junctions

27. Prediction of Splice Junction Sites • Splice junctions: • donor site: coding region | GT • acceptor: AG | coding region • Like translation starts, the flanks of splice junctions (acceptors and donors) show “biased” distributions of nucleotides in certain positions • These biased distributions of nucleotides are the basis for prediction of splice junctions

28. Prediction of Acceptor Sites • Nucleotide distribution in the flanks of acceptors Multiple positions have high “information content” Information content:  log (F(X)/0.25) If every nucleotide has 0.25 frequency in a position, then the position’s information content is ZERO. Use “information content” as a criterion for determining the length of flanks

29. Prediction of Acceptor Sites • Mathematical model: Fi (X): frequency of X (A, C, G, T) in position i • Score a segment as a candidate acceptor site by ilog (Fi (X)/0.25) • For each candidate acceptor sequence, apply the model and get a score • If the score if larger than zero, predict it is an “acceptor”; the higher score, the higher the probability the prediction is true AG

30. Prediction of Donor Sites • Nucleotide distribution in the flanks of donors • Mathematical model: Fi (X): frequency of X (A, C, G, T) in position i • Score a segment as a possible donor site by i log (Fi (X)/0.25)

31. Prediction of Donor Sites • For each candidate donor sequence, apply the model and get a score • If the score is larger than zero, predict it is a “donor”; the higher score, the higher the probability the prediction is true GT

32. Prediction of Donors/Acceptors • Position specific weight matrix model • Build a “position specific weight matrix model” • collect known {donor, acceptor} sequences and align them so that the GT or AG are aligned • Calculate the percentage of each type of nucleotide at each position • There are more sophisticated models for capturing higher order relationships between positions

33. Prediction of Exons • For each orf, find all donor and acceptor candidates by finding GT and YAG motifs • Score each donor and acceptor candidate using our position-specific weight matrix models • Find all pairs of (acceptor, donor) above some thresholds • Score the coding potential of the segment [donor, acceptor], using the hexmer model CAG GT CAG GT GT

34. Step 3: Classify Candidate Exons

35. Testing Exons: Codon Usage • Create a 64-element hash table and count the frequencies of codons in a candidate exon • Amino acids typically have more than one codon, but in nature certain codons are more in use • Uneven use of the codons may characterize a real gene

36. Codon Usage and Likelihood Ratio • An ORF is more “believable” than another if it has more “likely” codons • Do sliding window calculations to find ORFs that have the “likely” codon usage • Allows for higher precision in identifying true ORFs; much better than merely testing for length. • However, average vertebrate exon length is 130 nucleotides, which is often too small to produce reliable peaks in the likelihood ratio • Further improvement: in-frame hexamer count (frequencies of pairs of consecutive codons)

37. Codon Usage in Mouse Genome AA codon /1000 frac Ser TCG 4.31 0.05 Ser TCA 11.44 0.14 Ser TCT 15.70 0.19 Ser TCC 17.92 0.22 Ser AGT 12.25 0.15 Ser AGC 19.54 0.24 Pro CCG 6.33 0.11 Pro CCA 17.10 0.28 Pro CCT 18.31 0.30 Pro CCC 18.42 0.31 AA codon /1000 frac Leu CTG 39.95 0.40 Leu CTA 7.89 0.08 Leu CTT 12.97 0.13 Leu CTC 20.04 0.20 Ala GCG 6.72 0.10 Ala GCA 15.80 0.23 Ala GCT 20.12 0.29 Ala GCC 26.51 0.38 Gln CAG 34.18 0.75 Gln CAA 11.51 0.25

38. Exon Prediction Method 1:TestCode

39. TestCode • Statistical test described by James Fickett in 1982: tendency for nucleotides in coding regions to be repeated with periodicity of 3 • Judges randomness instead of codon frequency • Finds “putative” coding regions, not introns, exons, or splice sites • TestCode finds ORFs/Exons based on compositional bias with a periodicity of three

40. TestCode Statistics • Define a window size no less than 200 bp, slide the window the sequence down 3 bases. In each window: • Calculate for each base {A, T, G, C} • max (n3k+1, n3k+2, n3k) / min ( n3k+1, n3k+2, n3k) • Use these values to obtain a probability from a lookup table (which was a previously defined and determined experimentally with known coding and noncoding sequences

41. TestCode Statistics (cont’d) • Probabilities can be classified as indicative of " coding” or “noncoding” regions, or “no opinion” when it is unclear what level of randomization tolerance a sequence carries • The resulting sequence of probabilities can be plotted

42. TestCode Sample Output Coding No opinion Non-coding

43. Exon Prediction Method 2:Likelihood Ratio/Suprevised Learning

44. Prediction of Coding Regions • Build a coding-region predictor • Collect coding and non-coding sequences from GenBank (NCBI) • Calculate hexmer frequencies in coding and noncoding sequences (hexmer = 6-mer, or dimer of amino acids) • Application – coding-region prediction • consider both strands: forward and reverse • for each strand, identify all the orfs • for prokaryotic genome, for each possible translation start ATG, evaluate the coding of [ATG, STOP] • for eukaryotic genome, for the segment defined by each pair of possible start/acceptor & donor/stop, evaluate the coding potential

45. Prediction of Exons • For each segment [acceptor, donor], we get three scores (coding potential, donor score, acceptor score) • Various possibilities • all three scores are high – probably true exon • all three scores are low – probably not a real exon • all in the middle -- ?? • some scores are high and some are low -- ?? • What are the rules for exon prediction?

46. coding: noncoding: Prediction of Exons • Learning to classify exons from nonexons • Collect a set of exons and non-exons • Score them using our scoring schemes • Plot them as follows • “draw” a separating line between exons and non-exons • Making a prediction based on which side of the separating line a new point falls

47. Prediction of Exons • A “classifier” can be trained to separate exons from non-exons, based on the three scores • Closer to reality – other factors could also help to distinguish exons from non-exons exon length distribution 150 bp coding density is different in regions with different G+C contents 50% G+C A practical gene finding software may use many features to distinguish exons from non-exons

48. Prediction of Exons Each box represents a predicted exon A true exon typically has more than one predicted candidates, overlapping with each other

49. Gene Prediction in a New Genome • Dicodon (hexmer) frequencies are different from genome to genome – gene finder for one genome cannot be directly applied to another genome bovine shewanella

50. Popular Gene Prediction Algorithms • GENSCAN: uses Hidden Markov Models (HMMs) • TWINSCAN • Uses both HMM and similarity (e.g., between human and mouse genomes) HMM will be covered later in the course