CS 5263 Bioinformatics Lecture 18 Motif finding What

CS 5263 Bioinformatics Lecture 18 Motif finding

What is a (biological) motif? • A motif is a recurring fragment, theme or pattern • Sequence motif: a sequence pattern of nucleotides in a DNA sequence or amino acids in a protein • Structural motif: a pattern in a protein structure formed by the spatial arrangement of amino acids. • Network motif: patterns that occur in different parts of a network at frequencies much higher than those found in randomized network • Commonality: – higher frequency than would be expected by chance – Has, or is conjectured to have, a biological significance

(Sequence) motif finding • Given a set of sequences • Goal: find sequence motifs that appear in all or the majority of the sequences, and are likely associated with some functions – In DNA: regulatory sequences – In protein: functional/structural domains

Roadmap • • Biological background Representation of motifs Algorithms for finding motifs Other issues – Distinguish functional vs non-functional motifs – Search for instances of given motifs – Interpretation of motifs

• In motif finding, understanding the motivations, significance of the problems, difficulties, and ideas that have been explored are more important than knowing the details of the existing algorithms! • Most algorithms often perform poorly in real challenges! – Not necessarily a fault of algorithm designers • Algorithms will be improved

Biological background for motif finding

Cells respond to environment Various external messages Heat Responds to environmental conditions Food Supply

Genome is fixed – Cells are dynamic • A genome is static – Every cell in our body has a copy of same genome • A cell is dynamic – Responds to external conditions – Most cells follow a cell cycle of division • Cells differentiate during development

Gene regulation • … is responsible for the dynamic cell • Gene expression (production of protein) varies according to: – – Cell type Cell cycle External conditions Location

Where gene regulation takes place • Opening of chromatin • Transcription • Translation • Protein stability • Protein modifications

Transcriptional Regulation • Strongest regulation happens during transcription • Best place to regulate: No energy wasted making intermediate products • However, slowest response time After a receptor notices a change: 1. Cascade message to nucleus 2. Open chromatin & bind transcription factors 3. Recruit RNA polymerase and transcribe 4. Splice m. RNA and send to cytoplasm 5. Translate into protein

Transcription Factors Binding to DNA Transcriptional regulation: • Certain transcription factors bind to DNA Binding recognizes DNA substrings: • Regulatory motifs

Regulation of Genes Transcription Factor (TF) (Protein) RNA polymerase (Protein) DNA Promoter Gene

Regulation of Genes Transcription Factor (TF) (Protein) RNA polymerase (Protein) DNA Regulatory Element, TF binding site, TF binding motif, cis-regulatory motif (element) Gene

Regulation of Genes Transcription Factor (Protein) RNA polymerase DNA Regulatory Element Gene

Regulation of Genes New protein Transcription Factor DNA Regulatory Element Gene RNA polymerase

The Cell as a Regulatory Network If C then D gene D A B C Make D If B then NOT D If A and B then D D gene B D C Make B If D then B

Code for protein-DNA binding? Some knowledge exists

However, overall code still missing

Experimental methods • DNase footprinting

Experimental methods • To determine protein-DNA binding site is tedious and time-consuming • To determine the binding specificity is even harder – Involves mutating different combinations of nucleic acids in promoter region and observe the biological effects • Computational methods can help

Finding Regulatory Motifs . . . Given a collection of genes that are believed to be regulated by the same protein, Find the common TF-binding motif from promoters

Essentially a Multiple Local Alignment . . . • Find “best” multiple local alignment

• Then why don’t we just use multiple sequence alignment algorithms like the Multidimensional Dynamic Programming?

Characteristics of Regulatory Motifs • Tiny (6 -12 bp) • Intergenic regions are very long • Highly Variable • ~Constant Size – Because a constant-size transcription factor binds • Often repeated • Often conserved

Motif Representation

Motif representation • Collection of exact words – {ACGTTAC, ACGCTAC, AGGTGAC, …} • Consensus sequence (with wild cards) – {Ac. GTg. Tt. AC} – {ASGTKTKAC} S=C/G, K=G/T (IUPAC code) • Position specific weight matrices

