QPS -- quadratic programming sampler, a motif finder using biophysical modeling

We present a Markov chain Monte Carlo algorithm for local alignments of nucleotide sequences aiming to infer putative transcription factor binding sites, referred to as the quadratic programming sampler. The new motif finder incorporates detailed biophysical modeling of the transcription factor binding site recognition which arises an intrinsic threshold discriminating putative binding sites from other/background sequences. We validate the principal functioning of the algorithm on a sample of four promoter regions from Escherichia coli. The resulting description of the motif can be readily evaluated on the whole genome to identify new putative binding sites.


I. INTRODUCTION
Transcription factors (TFs) are DNA binding proteins with regulatory effects. They may either independently or in an interplay with other proteins activate or repress the expression of genes related to the sequence they bind to, as illustrated in Figure 1. In this simplified picture, they act by either facilitating or impeding the recruitment of RNA polymerase holoenzymes, protein complexes responsible for the transcription of DNA to RNA. Information on exact locations and sequences of TF binding sites, typically 8-15 nucleotides, not only reveals which genes may or may not be controlled by a specific TF, thus permitting the construction of networks of genetic interaction, but is also indispensable when predicting novel putative binding sites with sequence motifs defined by the known examples. Exact binding sequences are often still unknown as typically in Escherichia coli, where roughly 70 of a total of 231 activating and repressing TFs have experimentally verified binding motifs [1,2]. In eukaryotes the picture is, as expected, worse: the commercial database TRANSFAC R contains in its most recent version 10018 entries for eukaryotic TFs, of which just 834 have reported binding motifs [3].
In principle, sequences associated to coregulated or even homologous genes can be used to infer putative TF binding sites by merely aligning those sequences. Several methods have been proposed to perform this task efficiently, see [4] for a summary of the most popular approaches. In practice, however, each of these methods has its flaws [5,6].
The approach we present here is grounded on the representation of TFs by free energy matrices as developed * Electronic address: afd@kth.se in QPMEME by Djordjevic et al. [7]. This representation yields the contributions of specific nucleotides to the total free energy of interaction between the TF and a sequence of DNA. Energies are scaled in terms of the chemical potential, rendering an intrinsic binding threshold which simplifies the task of distinguishing possible binding sites from the background on genomic sequences.

II. METHOD
We proceed by reviewing a simple model of TF-DNA binding and the representation of binding site motifs before discussing the probabilistic model on which our approach is based, as well as the details of the algorithm.

A. TF-DNA Binding
The interaction between a TF and a specific sequence of DNA can be written as a pair of ordinary differen-arXiv:0802.0258v1 [q-bio.QM] 2 Feb 2008 tial equations, describing the variation of bound and free concentrations of both reactants. Such a model depends on the reaction rates for binding and dissociation of TFs with DNA, symbolically stating with equilibrium constants for binding and dissociation K bind and K diss , respectively. In a system with many particles, the equilibrium concentrations of free TFs, specific DNA sequences and bound TF • DNA complexes can be related by the Arrhenius equation where β is an inverse temperature and E(S) stands for the free energy of binding the TF to DNA with a specific sequence S. Accepting the probability for the TF binding the sequence S to be be given by it follows with the chemical potential µ, relating abundance of the TF and affinity to its binding sites. The chemical potential is given by up to the additive constant C.

B. Binding Site Motifs
To define a binding motif from a collection of known binding sites S 1 , . . . , S N , of length L each, the construction of a matrix w ηi containing statistical weights for the occurrence of a nucleotide η at position i in the motif is usually adopted [8].
are constructed by counting the occurrences c ηi of nucleotide η at position i in each of the binding sites and comparing the thus defined frequencies f ηi to the probabilities p η with which to expect η in the sequence. Those probabilities can be deduced from the whole genome in question, from shorter regions containing the binding sites or even just from the latter. f ηi usually takes into account the error due to the finite amount of sequences in the collection by adding pseudocounts [8,9] to the occurrence counts. The weight matrix construction can then be used to evaluate the information content of a motif denoted by w a measure for the dissimilarity between the motif and random sequences stemming from the probabilistic model defined by p η . Further, w ηi can be applied to find putative binding sites in a genome. Each subsequence of length L is thus associated to an information score [10], describing the likelihood of that sequence to belong to the set of binding sites. How to chose a threshold score discriminating putative binders from non-binders, however, remains an open question in the weight matrix approach [threshold]. A more subtle description of binding motifs by free energy matrices [7] addresses the problem of finding a threshold by inverting the interpretation of a binding motif. Instead of describing similarities in a set of sequences, one attempts to model the requirements of a sequence to be able to bind a specific TF, now itself represented by the motif. The construction of such energy matrices ε ηi is based upon the assumption that binding motifs represented by ε should maximize the probability of recovering the set M = {M 1 , M 2 , . . . , M N } of known binding sites from an ensemble of random sequences, while the probability of identifying binding sites on unrelated random sequences is minimized. This can be performed maximizing the likelihood with probabilities Pr(M ) of generating a binding site sequence M and binding probabilities P ε (M ) for a TF to bind to this sequence. Maximizing L can be shown to be equivalent to minimizing the variance σ 2 ε of free energies resulting from the TF (ε) binding to random sequences [7], leading to an optimal ε by solving arg min where R 0 is a threshold free energy defining binding sites. In the inference method described by Djordjevic et al. [7], ε is efficiently approximated in the low temperature limit β → ∞ by quadratic programming. The elements of the energy matrix are shifted by the mean free energy of the TF being bound to random sequences E ε and rescaled by the absolute value of its chemical potential µ. The evaluation of ε ηi on a specific nucleotide sequence M = (α 1 , α 2 , . . . , α w ) gives the at first sight somewhat cumbersome result where E ε (M ) is the free energy of a TF associated to a TF ε binding to M . Yet this representation has a major advantage as compared to weight matrices: the chemical potential discriminates between strong and weak binding sites and since ε ηi is directly inferred in terms of µ, the threshold is implicitly given as R 0 = −1. All sequences with R ε (M ) ≤ R 0 are thus presumably strong binding sites, while R ε (M ) > R 0 denotes weak and nonbinders. More explicitly, the probability to find a specific TF bound to motif sequence M is with the rescaled inverse temperaturê For further details we also refer to [7,11,12] and references therein. Note, however, thatβ remains a free parameter as long as estimates for the average energy and the chemical potential are missing. Varyingβ does obviously not change the qualitative result in (11) stating P ε > 0.5 if R ε > R 0 , but will lead to a sharper discrimination of binding sequences from non-binding ones.

