Genome evolution: a sequence-centric approach

1. Genome evolution: a sequence-centric approach Lecture 11: Transcription factor binding sites

2. (Probability, Calculus/Matrix theory, some graph theory, some statistics) Tree of life Genome Size Elements of genome structure Elements of genomic information Simple Tree Models HMMs and variants PhyloHMM,DBN Context-aware MM Factor Graphs Probabilistic models Genome structure Inference Mutations DP Sampling Variational apx. LBP Models for populations Drift Selection and fixation Draft Parameter estimation Population EM Generalized EM (optimize free energy) Protein coding genes Inferring Selection Today refs: Papers cited

3. Molecular clocks and lineage acceleration • How universal is the rate of the evolutionary process? • Mutations may depend on the number of cell division and thus in the length of generation • Mutations depends on the genomic machinery to prevent them (Lynch?) • Mutations may also depend on the environment • The molecular clock (MC) hypothesis state that evolution is working in a similar rate for all lineages Relative rate test: KOA – KOB = 0 ? Test: KCA – KCB O A B C

4. Different molecular clocks in apes and primates Kim et al., 2006 PLoS genet

5. Sequence specific transcription factors • Sequence specific transcription factors (TFs) are a critical part of any gene activation or gene repression machinary • TFs include a DNA binding domain that recognize specifically “regulatory elements” in the genome. • The TF-DNA duplex is then used to target larger transcriptional structure to the genomic locus.

6. Sequence specificity is represented using consensus sequences or weight matrices • The specificity of the TF binding is central to the understanding of the regulatory relations it can form. • We are therefore interested in defining the DNA motifs that can be recognize by each TF. • A simple representation of the binding motif is the consensus site, usually derived by studying a set of confirmed TF targets and identifying a (partial) consensus. Degeneracy can be introduced into the consensus by using N letters (matching any nucleotide) or IUPAC characters (erpresenting pairs of nucleotides, for exampe W=[A|T], S=[C|G] • A more flexible representation is using weight matrices (PWM/PSSM): • PWMs are frequently plotted using motif logos, in which the height of the character correspond to its probability, scaled by the position entropy ACGCGT ACGCGA ACGCAT TCGCGA TAGCGT

7. TF binding energy is approximated by weight matrices We can interpret weight matrices as energy functions: This linear approximation is reasonable for most TFs. Leu3 data (Liu and Clarke, JMB 2002)

8. TF binding affinity is kinetically important, with possible functional implications • s Ume6 11.5 Average PWM energy Stronger prediction 5.5 ChIP ranges Stronger binding Tanay. Genome Res 2006 Kalir et al. Science 2001

9. TFs are present at only a fraction of their optimal sequence tragets. Binding is combinatorially regulated by co-factors, nucleosomes and histone modifications Re TSS Re ATG Lee et al. Nat Gen 2007

10. TFs are present at only a fraction of their optimal sequence tragets. Binding is combinatorially regulated by co-factors, nucleosomes and histone modifications Active Inactive Barski et al. Cell 2007

11. TFBSs are clustered in promoters or in “sequence modules” • The distribution of binding sites in the genome is non uniform • In small genomes, most sites are in promoters, and there is a bias toward nucleosome free region near the TSS • In larger genomes (fly) we observe CRM (cis-regulatory-modules) which are frequently away from the TSS. These represent enhancers. • A single binding site, without the context of other co-sites, is unlikely to represent a functional loci

12. Constructing a weight matrix from aligned TFBSs is trivial • This is done by counting (or “voting”) • Several databases (e.g., TRANSFAC, JASPAR) contain matrices that were constructed from a set of curated and validated binding site • Validated site: usually using “promoter bashing” – testing reported constructs with and without the putative site Transfac 7.0/11.3 have 400/830 different PWMs, based on more than 11,000 papers However, there are no real different 830 matrices outthere – the real binding repertoire in nature is still somewhat unclear

13. Probabilistic interpretation of weight matrices and a generative model • One can think of a weight matrix as a probabilistic model for binding sites: • This is the site independent model, defining a probability space over k-mers • Given a set of aligned k-mers, we know that the ML motif model is derived by voting (a set of independent multinomial variables – like the dice case) • Now assume we are given a set of sequences that are supposed to include binding sites (one for each), but that we don’t know where the binding sites are. • In other words the position of the binding site is a hidden variable h. • We introduce a background model Pb that describes the sequence outside of the binding site (usually a d-order Markov model) • Given complete data we can write down the likelihood of a sequence s as:

14. Using EM to discover PWMs de-novo • Inference of the binding site location posterior: • Note that only k-factors should be computed for each location (Pb(s) is constant)) • Inference of the binding site location posterior: • Note that only k factors should be computed for each location (Pb(s) is constant)) • Starting with an initial motif model, we can apply a standard EM: E: M: • As always with the EM, initializing to reasonable PWM would be critical Following Baily and Elkan, MEME 1995