Position Specific Weight Matrix 1 2 3 4 5 6 7 8 9 A . 97 . 10 . 02 . 03 . 10 . 01 . 05 . 85 . 03 C . 01 . 40 . 01 . 04 . 05 . 01 . 05 . 03 G . 01 . 40 . 95 . 03 . 40 . 01 . 3 . 05 . 03 T . 01 . 10 . 02 . 90 . 45 . 97 . 6 . 05 . 91 S G T K A A C

frequency Sequence Logo A C G 1 2 3 4 5 6 7 8 9. 97. 10. 02. 03. 10. 01. 05. 85. 03. 01. 40. 01. 04. 05. 01. 05. 03. 01. 40. 95. 03. 40. 01. 3. 05. 03 T . 01. 10. 02. 90. 45. 97 . 6 . 05. 91

Sequence Logo A C G 1 2 3 4 5 6 7 8 9. 97. 10. 02. 03. 10. 01. 05. 85. 03. 01. 40. 01. 04. 05. 01. 05. 03. 01. 40. 95. 03. 40. 01. 3. 05. 03 T . 01. 10. 02. 90. 45. 97 . 6 . 05. 91

Entropy and information content • Entropy: a measure of uncertainty • The entropy of a random variable X that can assume the n different values x 1, x 2, . . . , xn with the respective probabilities p 1, p 2, . . . , pn is defined as

Entropy and information content • Example: A, C, G, T with equal probability § H = 4 * (-0. 25 log 2 0. 25) = log 2 4 = 2 bits § Need 2 bits to encode (e. g. 00 = A, 01 = C, 10 = G, 11 = T) § Maximum uncertainty • 50% A and 50% C: § H = 2 * (-0. 5 log 2 0. 5) = log 2 2 = 1 bit • 100% A § H = 1 * (-1 log 2 1) = 0 bit § Minimum uncertainty • Information: the opposite of uncertainty § I = maximum uncertainty – entropy § The above examples provide 0, 1, and 2 bits of information, respectively

Entropy and information content A C G 1. 97. 01 2. 10. 40 3. 02. 01. 95 4. 03. 04. 03 5. 10. 05. 40 6. 01. 01 7. 05. 3 T H I 8. 85. 05 9. 03. 03 . 01. 10. 02. 90. 45. 97. 6. 05. 91. 24 1. 72. 36. 63 1. 60 0. 24 1. 40 0. 85 0. 58 1. 76 0. 28 1. 64 1. 37 0. 40 1. 76 0. 60 1. 15 1. 42 Mean 1. 15 Total 10. 4 Expected occurrence in random DNA: 1 / 210. 4 = 1 / 1340 Expected occurrence of an exact 5 -mer: 1 / 210 = 1 / 1024

Sequence Logo 1 2 3 4 5 6 7 8 9 A C G . 97 . 10 . 02 . 03 . 10 . 01 . 05 . 85 . 03 . 01 . 40 . 01 . 04 . 05 . 01 . 05 . 03 . 01 . 40 . 95 . 03 . 40 . 01 . 3 . 05 . 03 T I . 01 . 10 . 02 . 90 . 45 . 97 . 6 . 05 . 91 1. 76 0. 28 1. 64 1. 37 0. 40 1. 76 0. 60 1. 15 1. 42

Background-normalized Seq Logo • Many genomes have skewed base distribution • In a thermophilic bacteria (i. e. living in a hot spring), GC content can be as high as 70%. • Thus a motif ATAT in the genome of a thermophilic bacteria would contain more information than a motif GCGC

Relative Entropy • Definition 6. 1. Let P and Q be two probability measures on the same alphabet X. Then the relative entropy (information divergence, Kullback-Leibler distance, discrimination) from P to Q is defined as • Easy to prove that if Q is a uniform distribution, D(P || Q) is equal to the Information content of P

