Estimation and efficient computation of the true probability of recurrence of short linear protein sequence motifs in unrelated proteins

Background Large datasets of protein interactions provide a rich resource for the discovery of Short Linear Motifs (SLiMs) that recur in unrelated proteins. However, existing methods for estimating the probability of motif recurrence may be biased by the size and composition of the search dataset, such that p-value estimates from different datasets, or from motifs containing different numbers of non-wildcard positions, are not strictly comparable. Here, we develop more exact methods and explore the potential biases of computationally efficient approximations. Results A widely used heuristic for the calculation of motif over-representation approximates motif probability by assuming that all proteins have the same length and composition. We introduce pv, which calculates the probability exactly. Secondly, the recently introduced SLiMFinder statistic Sig, accounts for multiple testing (across all possible motifs) in motif discovery. However, it approximates the probability of all other possible motifs, occurring with a score of p or less, as being equal to p. Here, we show that the exhaustive calculation of the probability of all possible motif occurrences that are as rare or rarer than the motif of interest, Sig', may be carried out efficiently by grouping motifs of a common probability (i.e. those which have permuted orders of the same residues). Sig'v, which corrects both approximations, is shown to be uniformly distributed in a random dataset when searching for non-ambiguous motifs, indicating that it is a robust significance measure. Conclusions A method is presented to compute exactly the true probability of a non-ambiguous short protein sequence motif, and the utility of an approximate approach for novel motif discovery across a large number of datasets is demonstrated.


Background
Short Linear Motifs (SLiMs) are abundant and ubiquitous recurring protein subsequences that function largely independently of the tertiary structure of the containing protein. SLiMs mediate tasks such as cell signalling, act as recognition sites for protein modification and target proteins for subcellular localisation [1,2]. SLiMs are short (typically between three and ten amino acids in length) and degenerate (positions are often flexible in terms of possible amino acids) making motif context important for specificity due to the limited number of residues in the interaction interface [3]. This simplicity gives them an evolutionary plasticity that is unavailable to globular protein domains, enabling de novo motifs to evolve convergently through point mutations, thereby adding new functions to proteins [4]. SLiMs preferentially occur in natively disordered regions of proteins [5] challenging the long standing dogma that structure equals function [6]. With 25% -30% of the human proteome natively disordered and large portions of the interactome unexplained by known domaindomain interactions [7] there is vast potential for novel SLiM discovery [8].
Increased knowledge of SLiM attributes, through the study of known functional motifs, has enabled advancements in computational methods for SLiM discovery. This culminated in the discovery of a Translin-binding motif, the first novel functional motif to be discovered by computational methods [9]. These methods can be split into two distinct classes; rediscovery of known SLiMs and prediction of novel SLiMs. SLiM rediscovery, which was pioneered by PROSITE [10], uses regular expression or profile matching to search for novel instances of previously known SLIMs. Tools such as the ELM [1] and MnM [11] servers use context information such as globularity, conservation or residue accessibility to improve confidence in returned motifs. Motif rediscovery techniques have also been used to search for enrichment of known SLiMs in proteins that are associated with particular protein functions, using statistical enrichment to define motifs important for these functions. These methods have been used to discover novel instance of both KEN box and EH1 transcriptional repressor motifs [12,13].
Several methods have also been suggested for the discovery of novel SLiMs. The most successful of these methods seek to find motifs that are over-represented in evolutionarily unrelated proteins sharing a common attribute (such as interacting with a common proteinbinding domain [14], those which have a common posttranslational modification [15] or which localise to the same sub-cellular location [16]). Dilimot [17] and SLiM-Disc [18] use this technique, to build upon the scoring schemes of ASSET [19] and PRATT [20] respectively, to successfully discover new, and rediscover known, functional motifs. However, these first generation methods have two important deficiencies: (a) their scores are biased on motif length and dataset size, making them incomparable across multiple datasets and between motifs of different lengths; and (b) they do not offer an indication of the likelihood that any motif could reach such a score by chance.
In previous work, we introduced SLiMFinder [21], a probabilistic method for SLiM discovery that heuristically accounts for these shortcomings with a two-step scoring scheme. Both steps of this scoring scheme use simplifying assumptions that are known to be violated in real data. Given that motif probability calculations are relatively approximate (it is difficult to model exactly the amino acid or word count background distributions in a set of proteins or protein regions under study) these approximations seem at first glance to be satisfactory. However, in practise we noted in studies searching for motifs across many datasets that consistent biases were emerging, that related to the complexity of the motif, and to the size of the dataset searched. Since computational motif discovery is most usually struggling to identify a weak signal against a background of noise, we were motivated to try and eliminate these biases by improving the statistical calculations.
The first step of the SLiMFinder scoring scheme applies the binomial probability mass function to mean success probabilities across all proteins, denoted here as p μ , to calculate the probability of a given motif occurring by chance in the number of proteins it occurs in, or greater. Therefore, this assumes that the proteins in a search dataset all have the same length and composition. Here, we introduce a more accurate calculation, p v which incorporates unique success probabilities for each protein, allowing for its length and composition. The advantage of correcting for this effect is that the p-value cut-offs from datasets that have varying sizes of proteins become more comparable, as well as being more realistic.
The second step of the SLiMFinder scoring scheme computes an approximate significance, Sig. This is a dataset-based score, which accounts for the multiple testing inherent in motif discovery, by calculating the probability that any motif in the dataset will be returned with that p-value or less by chance. However, Sig approximates the probability of all other possible motifs occurring with a score of p or less as p. Here, we introduce the true calculation, Sig', that replaces this approximation, and present an efficient algorithm for its calculation. The advantage of correcting for this effect is that motifs of different length, or containing amino acids of very different frequency, are treated equally, so that there is not a bias towards any class of motif.
We investigate the impact of replacing both approximations with the more exact calculations on the overall accuracy of the scoring scheme. Each step is independent allowing them to be implemented and tested separately, four significance scoring schemes are discussed (Table 1), investigating the trade-off between accuracy and efficiency. We discuss the optimal approaches for large-scale studies of multiple datasets.