15. Allowing false positive sequences • If we assume some of the sequences may lack a binding site, this should be incorporated into the model: • This is sometime called the ZOOPS model (Zero or one positions) l q hit s • In Bayesian terms: • Probability of sequence hit P(hit | S) • Probability of hit at position l = Pr(l|S) • We can consider the PWM parameters as variables in the model • Learning the parameters is then equivalent to inference

16. Using Gibbs sampling to discover PWMs de-novo • We can use Bibbs sampling to sample the hidden sites and estimate the PWM l q hit s • This is done by estimating the PWM from all locations except for the one we sample, and computing the hit probabilities as shown before l hit • Note that we are working with the MAP (Maximum a-posteriori) q to do the sampling: s l hit s • But this can be shown to approximate: Gibbs: Lawrence et al. Science 1993

17. Generalizing PWMs to allow site dependencies: mixture of PWMs and Trees Tree motif We only change the motif component of the likelihood model Mixture of PWMs Learning the model can become more difficult This is because computing the ML model parameter from complete data may be challenging Barash et al., RECOMB 2003

18. Discriminative scores for motifs • So far we used a generative probabilistic model to learn PWMs • The model was designed to generate the data from parameters • We assumed that TFBSs are distributed differently than some fixed background model • If our background model is wrong, we will get the wrong motifs.. • A different scoring approach try to maximize the discriminative power of the motif model. • We will not go here into the details of discriminative vs. generative models, but we shall exemplify the discriminative approach for PWMs. High specificity discriminator High sensitivity discriminator Lousy discriminator

19. Hypergeometric scores and thresholding PWMs Hyper geometric probability (sum for j>=k is the hg p-value) Positive Number of sequences True positive PWM score threshold For a discriminative score, we need to decide on both the PWM model and the threshold.

20. Exhaustive k-mer search • A very common strategy for motif finding is to do exhustive k-mer search. • Given a set of hits and a set of non hits, we will compute the number of occurrences of each k-mer in the two sets and report all cases that have a discriminative score higher than some threshold • Since k-mers either match or do not match, there is no issue with the threshold • For DNA, we will typically scan k=5-8. • This can be done efficiently using a map/hash: • Iterate on short sequence windows (of the desired k length) • For each window, mark the appearance of the k-mer in a table • Avoid double counting using a second map • It is easy to generalize such exhaustive approaches to include gaps or other types of degeneracy.

21. Refining k-mers to PWMs using heuristic “EM” • K-mer scan is an excellent intial step for finding refined weight matrices. For example, we can use them to initialize an EM. • If we want to find a weight matrix, but want to stick to the discriminative setting, we can heuristically use and “EM-like” algorithm: • Start with a k-mer seed • Add uniform prior to generate a PWM • Compute the optimal PWM threshold (maximal hyper-geometric score) • Restimate the PWM by voting from all PWM true positives • Consider additional PWM positions • Bound the position entropies to avoid over-fitting • Repeat two last steps until fail to improve score • There are of course no guarantees for improving the scores, but empirically this approach works very well.