Relative Entropy • Background: p. A = p. T = 0. 2, p. C = p. G = 0. 3 • Distribution on some column of a PWM: Case 1: p. A = 0. 85, p. C = p. G = p. T = 0. 05 Case 2: p. G = 0. 85 p. C = p. A = p. T = 0. 05 • Assuming uniform background distribution: I 1 = I 2 = 1. 15 • With the non-uniform background distribution: – D 1 = 1. 42 – D 2 = 0. 95

Background-normalized Seq Logo 1 2 3 4 5 6 7 8 9 A C G . 97 . 10 . 02 . 03 . 10 . 01 . 05 . 85 . 03 . 01 . 40 . 01 . 04 . 05 . 01 . 05 . 03 . 01 . 40 . 95 . 03 . 40 . 01 . 3 . 05 . 03 T I I’ . 01 . 10 . 02 . 90 . 45 . 97 . 6 . 05 . 91 1. 76 0. 28 1. 64 1. 37 0. 40 1. 76 0. 60 1. 15 1. 42 2 . 13 1. 35 1. 6 0. 45 2 . 70 1. 37 1. 65

Physical interpretation • Information content is reversely proportional to the binding energy – High information content => lower energy => high affinity of binding • Relative entropy represents the specificity of the binding sites compared to random DNA sequences

Real example • E. coli. Promoter • “TATA-Box” ~ 10 bp upstream of transcription start • TACGAT • TAAAAT • TATACT Consensus: TATAAT • GATAAT • TATGAT Note: none of the instances • TATGTT matches the consensus perfectly

Finding Motifs

Definitions of terms • Motif: a consensus sequence or a PWM • Pattern: alias for motif (used in combinatorial motif finding) • Instance of a motif: a substring of a sequence that “matches” to the motif – How to define “match” will be shown later

Motif finding schemes Conservation Yes No Whole Yes Genome 1 & 2 & 3 Genome 1 Dictionary building genome Phylogenetic footprinting No Gene 1 A & 1 B & 1 C or “Motif finding” Gene Set 1 & 2 & 3 1 A 1 B 1 C Gene set 1 Gene set 2 Gene set 3 Genome 2 Genome 3 Genome 1 Ideally, all information should be used, at some stage. i. e. , inside algorithm vs pre- or post-processing.

Classification of approaches • Combinatorial search – Based on enumeration of words and computing word similarities – Analogy to DP for sequence alignment • Probabilistic modeling – Construct models to distinguish motifs vs nonmotifs – Analogy to HMM for sequence alignment

Combinatorial motif finding • Idea 1: find all k-mers that appeared at least m times • Idea 2: find all k-mers that are statistically significant • Problem: most motifs allow divergence. Each variation may only appear once. • Idea 3: find all k-mers, considering IUPAC code – e. g. ASGTKTKAC, S = C/G, K = G/T – Still inflexible • Idea 4: find k-mers that approximately appeared at least m times – i. e. allow some mismatches

Combinatorial motif finding Given a set of sequences S = {x 1, …, xn} • A motif W is a consensus string w 1…w. K • Find motif W* with “best” match to x 1, …, xn Definition of “best”: d(W, xi) = min hamming dist. between W and a word in xi d(W, S) = i d(W, xi) W* = argmin( d(W, S) )

Exhaustive searches 1. Pattern-driven algorithm: For W = AA…A to TT…T (4 K possibilities) Find d( W, S ) Report W* = argmin( d(W, S) ) Running time: O( K N 4 K ) (where N = i |xi|) Guaranteed to find the optimal solution.

Exhaustive searches 2. Sample-driven algorithm: For W = a K-long word in some xi Find d( W, S ) Report W* = argmin( d( W, S ) ) OR Report a local improvement of W* Running time: O( K N 2 )

Exhaustive searches • Problem with sample-driven approach: • If: – True motif does not occur in data, and – True motif is “weak” • Then, – random strings may score better than any instance of true motif

Example • E. coli. Promoter • “TATA-Box” ~ 10 bp upstream of transcription start • TACGAT • TAAAAT • TATACT Consensus: TATAAT • GATAAT Each instance differs at most 2 • TATGAT bases from the consensus • TATGTT None of the instances matches the consensus perfectly