Methods
Motif discovery statistics can be defined by 3 basic probabilities; p 1+ , the probability that there are 1 or more occurrences of a motif in a protein; p, the probability in a dataset of N proteins that a given motif would occur by chance with its observed support or higher; and Sig, the probability that any motif would reach p or less by chance. In this section, each probability will be introduced and discussed, focussing on the speed and accuracy of the statistics, and an efficient and accurate method of calculating Sig is described.
Exact probability of a given motif support considering unique success probabilities Previous work: Mean Success Probability Heuristic The probability, p 1+ , that there are one or more occurrences of a motif in a protein is calculated differently by SLiMFinder and Dilimot. The former calculates p 1+ based on amino acid frequencies, whereas the later calculates p 1+ based on word counts. The merits of each method have been discussed previously [21]. Here, the calculation of p 1+ , is based on amino acid frequencies and allows ambiguous/degenerate positions to be evaluated. where the success probability, p motif , is the probability that the motif will occur at any position in the protein, n is the number of positions in the protein that the motif can occur, m is the motif, l is the number of nonwildcard positions in the motif, m i is position i in motif m, x is the number of degenerate/ambiguous possibilities at position i, m ik is the k th ambiguous possibility at position i, f(m ik ) is the background frequency of the amino acid m ik and B is the binomial function. The calculation is independent of alphabet and allows the use of terminal characters "^" and "$" denoting N-and Ctermini respectively, which is taken into account in calculating motif significance. Equation 1 considers both amino acid composition and the length of the protein to calculate the probability. By considering the length (the variable n accounts for the number of positions in a protein at which a motif can occur) the calculations allow short proteins/ peptides (for example peptides returned by phage display screens) to be scored without bias. The amino acid composition allows for any background amino acid frequency to be incorporated, including the distribution from within a protein, adding flexibility to the scoring scheme.
Generally, each protein (or cluster of proteins whenever several homologous proteins are grouped or weighted to correct dependencies introduced by divergent evolution [22]) will have a unique success probability, p 1+ . Previously, the mean success probability for the motif to occur in any given protein in the dataset, p 1+μ , was approximated as the mean of the p 1+ values.
. Approximate p (1 or more instances of a given motif in any given protein) The probability of k or more proteins in a dataset containing a given motif one or more times, p μ , may be estimated using I, the incomplete beta function, the cumulative distribution function of the binomial probability mass function.
Equation 3. Approximate p(k or more proteins with 1 or more occurrences of a given motif) where k is the support of the motif (i.e. the number of proteins containing it), N is the number of proteins in the dataset and p 1+μ is the mean success probability of a motif occurring in any protein in the dataset.

