Probabilities and Probabilistic Models

1
Probabilities and Probabilistic Models
2
Probabilistic models

• A model means a system that simulates an object
  under consideration.
• A probabilistic model is a model that produces
  different outcomes with different probabilities
  it can simulate a whole class of objects,
  assigning each an associated probability.
• In bioinformatics, the objects usually are DNA or
  protein sequences and a model might describe a
  family of related sequences.

3
Examples

1. The roll of a six-sided die six parameters p1,
   p2, , p6, where pi is the probability of rolling
   the number i. For probabilities, pi gt 0 and
2. Three rolls of a die the model might be that the
   rolls are independent, so that the probability of
   a sequence such as 2, 4, 6 would be p2 p4 p6.
3. An extremely simple model of any DNA or protein
   sequence is a string over a 4 (nucleotide) or 20
   (amino acid) letter alphabet. Let qa denote the
   probability, that residue a occurs at a given
   position, at random, independent of all other
   residues in the sequence. Then, for a given
   length n, the probability of the sequence
   x1,x2,,xn is

4
Conditional, joint, and marginal probabilities

• two dice D1 and D2. For j 1,2, assume that the
  probabi-lity of using die Dj is P(Dj ), and for
  i 1,2,?, 6, assume that the probability of
  rolling an i with dice Dj is
• In this simple two dice model, the conditional
  probability of rolling an i with dice Dj is
• The joint probability of picking die Dj and
  rolling an i is
• The probability of rolling i marginal
  probability

5
Maximum likelihood estimation

• Probabilistic models have parameters that are
  usually estimated from large sets of trusted
  examples, called a training set.
• For example, the probability qa for seeing amino
  acid a in a protein sequence can be estimated as
  the observed frequency fa of a in a database of
  known protein sequences, such as SWISS-PROT.
• This way of estimating models is called Maximum
  likelihood estimation, because it can be shown
  that using the observed frequencies maximizes the
  total probability of the training set, given the
  model.
• In general, given a model with parameters and a
  set of data D, the maximum likelihood estimate
  (MLE) for ? is the value which maximizes P(D ?).

6
Model comparison problem

• An occasionally dishonest casino uses two kinds
  of dice, of which 99 are fair, but 1 are
  loaded, so that a 6 appears 50 of the time.
• We pick up a dice and roll 6, 6, 6. This looks
  like a loaded die, is it? This is an example of a
  model comparison problem.
• I.e., our hypothesis Dloaded is that the die is
  loaded. The other alternative is Dfair. Which
  model fits the observed data better? We want to
  calculate
• P(Dloaded 6, 6, 6)

7
Prior and posterior probability

• P(Dloaded 6, 6, 6) is the posterior
  probability that the dice is loaded, given the
  observed data.
• Note that the prior probability of this
  hypothesis is 1/100 prior because it is our
  best guess about the dice before having seen any
  information about the it.
• The likelihood of the hypothesis Dloaded
• Posterior probability using Bayes theorem

8
Comparing models using Bayes theorem

• We set X Dloaded and Y 6, 6, 6, thus
  obtaining

• The probability P(Dloaded) of picking a loaded
  die is 0.01.
• The probability P(6, 6, 6 Dloaded) of rolling
  three sixes using a loaded die is 0.53 0.125.
• The total probability P(6, 6, 6) of three sixes
  is
• P(6, 6, 6 Dloaded) P(Dloaded) P(6, 6,
  6 Dfair)P(Dfair).
• Now,
• Thus, the die is probably fair.

9
Biological example

• Lets assume that extracellular (ext) proteins
  have a slightly different composition than
  intercellular (int) ones. We want to use this to
  judge whether a new protein sequence x1,, xn is
  ext or int.
• To obtain training data, classify all proteins in
  SWISS-PROT into ext, int and unclassifiable ones.
• Determine the frequencies and of
  each amino acid a in ext and int proteins,
  respectively.
• To be able to apply Bayes theorem, we need to
  determine the priors pint and pext, i.e. the
  probability that a new (unexamined) sequence is
  extracellular or intercellular, respectively.

10
Biological example - cont.

• We have and
• If we assume that any sequence is either
  extracellular or intercellular, then we have
• P(x) pext P(x ext) pintP(x int).
• By Bayes theorem, we obtain
• the posterior probability that a sequence is
  extracellular.
• (In reality, many transmembrane proteins have
  both intra- and extracellular components and more
  complex models such as HMMs are appropriate.)

11
Probability vs. likelihoodpravdepodobnost vs.
vierohodnost

• If we consider P( X Y ) as a function of X,
  then this is called a probability.
• If we consider P( X Y ) as a function of Y ,
  then this is called a likelihood.

12
Sequence comparison by compression
13
Motivation

• similarity as a marker for homology. And homology
  is used to infer function.
• Sometimes, we are only interested in a numerical
  distance between two sequences. For example, to
  infer a phylogeny.

Figure adapted from http//www.inf.ethz.ch/persona
l/gonnet/acat2000/side2.html
14
Text- vs DNA compression

• compress, gzip or zip routinely used to
  compress text files. They can be applied also to
  a text file containing DNA.
• E.g., a text file F containing chromosome 19 of
  human in fasta format F 61 MB, but
  compress(F) 8.5 MB.
• 8 bits are used for each character. However, DNA
  consists of only 4 different bases 2 bits per
  base are enough A 00, C 01, G 10, and T
  11.
• Applying a standard compression algorithm to a
  file containing DNA encoded using two bits per
  base will usually not be able to compress the
  file further.

15
The repetitive nature of DNA

• Take advantage of the repetitive nature of DNA!!
• LINEs (Long Interspersed Nuclear Elements),
  SINEs.
• UCSC Genome Browser http//genome.ucsc.edu

16
DNA compression

• DNA sequences are very compressible, especially
  for higher eukaryotes they contain many repeats
  of different size, with different numbers of
  instances and different amounts of identity.
• A first idea While processing the DNA string
  from left to right, detect exact repeats and/or
  palindromes (reverse-complemented repeats) that
  possess previous instances in the already
  processed text and encode them by the length and
  position of an earlier occurrence. For stretches
  of sequence for which no significant repeat is
  found, use two-bit encoding. (The program
  Biocompress is based on this idea.)
• Data structure for fast access to sequence
  patterns already encountered.
• Sliding window along unprocessed sequence.

17
DNA compression

• A second idea Build a suffix tree for the whole
  sequence and use it to detect maximal repeats of
  some fixed minimum size. Then code all repeats as
  above and use two-bit encoding for bases not
  contained inside repeats. (The program Cfact is
  based on this idea.)
• Both of these algorithms are lossless, meaning
  that the original sequences can be precisely
  reconstructed from their encodings. An number of
  lossy algorithms exist, which we will not discuss
  here.