Heuristic methods • Cannot afford exhaustive search on all patterns • Sample-driven approaches may miss real patterns • However, a real pattern should not differ too much from its instances in S • Start from the space of all words in S, extend to the space with real patterns

Some of the popular tools • Consensus (Hertz & Stormo, 1999) • WINNOWER (Pevzner & Sze, 2000) • MULTIPROFILER (Keich & Pevzner, 2002) • PROJECTION (Buhler & Tompa, 2001) • WEEDER (Pavesi et. al. 2001) • And a dozen of others

Consensus Algorithm: Cycle 1: For each word W in S For each word W’ in S Create alignment (gap free) of W, W’ Keep the C 1 best alignments, A 1, …, AC 1 ACGGTTG , ACGCCTG , CGAACTT , AGAACTA , GGGCTCT … GGGGTGT …

Algorithm (cont’d): Cycle i: For each word W in S For each alignment Aj from cycle i-1 Create alignment (gap free) of W, Aj Keep the Ci best alignments A 1, …, ACi

• C 1, …, Cn are user-defined heuristic constants Running time: O(k. N 2) + O(k. N C 1) + O(k. N C 2) + … + O(k. N Cn) = O(k. N 2 + NCtotal) Where Ctotal = i Ci, typically O(n. C), where C is a big constant

Extended sample-driven (ESD) approaches • Hybrid between pattern-driven and sample-driven • Assume each instance does not differ by more than α bases to the motif ( usually depends on k) motif instance α-neighborhood The real motif will reside in the neighborhood of some words in S. Instead of searching all 4 K patterns, we can search the -neighborhood of every word in S.

WEEDER • Naïve: N Kα 3α NK # of patterns to test # of words in sequences

Better idea • Using a joint suffix tree, find all patterns that: – Have length K – Appeared in at least m sequences with at most α mismatches • Post-processing

WEEDER: algorithm sketch Current pattern P, |P| < K # mismatches (e, B) Seq occ A C G T T • A list containing all eligible nodes: with at most α mismatches to P • For each node, remember #mismatches accumulated (e), and bit vector (B) for seq occ, e. g. [011100010] • Bit OR all B’s to get seq occurrence for P • Suppose #occ >= m – Pattern still valid • Now add a letter

WEEDER: algorithm sketch Current pattern P (e, B) A C G T T A • Simple extension: no branches. – No change to B – e may increase by 1 or no change – Drop node if e > α • Branches: replace a node with its child nodes – Drop if e > α – B may change • Re-do Bit OR using all B’s • Try a different char if #occ < m • Report P when |P| = K

WEEDER: complexity • Can get all D(P, S) in time O(n. N (K choose α) 3α) ~ O(n. N Kα 3α). n: # sequences. Needed for Bit OR. • Better than O(KN 4 K) since usually α << K • Kα 3α may still be expensive for large K – E. g. K = 20, α = 6

WEEDER: More tricks Current pattern P A C G T T A • Eligible nodes: with at most α mismatches to P • Eligible nodes: with at most min( L, α) mismatches to P – L: current pattern length – : error ratio – Require that mismatches to be somewhat evenly distributed among positions • Prune tree at length K

MULTIPROFILER W* W*: ACGTACG W: ATGTAAG W W differs from W* at positions. The consensus sequence for the words in the -neighborhood of W is similar to W. If we ignore all the chars that are similar to W, the rest may suggest the difference between W and W*

MULTIPROFILER: alg sketch • For each word P in S W* – Find its α-neighborhood in S – List of words: C W • For each position j from 1. . K of the words in C – Find the most popular char that differ from P[j] • Replace α positions in P with the chars found above W*: ACGTACG W: ATGTAAG – Call the new word P’ • W* = argmin D(P’, S)

MULTIPROFILER W* W*: ACGTACG W: ATGTAAG W • No complexity provided in the paper • More efficient than WEEDER for longer patterns: N < Kα 3α • How to choose α is an issue: – Large α: too many noises in neighborhood – Small α: few true instances in neighborhood

Probabilistic modeling approaches for motif finding