More exact calculation
The binomial function calculation in Equation 3, assumes all success probabilities are equal to p 1+μ (Equation 2). The binomial equation can be expanded to allow unique success probabilities for each protein, p v , by summing the product of success and failure probabilities for each possible combination of k successes, t, and N-k failure probabilities, f. This computation increases the number of calculations necessary N C k fold.
. Set of all possible combinations of k successes from N trials . Exact p(k or more proteins with 1 or more occurrences of a given motif)

Example
The error associated with p μ is most easily explained in terms of occasions where p 1+ for a protein is 0, a frequent occurrence when masking techniques are used to improve motif discovery (e.g., a protein entirely lacking W has zero probability of having the motif WxxS). Consider a set of 6 proteins, with success probabilities for a motif p 1+ of {0.1, 0.1, 0.1, 0.0, 0.0, 0.0}. 3 proteins contain an occurrence of the motif. The mean success probability p 1+ is 0.05 and the motif has a p μ of 1.64 e-3 . When unique success probabilities are used a p v of 7.29 e-4 is calculated. In this case, the motif occurring k times or more is roughly twice as unlikely to occur as approximated by p μ . The bias introduced by p μ can therefore be considered conservative estimate.
Exact probability that any motif will be returned with a given p-value or less Previous work: Motif Binomial p-value Significance Heuristic SLiMFinder introduced Sig, the probability of any motif, with the same number of non-wildcard positions, occurring with that probability p or less by chance.
. Approximate probability (any motif will occur with a binomial p-value of p or less) where R is calculated as 20 l (x+1) l-1 , l is the number of non-wildcard positions in the motif and x is the maximum length of a wildcard region allowed. Effectively, this score assumes that all other possible motifs are equally likely to give rise to a score equal to, or lower than, the motif of interest. However, this is not always true. For example, more likely motifs may have, for all supports, a p-value greater than the observed motif's pvalue, hence such motifs may never actually occur with a probability as low as p. Sig will always score a given p as being as likely or more likely than it is, making Sig a conservative approximation.

More Accurate Significance
To calculate Sig', the true probability of 1 or more fixed position motifs (i.e. with the same number of non-wildcard positions) in a dataset occurring with the observed p or lower, it is necessary to calculate, for each possible motif p', the probability that the motif will occur with a score of p or less. An efficient algorithm for the calculation of Sig' is available in the Appendix.
. probability (any given motif will occur with a p-value less than or equal to the p of the observed motif) . probability (any motif will occur with a pvalue of p or less) , p motif ). The probability that a protein will contain one or more occurrences of a motif is p 1+ (LL) = 0.9987, p 1+ (WW) = 0.0394, p 1+ (WL) = p 1+ (LW) = 0.3143 (see Equation 1). Using the binomial and cumulative binomial we can calculate, using p 1+ as the mean success probability (all protein sequences are the same length and use the same global amino acid frequencies so the success probabilities for each protein for a given motif will be the same), the probability that each motif will occur k times (Binomial in Table 2) and the probability that each motif will occur k or more times (p μ in Table 2).  Searching the dataset for over-represented motifs, the support k for the motifs are LL = 4, WW = 1, WL = 1 and LW = 1. p μ of each motif in the motif space is p μ LL = 0.9948, p μ WW = 0.1485, p μ LW = p μ WL = 0.7789 (Equation 3), making WW the top ranking motif occurring in 1 protein with a probability, p, of 0.1485. From Equation 6, the significance, Sig, of p μ WW is 0.474. For the example in Table 2

