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
Seq Db
--In the next iteration, the red sequence will
be thrown out. --It matches the query in
non-essential residues

• Iterate
• Find homologs using Blast on query
• Discard very similar homologs
• Align, make a profile, search with profile.
• Why is this more sensitive?

4
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
PrF0.71 PrY0.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
  ProbHMM 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
• PrS,PM PrSP,M PrPM
• PrALIL AND M1I1M2M3 M
• PrALIL M1I1M2M3,M PrM1I1M2M3 M
• PrALIL M ?

12
Formally

• The emitted sequence is SS1S2Sm
• The path traversed id P1P2P3..
• ej(s) emission probability of symbol s in state
  Pj
• Transition probability Tj,k Probability of
  transitioning from state j to state k.
• Pr(P,SM) eP1(S1) TP1,P2 eP2(S2)
• What is Pr(SM)?

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(SM) ?P Pr(S,PM)
• OR, it is the probability of the most likely path
• Pr(SM) maxP Pr(S,PM)
• 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,jM) be the probability of the most
  likely solution that emits S1Si, and ends in
  state j (is it sufficient to compute this?)
• Pmax(i,jM) max k Pmax(i-1,k) Tk,j ej(Si)
  (Viterbi)
• Psum(i,jM) ? k (Psum(i-1,k) Tk,j) ej(Si)

15
Profile HMM membership
A L - L A I V L A I - L
A L I L
Path M1 M2 I2 M3

• 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

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
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/expresse
  d 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
Output

• 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

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
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

41
(No Transcript)
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-i1 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.

51
DNA Signals

• Coding versus non-coding
• Splice Signals
• Translation start

52
Splice signals

• GT is a Donor signal, and AG is the acceptor
  signal

GT
AG
53
PWMs
321123456 AAGGTGAGT CCGGTAAGT GAGGTGAGG TAGGTAAGG

• Fixed length for the splice signal.
• Each position is generated independently
  according to a distribution
• Figure shows data from gt 1200 donor sites

54
MDD

• PWMs do not capture correlations between
  positions
• Many position pairs in the Donor signal are
  correlated

55

• Choose the position which has the highest
  correlation score.
• Split sequences into two those which have the
  consensus at position I, and the remaining.
• Recurse until ltTerminating conditionsgt

56
MDD for Donor sites
57
De novo Gene prediction Sumary

• Various signals distinguish coding regions from
  non-coding
• HMMs are a reasonable model for Gene structures,
  and provide a uniform method for combining
  various signals.
• Further improvement may come from improved signal
  detection

58
How many genes do we have?
Nature
Science
59
Alternative splicing
60
Comparative methods

• Gene prediction is harder with alternative
  splicing.
• One approach might be to use comparative methods
  to detect genes
• Given a similar mRNA/protein (from another
  species, perhaps?), can you find the best parse
  of a genomic sequence that matches that target
  sequence
• Yes, with a variant on alignment algorithms that
  penalize separately for introns, versus other
  gaps.

61
Gene Features
ATG
5 UTR
3 UTR
exon
intron
Translation start
Acceptor
Donor splice site
Transcription start
62
Gene Finding

• Given genomic DNA

63
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?

64
Generalizing
I
E
X
AAAAAA AAAAAC AAAAAG AAAAAT
10
10
5
5
20
10
Compute a frequency count for all hexamers. Use
this to decide whether a sequence X is an
exon/intron.
65
A geometric approach

• Plot the following vectors
• E 10, 20
• I 10, 5
• V3 5, 10
• V4 9, 15
• Is V3 more like E or more like I?

20
15
10
5
15
10
5
66
Choosing between Introns and Exons

• V V/V
• All vectors have the same length (lie on the unit
  circle)