22. High density arrays quantify TF binding preferences and identify binding sites in high throughput • Using microarrays (high resolution tiling arrays) we can now map binding sites in a genome-wide fashion for any genome • The problem is shifting from identifying binding sites to understanding their function and determining how sequences define them Harbison et al., Nature 2004

23. If only biology was that simple… Discrete and deterministic “binding sites” in yeast as identified by Young, Fraenkel and colleuges In fact, binding is rarely deterministic and discrete, and simple wiring is something you should treat with extreme caution.

24. r = 0.11 r = 0.74 r = 0.28 r = 0.8 r = 0.21 r = 0.72 ABF1 GCN4 MBP1 PWM sequence energy PWM sequence energy PWM sequence energy r = 0.42 r = 0.26 r = 0.42 r = 0.28 r = 0.42 r = 0.20 ChIP log(binding ratio) ChIP log(binding ratio) ChIP log(binding ratio) Correlation between PWM predicted binding and ChIP experiments spans high, medium and low affinity sites PWM regression exploits variable levels of binding affinity to robustly recover binding preferences. Motif regression optimizes the PWM given the overall correlation of the predicted binding energies and the measured ChIP values vs Tanay, GR 2004

25. TF1 Altered affinity CACGCGTT CACACGTT Rate? Selection? TFBS evolution: purifying selection and conservation TF1 TF1 Similar function Disrupted function CACGCGTT CACGCGTA CACGCGTT CACGAGTT Neutral evolution Low rate purifying selection TF1 TF2 Altered function CACGCGTT CACACGTT Low rate purifying selection

26. Binding sites conservation Kellis et al., 2003

27. Binding sites conservation: heuristic motif identification Kellis et al., 2003

28. Analyzing k-mer evolutionary dynamics • Instead of trying to identify conserved motifs try to infer the evolutionary rate of substitution between pairs of k-mers • Start from a multiple alignment and reconstruct ancestral sequences (assuming site independence, or even max parsimony) • Now estimate the number of substitution between pairs of 8-mers, compare this number to the number expected by the background model • Do it for a lot of sequence, so that statistics on the difference between observed and expected substitutions can be derived

29. TFBS evolution (so far): • Background: what are TFBSs doing? • Finding TFBS motifs • MEME/Gibbs • discriminative algorithms • affinity regression • Experiments • Comparative genomics of TFBSs: • kmer conservation • Kmer substituions • Evolutionary models • Halpren Bruno • Energy conservation • Intra-site epistasis and how to handle it • Using evolutionary models to define TFBSs

30. Inter-island organization in the Reb1 cluster: selection hints toward multi modality of Reb1 Saccharomyces TFBS Selection Network Nodes: octamers node conservation conserved @ 2SD conserved @ 3SD otherwise Arcs: 1nt substitution arc Rate Selection Normal neutral Low negative not enough stat Tanay et al., 2004

31. Substitution changing high affinity to high affinity motifs 0.3 0.2 0.1 0 TF1 -5 -4 -3 -2 -1 0 1 2 3 log delta affinity Altered affinity Motif 1 Motif 2 High Affinity (Kd < 60) Rate? Selection? Substitution changing high affinity to low affinity motifs High rate subs. Meidum Affinity (400 > Kd > 60) Leu3 selection network

32. TF5 TF4 AAATTT AATTTT AAAATT ACGCGT TCGCGT ACGCGT TF3 TF1 GATGAG GATGCG GATGAT CACGTG CACTTG All the rest TF2 TGACTG TGAGTG TGACTT A simple transcriptional code and its evolutionary implications

