MotifAdjuster: a tool for computational reassessment of transcription factor binding site annotations

MotifAdjuster helps to detect errors in binding site annotations.


Rationale
The regulation of gene expression involves a complex system of interacting components in all living organisms [1] and is of fundamental interest, for instance, for cell maintenance and development. One level of regulation is realized by DNAbinding transcription factors (TFs). The DNA-binding domain of a TF is capable of recognizing specific binding sites (BSs) in the promoter regions of its target genes [2]. Binding of a TF can induce (activator) or inhibit (repressor) the transcription of its target genes. The general ability to control a target gene may depend on the BS itself, its strand orientation, and its position with respect to the transcription start site. If other BSs are present, the ability of a TF to bind the DNA may additionally depend on strand orientations and positions of these BSs.
One important prerequisite for research on gene regulation is the reliable annotation of BSs. The approximate regions on the double-stranded DNA sequence bound by TFs can be determined by wet-lab experiments such as electrophoretic mobility shift assays (EMSAs) [3], DNAse footprinting [4], enzyme-linked immunosorbent assay (ELISA) [5,6], ChIPchip [7], or mutations of the putative BS and subsequent expression studies. Because TFs bind to double-stranded DNA, the strand annotations of nonpalindromic BSs in the databases are either missing or added, based on manual inspection or predictions from bioinformatics tools such as MEME [8], Gibbs Sampler [9,10], Improbizer [11], SeSiM-CMC [12], or A-GLAM [13].
First, erroneously annotated BS: This error may occur in the original study or during the transfer process from the scientific literature to the databases. A sequence is declared to contain a BS, although, in reality, it does not.
Second, shift of the BS: The BS may be erroneously shifted by one or a few base pairs. This typically happens during the transfer process from the scientific literature to the databases.
Third, missing or wrong strand orientation of the BS: The strand orientation of a BS is often not or incorrectly annotated. For example, all BS orientations are arbitrarily declared to be in 5'3' direction relative to the target gene in CoryneRegNet and in RegulonDB [14,16].
These problems can strongly affect any of the subsequent analysis steps, such as the inference of sequence motifs from "experimentally verified" data, the calculation of P values for the occurrence of BSs, the detection of putative BSs in genome-wide scans and their experimental validation, or the reconstruction of transcriptional gene-regulatory networks.
Here, we introduce MotifAdjuster, a software tool for detecting potential BS annotation errors and for proposing possible corrections. Existing bioinformatics tools [8][9][10][11][12][13] are not optimized for this task (Additional data file 1), because they do not allow shifting the BS by using a nonuniform distribution and considering both strands with unequal weights. In contrast, MotifAdjuster allows the user to incorporate prior knowledge about (i) the probability of erroneously annotated BSs, (ii) the distribution of possible shifts, and (iii) the strand preference.
One widely-used model for the representation of BSs is the position weight matrix (PWM) model [8][9][10][11][12][13]26,27], and many software tools for genome-wide scans of sequence motifs are based on PWM models [26,28,29]. MotifAdjuster is based on a simple mixture model using a PWM model on both strands for the motif sequences and a homogeneous Markov model of order 0 for the flanking sequences similar to MEME, Gibbs Sampler, Improbizer, SeSiMCMC, or A-GLAM. For a given set of BSs, MotifAdjuster tests whether each sequence contains a BS, and it refines the annotations of position and strand for each BS, if necessary, by maximizing the posterior of the mixture model by using a simple expectation maximization (EM) algorithm.
To test the efficacy of MotifAdjuster, we apply it to seven data sets from CoryneRegNet, and we record for each of them the set of potential annotation errors. For one example, the nitrate regulator NarL, we compare the proposed adjustments with the original literature, with a manual strand rean-notation of the BS strands, and with an independent and hand-curated reannotation provided by PRODORIC. Finally, we test whether the PWM estimated from the adjusted NarL BSs can help to detect unknown BSs in those promoter regions that are known to be bound by NarL, but for which no BS could be predicted in the past.

Algorithm
In this section, we present the MotifAdjuster algorithm including the mixture model, the prior, and the maximum a posteriori (MAP) estimation of the model parameters given the data.

