CSMET: Comparative Genomic Motif Detection via Multi-Resolution Phylogenetic Shadowing
We concern ourselves with uncovering motifs in eukaryotic cis-regulatory modules ( CRM ) from multiple evolutionarily related species , such as the members from the Drosophila clade ., Due to high degeneracy of motif instances , and complex motif organization within the CRMs , pattern-matching-based motif search in higher eukaryotes remains a difficult problem , even when representations such as the position weight matrices ( PWMs ) of the motifs are given ., Extant methods that operate on a single genome or simpler organisms such as yeast often yield a large number of false positives , especially when the sequence to be examined spans a long region ( e . g . , tens of thousands of bps ) beyond the basal promoters , where possible CRMs could be located ., As in gene finding , having orthologous sequences from multiple evolutionarily related taxa can potentially benefit motif detection because a reasonable alignment of these sequences could enhance the contrast of sequence conservation in motifs with respect to that of the non-motif regions , However , the alignment quality of non-coding regions is usually significantly worse than that of the coding regions , so that the aligned motif sequences are not reliably orthologous ., This is often unavoidable even for the best possible local alignment software because of the short lengths and weak conservation of TFBSs ., When applying a standard shadowing model on such alignments , motif instances aligned with non-orthologous sequences or gaps can be hard to identify due to low overall shadowing score of the aligned sequences ( Figure 1A ) ., In addition to the incomplete orthology due to imperfect alignment , a more serious concern comes from a legitimate uncertainty over the actual functional orthology of regions that are alignment-wise orthologous ., A number of recent investigations have shown that TFBS loss and gain are fairly common events during genome evolution 8 , 12 ., For example , Patel et al 13 showed that aligned “motif sites” in orthologous CRMs in the Drosophila clade may have varying functionality in different taxa ., Such cases usually occur in regions with reduced evolutionary constraints , such as regions where motifs are abundant , or near a duplication event ., The sequence dissimilarities of CRMs across taxa include indel events in the spacers , as well as gains and losses of binding sites for TFs such as the bcd-3 and hb-1 motifs in the evenskipped stripe 2 ( eve2 ) ( Figure 1B ) ., A recent statistical analysis of the Zeste binding sites in several Drosophila taxa also revealed existence of large-scale functional turnover 12 ., Nevertheless , the fact that sequence similarity is absent does not necessarily mean that the overall functional effect of the CRM as a whole is vastly different ., In fact , for the Drosophila clade , despite the substantial sequence dissimilarity in gap-gene CRMs such as eve2 , the expression of these gap genes shows similar spatio-temporal stripe patterns across the taxa 8 , 13 ., Although a clear understanding of the evolutionary dynamics underlying such inter- and intra-taxa diversity is still lacking , it is hypothesized that regulatory sequences such as TFBSs and CRMs may undergo adaptive evolution via stabilizing selections acting synergistically on different loci within the sequence elements 8 , 12 , which causes site evolution to be non-iid and non-isotropic across all taxa ., In such a scenario , it is crucial to be able to model the evolution of biological entities not only at the resolution of individual nucleotides , but also at more macroscopic levels , such as the functionality of whole sequence elements such as TFBSs over lineages ., To our knowledge , so far there have been few attempts along this line , especially in the context of motif detection ., The CSMET model presented in this paper intends to address this issue ., Orthology-based motif detection methods developed so far are mainly based on nucleotide-level conservation ., Some of the methods do not resort to a formal evolutionary model 14 , but are guided by either empirical conservation measures 15–17 , such as parsimonious substitution events or window-based nucleotide identity , or by empirical likelihood functions not explicitly modeling sequence evolution 4 , 18 , 19 ., The advantage of these non-phylogeny based methods lies in the simplicity of their design , and their non-reliance on strong evolutionary assumptions ., However , since they do not correspond to explicit evolutionary models , their utility is restricted to purely pattern search , and not for analytical tasks such as ancestral inference or evolutionary parameter estimation ., Some of these methods employ specialized heuristic search algorithms that are difficult to scale up to multiple species , or generalize to aligned sequences with high divergence ., Phylogenetic methods such as EMnEM 20 , MONKEY 21 , and our in-house implementation of PhyloHMM ( originally implemented in 1 for gene finding , but in our own version tailored for motif search ) explicitly adopt a complete and independent shadowing model at the nucleotide level ., These methods are all based on the assumption of homogeneity of functionality across orthologous nucleotides , which is not always true even among relatively closely related species ( e . g . , of divergence less than 50 mya in Drosophila ) ., Empirical estimation and simulation of turnover events is an emerging subject in the literature 12 , 22 , but to our knowledge , no explicit evolutionary model for functional turnover has been proposed and brought to bear in comparative genomic search of non-conserved motifs ., Thus our CSMET model represents an initial foray in this direction ., Closely related to our work , two recent algorithms , rMonkey 12—an extension over the MONKEY program , and PhyloGibbs 9—a Gibbs sampling based motif detection algorithm , can also explicitly account for differential functionality among orthologs , both using the technique of shuffling or reducing the input alignment to create well conserved local subalignments ., But in both methods , no explicit functional turnover model has been used to infer the turnover events ., Another recent program , PhyME 10 , partially addresses the incomplete orthology issue via a heuristic that allows motifs only present in a pre-chosen reference taxon to be also detectable , but it is not clear how to generalize this ability to motifs present in arbitrary combination of other taxa , and so far no well-founded evolutionary hypothesis and model is provided to explain the heuristic ., Non-homogeneous conservation due to selection across aligned sites has also been studied in DLESS 23 and PhastCons 24 , but unlike in CSMET , no explicit substitution model for lineage-specific functional evolution was used in these algorithms , and the HMM-based model employed there makes it computationally much more expensive than CSMET to systematically explore all possible evolutionary hypotheses ., A notable work in the context of protein classification proposed a phylogenomic model over protein functions , which employs a regression-like functional to model the evolution of protein functions represented as feature vectors along lineages in a complete phylogeny 25 , but such ideas have not been explored so far for comparative genomic motif search ., Various nucleotide substitution models , including the Jukes-Cantor 69 ( JC69 ) model 26 , and the Felsenstein 81 ( F81 ) model 27 , have been employed in current phylogenetic shadowing or footprinting algorithms ., PhyloGibbs and PhyME use an analogue of F81 proposed in 28 , which is one of the simplest models to handle arbitrary stationary distributions , necessary to model various specific PWMs of motifs ., Both PhyME and PhyloGibbs also offer an alternative to use a simplified star-phylogeny to replace the phylogenetic tree when dealing with a large number of taxa , which corresponds to an even simpler substitution process ., Our CSMET model differs from these existing methods in several important ways ., First , it uses a different evolutionary model based on a coupled-set of both functional and nucleotide substitution processes , rather than a single nucleotide substitution model to score every alignment block ., Second , it uses a more sophisticated and popular nucleotide substitution process based on the Felsenstein84 ( F84 ) model 29 , which captures the transition/transversion bias ., Third , it employs a hidden Markov model that explicitly models autocorrelation of evolutionary rates on successive sites in the genome ., Fourth , it uses an efficient deterministic inference algorithm that is linear to the length of the input sequence and either exponential ( under a full functional phylogeny ) or linear ( under a star-shaped functional phylogeny ) to the number of the aligned taxa , rather than the Monte Carlo or heuristic search algorithms that require long convergence times ., Essentially , CSMET is a context-dependent probabilistic graphical model that allows a single column in a multiple alignment to be modeled by multiple evolutionary trees conditioned on the functional specifications of each row ( i . e . , the functional identity of a substring in the corresponding taxon ) ( Figure 2 ) ., When conjoined with a hidden Markov model that auto-correlates the choices of different evolutionary rates on the phylogenetic trees at different sites , we have a stochastic generative model of phylogenetically related CRM sequences that allows both binding site turnover in arbitrary subsets of taxa , and coupling of evolutionary forces at different sites based on the motif organizations within CRMs ., Overall , CSMET offers an elegant and efficient way to take into consideration complex evolutionary mechanisms of regulatory sequences during motif detection ., When such a model is properly trained on annotated sequences , it can be used for comparative genomic motif search in all aligned taxa based on a posterior probabilistic inference algorithm ., This model can be also used for de novo motif finding as programs such as PhyloGibbs and PhyME , with a straightforward extension of the inference procedure that couples the training and prediction routines in an expectation-maximization ( EM ) iteration on unannotated sequence alignments ., In this paper , we focus on supervised motif search in higher eukaryotic genomes ., We compare CSMET with representative competing algorithms , including EMnEm , PhyloHMM , PhyloGibbs , and a mono-genomic baseline Stubb ( which uses an HMM on single species ) on both simulated data , and a pre-aligned Drosophila dataset containing 14 developmental CRMs for 11 aligned Drosophila species ., Annotations for motif occurrences in D . melanogaster of 5 gap-gene TFs - Bicoid , Caudal , Hunchback , Kruppel and Knirps - were obtained from the literature ., We show that CSMET outperforms the other methods on both synthetic and real data , and identifies a number of previously unknown occurrences of motifs within and near the study CRMs ., The CSMET program , the data used in this analysis , and the predicted TFBS in Drosophila sequences , are available for download at http://www . sailing . cs . cmu . edu/csmet/ ., At present , biologically validated orthologous motifs and CRMs across multiple taxa are extremely rare in the literature ., In most cases , motifs and CRMs are only known in some well-studied reference taxa such as the Drosophila melanogaster; and their orthologs in other species are deduced from multiple alignments of the corresponding regulatory sequences from these species according to the positions and PWMs of the “reference motifs” in the reference taxon ., This is a process that demands substantial manual curation and biological expertise; rarely are the outcomes from such analysis validated in vivo ( but see 8 for a few such validations in some selected Drosophila species where the transgenic platforms have been successfully developed ) ., At best , these real annotations would give us a limited number of true positives across taxa , but they are not suitable for a systematic performance evaluation based on precision and recall over true motif instances ., Thus we first compare CSMET with a carefully chosen collection of competing methods on simulated CRM sequences , where the motif profiles across all taxa are completely known ., We choose to compare CSMET with 3 representative algorithms for comparative genomic motif search , PhyloGibbs , EMnEM , PhyloHMM; and the program Stubb , which is specialized for motif search in eukaryotic CRMs , and in our paper , set to operate in mono-genomic mode ., The rationale for choosing these 4 benchmarks is detailed in the Material and Methods ., We applied CSMET and competing methods to a multi-specific dataset of Drosophila early developmental CRMs and motifs compiled from the literature 38 ., However , in this situation , we score accuracy only on the motifs annotated in Drosophila melanogaster ( rather than in all taxa ) , because they are the only available gold-standard ., Upon concluding this section , we also report some interesting findings by CSMET of putative motifs , some of which only exist in other taxa and do not have known counterparts in melanogaster ., CSMET is a novel phylogenetic shadowing method that can model biological sequence evolution at both nucleotide level at each individual site , and functional level of a whole TFBS ., It offers a principled way of addressing the problem that can seriously compromise the performance of many extant conservation-based motif finding algorithms: motif turnover in aligned CRM sequences from different species , an evolutionary event that results in functional heterogeneity across aligned sequence entities and shatters the basis of conventional alignment scoring methods based on a single function-specific phylogeny ., CSMET defines a new evolution-based score that explicitly models functional substitution along the phylogeny that causes motif turnover , and nucleotide divergence of aligned sites in each taxa under possibly different function-specific phylogenies conditioning on the turnover status of the site in each taxon ., In principle , CSMET can be used to estimate the rate of turnover of different motifs , which can elucidate the history and dynamics of functional diversification of regulatory binding sites ., But we notice that experimentally validated multi-species CRM/TFBS annotations that support an unbiased estimate of turnover rates are yet to be generated , as currently almost all biologically validated motifs only exist in a small number of representative species in each clade of the tree of life , such as melanogaster in the Drosophila clade ., Manual annotation on CRM alignments , as we used in this paper , tends to bias the model toward conserved motifs ., Thus , at this time , the biological interpretation of evolutionary parameters on the functional phylogeny remains preliminary ., Nevertheless , these estimated parameters do offer important utility from a statistical and algorithmic point of view , by elegantly controlling the trade-off between two competing molecular substitution processes—that of the motif sequence and of the background sequence—at every aligned site across all taxa beyond what is offered in any existing motif evolution model ., Empirically , we find that such modelling is useful in motif detection ., On both synthetic data and 14 CRMs from 11 Drosophila taxa , we find that the CSMET performs competitively against the state-of-the-art comparative genomic motif finding algorithm , PhyloGibbs , and significantly outperforms other methods such as EMnEM , PhyloHMM and Stubb ., In particular , CSMET demonstrates superior performance in certain important scenarios , such as cases where aligned sequences display significant divergence and motif functionalities are apparently not conserved across taxa or over multiple adjacent sites ., We also find that both CSMET and PhyloGibbs significantly outperform Stubb when the latter is naively applied to sequences of all taxa without exploiting their evolutionary relationships ., Our results suggest that a careful exploration of various levels of biological sequence evolution can significantly improve the performance of comparative genomic motif detection ., Recently , some alignment-free methods 19 have emerged which search for conserved TFBS rich regions across species based on a common scoring function , e . g . , distribution of word frequencies ( which in some ways mirrors the PWM of a reference species ) ., One may ask , given perhaps in the future a perfect search algorithm ( in terms of only computational efficiency ) , do we still need explicit model-based methods such as CSMET ?, We believe that even if exhaustive search of arbitrary string patterns becomes possible , models such as CSMET still offer important advantage not only in terms of interpretability and evolutionary insight as discussed above , but possibly also in terms of performance because of the more plausible scoring schemes they use ., This is because it is impractical to obtain the PWM of a motif in species other than a few reference taxa , thus the scores of putative motif instances in species where their own versions of the PWM are not available can be highly inaccurate under the PWM from the reference species due to evolution of the PWM itself in these study species with respect to the PWM in the reference species ., The CSMET places the reference PWM only at the tree root as an equilibrium distribution; for the tree leaves where all study species are placed , the nucleotide substitution model along tree branches allows sequences in each species to be appropriately scored under a species-specific distribution that is different from the reference PWM , thereby increasing its sensitivity to species-specific instantiations of motifs ., A possible future direction for this work lies in developing better approximate inference techniques for posterior inference under the CSMET model , especially under the scenarios of studying sequences from a large clade with many taxa , and/or searching for multiple motifs simultaneously ., It is noteworthy that our methods can be readily extended for de novo motif detection , for which an EM or a Monte Carlo algorithm can be applied for model-estimation based on the maximum likelihood principle ., Currently we are exploring such extensions ., Also we intend to develop a semi-supervised training algorithm that does not need manual annotation of motifs in other species on the training CRM alignment , so that we can obtain a less biased estimate of the evolutionary parameters of the CSMET model ., A problem with most of the extant motif finders , including the proposed CSMET , is that the length variation of aligned motifs ( e . g . , alignments with gaps ) cannot be accommodated ., In our model , while deletion events may be captured as gaps in the motif alignment , insertion events cannot be captured as the length of the motif is fixed ., This is because in a typical HMM sequence model the state transitions between sites within motifs are designed to be deterministic ., Thus stochastically accommodating gaps ( insertion events ) within motifs is not feasible ., Hence , some of the actual motifs missed by the competing algorithms were “gapped” motifs ., These issues deserve further investigation ., We use the Felsenstein 1984 model ( F84 ) 29 , which is similar to the Hasegawa–Kishino–Yanos 1985 model ( HKY85 ) 44 and widely used in the phylogenetic inference and footprinting literature 5 , 29 , for nucleotide substitution in our motif and background phylogeny ., Formally , F84 is a five-parameter model , based on a stationary distribution π ≡ πA , πT , πG , πC′ ( which constitutes three free parameters as the equilibrium frequencies sum to ) and the additional parameters κ and ι which impose the transition/transversion bias ., According to this model , the nucleotide-substitution probability from an internal node c to its descendent c′ along a tree branch of length b can be expressed as follows: ( 3 ) where i and j denote nucleotides , δij represents the Kronecker delta function , and εij is a function similar to the Kronecker delta function which is 1 if i and j are both pyrimidines or both purines , but 0 otherwise ., The summation in the denominator concisely computes purine frequency or pyrimidine frequency ., A more intuitive parameterization for F84 involves the overall substitution rate per site μ and the transition/transversion ratio ρ , which can be easily estimated or specified ., We can compute the transition matrix PN from μ and ρ using Equation 3 based on the following relationship between ( κ , ι ) and ( μ , ρ ) :To model functional turnover of aligned substrings along functional phylogeny Tf , we additionally define a substitution process over two characters ( 0 and 1 ) corresponding to presence or absence of functionality ., Now we use the single parameter JC69 model 26 for functional turnover due to its simplicity and straightforward adaptability to an alphabet of size, 2 . The transition probability along a tree branch of length β ( which represents the product of substitution rate μ and evolution time t , which are not identifiable independently , ) is defined by: ( 4 ) We estimate the evolutionary parameters from training data based on maximum likelihood , details are available in the Text S1 ., A complete phylogenetic tree T ≡ {τ , π , β , λ} with internal nodes {Vi; i\u200a=\u200a1:K′} and leaf nodes {Vi; i\u200a=\u200aK′+1:K} , where K denotes the total number of nodes ( i . e . , current and ancestral species ) instantiated in the tree and the node indexing follows a breath-first traversal from the root , defines a joint probability distribution of all-node configurations ( i . e . , the nucleotide contents at an aligned site in all species instantiated in the tree ) , which can be written as the following product of nt-substitution probabilities along tree branches: ( 5 ) where Vpa ( i ) denotes the parent-node of the node i in the tree , and the substitution probability PN ( ) is defined by Equation, 3 . For each position l of the multiple alignment , computing the probability of the entire column denoted by Al of aligned nucleotides from species corresponding to the leaves of a phylogenetic tree T ( l ) defined on position l , i . e . , P ( Al|T ( l ) ) , where Al correspond to an instantiation of the leaf nodes {Vi; i\u200a=\u200aK′+1:K} , takes exponential time if performed naively , since it involves the marginalization of all the internal nodes in the tree , i . e . , ( 6 ) We use the Felsenstein pruning algorithm 30 , which is a dynamic programming method that computes the probability of a leaf-configuration under a tree from the bottom up ., At each node of the tree , we store the probability of the subtree rooted at that node , for each possible nucleotide at that node ., At the leaves , only the probability for the particular nucleotide instantiated in the corresponding taxon is non-zero , and for all the other nucleotides , it is zero ., Unlike the naive algorithm , the pruning algorithm requires an amount of time that is proportional to the number of leaves in the tree ., We use a simple extension of this algorithm to compute the probabilities of a partial-alignment defined earlier under a marginal phylogeny , which is required in the coupled-pruning algorithm for CSMET , by considering only the leaves instantiated in ( but not in ) that is under a subtree T′ ( l ) that forms the marginal phylogeny we are interested in ., Specifically , let correspond to possible instantiations of the subset of nodes we need to marginalized out ., Since we already how to compute P ( Al|T ( l ) ) via marginalization over internal nodes , we simply further this marginalization over leaf nodes that corresponds to taxa instantiated in , i . e . , ( 7 ) where denotes the leaves instantiated in ., This amounts to replacing the leaf-instantiation step , which was originally operated on all leaves in the Felsenstein pruning algorithm , by a node-summation step over those leaves in ., In fact , in can be easily shown that this is equivalent to performing the Felsenstein pruning only on the partial tree T′ ( l ) that directly shadows , which is a smaller tree than the original T ( l ) , and only requires time ., Under the CSMET model , to perform the forward-backward algorithm for either motif prediction or unsupervised model training , we need to compute the emission probability given each functional state at every alignment site ., This is nontrivial because a CSMET is defined on an alignment block containing whole motifs across taxa rather than on a single alignment-column ., We adopt a “block-approximation” scheme , where the emission probability of each state at a sequence position , say , t , is defined on an alignment block of length L started at t , i . e . , , where At≡ ( A1 ( t ) , A2 ( t ) , … , AL ( t ) ) , and Al ( t ) denotes the lth column in an alignment block started from position t ., The conditional likelihood At given the nucleotide-evolutionary trees T and Tb coupled by the annotation tree Ta under a particular HMM state st is also hard to calculate directly , because the leaves of the two nucleotide trees are connected by the leaves of the annotation tree ( Figure 2B ) ., However , if the leaf-states of the annotation tree are known , the probability components coming from the two trees become conditionally independent and factor out ( see Equation 2 ) ., Recall that for a motif of length L , the motif tree actually contains L site-specific trees , i . e . , , and the the choice of these trees for every site in the same row ( i . e . , taxon ) , say , in the alignment block At , is coupled by a common annotation state ., Hence , given an annotation vector Zt for all rows of At , we actually calculate the probability of two subset of the rows given two subtrees ( i . e . , marginal phylogenies ) of the original phylogenetic trees for motif and backgrounds , respectively ( Figure 2B ) ., The subset is constructed by simply stacking the DNA bases of those taxon for which the annotation variables indicate that they were generated from the motif tree ., The subtree is constructed by simply retaining the set of nodes which correspond to the chosen subset , and the ancestors thereof ., Similarly we have and ., Hence , we obtain ( 8 ) The probability of a particular leaf-configuration of a tree , be it a partial or complete nucleotide tree , or an annotation tree , can be computed efficiently using the pruning algorithm ., Thus for each configuration of zt , we can readily compute and ., The block emission probability under CSMET can be expressed as: ( 9 ) where we use , , and to make explicit the dependence of the partial blocks and marginal trees on functional indicator vector zt ., We call this algorithm a coupled-pruning algorithm ., Note that in this algorithm we need to sum over a total number of 2M configurations of zt where M is the total number of taxa ( i . e . , rows ) in matrix At ., It is possible to reduce the computational complexity using a full junction tree algorithm on CSMET , which will turn the graphical model underlying CSMET into a clique tree of width ( i . e . , maximum clique size ) possibly smaller than M . But this algorithm is complicated and breaks the modularity of the tree-likelihood calculation by the coupled-pruning algorithm ., In typical comparative genomic analysis , we expect that M will not be prohibitively large , so our algorithm may still be a convenient and easy-to-implement alternative to the junction-tree algorithm ., Also this computation can be done off-line and in parallel ., Given the emission probabilities for each ancestral functional state at each site , we use the forward-backward algorithm for posterior decoding of the sequence of ancestral functional states along the input CRM alignment of length N . The procedure is the same as in a standard HMM applied to a single sequence , except that now the emission probability at each site , say with index t , is defined by the CSMET probability over an alignment block At at that position under an ancestral functional state , rather than the conditional probability of a single nucleotide observed at position t as in the standard HMM ., The complexity of this FB-algorithm is O ( Nk2 ) where k denotes the total number of functional states ., In this paper , we only implemented a simple HMM with one type motif allowed on either strand , so that k\u200a=\u200a3 . We defer a more elaborate implementation that allows multiple motifs and encodes sophisticated CRM architecture as in LOGOS 33 to a future extension ., Given an estimate of , we can infer the MAP estimates of —the functional annotation of every site t in every taxon i of the alignment ., Specifically , the posterior probability of a column of functional states Zt under ancestral functional state can be expressed as: ( 10 ) Recall that in the coupled-pruning algorithm , we can readily compute all the three conditional probability terms in the above equation ., Performing posterior inference allows us to make motif predictions in two ways ., A simple way is look at blocks in the alignment at which the posterior inference produces ones , and predict those to be motifs ., Alternatively , we can also use the inferred state of the alignment block together with the inferred ancestral state to compute a probability score ( as a heuristic ) based on the functional annotation tree ., The score for the block is the sum of probabilities of each block element being one ., Given blocks of aligned substrings {At} containing motif instances in at least one of the aligned taxa , in principle we can estimate both the annotation tree Tf ≡ {α , τf , βf} and the motif trees Tm ≡ {θ , τm , βm , λm} based on a maximum likelihood principle ., But since in our case most training CRM sequences do not have enough motif data to warrant correct estimation of the motif and function tree , we use the topology and branch lengths of a tree estimated by fastDNAml 36 from the entire CRM sequence alignment ( containing both motif and background ) as the common basis to build the Tf and Tm ., Specifically , fastDNAml estimates a maximum likelihood tree under the F84 model from the entire CRM alignment; we then scale the branch lengths of this tree to get the sets of branch lengths for Tf and Tm by doing a simple linear search ( see below ) of the scaling coefficient that maximize the likelihood of aligned motif sequences and aligned annotation sequences , under the Tm and Tf ( scaled based on the coefficients ) respectively ., For simplicity , we estimate the background tree Tb ≡ {θ , τb , βb , λb} separately from only aligned background sequences that are completely orthologous ( i . e . , containing no motifs in any taxon ) ., For both motifs and background phylogenies , the Felsenstein rate parameter μ for the corresponding nucleotide substitution models must also be estimated from the training data ., More technically , note that for Tm the scaling coefficient β and the rate parameter μ form a product in the expression of the substitution probability ( see Equation 3 ) and are not identifiable independently ., Thus we only need to estimate the compound rate parameter μ′\u200a=\u200aμβ ., Ideally , the optimal value of the μ′ should be obtained by performing a gradient descent on the likelihood under the corresponding phylogeny with respect to μ′ ., However , due to the phylogenetic tree probability terms involved in the likelihood computation , there is no closed form expression for the gradient that can be evaluated for a specific value of the compound rate parameter to determine the direction to choose for optimization ., Therefore , to find an approximation to the optimal value of μ′ , we perform a simple linear search in the space of μ′ as follows: and are lower and upper bounds respectively on the space of μ′ that is searched , and are heuristically chosen based on observation ., The step δ can be chosen to be as small as desired or is allowabl
Functional turnover of transcription factor binding sites ( TFBSs ) , such as whole-motif loss or gain , are common events during genome evolution ., Conventional probabilistic phylogenetic shadowing methods model the evolution of genomes only at nucleotide level , and lack the ability to capture the evolutionary dynamics of functional turnover of aligned sequence entities ., As a result , comparative genomic search of non-conserved motifs across evolutionarily related taxa remains a difficult challenge , especially in higher eukaryotes , where the cis-regulatory regions containing motifs can be long and divergent; existing methods rely heavily on specialized pattern-driven heuristic search or sampling algorithms , which can be difficult to generalize and hard to interpret based on phylogenetic principles ., We propose a new method: Conditional Shadowing via Multi-resolution Evolutionary Trees , or CSMET , which uses a context-dependent probabilistic graphical model that allows aligned sites from different taxa in a multiple alignment to be modeled by either a background or an appropriate motif phylogeny conditioning on the functional specifications of each taxon ., The functional specifications themselves are the output of a phylogeny which models the evolution not of individual nucleotides , but of the overall functionality ( e . g . , functional retention or loss ) of the aligned sequence segments over lineages ., Combining this method with a hidden Markov model that autocorrelates evolutionary rates on successive sites in the genome , CSMET offers a principled way to take into consideration lineage-specific evolution of TFBSs during motif detection , and a readily computable analytical form of the posterior distribution of motifs under TFBS turnover ., On both simulated and real Drosophila cis-regulatory modules , CSMET outperforms other state-of-the-art comparative genomic motif finders .
Functional turnover of transcription factor binding sites ( TFBSs ) , such as whole-motif loss or gain , are common events during genome evolution , and play a major role in shaping the genome and regulatory circuitry of contemporary species ., Conventional methods for searching non-conserved motifs across evolutionarily related species have little or no probabilistic machinery to explicitly model this important evolutionary process; therefore , they offer little insight into the mechanism and dynamics of TFBS turnover and have limited power in finding motif patterns shaped by such processes ., In this paper , we propose a new method: Conditional Shadowing via Multi-resolution Evolutionary Trees , or CSMET , which uses a mathematically elegant and computationally efficient way to model biological sequence evolution at both nucleotide level at each individual site , and functional level of a whole TFBS ., CSMET offers the first principled way to take into consideration lineage-specific evolution of TFBSs and CRMs during motif detection , and offers a readily computable analytical form of the posterior distribution of motifs under TFBS turnover ., Its performance improves upon current state-of-the-art programs ., It represents an initial foray into the problem of statistical inference of functional evolution of TFBS , and offers a well-founded mathematical basis for the development of more realistic and informative models .
computational biology/evolutionary modeling, computational biology/comparative sequence analysis, computational biology/sequence motif analysis