• In the following we will discuss the GenCompress
  algorithm due to Xin Chen, Sam Kwong and Ming Li.

18
Edit operations

• The main idea of GenCompress is to use inexact
  matching, followed by edit operations. In other
  words, instances of inexact repeats are encoded
  by a reference to an earlier instance of the
  repeat, followed by some edit operations that
  modify the earlier instance to obtain the current
  instance.
• Three standard edit operations
• Replace (R, i, char) replace the character at
  position i by character char.
• Insert (I, i, char) insert character char at
  position i.
• Delete (D, i) delete the character at position
  i.

19
Edit operations

• different edit operation sequences
• (a) C C C C R C C C C C or (b) C C C C D
  C I C C C C
• g a c c g t c a t t g a c c g t c a
  t t
• g a c c t t c a t t g a c c
  t t c a t t
• infinite number of ways to convert one string
  into another.
• Given two strings q and p. An edit transcript
  ?(q, p) is a list of edit operations that
  transforms q into p.
• E.g., in case (a) the edit transcript is
  ?(gaccgtcatt,gaccttcatt) (R, 4, t),
• whereas in case (b) it is ?
  ?(gaccgtcatt,gaccttcatt) (D, 4), (I, 4, g).
• (positions start at 0 and are given relative to
  current state of the
• string, as obtained by application of any
  preceding edit operations.)

20
Encoding DNA

• Using the two-bit encoding method, gaccttcatt can
  be encoded in 20 bits
• 10 00 01 01 11 11 01 00 11 11
• The following three encode gaccttcatt relative
  to gaccgtcatt
• In the exact matching method we use a pair of
  numbers (repeat - position, repeat - length) to
  represent an exact repeat. We can encode
  gaccttcatt as (0, 4), t, (5, 5), relative to
  gaccgtcatt. Let 4 bits encode an integer, 2 bits
  encode a base and one bit to indicate whether the
  next part is a pair (indicating a repeat) or a
  base. We obtain an encoding in 21 bits
• 0 0000 0100 1 11 0 0101 0101.

21
Encoding DNA

• In the approximate matching method we can encode
  gaccttcatt as (0, 10), (R, 4, t) relative to
  gaccgtcatt. Let use encode R by 00, I by 01, D
  by 11 and use a single bit to indicate whether
  the next part is a pair or a triplet. We obtain
  an encoding in 18 bits
• 0 0000 1010 1 00 0100 11
• For approximate matching, we could also use the
  edit sequence (D, 4), (I, 4, t), for example,
  yielding the relative encoding (0, 10), (D,
  4), (I, 4, t), which uses 25 bits
• 0 0000 1010 1 11 0100 1 01 0100 11.

22
GenCompress

• a one-pass algorithm based on approximate
  matching
• For a given input string w, assume that a part v
  has already been compressed and the remaining
  part is u, with w vu. The algorithm finds an
  optimal prefix p of u that approximately
  matches some substring of v such that p can be
  encoded economically. After outputting the
  encoding of p, remove the prefix p from u and
  append it to v. If no optimal prefix is found,
  output the next base and then move it from u to
  v. Repeat until u ?.

w
v
u
p
p
23
The condition C

• How do we find an optimal prefix p? The
  following condition will be used to limit the
  search.
• Given two numbers k and b. Let p be a prefix of
  the unprocessed part u and q a substring of the
  processed part v. If q gt k, then any transcript
  (q, p) is said to satisfy the condition C (k,
  b) for compression, if its number of edit
  operations is ? b.
• Experimental studies indicate that C (k, b)
  (12, 3) gives good results.
• In other words, when attempting to determine the
  optimal prefix for compression, we will only
  consider repeats of length k that require at most
  b edit operations.

24
The compression gain function

• We define a compression gain function G to
  determine whether a particular approximate repeat
  q, p and edit transcript are beneficial for the
  encoding
• G(q, p, ?) max 2p - (i, q) - wl (q,
  p) - c, 0.
• where
• p is a prefix of the unprocessed part u,
• q is a substring of the processed part v of
  length q that starts at position i,
• 2p is the number of bits that the two-bit
  encoding would use,
• (i, q) is the encoding size of (i, q),
• w is the cost of encoding an edit operation,
• (q, p) is the number of edit operations in (q,
  p),
• and c is the overhead proportional to the size of
  control bits.

25
The GenCompress algorithm

• Input A DNA sequence w, parameter C (k, b)
• Output A lossless compressed encoding of w
• Initialization u w and v e
• while u ¹ e do
• Search for an optimal prefix p of u
• if an optimal prefix p with repeat q in v is
  found then
• Encode the repeat representation (i, q), where
  i is
• the position of q in v, together with the
  shortest edit
• transcript (q, p). Output the code.
• else Set p equal to the first character of u,
• encode and output it.
• Remove the prefix p from u and append it to v
• end

26
Implementing the optimal prefix search

• search for the optimal prefix too slow
• Lemma An optimal prefix p always ends right
  before a mismatch.
• Lemma Let l(q, p) be an optimal edit sequence
  from q to p. If qi is copied onto pj in l, then l
  restricted to (q0i, p0j) is an optimal
  transcript from q0i q0q1 . . . qi to p0j
  p0p1 . . . , pj .
• simplified as follows
• to find an approximate match for p in v, we
  look for an exact match of the first l bases in
  p, where l is a fixed small number
• an integer is stored at each position i of v that
  is determined by the the word of length l
  starting at i.

27
Implementing the optimal prefix search

1. Let w vu where v has already been compressed.
2. Find all occurrences u0l in v, for some small l.
   For each such occurrence, try to extend it,
   allowing mismatches, limited by the above
   observations and condition C. Return the prefix p
   with the largest compression gain value G.

28
Performance of GenCompress

• any nucleotide can be encoded canonically using 2
  bits, we define the compression ratio of a
  compression algorithm as
• I is the number of bases in the input DNA
  sequence
• O is the length (number of bits) of the output
  sequence
• Alternatively, if our DNA string is already
  encoded canonically, we can define the
  compression ratio of a compression algorithm as
• I is the number of bits in the canonical
  encoding of the input DNA sequence and O is the
• length (number of bits) of the output sequence.

29
Performance of GenCompress
30
Recent approaches Encoding of non-repeat regions

• Order-2 arithmetic encoding the adaptive
  probability of a symbol is computed from the
  context (the last 2 symbols) after which it
  appears 3 symbols code one amino-acid???
• Context tree weighting coding (CTW) a tree
  containing all processed substrings of length k
  is built dynamically and each path (string) in
  the tree is weighted by its probability these
  probabilities are used in an arithmetic encoder

Encoded data
input
Arithemtic encoder
CTW model estimator
31
Recent approachesEncoding of numbers