Probabilistic modeling approaches • A motif model – Usually a PWM – M = (Pij), i = 1. . 4, j = 1. . k, k: motif length • A background model – Usually the distribution of base frequencies in the genome (or other selected subsets of sequences) – B = (bi), i = 1. . 4 • A word can be generated by M or B

Expectation-Maximization • For any word W, § P(W | M) = PW[1] 1 PW[2] 2…PW[K] K § P(W | B) = b. W[1] b. W[2] …b. W[K] • Let = P(M), i. e. , the probability for any word to be generated by M. • Then P(B) = 1 - • Can compute the posterior probability P(M|W) and P(B|W) § P(M|W) ~ P(W|M) * § P(B|W) ~ P(W|B) * (1 - )

Expectation-Maximization Initialize: Randomly assign each word to M or B • Let Zxy = 1 if position y in sequence x is a motif, and 0 otherwise • Estimate parameters M, , B Iterate until converge: • E-step: Zxy = P(M | X[y. . y+k-1]) for all x and y • M-step: re-estimate M, given Z (B usually fixed)

Expectation-Maximization position 5 1 Initialize E-step probability 1 5 9 9 M-step • E-step: Zxy = P(M | X[y. . y+k-1]) for all x and y • M-step: re-estimate M, given Z

MEME • • • Multiple EM for Motif Elicitation Bailey and Elkan, UCSD http: //meme. sdsc. edu/ Multiple starting points Multiple modes: ZOOPS, TCM

Gibbs Sampling • Another very useful technique for estimating missing parameters • EM is deterministic – Often trapped by local optima • Gibbs sampling: stochastic behavior to avoid local optima

Gibbs sampling Initialize: Randomly assign each word to M or B • Let Zxy = 1 if position y in sequence x is a motif, and 0 otherwise • Estimate parameters M, B, Iterate: • • • Randomly remove a sequence X* from S Recalculate model parameters using S X* Compute Zx*y for X* Sample a y* from Zx*y. Let Zx*y = 1 for y = y* and 0 otherwise

Gibbs Sampling probability position Sampling • Gibbs sampling: sample one position according to probability • • – Update prediction of one training sequence at a time Viterbi: always take the highest Simultaneously update predictions of all sequences EM: take weighted average

Gibbs sampling motif finders • Gibbs Sampler, based on C. Larence et. al. Science, 1993 • Align. ACE, Nat Biotech 1998, developed in Church lab, Harvard Univ • Bio. Prospector, X. Liu et. al. PSB 2001 , an improvement of Align. ACE

Better background model • Repeat DNA can be confused as motif – Especially low-complexity CACACA… AAAAA, etc. • Solution: more elaborate background model – Higher-order Markov model 0 th order: B = { p. A, p. C, p. G, p. T } 1 st order: B = { P(A|A), P(A|C), …, P(T|T) } … Kth order: B = { P(X | b 1…b. K); X, bi {A, C, G, T} } Has been applied to EM and Gibbs (up to 3 rd order)

Limits of Motif Finders 0 ? ? ? gene • Given upstream regions of coregulated genes: – Increasing length makes motif finding harder – random motifs clutter the true ones – Decreasing length makes motif finding harder – true motif missing in some sequences

Challenging problem d mutations n = 20 k L = 1000 • (k, d)-motif challenge problem • Many algorithms fail at (15, 4)-motif for n = 20 and L = 1000 • Combinatorial algorithms usually work better on challenge problem – However, they are usually designed to find (k, d)-motifs – Performance in real data varies

Motif finding in practice • Now we’ve found some good looking motifs – Easiest step? • What to do next? – Are they real? – How do we find more instances in the rest of the genome? – What are their functional meaning? • Motifs => regulatory networks

To make sense about the motifs • Each program usually reports a number of motifs (tens to hundreds) – Many motifs are variations of each other – Each program also report some different ones • Each program has its own way of scoring motifs – – Best scored motifs often not interesting AAAA ACAC TATAT