Mixture model
We denote a DNA sequence of length L by x:= (x 1 , x 2 , ..., x L ), the nucleotide at position ᐍ  [1, L] by x ᐍ  {A, C, G, T}, and the reverse complement of x by x RC . For modeling a BS x of length w, we use a PWM model, which assumes that the nucleotides at all positions are statistically independent of each other, resulting in an additive log-likelihood of sequence x given the model parameters  [30,31], where the subscript f stands for foreground. Here, denotes the logarithm of the probability of finding nucleotide a  {A, C, G, T} at position ᐍ,  ᐍ denotes the four-dimensional vector , and  denotes the (4 × w) matrix, that is,  denotes the PWM [32][33][34][35][36].
For modeling the flanking sequences, we use a homogeneous Markov model of order 0, which assumes that all nucleotides are statistically independent, resulting in an additive log-likelihood of sequence x given model parameters  [32][33][34][35][36], where the subscript b stands for background. Here,  a denotes the logarithm of the probability of nucleotide a, and  denotes the vector ( A , ...,  T ) T .
For the detection of sequences (i) erroneously annotated as containing BSs, (ii) with shifted BSs, or (iii) with missing or wrong strand annotations, we introduce the three random variables u 1 , u 2 , and u 3 .
The variable u 1 handles the possibility that a sequence annotated as containing a BS does not contain a BS. u 1 = 0 denotes the case that the sequence contains no BS, and u 1 = 1 denotes the case that the sequence contains exactly one BS. If the sequence contains one BS, it can be located at different positions and on both strands.
The variable u 2 handles the possibility of shifts of a BS caused by annotation errors. u 2 models the start position of the BS in the sequence with respect to the annotated start position. This variable can assume the integer values {-s, -(s-1), ..., s-1, s}, where s is the maximal shift of the BS upstream or downstream of the annotated position.
The variable u 3 handles the possibility that a BS can have two orientations in the double-stranded upstream region of the target gene. According to the notation of CoryneRegNet, u 3 = 0 denotes the forward strand defined as the strand in 5'3' direction relative to the target gene, and u 3 = 1 denotes the reverse complementary strand.
For shortness of notation, we define u := (u 1 , u 2 , u 3 ). Because we do not know the values of u, these variables are modeled as hidden variables. We assume that u 2 and u 3 are conditionally independent of each other given u 1 ; that is, we assume that annotation errors of position and strand are conditionally independent given the occurrence of the BS. We define where the subscript h stands for hidden, and where f:= (f 1 , f 2 , f 3 ) denotes the vector of parameters of this distribution. MotifAdjuster allows the user to specify the probability P h (u 1 |f 1 ) that a sequence contains (or does not contain) a BS and the probability distribution P h (u 2 |u 1 , f 2 ) for the length of the erroneous shift. In addition, MotifAdjuster estimates the logarithm of the probability that the BS is located on the forward (v = 0) or the reverse complementary (v = 1) strand, , from the user-provided data as described in subsection Expectation maximization algorithm.
The hidden values of u lead to the likelihood of the data x given the model parameters (, , f), where the sum runs over all possible values of u. Here, the subscript a stands for accumulated, and the subscript c stands for composite. In the following, we define the likelihood in close analogy to [8,37]. If sequence x contains no BS, we assume that x is generated by a homogeneous Markov model of order 0; that is, If the sequence x contains a BS, then u 2 encodes its start position, u 3 encodes its strand, and we assume that the nucleotides upstream and downstream of the BS are generated by a homogeneous Markov model of order 0, yielding and where the subscript m stands for motif.

Prior
As prior of the parameters of the PWM model, we use the "common choice" [34][35][36]  The choice of this prior is pragmatic rather than biologically motivated. This prior is conjugate to the likelihood, allowing to write the posterior as a product of transformed Dirichlets. As PWM models are special cases of Bayesian networks, the chosen prior can be understood as a special case of the Bayesian Dirichlet (BD) prior [38].
Analogously, for homogeneous Markov models of order 0, we choose a transformed Dirichlet P(|) := D(|), where  a denotes the positive hyperparameter of  a .
MotifAdjuster allows the user to specify P(u 1 |f 1 ) and P(u 2 |u 1 , f 2 ). In principle, MotifAdjuster allows the user to specify any probability distribution P(u 2 |u 1 , f 2 ) for the length of the erroneous shift, allowing also asymmetric or bimodal distributions, if needed. For an easy and user-friendly execution, MotifAdjuster also offers a discrete and symmetrically truncated Gaussian distribution defined by where z is an integer value ranging from -s to s. The real-valued parameter  is similar to the standard deviation of a Gaussian distribution and can be specified by the user, and we denote f 2 := (s, ).
We expect that some sequences are annotated to contain a BS, although they do not contain a BS in reality, but we believe that the fraction of such incorrectly annotated sequences is small. Hence, we choose P(u 1 = 0|f 1 )=0.2 for the studies presented in this article; that is, we assume that only 20% of the sequences annotated to contain a BS do not contain a BS in reality. We further expect that the annotated position of the BS might be shifted accidentally by a few base pairs, so we choose s = 5 and a discrete and symmetrically truncated Gaussian distribution with  = 1. This choice results in a conditional probability of approximately 40% that the BS is not shifted, of approximately 25% that it is shifted 1 bp, and of approximately 5% that it is shifted by more than 1 bp upstream or downstream of the annotated start position, respectively, given that a BS is present in sequence x.
As prior of the parameter f 3 , we choose a transformed Dirich- Putting all pieces together, we define the prior of the parameters of the mixture model of Equation (4) by: stating that we assume , , and f 3 to be statistically independent.
We denote the ESS of the mixture model chosen before inspecting any database by , and we set the ESS of the PWM model to P(u 1 = 1|f 1 )·, the positive hyperparameters of the strand parameters to , and the ESS of the homogeneous Markov model of order 0 to (L -P(u 1 = 1|f 1 )·w)·. For the reassessment of BSs presented in this article, we choose an ESS of  = 5, yielding an ESS of 4 for the PWM model,  0 =  1 = 2, and an ESS of 57 for the homogeneous Markov model of order 0. This choice yields for every a  {A, C, G, T} and every ᐍ  [1, w], stating that the chosen prior of the PWM model can be understood as a special case of the BDeu prior [39,40], which in turn is a special case of the BD prior.