• Fibonacci encoding
• any positive integer can be uniquely expressed as
  the sum of distinct Fibonacci numbers so, that
  no two consecutive Fibonacci numbers are used
• by adding a 1 after the bit corresponding the
  largest Fibonacci number used in the sum the
  representation becomes self-delimited

1 2 3 4 8 18
Fibonacci 11 011 0011 1011 000011 0001011
1-shifted Fib. 1 011 0011 00011 001011 01010011
3-shifted Fib. 001 010 011 100 00011 00001011
1,2,3,5,8,13,21,
32
Recent approachesEncoding of numbers

• k- Shifted Fibonacci encoding
• usually there are many small numbers and few
  large numbers to encode
• n Î 1,,2k 1 normal binary encoding
• n ³ 2k as 0k followed by Fibonacci encoding of
  n (2k 1)

1 2 3 4 8 18
Fibonacci 11 011 0011 1011 000011 0001011
1-shited Fib. 1 011 0011 00011 001011 01010011
3-shited Fib. 001 010 011 100 00011 00001011
1,2,3,5,8,13,21,
33
Recent approaches

• E.g. DNAPack
• uses dynamic programming for selection of
  segments (copied, reversed and/or modified)
• O(n3) still too slow, hence heuristics used

34
Conditional compression

• Given sequence z, compress a sequence w relative
  to z
• Let Compress(w z) denote the length of the
  compression of w, given z. Similarly, let
  Compress(w) Compress(w e), where e denotes
  the empty word. Compress is not the unix
  compression program compress.
• In general, Compress(w z) ¹ Compress(z w).
• For example, the Biocompress-2 program produces
• CompressRatio (brucella rochalima) 55.95,
  and
• CompressRatio (rochalima brucella) 34.56.
• If zw, then Compress(w z) is very small.
• If z and w are completely independent, then
• Compress(w z) Compress(w)

35
Evolutionary distance

• How to define evolutionary distance between
  strings based on conditional compression?
• E.g.
• is used in literature, but it has no good
  reason.
• We need a symmetric measure!
• Þ we can use Kolmogorov complexity

36
Kolmogorov complexity

• Let K(w z) denote the Kolmogorov complexity of
  string w, given z. Informally, this is the length
  of the shortest program that outputs w given z.
  Similarly, set
• K(w) K(w e).
• The following result is due to Kolmogorov and
  Levin
• Theorem Within an additive logarithmic factor,
• K(w z) K(z) K(z w) K(w).
• This implies
• K(w) - K(w z) K(z) - K(z w).
• Normalizing by the sequence length we obtain the
  following symmetric map relatedness

37
A symmetric measure of similarity

• A distance between two sequences w and z
• Unfortunately, K(w) and K(w z) are not
  computable!
• Approximation
• K(w) Compress(w) GenCompress(w)
• K(w z) Compress(w z) Compress(zw) -
  Compress(z)
• GenCompress(zw) - GenCompress(z)

38
Application to genomic sequences
39
Application to genomic sequences
40
Application to genomic sequences

• 16S rRNA for procaryotes corresponds to 18S rRNA
  for eucaryotes

H. butylicus
H. gomorrense
A. urina
H. glauca
R. globiformis
L.sp. nakagiri
U.crescens
41
Finding Regulatory Motifs in DNA Sequences
Micro-array experiments indicate that sets of
genes are regulated by common transcription
factors (TFs). These attach to the DNA upstream
of the coding sequence, at certain binding sites.
Such a site displays a short motif of DNA that is
specific to a given type of TF.
42
Random Sample

• atgaccgggatactgataccgtatttggcctaggcgtacacattagataa
  acgtatgaagtacgttagactcggcgccgccgacccctattttttgag
  cagatttagtgacctggaaaaaaaatttgagtacaaaacttttccgaata
  ctgggcataaggtacatgagtatccctgggatgacttttgggaacact
  atagtgctctcccgatttttgaatatgtaggatcattcgccagggtccga
  gctgagaattggatgaccttgtaagtgttttccacgcaatcgcgaacc
  aacgcggacccaaaggcaagaccgataaaggagatcccttttgcggta
  atgtgccgggaggctggttacgtagggaagccctaacggacttaatggcc
  cacttagtccacttataggtcaatcatgttcttgtgaatggattttta
  actgagggcatagaccgcttggcgcacccaaattcagtgtgggcgagcgc
  aacggttttggcccttgttagaggcccccgtactgatggaaactttca
  attatgagagagctaatctatcgcgtgcgtgttcataacttgagttgg
  tttcgaaaatgctctggggcacatacaagaggagtcttccttatcagtta
  atgctgtatgacactatgtattggcccattggctaaaagcccaacttg
  acaaatggaagatagaatccttgcatttcaacgtatgccgaaccgaaagg
  gaagctggtgagcaacgacagattcttacgtgcattagctcgcttccg
  gggatctaatagcacgaagcttctgggtactgatagca

43
Implanting Motif AAAAAAAGGGGGGG

• atgaccgggatactgatAAAAAAAAGGGGGGGggcgtacacattagataa
  acgtatgaagtacgttagactcggcgccgccgacccctattttttgag
  cagatttagtgacctggaaaaaaaatttgagtacaaaacttttccgaata
  AAAAAAAAGGGGGGGatgagtatccctgggatgacttAAAAAAAAGGG
  GGGGtgctctcccgatttttgaatatgtaggatcattcgccagggtccga
  gctgagaattggatgAAAAAAAAGGGGGGGtccacgcaatcgcgaacc
  aacgcggacccaaaggcaagaccgataaaggagatcccttttgcggta
  atgtgccgggaggctggttacgtagggaagccctaacggacttaatAAAA
  AAAAGGGGGGGcttataggtcaatcatgttcttgtgaatggatttAAA
  AAAAAGGGGGGGgaccgcttggcgcacccaaattcagtgtgggcgagcgc
  aacggttttggcccttgttagaggcccccgtAAAAAAAAGGGGGGGca
  attatgagagagctaatctatcgcgtgcgtgttcataacttgagttAA
  AAAAAAGGGGGGGctggggcacatacaagaggagtcttccttatcagtta
  atgctgtatgacactatgtattggcccattggctaaaagcccaacttg
  acaaatggaagatagaatccttgcatAAAAAAAAGGGGGGGaccgaaagg
  gaagctggtgagcaacgacagattcttacgtgcattagctcgcttccg
  gggatctaatagcacgaagcttAAAAAAAAGGGGGGGa

44
Where is the Implanted Motif?