C. Probabilistic Model
Let us introduce a first order Markov model for the genomic background with conditional probabilities for the generation of a sequence M = (α 1 , α 2 , . . . , α w ) written as Pr(η | ν) Pr(find η preceded by ν) where η and ν represent single nucleotides. The probability of finding M among random sequences is thus given by understanding the boundary condition Adopting a more compact notation, we introduce the passage matrices W i (β) of probabilities for a TF to be bound to a site featuring the pair of nucleotides νη at position i Products W i (β) of passage matrices apparently yield the probabilities of coming across TF-DNA hybrids of corresponding length. Consequently, the generating function for a motif sequence bound by ε can be written as trace of the matrix product from which common statistical quantities describing binding of the TF to DNA can be derived. Of special interest is clearly the variance of free energies of a TF binding to random sequences, expressed as second derivative of the generating function Evaluating this expression, Djordjevic et al. show how to solve for ε by minimizing the variance [7].

D. Monte Carlo Sampling
In the set of N sequences {S 1 , . . . , S N }, e.g. promoter regions of co-regulated genes or upstream regions of homologous genes, we want to identify locally conserved subsequences of length w, supposedly sharing common TF binding sites. The alignment of the motif sequences M i is represented by their positions a i on the respective sequence S i . Let us first introduce the alignment probability distributionP ε (x) for a binding motif at position x on sequence S = (α 1 , α 2 , . . . , α L ), which is constructed from the binding probabilities P ε by settinĝ FIG. 2: Schematics of the QPS algorithm taking N input sequences S1, . . . , SN with initialisation (i), motif extraction (m), matrix computation (c) and evaluation (e) steps. S2 is highlighted being retained for updating.
The inference of an optimal local alignment is accomplished by a standard Monte Carlo method following the procedure (m') draw new alignment position a i from the alignment probability distributionP ε (a) and iterate with (c), Figure 2 illustrates the procedure in which each iteration draws a new alignment on the skipped sequence.
Iteratively updating the alignments, we sample the distribution of alignment positions until reaching stationarity on all of the sequences S i . The evolution of the distribution on biological sequences when inferring a binding motif of length 15 is shown in figure 3.
The sampling is performed at finite temperature (β > 0) in equation (11). Varying the rescaled inverse temperature during sampling allows to define a simple annealing schedule with stronger discrimination of possible binding from non-binders in the actual model asβ grows.

III. RESULTS AND DISCUSSION
We validated the functionality of QPS on a small set of coregulated promotors in in Escherichia coli consisting of aceBAKp, icdAp, pckAp, and ptsHp. Each region contains an experimentally known binding site for the fructose repressor protein FruR, which we tried to infer. The heatmaps in figure 3   the alignment probability distributionP ε when iterating with QPS and one of the obtained alignments is presented in table I. We settled for initiating the algorithm witĥ β = 20 and allowing for subsequent linear augmentation as the sampling proceeds,β with sampling iteration k. This corresponds approximately to an annealing schedule ∼ T −1 for the temperature. Exact binding sites remaining unknown for a vast number of TFs, and it is an interesting problem to try to infer a binding motif by aligning a set of sequences which are supposed to share sites for a specific TF. A wide range of different approaches have been developed [4,13], greedy pattern search algorithms [14], context free grammar constructors [15], and several statistical methods [16,17], to cite but just a small selection. Still it appears that no single method is capable of identifying motifs in a reliable way [5,18] and more recent approaches tend to combine several algorithms [4,19] to get a certain degree of cross-validation between individual methods. The method we present has been conceptually verified on a small sample of Escherichia coli promoter regions and might prove useful in combination with other approaches. The advantage of our algorithm is that it makes direct use of a biophysical representation of the TF. This representation is provided as result and can be readily applied to predict yet unknown binding sites elsewhere on the genome.
The here described algorithm has been implemented in C++ and is publicly available under the GPL on http://www.csc.kth.se/ ∼ afd/qps/.