Expectation maximization algorithm
The model parameters of the mixture model defined by Equation (4) cannot be estimated analytically, but any numeric optimization algorithm can be used for maximizing the posterior. One popular optimization algorithm for maximizing the likelihood P(S|, , f) is the EM algorithm [41]. The EM algorithm can be easily modified for maximizing the posterior P(, , f|S, , , ) of the data set S by iteratively maximizing: can be maximized analytically with respect to , , and f 3 , yielding the familiar expressions provided in Additional data file 2. The posterior P(, , f|S, , , ) increases monotonically with each iteration, implying that the modified EM algorithm converges to the global maximum, a local maximum, or a saddle point. We stop the algorithm if the logarithmic increase of the posterior between two subsequent iterations becomes smaller than 10 -6 , restart the algorithm 10 times with randomly chosen initial values of , and choose the parameters of that start with the highest posterior, similar to [8,37]. If we restrict P h (u 2 |u 1 , f 2 ) to a uniform distribution over all possible start positions, if we set P h (u 3 |u 1 = 1) = 0.5, and if we restrict the background model to be strand symmetric, then we obtain the probabilistic model that is the basis of [8,37].
The flexibility allowed by MotifAdjuster is important for its practical applicability. Typically, the user has prior knowledge about (i) the expected motif occurrence and (ii) the shift distribution, but (iii) no or only limited prior knowledge about the distribution of the BS strand orientation. Hence, we allow the user to specify the logarithm of the probability that a sequence contains a BS f 1,0 , a nonuniform distribution to incorporate the prior knowledge of the shift distribution, and we estimate the logarithm of the probability that the BS is located on the forward strand f 3,0 from the data. This setting allows MotifAdjuster to work, without additional intervention, also in the two extreme cases that the BSs lie predominantly either on the forward or on the reverse complementary strand.

Case studies
In this section we present the results of MotifAdjuster applied to seven data sets of Escherichia coli, the validation of Motif-Adjuster results for NarL BSs, and the prediction of a novel NarL BS. For each of the seven TFs, we analyze whether the adjustments proposed by MotifAdjuster result in an improved motif of the BSs (Figure 1). We compute the sequence logos [54,55] of the original BSs obtained from CoryneRegNet and those of the BSs proposed by MotifAdjuster, which we call original sequence logos and adjusted sequence logos, respectively. Comparing these sequence logos, we find that the adjusted sequence logos show a higher conservation than the original sequence logos in all seven cases. We also compare the sequence logos with consensus sequences obtained from the literature [56][57][58][59][60][61], and we find that the adjusted sequence logos are more similar to the consensus sequences than the original sequence logos. In addition, we find, for the TFs CpxR, Fur, and NarL, that the adjusted sequence logos allow us to recognize clear motifs that could not be recognized in the original sequence logos obtained from CoryneRegNet.

Results for seven data sets of Escherichia coli
We investigate whether there exists any systematic dependence of the observed rate of proposed adjustments exists on the number of BSs, the BS length, and the GC content of the BSs. We find no obvious dependence of the error rate on the number of BSs and on the BS length. Comparing the GC content of the BSs, we find that the GC content of the BSs of all but one TF ranges from 30% to 40%. However, the GC content of the Fur BSs is only 20%. This low GC content might be the reason for the unexpectedly high percentage of shifts in this data set, because it is more likely to shift a BS accidentally in a sequence composed of a virtually binary alphabet.