• atgaccgggatactgataaaaaaaagggggggggcgtacacattagataa
  acgtatgaagtacgttagactcggcgccgccgacccctattttttgag
  cagatttagtgacctggaaaaaaaatttgagtacaaaacttttccgaata
  aaaaaaaagggggggatgagtatccctgggatgacttaaaaaaaaggg
  ggggtgctctcccgatttttgaatatgtaggatcattcgccagggtccga
  gctgagaattggatgaaaaaaaagggggggtccacgcaatcgcgaacc
  aacgcggacccaaaggcaagaccgataaaggagatcccttttgcggta
  atgtgccgggaggctggttacgtagggaagccctaacggacttaataaaa
  aaaagggggggcttataggtcaatcatgttcttgtgaatggatttaaa
  aaaaaggggggggaccgcttggcgcacccaaattcagtgtgggcgagcgc
  aacggttttggcccttgttagaggcccccgtaaaaaaaagggggggca
  attatgagagagctaatctatcgcgtgcgtgttcataacttgagttaa
  aaaaaagggggggctggggcacatacaagaggagtcttccttatcagtta
  atgctgtatgacactatgtattggcccattggctaaaagcccaacttg
  acaaatggaagatagaatccttgcataaaaaaaagggggggaccgaaagg
  gaagctggtgagcaacgacagattcttacgtgcattagctcgcttccg
  gggatctaatagcacgaagcttaaaaaaaaggggggga

45
Implanting Motif AAAAAAGGGGGGG with Four
Mutations

• atgaccgggatactgatAgAAgAAAGGttGGGggcgtacacattagataa
  acgtatgaagtacgttagactcggcgccgccgacccctattttttgag
  cagatttagtgacctggaaaaaaaatttgagtacaaaacttttccgaata
  cAAtAAAAcGGcGGGatgagtatccctgggatgacttAAAAtAAtGGa
  GtGGtgctctcccgatttttgaatatgtaggatcattcgccagggtccga
  gctgagaattggatgcAAAAAAAGGGattGtccacgcaatcgcgaacc
  aacgcggacccaaaggcaagaccgataaaggagatcccttttgcggta
  atgtgccgggaggctggttacgtagggaagccctaacggacttaatAtAA
  tAAAGGaaGGGcttataggtcaatcatgttcttgtgaatggatttAAc
  AAtAAGGGctGGgaccgcttggcgcacccaaattcagtgtgggcgagcgc
  aacggttttggcccttgttagaggcccccgtAtAAAcAAGGaGGGcca
  attatgagagagctaatctatcgcgtgcgtgttcataacttgagttAA
  AAAAtAGGGaGccctggggcacatacaagaggagtcttccttatcagtta
  atgctgtatgacactatgtattggcccattggctaaaagcccaacttg
  acaaatggaagatagaatccttgcatActAAAAAGGaGcGGaccgaaagg
  gaagctggtgagcaacgacagattcttacgtgcattagctcgcttccg
  gggatctaatagcacgaagcttActAAAAAGGaGcGGa

46
Where is the Motif???

• atgaccgggatactgatagaagaaaggttgggggcgtacacattagataa
  acgtatgaagtacgttagactcggcgccgccgacccctattttttgag
  cagatttagtgacctggaaaaaaaatttgagtacaaaacttttccgaata
  caataaaacggcgggatgagtatccctgggatgacttaaaataatgga
  gtggtgctctcccgatttttgaatatgtaggatcattcgccagggtccga
  gctgagaattggatgcaaaaaaagggattgtccacgcaatcgcgaacc
  aacgcggacccaaaggcaagaccgataaaggagatcccttttgcggta
  atgtgccgggaggctggttacgtagggaagccctaacggacttaatataa
  taaaggaagggcttataggtcaatcatgttcttgtgaatggatttaac
  aataagggctgggaccgcttggcgcacccaaattcagtgtgggcgagcgc
  aacggttttggcccttgttagaggcccccgtataaacaaggagggcca
  attatgagagagctaatctatcgcgtgcgtgttcataacttgagttaa
  aaaatagggagccctggggcacatacaagaggagtcttccttatcagtta
  atgctgtatgacactatgtattggcccattggctaaaagcccaacttg
  acaaatggaagatagaatccttgcatactaaaaaggagcggaccgaaagg
  gaagctggtgagcaacgacagattcttacgtgcattagctcgcttccg
  gggatctaatagcacgaagcttactaaaaaggagcgga

47
Why Finding (15,4) Motif is Difficult?

• atgaccgggatactgatAgAAgAAAGGttGGGggcgtacacattagataa
  acgtatgaagtacgttagactcggcgccgccgacccctattttttgag
  cagatttagtgacctggaaaaaaaatttgagtacaaaacttttccgaata
  cAAtAAAAcGGcGGGatgagtatccctgggatgacttAAAAtAAtGGa
  GtGGtgctctcccgatttttgaatatgtaggatcattcgccagggtccga
  gctgagaattggatgcAAAAAAAGGGattGtccacgcaatcgcgaacc
  aacgcggacccaaaggcaagaccgataaaggagatcccttttgcggta
  atgtgccgggaggctggttacgtagggaagccctaacggacttaatAtAA
  tAAAGGaaGGGcttataggtcaatcatgttcttgtgaatggatttAAc
  AAtAAGGGctGGgaccgcttggcgcacccaaattcagtgtgggcgagcgc
  aacggttttggcccttgttagaggcccccgtAtAAAcAAGGaGGGcca
  attatgagagagctaatctatcgcgtgcgtgttcataacttgagttAA
  AAAAtAGGGaGccctggggcacatacaagaggagtcttccttatcagtta
  atgctgtatgacactatgtattggcccattggctaaaagcccaacttg
  acaaatggaagatagaatccttgcatActAAAAAGGaGcGGaccgaaagg
  gaagctggtgagcaacgacagattcttacgtgcattagctcgcttccg
  gggatctaatagcacgaagcttActAAAAAGGaGcGGa

AgAAgAAAGGttGGG
.......
cAAtAAAAcGGcGGG
48
Challenge Problem(Pevzner and Sze)

• Find a motif in a sample of
• - 20 random sequences (e.g. 600 nt
  long)
• - each sequence containing an implanted
• pattern of length 15,
• - each pattern appearing with 4
  mismatches
• as (15,4)-motif.

49
Combinatorial Gene Regulation

• A microarray experiment showed that when gene X
  is knocked out, 20 other genes are not expressed
• How can one gene have such drastic effects?

50
Regulatory Proteins

• Gene X encodes regulatory protein, a.k.a. a
  transcription factor (TF)
• The 20 unexpressed genes rely on gene Xs TF to
  induce transcription
• A single TF may regulate multiple genes

51
Regulatory Regions

• Every gene contains a regulatory region (RR)
  typically stretching 100-1000 bp upstream of the
  transcriptional start site
• Located within the RR are the Transcription
  Factor Binding Sites (TFBS), also known as
  motifs, specific for a given transcription factor
• TFs influence gene expression by binding to a
  specific location in the respective genes
  regulatory region - TFBS