• Next, compute the angle to E, and I.
• Choose the feature that is closer (smaller
  angle.

V3
E
I
67
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

68
Coding region can be detected

• Plot the E-score using a sliding window of fixed
  length.
• The (large) exons will show up reliably.
• Not enough to predict gene boundaries reliably

E-score
69
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
70
Combining Signals

• We can compute the following
• E-scorei,j
• I-scorei,j
• D-scorei
• I-scorei
• Goal is to find coordinates that maximize the
  total score

i
j
71
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

72
Combining signals using D.P.

• An HMM is the best way to model and optimize the
  combination of signals
• Here, we will use a simpler approach which is
  essentially the same as the Viterbi algorithm for
  HMMs, but without the formalism.

73
Gene finding reformulated
IIIIIEEEEEEIIIIIIEEEEEEIIIIEEEEEEEIIIII

• Recall that our goal was to identify the
  coordinates of the exons.
• Instead, we label every nucleotide as I
  (Intron/Intergenic) or E (Exon). For simplicity,
  we treat intergenic and introns as identical.

74
Gene finding reformulated
i1
i2
i3
i4
IIIIIEEEEEEIIIIIIEEEEEEIIIIEEEEEE IIIII

• Given a labeling L, we can score it as
• I-score0..i1 E-scorei1..i2 D-scorei21
  I-scorei21..i3-1 A-scorei3-1
  E-scorei3..i4 .
• Goal is to compute a labeling with maximum score.

75
Optimum labeling using D.P. (Viterbi)

• Define VE(i) Best score of a labeling of the
  prefix 1..i such that the i-th position is
  labeled E
• Define VI(i) Best score of a labeling of the
  prefix 1..i such that the i-th position is
  labeled I
• Why is it enough to compute VE(i) VI(i) ?

76
Optimum parse of the gene
j
i
j
i
77
Generalizing

• Note that we deal with two states, and consider
  all paths that move between the two states.

E
I
i
78
Generalizing

• We did not deal with the boundary cases in the
  recurrence.
• Instead of labeling with two states, we can label
  with multiple states,
• Einit, Efin, Emid,
• I, IG (intergenic)

IG
I
Efin
Emid
Note all links are not shown here
Einit
79
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.
• They allow an automated approach to weighting
  features.

80
An HMM for Gene structure
81
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.

82
Length distributions of Introns Exons
83
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.

84
Forward algorithm for gene finding
qk
j
i
Duration Prob. Probability that you stayed in
state qk for j-i1 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
85
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.

86
DNA Signals

• Coding versus non-coding
• Splice Signals
• Translation start

87
Splice signals

• GT is a Donor signal, and AG is the acceptor
  signal

GT
AG
88
PWMs
321123456 AAGGTGAGT CCGGTAAGT GAGGTGAGG TAGGTAAGG

• Fixed length for the splice signal.
• Each position is generated independently
  according to a distribution
• Figure shows data from gt 1200 donor sites

89
MDD

• PWMs do not capture correlations between
  positions
• Many position pairs in the Donor signal are
  correlated

90

• Choose the position which has the highest
  correlation score.
• Split sequences into two those which have the
  consensus at position I, and the remaining.
• Recurse until ltTerminating conditionsgt

91
MDD for Donor sites
92
Gene prediction Summary

• Various signals distinguish coding regions from
  non-coding
• HMMs are a reasonable model for Gene structures,
  and provide a uniform method for combining
  various signals.
• Further improvement may come from improved signal
  detection

93
How many genes do we have?
Nature
Science
94
Alternative splicing
95
Comparative methods

• Gene prediction is harder with alternative
  splicing.
• One approach might be to use comparative methods
  to detect genes
• Given a similar mRNA/protein (from another
  species, perhaps?), can you find the best parse
  of a genomic sequence that matches that target
  sequence
• Yes, with a variant on alignment algorithms that
  penalize separately for introns, versus other
  gaps.

96
Comparative gene finding tools

• Procrustes/Sim4 mRNA vs. genomic
• Genewise proteins versus genomic
• CEM genomic versus genomic
• Twinscan Combines comparative and de novo
  approach.

97
Databases

• RefSeq and other databases maintain sequences of
  full-length transcripts.
• We can query using sequence.