Validation of MotifAdjuster results for NarL
To evaluate the previous results, we choose NarL as example and scrutinize the proposed reannotations of MotifAdjuster for this case. The nitrate regulator NarL of Escherichia coli is one of the key factors controlling the upregulation of the nitrate respiratory pathway and the downregulation of other respiratory chains. In the absence of oxygen, the energetically most efficient anaerobic respiratory chain uses nitrate and nitrite as electron acceptors [62]. Detection of and adaptation to extracellular nitrate levels are accomplished by complex interactions of a double two-component regulatory system, which consists of the homologous sensory proteins NarQ and NarX, and the homologous TFs NarL and NarP. Depending on the BS arrangement and localization relative to the transcription start site, NarL and NarP act as activators or repressors, thereby enabling a flexible control of the expression of nearly 100 genes.
CoryneRegNet stores 74 NarL BSs, each of length 7 bp (Table  1). Of these 74 BSs, only 36 are considered accurate by Motif-Adjuster, whereas 38 are considered to be questionable. In 25 cases, MotifAdjuster proposes to switch the strand orientation of the BS; in five cases, it proposes to shift the location of the BS, and for six BSs, it proposes both a switch of strand ori-entation and a shift of position. In addition, two BSs are proposed for removal. We present a summary of these results in Table 2, we provide a complete list of the results in Additional data file 4, and we summarize in Table 3 those 13 BSs of the regulator NarL where MotifAdjuster proposes to shift the location of the BS or to remove it from the databases.
Comparison of binding-site conservation, showing the original sequence logos, the consensus sequences for the TFs obtained from the literature [56][57][58][59][60][61], and the adjusted sequence logos for the data sets of the TFs CpxR, Crp, Fis, Fnr, Fur, Lrp, and NarL Figure 1 Comparison of binding-site conservation, showing the original sequence logos, the consensus sequences for the TFs obtained from the literature [56][57][58][59][60][61], and the adjusted sequence logos for the data sets of the TFs CpxR, Crp, Fis, Fnr, Fur, Lrp, and NarL. We find in all seven cases that (i) the adjusted sequence logos show a higher conservation than the original sequence logos, (ii) the adjusted sequence logos are more similar to the consensus sequences than to the original sequence logos; and (iii) clear motifs can be recognized in the adjusted sequence logos of the TFs CpxR, Fur, and NarL that could not be recognized in the original sequence logos. To evaluate the accuracy of MotifAdjuster, we check the original literature [63,37] for each of the 13 questionable BS candidates. Comparing both, we find that the proposed annotations agree with those in the literature in all cases but one (BS of gene b1224). That is, in 12 of 13 cases signaled by MotifAdjuster as being questionable, the detected error was indeed caused by an inaccurate transfer from the original literature into the gene-regulatory databases RegulonDB and CoryneRegNet. Of those 12 questionable BSs, 10 BSs are correctly proposed to be shifted, and two are correctly proposed to be removed.