• So finding the same motif in multiple genes
  regulatory regions suggests a regulatory
  relationship amongst those genes

52
Motifs and Transcriptional Start Sites

• A TFBS can be located anywhere within the
  Regulatory Region.
• TFBS may vary slightly across different
  regulatory regions since non-essential bases
  could mutate

53
Transcription Factors and Motifs
54
Motif Logo

• TGGGGGA
• TGAGAGA
• TGGGGGA
• TGAGAGA
• TGAGGGA

• Motifs can mutate on non important bases
• The five motifs in five different genes have
  mutations in position 3 and 5
• Representations called motif logos illustrate the
  conserved and variable regions of a motif

55
Motif Logos An Example
(http//www-lmmb.ncifcrf.gov/toms/sequencelogo.ht
ml)
56
Identifying Motifs Complications

• We do not know the motif sequence
• We do not know where it is located relative to
  the genes start
• Motifs can differ slightly from one gene to the
  next
• How to discern it from random motifs?

57
A Motif Finding Analogy

• The Motif Finding Problem is similar to the
  problem posed by Edgar Allan Poe (1809 1849) in
  his Gold Bug story

58
The Gold Bug Problem

• Given a secret message
• 53!305))64826)4.)4)80648!860))8588
  !83(88)5!
• 46(8896?8)(485)5!2(49562(5-4)88
  4069285))6
• !8)41(94808188148!854)485!52880681(94
  8(884(?3
• 448)4161188?
• Decipher the message encrypted in the fragment

59
Hints for The Gold Bug Problem

• Additional hints
• The encrypted message is in English
• Each symbol correspond to one letter in the
  English alphabet
• No punctuation marks are encoded

60
The Gold Bug Problem Symbol Counts

• Naive approach to solving the problem
• Count the frequency of each symbol in the
  encrypted message
• Find the frequency of each letter in the alphabet
  in the English language
• Compare the frequencies of the previous steps,
  try to find a correlation and map the symbols to
  a letter in the alphabet

61
Symbol Frequencies in the Gold Bug Message

• Gold Bug Message

Symbol 8 4 ) 5 6 ( ! 1 0 2 9 3 ? - .
Frequency 34 25 19 16 15 14 12 11 9 8 7 6 5 5 4 4 3 2 1 1 1

• English Language
• e t a o i n s r h l d c u m f p g w y b v k x j q
  z
• Most frequent
  Least frequent

62
The Gold Bug Message Decoding First Attempt

• By simply mapping the most frequent symbols to
  the most frequent letters of the alphabet
• sfiilfcsoorntaeuroaikoaiotecrntaeleyrcooestvenpin
  elefheeosnlt
• arhteenmrnwteonihtaesotsnlupnihtamsrnuhsnbaoeyent
  acrmuesotorl
• eoaiitdhimtaecedtepeidtaelestaoaeslsueecrnedhimta
  etheetahiwfa
• taeoaitdrdtpdeetiwt
• The result does not make sense

63
The Gold Bug Problem l-tuple count

• A better approach
• Examine frequencies of l-tuples, combinations of
  2 symbols, 3 symbols, etc.
• The is the most frequent 3-tuple in English and
  48 is the most frequent 3-tuple in the
  encrypted text
• Make inferences of unknown symbols by examining
  other frequent l-tuples

64
The Gold Bug Problem the 48 clue

• Mapping the to 48 and substituting all
  occurrences of the symbols
• 53!305))6the26)h.)h)te06the!e60))e5te
  e!e3(ee)5!t
• h6(tee96?te)(the5)t5!2(th9562(5-h)eet
  h0692e5)t)6!e
• )ht1(9the0e1tee1the!e5th)he5!52ee06e1(9the
  t(eeth(?3ht
• he)ht161t1eet?t

65
The Gold Bug Message Decoding Second Attempt

• Make inferences
• 53!305))6the26)h.)h)te06the!e60))e5te
  e!e3(ee)5!t
• h6(tee96?te)(the5)t5!2(th9562(5-h)eet
  h0692e5)t)6!e