33. The Halpren-Bruno model for selection on affinity The basic notion here is of the relations between sequence, binding and function/fitness Sequence Binding energy Function We argued that E(S) can be approximated by a PWM F(E) is a completely different story, for example: Is there any function at all to low affinity binding sites? Is there a difference between very high affinity and plain strong binding sites? Are all appearances of the site subject to the same fitness landscape? In the Halpern-Bruno fraemwork, we assume that F(E)~E~PWM(S) In other words: everything is linear, there are no binding energy threshold for function. Also, all sites are equivalent.

34. The Halpren-Bruno model for selection on affinity We work on deriving the substitution rate at each position of the binding site, given its observed stationary frequency. We are assuming that the fitness of the site is defined by multiplying the fitnesses of each locus. According to Kimura’s theory, an allele with fitness s and a homogeneous population would fixate with probability: Assuming slow mutation rate (which allow us to assume a homogenous population) and motifs a and b with relative fitness s the fixation probabilities (chance of fixation given that mutation occurred!) are: If p represent the mutation probability, and p the stationary distribution, and if we assume the process as a whole is reversible then: (Halpern and Bruno, MBE 1998)

35. The Halpren-Bruno model for selection on affinity The HB model is extremely limited for the study of general sequences. When restricting the analysis to relatively specific sites, HB is not completely off Moses et al., 2003