CpxR
Turning to the BS of the gene b1224, we find it is published as given in the databases [64], in contrast to the proposal of MotifAdjuster. However, Darwin et al.
[67] report that a mutation of this BS has little or no effect on the expression of b1224. Hence, the proposal could possibly be correct, and the BS could be shifted or even be deleted.
In addition, MotifAdjuster checks the strand annotation of BSs and proposes strand switches if needed. To validate these annotations, we cannot use the annotations from RegulonDB and CoryneRegNet, because these databases contain all BSs in 5'3' direction relative to the target gene. Hence, we consult annotation experts at the Center for Biotechnology in Bielefeld to reannotate the strand orientation of the BSs manually, and we compare the results with those of MotifAdjuster. Interestingly, we find that the strand orientations proposed by MotifAdjuster are in perfect (100%) agreement with the manually-curated strand orientations. As an inde-pendent test of the efficacy of MotifAdjuster for NarL BSs, we use the manually annotated BSs provided by the PRODORIC database [68]. Remarkably, we find also in this case that the results of MotifAdjuster perfectly agree with the annotations.
Another hint that the proposed adjustments of MotifAdjuster could be reasonable is based on the observation that NarL and NarP homodimers bind to a 7-2-7' BS arrangement [61], an inverted repeat structure consisting of a BS on the forward strand, a 2-bp spacer, and a BS on the reverse complementary strand. NarP exclusively binds as homodimer to this 7-2-7' structure. NarL homodimers bind at 7-2-7' sites with highaffinity, but NarL monomers can also bind to a variety of other heptamer arrangements. Instances of this 7-2-7' structure have been reported for four genes: fdnG, napF, nirB, and nrfA [61,65]. In contrast to this observation, all BSs in Cory-neRegNet as well as RegulonDB are annotated to be on the forward strand, including the second half of the inverted repeat. When applied to these four genes, MotifAdjuster proposes all heptamers of the second half of the 7-2-7' structure to be switched to the reverse strand, in agreement with [61,65]. In addition, MotifAdjuster proposes six additional 7-2-7' BS arrangements, located in the upstream regions of the genes adhE, aspA, dcuS, frdA, hcp, and norV. The positions and the orientations are presented in Additional data file 4. . We extract the upstream regions of these genes, where an upstream sequence is defined by CoryneRegNet as the sequence between positions -560 bp and +20 bp relative to the first position of the annotated start codon of the first gene of the target operon. In addition, we extract those upstream regions of Escherichia coli that belong to operons not annotated as being regulated by NarL (background data set).
We investigate whether we can now detect NarL BSs based on the adjusted data set that could not be detected based on the original data set from CoryneRegNet. For that purpose, we estimate the parameters  of the PWM model on the adjusted data set as proposed by MotifAdjuster and  of the homogeneous Markov model on the background data set. From the adjusted PWM, we build a mixture model over both strands with the same probability for each strand; that is, exp(f 3,0 ) = exp(f 3,1 ) = 0.5. For the classification of an unknown heptamer x, we build a simple likelihood-ratio classifier with these parameters , , f 3 and define the log-likelihood ratio by For an upstream region, we compute r max defined as the highest log-likelihood ratio of any heptamer x in this upstream region. We compute the P value of a potential BS x with value r(x) as fraction of the background sequences whose r max -values exceed r(x).
With this classifier, a significant NarL BS can now be detected in the upstream region of torC. Figure 2a shows the doublestranded DNA fragment with the predicted BS (TACCCCT) located on the forward strand starting at -209 bp relative to the start codon, and at -181 bp relative to the annotated transcription start site [76]. The distance of the predicted BS to the start codon agrees with the distance distribution of previously known NarL BS (Figure 2b), providing additional evidence for the predicted BS. This finding closes the gap between sequence-analysis and gene-expression studies, as the torCAD operon consists of three genes that are essential for the trimethylamine N-oxide (TMAO) respiratory pathway [76]. TMAO is present as an osmoprotector in tissues of invertebrates and can be used as respiratory electron acceptor by Escherichia coli.

Conclusions
Gene-regulatory databases, such as AGRIS, AthaMap, Cory-neRegNet, CTCFBSDB, JASPAR, ORegAnno, PRODORIC, RegulonDB, SCPD, TRANSFAC, TRED, or TRRD store valuable information about gene-regulatory networks, including TFs and their BSs. These BSs are usually manually extracted from the original literature and subsequently stored in databases. The whole pipeline of wet-lab BS identification and annotation, publication, and manual transfer from the scientific literature to data repositories is not just time consuming but also error prone, leading to many false annotations currently present in databases.
MotifAdjuster is a software tool that supports the (re-)annotation process of BSs in silico. It can be applied as a qualityassurance tool for monitoring putative errors in existing BS repositories and for assisting with a manual strand annotation. MotifAdjuster maximizes the posterior of the parameters of a simple mixture model by considering the possibilities that (i) a sequence being annotated as containing a BS in reality does not contain a BS; (ii) the annotated BS is erroneously shifted by a few base pairs; and (iii) the annotated BS is erro-Position of the predicted NarL binding site in the upstream region of torC neously located on the false strand and must be reverse complemented. In contrast to existing de-novo motif-discovery algorithms, MotifAdjuster allows the user to specify the probability of finding a BS in a sequence and to specify a nonuniform shift distribution.
We apply MotifAdjuster to seven data sets of BSs for the TFs CpxR, Crp, Fis, Fnr, Fur, Lrp, and NarL with a total of 536 BSs, and we find 51 BSs proposed for removal and 134 BSs proposed for shifts. In total, this results in 34.5% of the BSs being proposed for adjustments. We choose NarL as an example to scrutinize the proposed reannotations of MotifAdjuster. Checking the original literature for each of the 13 cases shows that the proposed deletions and shifts of MotifAdjuster are in agreement with the published data. Comparing the strand annotation of MotifAdjuster with independent information indicates that the proposals of MotifAdjuster are in accordance with human expertise. Furthermore, MotifAdjuster enables the detection of a novel BS responsible for the regulation of the torCAD operon, finally augmenting experimental evidence of its NarL regulation.
MotifAdjuster is an open-source software tool that can be downloaded, extended easily if needed, and used for computational reassessments of BS annotations.