• )ht1(9the0e1tee1the!e5th)he5!52ee06e1(9the
  t(eeth(?3ht
• he)ht161t1eet?t
• thet(ee most likely means the tree
• Infer ( r
• th(?3h becomes thr?3h
• Can we guess and ??

66
The Gold Bug Problem The Solution

• After figuring out all the mappings, the final
  message is
• agoodglassinthebishopshostelinthedevilsseatwenyon
  edegreesandt
• hirteenminutesnortheastandbynorthmainbranchseven
  thlimbeastside
• shootfromthelefteyeofthedeathsheadabeelinefromthe
  treethrought
• heshotfiftyfeetout
• Punctuation is important
• A GOOD GLASS IN THE BISHOPS HOSTEL IN THE
  DEVILS SEA,
• TWENY ONE DEGREES AND THIRTEEN MINUTES NORTHEAST
  AND BY NORTH,
• MAIN BRANCH SEVENTH LIMB, EAST SIDE, SHOOT FROM
  THE LEFT EYE OF
• THE DEATHS HEAD A BEE LINE FROM THE TREE
  THROUGH THE SHOT,
• FIFTY FEET OUT.

67
Solving the Gold Bug Problem

• Prerequisites to solve the problem
• Need to know the relative frequencies of single
  letters, and combinations of two and three
  letters in English
• Knowledge of all the words in the English
  dictionary is highly desired to make accurate
  inferences

68
Motif Finding and The Gold Bug Problem
Similarities

• Nucleotides in motifs encode for a message in the
  genetic language. Symbols in The Gold Bug
  encode for a message in English
• In order to solve the problem, we analyze the
  frequencies of patterns in DNA/Gold Bug message.
• Knowledge of established regulatory motifs makes
  the Motif Finding problem simpler. Knowledge of
  the words in the English dictionary helps to
  solve the Gold Bug problem.

69
Similarities (contd)

• Motif Finding
• In order to solve the problem, we analyze the
  frequencies of patterns in the nucleotide
  sequences
• Gold Bug Problem
• In order to solve the problem, we analyze the
  frequencies of patterns in the text written in
  English

70
Similarities (contd)

• Motif Finding
• Knowledge of established motifs reduces the
  complexity of the problem
• Gold Bug Problem
• Knowledge of the words in the dictionary is
  highly desirable

71
Motif Finding and The Gold Bug Problem
Differences

• Motif Finding is harder than Gold Bug problem
• We dont have the complete dictionary of motifs
• The genetic language does not have a standard
  grammar
• Only a small fraction of nucleotide sequences
  encode for motifs the size of data is enormous

72
The Motif Finding Problem

• Given a random sample of DNA sequences
• cctgatagacgctatctggctatccacgtacgtaggtcctctgtgcgaat
  ctatgcgtttccaaccat
• agtactggtgtacatttgatacgtacgtacaccggcaacctgaaacaaac
  gctcagaaccagaagtgc
• aaacgtacgtgcaccctctttcttcgtggctctggccaacgagggctgat
  gtataagacgaaaatttt
• agcctccgatgtaagtcatagctgtaactattacctgccacccctattac
  atcttacgtacgtataca
• ctgttatacaacgcgtcatggcggggtatgcgttttggtcgtcgtacgct
  cgatcgttaacgtacgtc
• Find the pattern that is implanted in each of the
  individual sequences, namely, the motif

73
The Motif Finding Problem (contd)

• Additional information
• The hidden sequence is of length 8
• The pattern is not exactly the same in each array
  because random point mutations may occur in the
  sequences

74
The Motif Finding Problem (contd)

• The patterns revealed with no mutations
• cctgatagacgctatctggctatccacgtacgtaggtcctctgtgcgaat
  ctatgcgtttccaaccat
• agtactggtgtacatttgatacgtacgtacaccggcaacctgaaacaaac
  gctcagaaccagaagtgc
• aaacgtacgtgcaccctctttcttcgtggctctggccaacgagggctgat
  gtataagacgaaaatttt
• agcctccgatgtaagtcatagctgtaactattacctgccacccctattac
  atcttacgtacgtataca
• ctgttatacaacgcgtcatggcggggtatgcgttttggtcgtcgtacgct
  cgatcgttaacgtacgtc
• acgtacgt
• Consensus String

75
The Motif Finding Problem (contd)

• The patterns with 2 point mutations
• cctgatagacgctatctggctatccaGgtacTtaggtcctctgtgcgaat
  ctatgcgtttccaaccat
• agtactggtgtacatttgatCcAtacgtacaccggcaacctgaaacaaac
  gctcagaaccagaagtgc
• aaacgtTAgtgcaccctctttcttcgtggctctggccaacgagggctgat
  gtataagacgaaaatttt
• agcctccgatgtaagtcatagctgtaactattacctgccacccctattac
  atcttacgtCcAtataca
• ctgttatacaacgcgtcatggcggggtatgcgttttggtcgtcgtacgct
  cgatcgttaCcgtacgGc

76
The Motif Finding Problem (contd)

• The patterns with 2 point mutations
• cctgatagacgctatctggctatccaGgtacTtaggtcctctgtgcgaat
  ctatgcgtttccaaccat
• agtactggtgtacatttgatCcAtacgtacaccggcaacctgaaacaaac
  gctcagaaccagaagtgc
• aaacgtTAgtgcaccctctttcttcgtggctctggccaacgagggctgat
  gtataagacgaaaatttt
• agcctccgatgtaagtcatagctgtaactattacctgccacccctattac
  atcttacgtCcAtataca
• ctgttatacaacgcgtcatggcggggtatgcgttttggtcgtcgtacgct
  cgatcgttaCcgtacgGc

Can we still find the motif, now that we have 2
mutations?
77
Defining Motifs

• To define a motif, let us say we know where the
  motif starts in the sequence
• The motif start positions in their sequences can
  be represented as s (s1,s2,s3,,st)

78
Motifs Profiles and Consensus

• a G g t a c T t
• C c A t a c g t
• Alignment a c g t T A g t
• a c g t C c A t
• C c g t a c g G
• 
  _________________
• A 3 0 1 0 3 1 1 0
• Profile C 2 4 0 0 1 4 0 0
• G 0 1 4 0 0 0 3 1
• T 0 0 0 5 1 0 1 4
• _________________
• Consensus A C G T A C G T

• Line up the patterns by their start indexes
• s (s1, s2, , st)
• Construct matrix profile with frequencies of each
  nucleotide in columns
• Consensus nucleotide in each position has the
  highest score in column

79
Consensus

• Think of consensus as an ancestor motif, from
  which mutated motifs emerged
• The distance between a real motif and the
  consensus sequence is generally less than that
  for two real motifs

80
Evaluating Motifs

• We have a guess about the consensus sequence, but
  how good is this consensus?
• Need to introduce a scoring function to compare
  different guesses and choose the best one.
• t number of sample DNA sequences
• n length of each DNA sequence
• DNA sample of DNA sequences (t x n array)
• l length of the motif (l-mer)
• si starting position of an l-mer in sequence i
• s(s1, s2, st) array of motifs starting
  positions

81
Parameters

• cctgatagacgctatctggctatccaGgtacTtaggtcctctgtgcgaa
  tctatgcgtttccaaccat
• agtactggtgtacatttgatCcAtacgtacaccggcaacctgaaacaaa
  cgctcagaaccagaagtgc
• aaacgtTAgtgcaccctctttcttcgtggctctggccaacgagggctga
  tgtataagacgaaaatttt
• agcctccgatgtaagtcatagctgtaactattacctgccacccctatta
  catcttacgtCcAtataca
• ctgttatacaacgcgtcatggcggggtatgcgttttggtcgtcgtacgc
  tcgatcgttaCcgtacgGc

DNA
t5
n 69
s1 26 s2 21 s3 3 s4 56 s5
60
s
82
Scoring Motifs
l

• Given s (s1, st ) and DNA
• Score(s,DNA)

• a G g t a c T t
• C c A t a c g t
• a c g t T A g t
• a c g t C c A t
• C c g t a c g G
• _________________
• A 3 0 1 0 3 1 1 0
• C 2 4 0 0 1 4 0 0
• G 0 1 4 0 0 0 3 1
• T 0 0 0 5 1 0 1 4
• _________________
• Consensus a c g t a c g t
• Score 3445343430

t
83
The Motif Finding Problem

• If starting positions s(s1, s2, st) are given,
  finding consensus is easy even with mutations in
  the sequences because we can simply construct the
  profile to find the motif (consensus)
• But the starting positions s are usually not
  given. How can we find the best profile matrix?

84
The Motif Finding Problem Formulation

• Goal Given a set of DNA sequences, find a set of
  l-mers, one from each sequence, that maximizes
  the consensus score
• Input A t x n matrix of DNA, and l, the length
  of the pattern to find
• Output An array of l starting positions s
  (s1, s2, st) maximizing Score(s,DNA)

85
The Motif Finding Problem Brute Force Solution

• Compute the scores for each possible combination
  of starting positions s
• The best score will determine the best profile
  and the consensus pattern in DNA
• The goal is to maximize Score(s,DNA) by varying
  the starting positions si, where

• si 1, , n - l 1
• i 1, , t

86
BruteForceMotifSearch

• BruteForceMotifSearch(DNA, t, n, l)
• bestScore ? 0
• for each s(s1,s2 , . . ., st) from (1,1 . . . 1)
  to (n- l 1, . . ., n- l 1)
• if (Score (s,DNA) gt bestScore )
• bestScore ? score(s, DNA )
• bestMotif ? (s1,s2 , . . . , st )
• return bestMotif

87
Running Time of BruteForceMotifSearch

• Varying (n - l 1) positions in each of t
  sequences, we are looking at (n - l 1) t sets
  of starting positions
• For each set of starting positions, the scoring
  function makes l operations, so complexity is l
  (n l 1) t O(l n t )
• That means that for t 8, n 1000, l 10 we
  must perform approximately 10 20 computations
  it will take billions years

88
The Median String Problem

• Given a set of t DNA sequences find a pattern
  that appears in all t sequences with the minimum
  number of mutations
• This pattern will be the motif

89
Hamming Distance

• Hamming distance
• dH(v,w) is the number of nucleotide pairs that do
  not match when v and w are aligned. For example
• dH(AAAAAA, ACAAAC) 2

90
Total Distance An Example

• Given v acgtacgt and s
• acgtacgt
• cctgatagacgctatctggctatccacgtacAtaggtcctctgtgcgaat
  ctatgcgtttccaaccat
• acgtacgt
• agtactggtgtacatttgatacgtacgtacaccggcaacctgaaacaaac
  gctcagaaccagaagtgc
• acgtacgt
• aaaAgtCcgtgcaccctctttcttcgtggctctggccaacgagggctgat
  gtataagacgaaaatttt
• 
  acgtacgt
• agcctccgatgtaagtcatagctgtaactattacctgccacccctattac
  atcttacgtacgtataca
• 
  acgtacgt
• ctgttatacaacgcgtcatggcggggtatgcgttttggtcgtcgtacgct
  cgatcgttaacgtaGgtc
• v is the sequence in red, x is the sequence in
  blue
• TotalDistance(v,DNA) 10201 4

dH(v, x) 1
dH(v, x) 0
dH(v, x) 0
dH(v, x) 2
dH(v, x) 1
91
Total Distance Definition

• For each DNA sequence i, compute all dH(v, x),
  where x is an l-mer with starting position si
• (1 lt si lt n l 1)
• Find minimum of dH(v, x) among all l -mers in
  sequence i
• TotalDistance(v,DNA) is the sum of the minimum
  Hamming distances for each DNA sequence i
• TotalDistance(v,DNA) mins dH(v, s), where s is
  the set of starting positions s1, s2, st

92
The Median String Problem Formulation

• Goal Given a set of DNA sequences, find a median
  string
• Input A t x n matrix DNA, and l, the length of
  the pattern to find
• Output A string v of l nucleotides that
  minimizes TotalDistance(v,DNA) over all strings
  of that length

93
Median String Search Algorithm

• MedianStringSearch (DNA, t, n, l)
• bestWord ? AAAA
• bestDistance ? 8
• for each l -mer s from AAAA to TTTT
• if TotalDistance(s,DNA) lt bestDistance
• bestDistance?TotalDistance(s,DNA)
• bestWord ? s
• return bestWord

94
Motif Finding Problem Median String Problem

• The Motif Finding is a maximization problem while
  Median String is a minimization problem
• However, the Motif Finding problem and Median
  String problem are computationally equivalent
• Need to show that minimizing TotalDistance is
  equivalent to maximizing Score

95
We are looking for the same thing
l

• At any column iScorei TotalDistancei t
• Because there are l columns
• Score TotalDistance l t
• Rearranging
• Score l t - TotalDistance
• l t is constant the minimization of the right
  side is equivalent to the maximization of the
  left side

• a G g t a c T t
• C c A t a c g t
• Alignment a c g t T A g t
• a c g t C c A t
• C c g t a c g G
• _________________
• A 3 0 1 0 3 1 1 0
• Profile C 2 4 0 0 1 4 0 0
• G 0 1 4 0 0 0 3 1
• T 0 0 0 5 1 0 1 4
• _________________
• Consensus a c g t a c g t
• Score 34453434
• TotalDistance 21102121

t
96
Motif Finding Problem vs. Median String Problem

• Why bother reformulating the Motif Finding
  problem into the Median String problem?
• The Motif Finding Problem needs to examine all
  the combinations for s. That is (n - l 1)t
  combinations!!!
• The Median String Problem needs to examine all 4l
  combinations for v. This number is relatively
  smaller

97
Motif Finding Improving the Running Time

• Recall the BruteForceMotifSearch
• BruteForceMotifSearch(DNA, t, n, l)
• bestScore ? 0
• for each s(s1,s2 , . . ., st) from (1,1 . . .
  1) to (n- l 1, . . ., n- l1)
• if (Score(s,DNA) gt bestScore)
• bestScore ? Score(s, DNA)
• bestMotif ? (s1,s2 , . . . , st)
• return bestMotif
• Branch-bound searching

98
Structuring the Search

• How can we perform the line
• for each s(s1,s2 , . . ., st) from (1,1 . . . 1)
  to (n-l1, . . ., n-l1) ?
• We need a method for efficiently structuring and
  navigating the many possible motifs
• This is not very different than exploring all
  t-digit numbers

99
Median String Improving the Running Time

• MedianStringSearch (DNA, t, n, l)
• bestWord ? AAAA
• bestDistance ? 8
• for each l-mer s from AAAA to TTTT
• if TotalDistance(s,DNA) lt bestDistance
• bestDistance ?TotalDistance(s,DNA)
• bestWord ? s
• return bestWord

100
Structuring the Search

• For the Median String Problem we need to consider
  all 4l possible l-mers
• aa aa
• aa ac
• aa ag
• aa at
• .
• .
• tt tt
• How to organize this search?

l
101
Alternative Representation of the Search Space

• Let A 1, C 2, G 3, T 4
• Then the sequences from AAA to TTT become
• 1111
• 1112
• 1113
• 1114
• .
• .
• 4444
• Notice that the sequences above simply list all
  numbers as if we were counting on base 4 without
  using 0 as a digit

l
102
Linked List

• Suppose l 2
• aa ac ag at ca cc cg ct ga gc gg gt
  ta tc tg tt
• Need to visit all the predecessors of a sequence
  before visiting the sequence itself

Start
103
Linked List (contd)

• Linked list is not the most efficient data
  structure for motif finding
• Lets try grouping the sequences by their
  prefixes
• aa ac ag at ca cc cg ct ga gc gg gt
  ta tc tg tt

104
Search Tree

• a- c- g-
  t-
• aa ac ag at ca cc cg ct ga gc gg gt
  ta tc tg tt

root
--
105
Analyzing Search Trees

• Characteristics of the search trees
• The sequences are contained in its leaves
• The parent of a node is the prefix of its
  children
• How can we move through the tree?

106
Moving through the Search Trees

• Four common moves in a search tree that we are
  about to explore
• Move to the next leaf
• Visit all the leaves
• Visit the next node
• Bypass the children of a node

107
Visit the Next Leaf
Given a current leaf a , we need to compute the
next leaf

• NextLeaf( a,L, k ) // a the array of
  digits
• for i ? L to 1 // L length of the
  array
• if ai lt k // k max digit
  value
• ai ? ai 1
• return a
• ai ? 1
• return a

108
NextLeaf (contd)

• The algorithm is common addition in radix k
• Increment the least significant digit
• Carry the one to the next digit position when
  the digit is at maximal value

109
NextLeaf Example

• Moving to the next leaf
• 1- 2- 3-
  4-
• 11 12 13 14 21 22 23 24 31 32 33 34
  41 42 43 44

--
Current Location
110
NextLeaf Example (contd)

• Moving to the next leaf
• 1- 2- 3-
  4-
• 11 12 13 14 21 22 23 24 31 32 33 34
  41 42 43 44

--
Next Location
111
Visit All Leaves

• Printing all permutations in ascending order
• AllLeaves(L,k) // L length of the sequence
• a ? (1,...,1) // k max digit value
• while forever // a array of digits
• output a
• a ? NextLeaf(a,L,k)
• if a (1,...,1)
• return

112
Visit All Leaves Example

• Moving through all the leaves in order
• 1- 2- 3-
  4-
• 11 12 13 14 21 22 23 24 31 32 33 34
  41 42 43 44
• 1 2 3 4 5 6 7 8 9
  10 11 12 13 14 15

--
Order of steps
113
Depth First Search

• So we can search leaves
• How about searching all vertices of the tree?
• We can do this with a depth first search

114
Visit the Next Vertex

• NextVertex(a,i,L,k) // a the array of
  digits
• if i lt L // i prefix
  length
• a i1 ? 1 // L max length
• return ( a,i 1) // k max digit value
• else
• for j ? l to 1
• if aj lt k
• aj ? aj 1
• return( a,j )
• return(a,0)

115
Example

• Moving to the next vertex
• 1- 2- 3-
  4-
• 11 12 13 14 21 22 23 24 31 32 33 34
  41 42 43 44

Current Location
--
116
Example

• Moving to the next vertices
• 1- 2- 3-
  4-
• 11 12 13 14 21 22 23 24 31 32 33 34
  41 42 43 44

Location after 5 next vertex moves
--
117
Bypass Move

• Given a prefix (internal vertex), find next
  vertex after skipping all its children
• Bypass(a,i,L,k) // a array of digits
• for j ? i to 1 // i prefix length
• if aj lt k // L maximum length
• aj ? aj 1 // k max digit value
• return(a,j)
• return(a,0)

118
Bypass Move Example

• Bypassing the descendants of 2-
• 1- 2- 3-
  4-
• 11 12 13 14 21 22 23 24 31 32 33 34
  41 42 43 44

Current Location
--
119
Example

• Bypassing the descendants of 2-
• 1- 2- 3-
  4-
• 11 12 13 14 21 22 23 24 31 32 33 34
  41 42 43 44

Next Location
--
120
Revisiting Brute Force Search

• Now that we have method for navigating the tree,
  lets look again at BruteForceMotifSearch

121
Brute Force Search Again

• BruteForceMotifSearchAgain(DNA, t, n, l)
• s ? (1,1,, 1)
• bestScore ? Score(s,DNA)
• while forever
• s ? NextLeaf (s, t, n- l 1)
• if (Score(s,DNA) gt bestScore)
• bestScore ? Score(s, DNA)
• bestMotif ? (s1,s2 , . . . , st)
• return bestMotif

122
Can We Do Better?

• Sets of s(s1, s2, ,st) may have a weak profile
  for the first i positions (s1, s2, ,si)
• Every row of alignment may add at most l to Score
• Optimism if all subsequent (t-i) positions
  (si1, st) add
• (t i ) l to Score(s,i,DNA)
• If Score(s,i,DNA) (t i ) l lt BestScore, it
  makes no sense to search in vertices of the
  current subtree
• Use ByPass()

123
Branch and Bound Algorithm for Motif Search

• Since each level of the tree goes deeper into
  search, discarding a prefix discards all
  following branches
• This saves us from looking at (n l 1)t-i
  leaves
• Use NextVertex() and ByPass() to navigate the
  tree

124
Pseudocode for Branch and Bound Motif Search

• BranchAndBoundMotifSearch(DNA,t,n,l)
• s ? (1,,1)
• bestScore ? 0
• i ? 1
• while i gt 0
• if i lt t
• optimisticScore ? Score(s, i, DNA) (t i )
  l
• if optimisticScore lt bestScore
• (s, i) ? Bypass(s,i, n- l 1)
• else
• (s, i) ? NextVertex(s, i, n- l 1)
• else
• if Score(s,DNA) gt bestScore
• bestScore ? Score(s)
• bestMotif ? (s1, s2, s3, , st)
• (s,i) ? NextVertex(s,i,t,n- l 1)
• return bestMotif

125
Median String Search Improvements

• Recall the computational differences between
  motif search and median string search
• The Motif Finding Problem needs to examine all
  (n-l1)t combinations for s.
• The Median String Problem needs to examine 4l
  combinations of v. This number is relatively
  small
• We want to use median string algorithm with the
  Branch and Bound trick!

126
Branch and Bound Applied to Median String Search

• Note that if the total distance for a prefix is
  greater than that for the best word so far
• TotalDistance (prefix, DNA ) gt BestDistance
• there is no use exploring the remaining part of
  the word
• We can eliminate that branch and BYPASS exploring
  that branch further

127
Bounded Median String Search

• BranchAndBoundMedianStringSearch(DNA,t,n,l )
• s ? (1,,1)
• bestDistance ? 8
• i ? 1
• while i gt 0
• if i lt l
• prefix ? string corresponding to the
  first i nucleotides of s
• optimisticDistance ? TotalDistance(prefix,D
  NA)
• if optimisticDistance gt bestDistance
• (s, i ) ? Bypass(s,i, l, 4)
• else
• (s, i ) ? NextVertex(s, i, l, 4)
• else
• word ? nucleotide string corresponding to s
• if TotalDistance(s,DNA) lt bestDistance
• bestDistance ? TotalDistance(word, DNA)
• bestWord ? word
• (s,i ) ? NextVertex(s,i, l, 4)
• return bestWord

128
Improving the Bounds

• Given an l-mer w, divided into two parts at point
  i
• u prefix w1, , wi,
• v suffix wi1, ..., wl
• Find minimum distance for u in a sequence
• No instances of u in the sequence have distance
  less than the minimum distance
• Note this doesnt tell us anything about whether
  u is part of any motif. We only get a minimum
  distance for prefix u

129
Improving the Bounds (contd)

• Repeating the process for the suffix v gives us a
  minimum distance for v