36. Testing the general binding energy – fitness correspondence Expected and observed energy distribution in E.Coli CRP sites (left) and background (right) • While E(S) is approximated by a PWM, F(E) is unlikely to be linear • Assume that the background probability of a motif a is P0(a). In detailed balance, and assuming the fitness of a at functional sites is F(a), the stationary distribution at sites can be shown to be: • If we collapse all sites with binding energy E (and hence the same F(a)=F(E(a)) • The entire genome should behave like a mixture of background sequance and functional loci: • So we can try and recover Q(E) and therefore F(E) from the maximum likelihood parameters fitting an empirical W(E) Comparison of CRP energies in E.coli and S. typhimurium Inferred F(E), is shown in Orange Mustonen and Lassig, PNAS 2005

37. KS statistics More tests for possible conservation of low binding energy sites Simulation (Neutral, context aware) S. mikitae S. cerevisiae High affinity ΔE ΔE .. .. ΔE ΔE .. .. Low affinity

38. S S S Binding site conservation Conservation of total energy More tests for possible conservation of low binding energy sites Reb1 Conservation score binding energy percentile Gcn4 Ume6 Cbf1 Mbp1 Conservation score binding energy percentile binding energy percentile binding energy percentile binding energy percentile Tanay, GR 2006

39. Algorithms for discovering conserved sites: • Assuming PWM models • Search for loci that behave as predicted by the model (P(s|TFBS)/P(s|back) > Threshold) • Search for genomic regions with surprisingly many conserved binding sites • Search for an ML PWM • Search for motifs and conserved sites in aligned sequences • Assume the phylogeny, alignment, look for a PWM that will optimize the likelihood of the data • Search for motifs and conserved sites in unaligned sequences • Assume the phylogeny, look for an ML PWM • Search for a general ML evolutionary model (many PWMs) • Search for a set of PWMs/Motif model that will maximize the likelihood of the data

40. An evolutionary model P(si>si’|si-1) Phylo Tree Neutral model Motif (PWM?) TFBS evo model Substitution probability = background *(fixation) Substitution probability = HB model

41. The PhyME/PhyloGibbs model (Sinha, Blanchette, Tompa 2004, Sidharthan,Siggia,Nimwegen 2005, based on evo model by Siggia and collegues)our presentation is a bit generalized and adapted.. • Earlier/Other similar approaches using EM but practicxallt • Evolutionary model: • A phylogeny • Neutral independence model (simple tree, probabilities on branches) • A PWM-induced fixation probabilities that are proportional to the matrix weights: • Recall the single species generative model • The generalization to alignments is similar, but should consider a phylogeny at each locus Back Motif Back

42. Phylogenetic generative model TFBS 1 2 3 4 Prob of emitting the alignment given the mode Joint probability – product over independent trees Prior of having a motif/background Inference – for loci modes (l(i)) and for ancestral sequences (hierachicaly as in Ex 2) Learning – find ML PWM matrix

43. Posterior of the missing data Joint probability Complete data! EM for the PWM parameters Key point: the positions in the tree are independent once we decide on their “mode” (background or a certain positions in a binding site) Complete data: the mode of each position (l(i)). The ancestral sequences (h) Joint probability: The EM target: We can use the usual trick and decompose this into a sum of independent terms, one for each PWM position. For position k we have:

44. EM for the PWM parameters Average number of times of we observed the variables Log of the optimized parameter • In simple words: • we are optimizing a weighted sum of log(x*w) or log(1-Sxw) • The coefficient for the mutations paxj -> xj at PWM positions k is the average number of times we observe it given the previous parameter set w’ • The optimization problem is solved using Lagrange multipliers (to satisfy the constraint on the wk). But because each parameter appears in terms of two forms (log(x*w) and log(1-x*w), the solution is not as trivial as the dice case.

45. Example for intra-site epistasis AAACGT ACGCGT Assume the following TFBS alignment (where ACGCGT is the motif) According to the simple TFBS evo model, we should “pay” twice for the loss of the ACGCGT site, since each of the two mutations would be multiplied by a very low fixation probability ACGCGT “Double Loss” AAACGT ACGCGT The more realistic scenario involve one mutation that was under pressure, but then neutrality ACGCGT “Loss” Multiple possible trajectories can have different loss/gain/neutrality dynamics AAGCGT ACGCGT “Neutrality” AAACGT Assuming Affinity=PWM=Fitness gives the Halpren Burno model, or the frequency selection approximation When fitness(PWM) is non linear, we have epistasis which means problems for the simple loci-independent model Affinity Affinity PWM PWM Fitness Fitness

46. *selection factor An evolutionary model (2) P(si>si’|si-1) Phylo Tree Neutral model TF Target Sets Selection factors (per TF) To fully capture the effect of binding sites, we need to study a Markov model over the entire sequence. Instantaneous rates: background * selection on TFBS changes Mutation rate = background *selection factor

47. Controlling context effects: approximately independent blocks • The sequence is decomposed into intervals, each containing one or more overlapping binding sites, or sequences that are one mutation away from a binding site. Such intervals are called epistatic intervals • The joint probability is written as a product over epistatic intervals • For each epistatic interval we should compute exp(Qt)(from, to), where the rate matrix is large (4d) when d is the interval size • This is approximated by (exp(Qt/N))N

48. Learning model parameters • Finding the optimal selection factor is solved by non-linear optimization of the likelihood • Looking for target set is done using a greedy algorithm. • The resulted target sets and their selection factors are analogous to motifs/PWMs with some additional evolutionary parameter indicating the strength of selection (or the sufficiency of the motif to determine a functional site) • As shown to the left, similar motifs are sometime separated by significant selection factors, suggesting functional partitioning.

49. RNA folds and the function of RNA moelcules • RNA molecular perform a wide variety of functions in the cell • They differ in length and class, from very short miRNA to much longer rRNA or other structural RNAs. • They are all affected strongly by base-pairing – which make their structural mostly planar (with many exceptions!!) and relatively easy to model Simple RNA folding energy: number of matching basepairs or sum over basepairing weights More complex energy (following Zucker): each feature have an empirically determined parameters stem stacking energy (adding a pair to a stem) bulge loop length interior loop length hairpin loop length dangling nucleotides and so on. Pseudoknots (breaking of the basepairing hierarchy) are typically forbidden:

50. Predicting fold structure Due to the hierarchical nature of the structure (assuming no pseudoknots), the situation can be analyzed efficiently using dynamic programming. We usually cannot be certain that there is a single, optimal fold, especially if we are not at all sure we are looking at a functional RNA. It would be better to have posterior probabilities for basepairing given the data and an energy model…This can be achieved using a generalization of HMM called Stochastic Context Free Grammar