Modeling the Evolution of Regulatory Elements by Simultaneous Detection and Alignment with Phylogenetic Pair HMMs

The computational detection of regulatory elements in DNA is a difficult but important problem impacting our progress in understanding the complex nature of eukaryotic gene regulation. Attempts to utilize cross-species conservation for this task have been hampered both by evolutionary changes of functional sites and poor performance of general-purpose alignment programs when applied to non-coding sequence. We describe a new and flexible framework for modeling binding site evolution in multiple related genomes, based on phylogenetic pair hidden Markov models which explicitly model the gain and loss of binding sites along a phylogeny. We demonstrate the value of this framework for both the alignment of regulatory regions and the inference of precise binding-site locations within those regions. As the underlying formalism is a stochastic, generative model, it can also be used to simulate the evolution of regulatory elements. Our implementation is scalable in terms of numbers of species and sequence lengths and can produce alignments and binding-site predictions with accuracy rivaling or exceeding current systems that specialize in only alignment or only binding-site prediction. We demonstrate the validity and power of various model components on extensive simulations of realistic sequence data and apply a specific model to study Drosophila enhancers in as many as ten related genomes and in the presence of gain and loss of binding sites. Different models and modeling assumptions can be easily specified, thus providing an invaluable tool for the exploration of biological hypotheses that can drive improvements in our understanding of the mechanisms and evolution of gene regulation.Author Summary: The computational detection of regulatory elements in DNA is a difficult but important problem for decoding eukaryotic gene regulation. Increasing sequence data has made it possible to utilize related genomes, but this is not as straightforward as it may seem, as the evolution of noncoding regulatory regions is relatively poorly understood. In this work we describe a modeling framework and software implementation for aligning multiple DNA sequences to each other while simultaneously predicting functional regions in that DNA (such as the locations where proteins bind to the DNA for the purpose of regulating genes). Those functional regions may or may not be evolutionarily conserved across the sequences. Our framework allows for explicit modeling of evolutionary change across sequences in both the individual nucleotides making up the sequences and in the functional significance of the sequences (functional versus nonfunctional). While most competing frameworks and implementations are limited to a maximum number of sequences and their lengths, ours is scalable. We demonstrate the value of our system by using it to align a set of complex regulatory regions across ten Drosophila species and to predict protein-binding sites in those sequences.