Strategies to improve results • Combine results from different algorithms usually helpful – Ones that appeared multiple times are probably more interesting • Except simple repeats like AAAAA or ATATA • Will talk about this later. – Cluster motifs into groups. Issues: • Measure similarities between two motifs (PWMs) • # of clusters

Strategies to improve results • Compare with known motifs in database – TRANSFAC – JASPAR • Issues: – Compute similarities among motifs – How similar is similar?

Strategies to improve results • Statistical test of significance – Enrichment in target sequences vs background sequences Target set T Assumed to contain a common motif, P Background set B Assumed to not contain P, or with very low frequency Ideal case: every sequence in T has P, no sequence in B has P

Statistical test for significance P Target set T Background set + target set B+T N Appeared in n sequences Appeared in m sequences • If n / N >> m / M – P is enriched (over-represented) in T – Statistical significance? • If we randomly draw N sequences from (B+T), how likely we will see at least n sequences with P? M

Hypergeometric distribution • A box with M balls, of which N are red, and the rest are blue. • We randomly draw m balls from the box • What’s the probability we’ll see n red balls? • Red ball: target sequences • Blue ball: background sequences • Total # of choices: (M choose m) • # of choices to have n red balls: (N choose n) x (M-N choose m-n)

Cumulative hypergeometric test for motif significance • We are interested in: if we randomly pick m balls, how likely that we’ll see at least n red balls? This can be interpreted as the p-value for the null hypothesis that we are randomly picking. Alternative hypothesis: our selection favors red balls. Equivalent: the target set T is enriched with motif P. Or: P is over-represented in T.

Examples • • Yeast genome has 6000 genes Select 50 genes believed to be co-regulated by a common TF Found a motif for these 50 genes It appeared in 20 out of these 50 genes In the whole genome, 100 genes have this motif M = 6000, N = 50, m = 100+20 = 120, n = 20 Intuition: – m/M = 120/6000. In Genome, 1 out 50 genes have the motif – N = 50, would expect only 1 gene in the target set to have the motif – 20 -fold enrichment • P-value = 6 x 10 -22 • n = 5. 5 -fold enrichment. P-value = 0. 003 • Normally a very low p-value is needed, e. g. 10 -10

ROC curve for motif significance • Motif is usually a PWM • Any word will have a score – – Typical scoring function: Log P(W | M) / P(W | B) W: a word. M: a PWM. B: background model • To determine whether a sequence contains a motif, a cutoff has to be decided – With different cutoffs, you get different number of genes with the motif – Hyper-geometric test first assumes a cutoff – It may be better to look at a range of cutoffs

ROC curve for motif significance P Target set T N Given a score cutoff Appeared in n sequences • • • Background set + target set B+T Appeared in m sequences With different score cutoff, will have different m and n Assume you want to use P to classify T and B Sensitivity: n / N Specificity: (M-N-m+n) / (M-N) False Positive Rate = 1 – specificity: (m – n) / (M-N) With decreasing cutoff, sensitivity , FPR M

ROC curve for motif significance A good cutoff Lowest cutoff. Every sequence has the motif. Sensitivity = 1. specificity = 0. sensitivity 1 0 ROC-AUC: area under curve. 1: the best. 0. 5: random. Motif 1 Motif 2 Random 0 1 -specificity Motif 1 is more enriched than motif 2. 1 Highest cutoff. No motif can pass the cutoff. Sensitivity = 0. specificity = 1.

Other strategies • Cross-validation – Randomly divide sequences into 10 sets, hold 1 set for test. – Do motif finding on 9 sets. Does the motif also appear in the testing set? • Phylogenetic conservation information – Does a motif also appears in the homologous genes of another species? – Strongest evidence – However, will not be able to find species-specific ones

Other strategies • Finding motif modules – Will two motifs always appear in the same gene? • Location preference – Some motifs appear to be in certain location • E. g. , within 50 -150 bp upstream to transcription start – If a detected motif has strong positional bias, may be a sign of its function • Evidence from other types of data sources – Do the genes having the motif always have similar activities (gene expression levels) across different conditions? – Interact with the same set of proteins? – Similar functions? – etc.