Randomised test datasets
To test each of the 4 significance scoring schemes, 3 sets of 100 datasets containing 10, 30 and 60 proteins were picked at random from the UNIPROT database [22]. These datasets model random interaction networks that are likely to have little or no enrichment for any functional motifs. To allow for the calculation of Sig' v , a process that is currently computationally intractable using the calculations described here, a similar set of proteins was created to artificially force p 1+ values to be equal for each protein thereby testing the scoring schemes accuracy when no bias is present. A human protein was chosen at random from the UNIPROT database and the residues where shuffled to create datasets of 10, 30 and 60 proteins, therefore each protein is same length and has the same amino acid frequency and hence, for any given motif, has the same p 1+ for each protein.
Datasets were analysed using SLiMFinder with default settings with the exception of statistical calculations that were completed as described in this paper. The statistics described in this paper are implemented in freely downloadable software, as a modification to version 4.0 of the previously described SLiMFinder software [21], with additional command line switches introduced for the two calculations (sigprime = T/F and sigv = T/F).
The programme is available at http://bioinformatics. ucd.ie/shields/software/slimfinder. The version of the SLiMSuite software package available at time of submission is also included as supplemental information with this paper (Additional File 1).

Results and Discussion
Comparing the four statistical scoring schemes to expectation When tested on randomly created datasets, significance scores for fixed motifs should ideally be uniformly distributed [23] (i.e. if 100 datasets are analysed we would expect to see 1 motif with a significance of 0.01 or less, and 50% of the datasets should have a significance of 0.5 or greater). Using datasets of proteins randomly selected from UniProt, we tested to what extent each of the 4 scoring schemes rejected the null hypothesis, that the top ranking significance scores from these random datasets are uniformly distributed.
For all motif lengths and dataset sizes tested, Sig' v did not reject the null hypothesis, consistent with it being a true significance measure (Fig. 1a, top panel). The other three statistics all reject the null hypothesis to differing degrees, with Sig' showing less departure from expectation compared to Sig and Sig v . Sig rejects the null hypothesis for all the datasets. The statistics showed similar effects when considering the divergence of the observed p-values (root mean squared error) from the uniform distribution (Fig. 1a, side panel). The degree of rejection of the null hypothesis varies considerably according to the dataset size and the motif length (Fig  1a, main panel). A heatmap illustrates these effects (Fig  1b), showing that the degree of rejection of the null hypothesis for each of the three approximate statistics depends on both the dataset size, and on the length of the motif under consideration. Overall, the greatest improvement over the Sig statistic is provided by considering the exact probability that any motif may have a p-value of the observed significance or less (implemented in Sig' and in Sig' v ); while the correction for varying protein lengths and composition in a dataset has a lesser effect when calculating the p-value of the observed motif (implemented in Sig v and in Sig' v ).
How comparable are significance values of motifs of differing length, or motifs discovered in different sized datasets?
One of the major biases with using a binomial p-value, p μ , which SLiMFinder attempted to improve upon, was the incomparability of p-values for motifs of different lengths. p μ is strikingly different for random datasets, depending on whether it is the value for a 3mer or a 5mer (SF2 in Additional File 2). We were interested to discover to what extent a given motif significance score might be comparable for different discovered motifs, and to what extent the score is biased, depending on dataset size and motif length. For each scoring scheme, we compared the returned significance scores obtained with each motif size in each dataset size.
No biases are evident with the Sig' v scoring scheme where the hypothesis that the significance scores were sampled from the same distribution was only rejected (p < 0.05) for one of the 36 comparisons (Fig. 2b), consistent with expectations. The Sig statistic showed some biases, but somewhat surprisingly, the Sig v statistic appeared more biased than the Sig statistic (Fig. 2b): in  particular, the larger datasets seemed to be biased in a way that the smaller datasets are not (Fig 2a). However, It can be seen that the majority of Sig correlation between datasets relates to p-values in the range of 0.8 to 1.0 (Fig 2a) and is a related to the schemes conservative bias. The box-plots in Fig 2a indicate that Sig' is reasonably similarly distributed for all motif lengths and dataset sizes, but it is noticeable that the mean p-value varies more across datasets and motif sizes than Sig' v does. What is not clear from these random datasets is the extent to which these two statistics are biased for motifs that are significantly over-enriched (i.e. in ranking true positives, at the other end of the p-value scale).

Ambiguity
SLiMFinder considers subsets of pre-defined residue groupings to allow ambiguity in motif positions; this further adds to the complications of significance calculations and the probability calculations have not been extended to consider these. In general, the true significance calculations described in this paper, when used to calculate the significance of an ambiguous motif, will over-estimate the significance of the motif. In SF3 in Additional File 2, this is illustrated. For these examples, ambiguous motifs were defined as those with one or more ambiguous position or with at least one variable length spacer, using the default set of ambiguous groupings (KRH, DE, ILMV, FYWH, AGS, ST). Sig' v overestimates the significance, ignoring the greatly increased search space of potential ambiguous motifs, in addition to the search space of fixed-position motifs (SF3 in Additional File 2). We note that under the conditions tested, the Sig calculation is also conservative in all of the datasets considered here, even in the presence of ambiguity. Thus, for ambiguous motifs the Sig statistic may provide a more robust heuristic for very roughly approximating absolute significance. However, it must be noted that in general ambiguous motifs will still have a relatively higher Sig score than non-ambiguous motifs, so that a single model for ranking of ambiguous and non-ambiguous motifs is not ideal. Accordingly, search results within a dataset or across many datasets for ambiguous and non-ambiguous motifs should be considered as two separate rankings, given the current available models.

Calculation redundancy
The level of redundancy in the brute force calculation of Sig' compared to the use of non-redundant motif groupings (Appendix, Equations 9-13) is provide in Table 3. For 5mer motifs there is a 75-fold decrease in the number of calculations necessary for the calculation of the same value. Since the level of redundancy increases as the length of the motif increases, true Sig' calculations are possible for motifs which would been computationally infeasible using brute force calculations.

Computational efficiency
For each motif length considered, the number of calculations necessary to calculate Sig', with the algorithmic speedup described (Appendix, Equations 9-13), increases the number of calculations a+l-1 C l fold, compared to the Sig calculation, where a is the alphabet size and l is the motif length (without the speedup the increase is a l fold; see ST2 in Additional File 2). While Sig' is in these examples 3-50 times more computationally expensive than Sig, its calculation is tractable, increasing approximately linearly as dataset size increases (ST2 in Additional File 2). For very large datasets, the computational costs seem to be much less. Sig v increases the number of calculations ∑ n C j fold, for all values of j between k and n, where n is the number of clusters in the dataset and k is the support of the motif. The rapid increase of Sig v calculation runtime as the dataset size increases makes it impractical for larger datasets (ST2 in Additional File 2). Sig' v performs poorly compared to Sig, it increases the number of calculations a+l-1 C l *∑ n C j fold, for all values of j between k and n, for each motif length considered. The computational cost of Sig' v calculation is unsuitable for high throughput analysis in spite of its superior accuracy; it took over a day on a standard workstation for datasets with more than 10 proteins.

Comparison of the four scoring schemes with the ELM benchmarking dataset
The ELM benchmarking dataset is a widely used standard by which motif discovery methods can be tested, compared and trained. Although the dataset is Table 3 Redundancy of motif probabilities (see . biologically unrealistic for most discovery situations (each of the proteins in each dataset is guaranteed to contain one or more instances of the functional motif), the data allows direct comparison of the attributes of scoring schemes. Table 4 compares the Sig' and Sig scoring schemes on the 16 ELM benchmarking datasets (version -8 April 2009). Residues with an IUPRED [24] score of less then 0.2 were masked, as were under-conserved residues by relative local conservation masking using metazoan orthologues alignments retrieved by GOPHER [24].
Sig' scored the true positive motif as more significant than the Sig scoring scheme in all datasets with the exception of the highly significant TRG_ER_KDEL motif which was scored as more significant by Sig. For one of these, LIG_CYCLIN_1, this resulted in the returned motifs becoming significant (cut-off of 0.01). Two false positives (LIG_14-3-3_3 and LIG_NRBOX) were also returned as significant by Sig' which were not returned by Sig, however in both cases these motifs were off target motifs, true functional motifs involved in the regulation of the proteins in the datasets (See ST1 in Additional File 2). One other example of an off target motifs was significantly over represented by both scoring schemes, two motifs with strong nuclear localisation, LIG_RB and LIG_PCNA, returned arginine/lysine rich nuclear localisation motifs.

Conclusions
We have shown how more exact calculations of short non-ambiguous protein motif probabilities may be calculated efficiently, and the computational cost associated with these calculations. We propose the use of Sig' calculations as the scoring scheme for motif discovery with current typical computational resources, as it offers the best trade off between speed and accuracy. Although not as accurate, the Sig score will, in circumstances where a motif is obviously enriched, return a close approximation of the true score. Thus, a useful technique for high-throughput analyses is to recover datasets of interest using the quicker Sig statistical framework at low stringency (it is a conservative score) before re-running the data using the more accurate statistics. However, for more general use the Sig' statistical framework is both efficient and accurate enough to suffice.
Improving these statistical measures has a strong bearing on how to efficiently rank potential novel motifs discovered in searches of many datasets derived from across the interactome. Future work will need to focus on developing better statistical measures for ambiguous motifs (ambiguous in terms of both alternative amino acids at certain positions, and in terms of variable length spacers). The problem is non-trivial as the combinatorial explosion introduced by ambiguity and the complex protein groupings that allow a motif to reach a given  support makes a true ambiguous motif significance computationally intractable. Advances in methodology [21], dataset design [25] and motif enrichment [26], in association with the statistical framework described here, have put computational SLiM discovery in a position where it can complement experimental methods for the discovery of novel functional SLiMs. As our knowledge of protein-protein interaction, cellular localisation and keyword classification of proteins improves, these methods should aid in the discovery of the vast number of SLiMs remaining to be discovered by pinpointing particular residues of strong biological interest within a robust statistical framework.

Appendix
Improving the computational efficiency of the Sig' calculation The major bottleneck in the calculation of Sig' for a motif p-value, p μ , is the calculation of p' of all motifs with a cumulative p-value equal to or less than p. For longer motifs the complexity can make search times unfeasable, since the number of possible motifs grows exponentially with motif length. However, methods can take advantage of the inherent redundancy in the motif space due to the commutative nature of the motif occurrence probability calculations, p motif (Equation 1). Motifs which contain the same number of each residue from the alphabet A (where each residue can have between 0 and l occurrences in the motif) will have the same p motif . This means that the order of the residues in the motif has no effect on a motif's occurrence probability (e.g. RGG, GRG, GGR all have the same p motif ). Exploiting the redundancy by calculating motif groupings with the same value p motif yields large increases in the speed of calculations. The value p' will be the same for each motif in a grouping and therefore it is necessary to only calculate this value once for each group, p' group , rather than for each motif. Equation 9. Sig' calculation for non-redundant motif groupings A non-redundant Sig' can be calculated using Equation 9. G is the set of all distinct unordered groupings of residues for each partition in P. A partition describes the number of unique residues occurring in a grouping. For example, there are 3 possible partitions for motifs of length 3, q = {(3),(1,2),(1,1,1)} meaning a 3-mer can have 3 residues which are the same, two residues which are the same and another different residue or 3 residues which are different respectively. The set q corresponds to the set P = {(0,0,1),(1,1,0),(3,0,0)}. For example, RGG, GRG, GGR can be considered as a grouping R|GG, which has a partition P = (1,1,0) and q = (1,2). RGD, RDG, GRD, GDR, DGR and DRG can be considered as a grouping R|G|D, which has a partition P = (3,0,0) and q = (1,1,1) and RRR, can be considered as a grouping RRR, which has a partition P = (0,0,1) and q = (3). The number of possible partitions for a given motif length can be described as the cardinality of set P and follows the series {1,2,3,5,7,11,15,22} as l goes from 1 to 8. Each partition has Cg i distinct unordered groupings of residues and each of these groupings has Cm i possible ordered motif combinations. For example, a grouping with P = (0,0,1), 3 of the same amino acids eg. RRR, can only be chosen 1 way while P = (3,0,0), 3 different amino acids, can be chosen 6 ways.

Equation 12 Number of ordered motif combinations in a grouping
To calculate the p' group value for each non-redundant motif group, it is first necessary to define all such groupings. This can be achieved efficiently using the recursive algorithm described in the pseudocode found in SF1 in Additional File 2. The total number of distinct non-redundant motifs of length l and an alphabet of size a, R nr , is calculated as: