Quantifying similarity between motifs
© Gupta et al.; licensee BioMed Central Ltd. 2007
Received: 13 September 2006
Accepted: 26 February 2007
Published: 26 February 2007
A common question within the context of de novo motif discovery is whether a newly discovered, putative motif resembles any previously discovered motif in an existing database. To answer this question, we define a statistical measure of motif-motif similarity, and we describe an algorithm, called Tomtom, for searching a database of motifs with a given query motif. Experimental simulations demonstrate the accuracy of Tomtom's E values and its effectiveness in finding similar motifs.
Discovering and characterizing DNA and protein sequence motifs are fundamental problems in computational biology. Here, we use the term 'motif' to refer to a position-specific probability matrix that describes a short sequence of amino acids or nucleotides that is important to the functioning of the cell. For example, the regulation of transcription requires sequence-specific binding of transcription factors to certain cis-acting motifs, which typically are located upstream of transcriptional start sites . On the other hand, protein sequence motifs might correspond to active sites in enzymes or to binding sites in receptors .
We are not the first to describe a method for quantifying the similarities between pairs of motifs. Pietrokovski  compared protein motifs using a straightforward algorithm based on the Pearson correlation coefficient (PCC). Subsequently, Hughes and coworkers  applied a similar method to DNA motifs. Wang and Stormo  introduced an alternate motif column comparison function, termed the average log-likelihood ratio (ALLR). More recently, Schones and coworkers  introduced two motif similarity functions, one based on the Pearson χ2 test and the other on the Fisher-Irwin exact test (FIET). They showed that these two new functions have better discriminative power than the PCC and ALLR similarity functions. In addition, multiple research groups have used Kullback-Leibler divergence (KLD) to compare motifs [11–13], and Choi and coworkers  used euclidean distance (ED) to compare protein profiles. Finally, Sandelin and Wasserman  used their own column comparison function (SW) within the context of a dynamic programming alignment approach to compare DNA motifs. This method differs significantly from all other DNA-motif based approaches in the sense that it allows gaps in the motif-motif alignments.
In this report we focus on ungapped alignments of motifs. We describe a general method for accurately modeling the empirical null distribution of scores from an arbitrary, additive column comparison function. We estimate the null distribution of scores for each column in a 'query' motif using the observed scores of aligning it with each motif column in a database of 'target' motifs. Using a dynamic programming algorithm inspired by earlier work on searching a sequence database with a motif [16–18], we estimate the null distribution of the sum of scores for any range of contiguous columns in the query motif. This makes it possible for the user to determine whether the motif comparison score between the query motif and a particular target motif is statistically significant. Previous methods begin by defining a score between two motif columns, and then they combine these scores either by summing (as we do) [7–9, 14] or by taking the mean [11–13] or geometric mean  of the column scores. Our scoring method differs in that it computes the P values of the match scores for the columns of the query motif aligned with a given target motif in all possible ways (without gaps). These 'offset' P values are computed using the cumulative density functions estimated from the target database, as described above. The minimum P value among these offset P values is used to compute the overall P value of the match between the query motif and the target motif, assuming independence of the offset P values. This is called the 'motif' P value. Finally, we apply a Bonferroni correction to the motif P values to derive an E value.
This algorithm is implemented in a software tool called Tomtom, which is publicly available as part of the MEME Suite of motif analysis tools [19–21]. Tomtom can compute E values based on any one of seven column comparison functions: PCC, ALLR, PCS, FIET, KLD, ED, or SW. In this work, we demonstrate the accuracy of Tomtom's statistical estimates. We also validate Tomtom'smotif retrieval accuracy via a simulation experiment. The results show that, in addition to providing formal semantics for motif similarity scores, Tomtom's P value estimation yields improved rankings relative to ad hoc normalization schemes.
In this section, we describe the motif-motif comparison problem and outline our solution. Say we are given two motifs, Q and T. Our goal is to define a motif comparison function S(·,·), such that S(Q,T) is small if and only if Q and T resemble one another in some biologically relevant way. For now, let us sidestep the issue of defining 'biologically relevant' and assume that someone has given us a function s(·,·) that compares two motif columns. Thus, we can compare, for example, the ith column of Q and the jth column of T using s(Q i ,T j ). Our problem is to use the column comparison function s(·,·) to define the motif similarity function S(·,·).
This problem can be further subdivided into two subproblems. One subproblem is that we do not know a priori whether the motifs Q and T should be offset with respect to one another. Indeed, in the case of DNA motifs, we often do not even know whether the motifs lie on the same DNA strand. Therefore, our motif similarity function must take into account all possible offsets and relative orientations. A second subproblem is that even if we knew the correct offset and relative orientation, we need a method for combining the column comparison scores into a single score. Below, we describe solutions to each of these problems.
Computing offset Pvalues
Initially, let us simplify our problem even further. Not only has someone told us the correct column-wise similarity function s(·,·), but they have also specified the correct relative offset and orientation of our motifs Q and T. For now, we assume that the motifs are of equal width w, that they lie in the same orientation, and that they have a relative offset of zero. Furthermore, we assume that columns of the motifs are independent and that our scores can be summed. Our problem is to compute a P value for this summation. Because the P value is relative to the given offset, we refer to this as the 'offset P value'. We adopt a dynamic programming method to calculate the null distribution of summed similarity scores with respect to the motif Q.
A similar method has been used to compute a P value for the match between a motif and a given sequence . Briefly, that method can be described as follows. Say that we have a motif Q of width w, and we have a score function ŝ(Q i ,a) that yields a positive integer score for the similarity of the ith column of Q and the letter a ∈ A. These integral scores correspond to indices x of an array A defining the desired probability density function (PDF). A is filled recursively by noting that, if we know the PDF A(i) for matches to the first i positions in Q, then we can calculate the PDF A(i+1) as follows:
Where P a is the null probability of letter a. The recursion is initialized with A(0)(0) = 1 and A(0)(x) = 0 for x ≥ 1. Iterating with i = 1 ... w yields the PDF for a random sequence matching the motif, which is used to calculate a cumulative probability distribution and thus P values. The challenge in generalizing the above algorithm to the motif-motif comparison problem arises because we do not have a fixed alphabet of amino acids or nucleotides for the summation in Eqn 1. Instead, we have an infinite 'alphabet' of motif columns. Our solution involves constructing an implicit alphabet of motif columns from the distribution of scores between all query motif columns versus all columns in a database of target motifs. This is an efficient solution because the matrix of query-versus-target motif column scores must be computed during the database search procedure.
In detail, the algorithm proceeds in five steps. First, for a given motif Q of width w Q and a given collection of target motifs T1 ... T n whose total width is w T , we compute a w Q -by-w T matrix Γ such that Γi,j= s(Q i ,T j ). This matrix constitutes the null distribution for our P value calculation. Second, we linearly rescale the values in Γ such that the minimum value is 1 and the maximum is t, where t is the (user-specified) number of letters in the motif column alphabet (in Tomtom, t = 100). We then round the values in Γ to integers. Third, for each column i in Q and for each possible scaled, integer score 1 ≤ x ≤ t, we compute the frequency of x in the null distribution of the ith column of Γ:
Where δ(·) is the Kronecker delta function. In the fourth step, we initialize a PDF A(0), as described above, and then perform the recursion as follows:
Computing motif Pvalues
The above procedure yields a P value for a query and target motif with a particular offset and relative orientation. In order to compute a motif P value, Tomtom identifies the offset and relative orientation for which the offset P value is minimal. The probability of observing a minimum P value of P* among a collection of N independent P values is 1 - (1 - P*) N . This value is the motif P value.
Tomtom searches a target database of motifs using a given motif as the query. The resulting motif P values must therefore be corrected for multiple tests. Tomtom uses a form of Bonferroni correction that assumes that the targets are independent of one another. The correction consists of multiplying the motif P value by the number of targets in the database. The result is an E value - the expected number of times that the given query would be expected to match a target as well or better than the observed match in a randomized target database of the given size.
We perform three separate experiments to assess the validity of Tomtom's statistical confidence estimates and measure Tomtom's ability to recognize related motifs.
Assessing Pvalue accuracy
Measuring retrieval accuracy
Next, we designed a simulation experiment to test Tomtom's ability to retrieve a related target motif from a database. The experiment is designed to simulate the following situation. Suppose that a researcher discovers a 'new' motif that is actually the same as one in a motif database. The new motif may contain some of the same sequences that were used to create the database motif plus some new sequences. Moreover, the exact boundaries of the novel motif may not exactly match the boundaries of the corresponding motif in the database. We simulate this situation, and then measure Tomtom's ability to identify the correct motif in the database.
In detail, the experiment proceeds as follows. We begin by selecting all 107 motifs in the JASPAR database (jaspar core). Then, we simulate a collection of 10 query and 10 target motifs for each of these JASPAR motifs by subsampling with replacement from the original sites of the JASPAR motif. The difficulty of the retrieval task can be modulated by reducing the number of sites sampled. In our first experiment, if the JASPAR motif has S associated sites, then the query and target motifs are simulated using S/2, S/4, S/8 or S/16 sites. In this step, we eliminate motifs that would yield fewer than two sampled sites, thus leaving 82 JASPAR motifs for the S/8 subsampling experiment. In the next step, we trim the edges of half of the motifs in each database. The number of columns to be deleted from a given motif is determined by selecting a random number uniformly from [-⌊0.8w⌋, ..., 0, 0, ..., ⌊0.8w⌋], where w is the motif width. The sign of the selected number determines which end of the motif is truncated. After this procedure, each motif in the JASPAR database has 10 corresponding motifs in the query and in the target database. Tomtom's task, given one of the query motifs, is to retrieve all 10 corresponding target motifs with smaller motif P values than any of the unrelated target motifs.
This result is encouraging; however, sorting the results from all of the queries into a single list is somewhat unrealistic. In practice, the user is only concerned with the quality of the ranking with respect to one motif at a time. We therefore compute ROC curves separately for each query motif. In order to quantify the extent to which the correct pairs appear near the top of the ranked list, we compute the area under each ROC curve (the ROC score). A perfect ranking would receive an ROC score of 1.0, whereas a random ranking would receive an ROC score of 0.5.
Mean ROC scores for various motif column comparison functions and score combination methods
Comparison of motif P values with other methods of combining column scores
P value versus sum
1.87 × e-03
P value versus AM
3.55 × e-13
2.43 × e-13
2.35 × e-13
6.93 × e-14
2.05 × e-13
1.38 × e-12
2.19 × e-13
P value versus GM
7.00 × e-11
2.52 × e-13
2.04 × e-12
2.06 × e-10
1.09 × e-13
1.58 × e-12
1.98 × e-13
E-value based retrieval rate
In a third experiment, we tested the utility of the E values computed by Tomtom. This experiment was conducted to determine whether, using a reasonable significance threshold, Tomtom can successfully retrieve a JASPAR motif from the database. In this experiment, it is not sufficient for the correct target motif to have the best score; the score must also be statistically significant.
Tomtom is a motif comparison algorithm that ranks the target motifs in a given database according to the estimated statistical significance of the match between the query and the target. In this work we show that the motif P values computed by Tomtom are accurate, in the sense that they are uniformly distributed when computed on randomized data. We also show that the P value calculation produces rankings that are significantly better than the rankings produced by ad hoc normalization schemes. It is important to emphasize, however, that even if the rankings produced by Tomtom were no better than ad hoc rankings produced, P value normalization would still be the preferred method because of the inherent advantages of having a measure of the statistical significance of query-target matches. Finally, we show that Tomtom correctly assigns E values less than 0.01 to a large percentage of positive matches. This result indicates that it is highly probable that Tomtom successfully retrieves a related motif with a significant E value. All of these properties make Tomtom a valuable tool for identifying truly related motifs.
During the course of our experiments, we compared seven different motif column comparison functions. Surprisingly, the simple ED between motif columns performs best. Consequently, Tomtom's default behavior is to compare columns using ED. However, for some types of motifs (for instance, protein motifs) other comparison functions may be more appropriate. Consequently, Tomtom provides an option to use any of the seven column comparison functions.
In terms of practical applicability, Tomtom is especially relevant in conjunction with MEME, an ab initio motif discovery tool. Novel motifs identified using MEME can be reliably searched against known motifs using Tomtom. Both Tomtom and MEME are currently available as part of the MEME Suite of motif analysis tools [19, 20], and a Tomtom website is under development.
Materials and methods
Motif column comparison functions
At Tomtom's core is a function that defines the similarity between one column of the query and one column of the target motif. Tomtom implements seven such functions, described below. The 'raw' score for an ungapped alignment of columns from a query motif and a target motif is computed by summing the column comparison scores computed using any of the following functions. Tomtom converts the raw scores into P values and E values, as described above.
In the following discussion, X refers to a column of the query motif, and is a multinomial probability vector. The quantity X a refers to the probability of letter a ∈ A in X. For some of these functions, these probabilities are multiplied by a motif-dependent constant to give the 'counts' of different letters in each column of the motif. We use N Xa to refer to the count of letter a in column X. Similar definitions apply for Y, a column from the target motif. The quantity |A| refers to the length of the motif alphabet (four for DNA, 20 for proteins).
Pearson correlation coefficient
The PCC was first introduced for computing motif-motif similarity by Pietrokovski . For two columns X and Y, PCC is computed using the following formula:
The latter two expressions reduce to , because
for multinomial probability vectors X and Y.
Average log-likelihood ratio
The ALLR formula described by Wang and Stormo  to quantify similarity between columns X and Y for position specific weight matrix motifs is as follows:
where P a is the background (prior) frequency of letter a.
Contingency table (with margins) for DNA motifs
Where N o ja = N ja is the 'observed' count of letter a in column j, and N e ja = N j N a /N is the 'expected' count of letter a at column j (Table 3 for definitions).
The P value is calculated from this χ2 score using |A| - 1 degrees of freedom. Because our null hypothesis is that these two columns are derived from the same multinomial distribution, a higher P value implies similarity. This P value is treated as an additive score.
Fisher-Irwin exact test
The FIET  is an analytical computation of the Pearson χ2 P value. In particular, this calculation is important when marginal frequencies are small, which is often the case in position frequency matrices. The marginal P value of the contingency table for DNA motifs (Table 3) follows the multiple hypergeometric distribution :
The formula for protein motifs is similar. The two-sided P value for the table is the sum of probabilities of all tables that are at least as extreme. This P value is computed using the algorithm described by Mehta and Patel . As with the χ2 test, this P value is used as an additive score.
Choi and coworkers  introduced the ED as a means to compare protein motifs. The ED for two DNA profile columns X and Y is computed using the following formula:
Sandelin-Wasserman similarity function
Sandelin and Wasserman  introduced their own motif column comparison function for the construction of familial binding profiles. The SW score for two columns X and Y is computed using the following formula:
Additional data files
The following additional data are available with the online version of this paper. Additional data file 1 is a figure showing the accuracy of motif comparison P values. Additional data file 2 is a figure showing the accuracy of offset P values. Additional data file 3 is a table summarizing the mean ROCs for various motif column comparison functions and score combination methods for various sampling rates. Additional data file 4 is a table comparing motif P values with other methods of combining column scores for various sampling rates. Additional data file 5 is a figure showing the motif retrieval accuracy for various column similarity functions at a sampling rate of S/8. Additional data file 6 is a figure whosing the E value based retrieval rate for two additional significance levels (E value less than 0.05 or 0.001).
- Maniatis T, Goodbourn S, Fischer JA: Regulation of inducible and tissue-specific gene expression. Science. 1987, 236: 1237-1245. 10.1126/science.3296191.PubMedView ArticleGoogle Scholar
- Pawson T, Nash P: Assembly of cell regulatory systems through protein interaction domains. Science. 2003, 300: 445-452. 10.1126/science.1083653.PubMedView ArticleGoogle Scholar
- Tompa M, Li N, Bailey T, Church G, Moor BD, Eskin E, Favorov A, Frith M, Fu Y, Kent W, et al: Assessing computational tools for the discovery of transcription factor binding sites. Nat Biotechnol. 2005, 23: 137-144. 10.1038/nbt1053.PubMedView ArticleGoogle Scholar
- Sandelin A, Alkema W, Engstrom P, Wasserman W, Lenhard B: JASPAR: an open-access database for eukaryotic transcription factor binding profiles. Nucliec Acids Res. 2004, 32: D91-D94. 10.1093/nar/gkh012.View ArticleGoogle Scholar
- Wingender E, Chen X, Hehl R, Karas H, Liebich I, Matys V, Meinhardt T, Pruss M, Reuter I, Schacherer F: TRANSFAC: an integrated system for gene expression regulation. Nucleic Acids Res. 2000, 28: 316-319. 10.1093/nar/28.1.316.PubMedPubMed CentralView ArticleGoogle Scholar
- Henikoff S, Henikoff JG: Protein family classification based on searching a database of blocks. Genomics. 1994, 19: 97-107. 10.1006/geno.1994.1018.PubMedView ArticleGoogle Scholar
- Pietrokovski S: Searching databases of conserved sequence regions by aligning protein multiple-alignments. Nucleic Acids Res. 1996, 24: 3836-3845. 10.1093/nar/24.19.3836.PubMedPubMed CentralView ArticleGoogle Scholar
- Hughes JD, Estep PW, Tavazoie S, Church GM: Computational identification of cis-regulatory elements associated with groups of functionally related genes in Saccharomyces cerevisiae. J Mol Biol. 2000, 296: 1205-1214. 10.1006/jmbi.2000.3519.PubMedView ArticleGoogle Scholar
- Wang T, Stormo GD: Combining phylogenetic data with co-regulated genes to to identify regulatory motifs. Bioinformatics. 2003, 19: 2369-2380. 10.1093/bioinformatics/btg329.PubMedView ArticleGoogle Scholar
- Schones DE, Sumazin P, Zhang MQ: Similarity of position frequency matrices for transcription factor binding sites. Bioinformatics. 2005, 21: 307-313. 10.1093/bioinformatics/bth480.PubMedView ArticleGoogle Scholar
- Roepcke S, Grossmann S, Rahmann S, Vingron M: T-Reg Comparator: an analysis tool for the comparison of position weight matrices. Nucleic Acids Res. 2005, 33: W438-W441. 10.1093/nar/gki590.PubMedPubMed CentralView ArticleGoogle Scholar
- Thijs G, Marchal K, Lescot M, Rombauts S, Moor BD, Rouze P, Moreau Y: A Gibbs sampling method to detect overrepresented motifs in the upstream regions of coexpressed genes. J Comput Biol. 2002, 9: 447-464. 10.1089/10665270252935566.PubMedView ArticleGoogle Scholar
- Aerts S, Loo PV, Thijs G, Moreau Y, Moor BD: Computational detection of cis-regulatory modules. Bioinformatics. 2003, 19: ii5-ii14. 10.1093/bioinformatics/btg1052.PubMedView ArticleGoogle Scholar
- Choi I, Kwon J, Kim S: Local feature frequency profile: a method to measure structural similarity in proteins. Proc Natl Acad Sci USA. 2004, 101: 3797-3802. 10.1073/pnas.0308656100.PubMedPubMed CentralView ArticleGoogle Scholar
- Sandelin A, Wasserman WW: Constrained binding site diversity within families of transcription factors enhances pattern discovery bioinformatics. J Mol Biol. 2004, 338: 207-215. 10.1016/j.jmb.2004.02.048.PubMedView ArticleGoogle Scholar
- Staden R: Searching for motifs in nucleic acid sequences. Methods Mol Biol. 1994, 25: 93-102.PubMedGoogle Scholar
- Tatusov RL, Altschul SF, Koonin EV: Detection of conserved segments in proteins: iterative scanning of sequence databases with alignment blocks. Proc Natl Acad Sci USA. 1994, 91: 12091-12095. 10.1073/pnas.91.25.12091.PubMedPubMed CentralView ArticleGoogle Scholar
- Bailey TL, Gribskov M: Methods and statistics for combining motif match scores. J Comput Biol. 1998, 5: 211-221.PubMedView ArticleGoogle Scholar
- Bailey TL, Elkan CP: Fitting a mixture model by expectation-maximization to discover motifs in biopolymers. Proceedings of the Second International Conference on Intelligent Systems for Molecular Biology. Edited by: Altman R, Brutlag D, Karp P, Lathrop R, Searls D. 1994, AAAI Press, Menlo Park, CA, 28-36.Google Scholar
- Grundy WN, Bailey TL, Elkan CP, Baker ME: Meta-MEME: motif-based hidden Markov models of protein families. Comp Appl Biosci. 1997, 13: 397-406.PubMedGoogle Scholar
- MEME: a suite of motif analysis tools. [http://meme.sdsc.edu]
- Tomtom: a program to quantify motif-motifsimilarity. [http://noble.gs.washington.edu/proj/tomtom]
- Hanley JA, McNeil BJ: The meaning and use of the area under a receiver operating characteristic (ROC) curve. Radiology. 1982, 143: 29-36.PubMedView ArticleGoogle Scholar
- Agresti A: A survey of exact inference for contingency tables. Stat Sci. 1992, 7: 131-177.View ArticleGoogle Scholar
- Mehta CR, Patel NR: Algorithm 643 FEXACT: a FORTRAN subrouting for Fisher's exact test on unordered rXc contingency tables. ACM Trans Mathematical Software. 1986, 12 (June):Google Scholar
- Crooks GE, Hon G, Chandonia JM, Brenner SE: WebLogo: a sequence logo generator. Genome Res. 2004, 14: 1188-1190. 10.1101/gr.849004.PubMedPubMed CentralView ArticleGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.