Mutual enrichment in ranked lists and the statistical assessment of position weight matrix motifs
 Limor Leibovich^{1} and
 Zohar Yakhini^{1, 2}Email author
https://doi.org/10.1186/17487188911
© Leibovich and Yakhini; licensee BioMed Central Ltd. 2014
 Received: 1 December 2013
 Accepted: 30 March 2014
 Published: 5 April 2014
Abstract
Background
Statistics in ranked lists is useful in analysing molecular biology measurement data, such as differential expression, resulting in ranked lists of genes, or ChIPSeq, which yields ranked lists of genomic sequences. State of the art methods study fixed motifs in ranked lists of sequences. More flexible models such as position weight matrix (PWM) motifs are more challenging in this context, partially because it is not clear how to avoid the use of arbitrary thresholds.
Results
To assess the enrichment of a PWM motif in a ranked list we use a second ranking on the same set of elements induced by the PWM. Possible orders of one ranked list relative to another can be modelled as permutations. Due to sample space complexity, it is difficult to accurately characterize tail distributions in the group of permutations. In this paper we develop tight upper bounds on tail distributions of the size of the intersection of the top parts of two uniformly and independently drawn permutations. We further demonstrate advantages of this approach using our software implementation, mmHGFinder, which is publicly available, to study PWM motifs in several datasets. In addition to validating known motifs, we found GCrich strings to be enriched amongst the promoter sequences of long noncoding RNAs that are specifically expressed in thyroid and prostate tissue samples and observed a statistical association with tissue specific CpG hypomethylation.
Conclusions
We develop tight bounds that can be calculated in polynomial time. We demonstrate utility of mutual enrichment in motif search and assess performance for synthetic and biological datasets. We suggest that thyroid and prostatespecific long noncoding RNAs are regulated by transcription factors that bind GCrich sequences, such as EGR1, SP1 and E2F3. We further suggest that this regulation is associated with DNA hypomethylation.
Keywords
 Statistical enrichment
 Position weight matrices
 Highthroughput sequencing data analysis
 Tissue specific methylation patterns
 lncRNA
Background
Modern data analysis often faces the task of extracting characteristic features from sets of elements singled out according to some measurement. In molecular biology, for example, an experiment may lead to measurement results pertaining to genes and then questions are asked about the properties of genes for which these were high or low. This is an example, of course, and the set of elements does not have to be genes. They can be genomic regions, proteins, structures, etc. A central technique for addressing the analysis of characteristic properties of sets of elements is statistical enrichment. More specifically – the experiment results are often representable as ranked lists of elements and we then seek enrichment of other properties of these elements at the top or bottom of the ranked list. GSEA [1], for example, is a tool that addresses characteristic features of genes that are found to be differentially expressed in a comparative transcriptomics study. GOrilla [2, 3] addresses GO terms enriched in ranked lists of genes where the ranking can be, for example, the result of processing differential expression data or of correlations computed between genomic DNA copy number and expression [4–6]. FATIGO [7] is also a tool that is useful in the context of analysing GO terms in ranked lists of genes. DRIMust [8–10] searches for sequence motifs that are enriched, in a statistically significant manner, in the top of a ranked list of sequences, which can be produced by techniques like ChIPSeq.
All the aforementioned tools utilize a statistical approach that is based on assessing enrichment of an input set in an input ranked list by quantifying the enrichment obtained at various cutoffs applied to the ranked list. It is often the case, however, that two quantitative properties need to be compared to each other. For example, when the elements are genes, we may have measured differential expression values, as well as measured ChIPSeq signals. We are therefore interested in assessing mutual enrichment in two ranked lists. Another example consists of one ranking according to differential expression and one according to prediction scores for miRNA targets. miTEA [11, 12] addresses this latter case by statistically assessing the enrichment of miRNA targets in a ranked list of genes (also see [13]). To address mutual enrichment in two ranked lists over the same set of N elements, miTEA [11] performs analysis on permutations. Mutual enrichment in the top of two ranked lists can be simplified, from a mathematical point of view, by arbitrarily setting the indices of one list to the identity permutation (1,2,…,N) and treating the other list as a permutation π = π(1), …, π(N) over these numbers. For the purpose of assessing the intersection of the top of the two ranked lists in a data driven manner, miTEA asks which prefix [1,…,n_{1}] is enriched in the first n_{2} elements of π, that is in the set π(1), …, π(n_{2}). The statistics introduced by miTEA is called mmHG (minimumminimumHyperGeometric). A slightly different variant of mmHG is described later in this section.
Statistics in the group of permutations S_{ N } is often difficult because the number of entities to be considered by any null model is N!. Direct exhaustive calculation of tail distributions over S_{ N } is therefore impractical for all but very small values of N. This difficulty is addressed by several heuristic techniques. Mapping into continuous spaces, such as in [14], has proven useful in certain cases but not for studying large deviations. In the case of enrichment statistics that focuses on the top of the permutation and seeks to assess extreme events, such as mmHG, we prefer to use bounds on tail probabilities. Tail probabilities are useful constructs when applied to analysing molecular biology measurement data as they enable statistical assessment of observed results.
In this work we derive tight bounds on the tail probabilities of mutual enrichment at the top of two random permutations uniformly drawn over S_{ N } and demonstrate the utility of this approach in the context of flexible motif discovery. Our bounds are computable in polynomial time and potentially add to the accuracy of reported position weight matrix (PWM) motifs for nucleic acid sequences.
Mutual enrichment in ranked lists – the mmHG statistics
The mmHG statistics [11] is a generalization of the mHG statistics [2, 15–17]. The mHG statistics quantifies the enrichment level of a set of elements at the top of a ranked list of elements of the same type, whereas the mmHG statistics assesses the level of mutual enrichment in two ranked lists over the same set of elements. While any parametric or nonparametric correlation statistics (e.g. Spearman’s correlation coefficient), that takes the same input, calculates the overall agreement between the two ranked lists, the mmHG statistic focuses only on agreement at the top of the two ranked lists. mmHG counts elements common to the top of both lists, without predefining what top is. Its intended output is the probability of observing an intersection at least as large in two randomly ranked lists (defined as the enrichment pvalue). In this section we describe the mmHG statistics and in later sections we suggest tight bounds for the pvalue. Our definition of the mmHG statistics varies slightly from that of Steinfeld et al.[11], which is used by miTEA.
Mutual enrichment in the top of two ranked lists can be simplified, from a mathematical point of view, by arbitrarily setting the indices of one list to the identity permutation (1,2,…,N) and treating the other list as a permutation. Details of this transform are given in the next section. We now define mmHG for the simple case of one permutation. Consider a permutation π = π(1), …, π(N) ∈ S_{ N }  the group of all permutations over the numbers 1,…,N. mmHG is a function that takes π and calculates two numbers 1 ≤ n_{1}, n_{2} ≤ N such that the observed intersection between the numbers 1,…,n_{1} and the first n_{2} elements of π – namely, π(1), …, π(n_{2})  is the most surprising in terms of the hypergeometric pvalue.
In the Results section we present significantly tighter bounds.
Position weight matrix motifs
Data produced by techniques such as ChIPSeq [18], ChIPexo [19], CLIP [20], PARCLIP [21] and others are readily representable as ranked lists of sequences, where the ranking is according to the measured binding affinity. Computational tools and approaches to motif discovery form part of the data analysis workflow that is used to extract knowledge and understanding from this type of studies. We are often interested in sequence motifs that are observed to be enriched in sequences where strong binding affinity is measured. A position weight matrix (PWM) is a commonly used representation of motifs in biological sequences [22–24]. This representation is more faithful to the underlying biology than representation by exact words, owing to the tendency of binding sites to be short and degenerate [25]. A PWM is a matrix of score values that gives a weighted match to any given substring of fixed length. It has one row for each symbol in the alphabet, and one column for each position in the pattern. Assuming an input sequence of length equal to the PWM width, we simply multiply the scores assigned to each letter in each of the positions in the input sequence to obtain the likelihood of the input string (alternatively, we can sum the logarithms of the probabilities). That is, the score assigned by a PWM to a substring S = S_{1}…S_{K} is defined as $\prod _{\mathit{j}=1}^{\mathit{K}}}{\mathit{p}}_{{\mathit{s}}_{\mathit{j}},\mathit{j}$, where j represents a position in the substring; s_{ j } is the symbol at position j in the substring; and p_{α,j} is the score in row α, column j of the matrix. In other words, a PWM score is the product of positionspecific scores for each symbol in the substring. This definition can be generalized to yield a score for a sequence S = S_{1}…S_{M} longer than the PWM by calculating $\mathit{ma}{\mathit{x}}_{1\le \mathit{i}\le \mathit{M}\mathit{K}+1}{\displaystyle \prod _{\mathit{j}=1}^{\mathit{K}}}{\mathit{p}}_{{\mathit{s}}_{\mathit{i}+\mathit{j}1},\mathit{j}}$. Alternatively, an enhanced model that takes into account multiple occurrences of the PWM in the sequence can be applied by summing over sufficiently strong occurrences of the PWM or by other more sophisticated approaches [26].
mmHG statistics for PWM motifs
Given a set of sequences that were tested in a high throughput experiment such as ChIPSeq [18], CLIP [20] and others, they can be ranked according to the measured binding affinities, yielding a ranked list L_{1}. Since usually we are interested in finding motifs amongst sequences having strong binding affinities, we actually search for motifs that are more prevalent at the top of this list. It is clear that any algorithm for denovo flexible motif search would need to evaluate candidate PWMs. Given a PWM which we want to assess, the sequences can also be ranked according to their PWM scores, yielding another ranked list L_{2}, different from L_{1}. A significant PWM motif would yield significant scores for sequences having strong binding affinities. Therefore, the question of PWM motif discovery from ranked experimental data can be formulated as quantifying the mutual enrichment level for the two ranked lists L_{1} and L_{2}. Given two ranked lists L_{1} and L_{2} over the universe of N sequences, they can be transformed into two respective permutations, π_{1} = (π_{1}(1), …, π_{1}(N)) and π_{2} = (π_{2}(1), …, π_{2}(N)). The relative permutation π, of π_{2} w.r.t. π_{1}, is defined by π(π_{1}(j)) = π_{2}(j), for every j = 1,…,N, or simply, using operations in the group S_{ N }: π = π_{2} ∙ π_{1}^{ 1}. Using the relative permutation π, we can represent the mutual enrichment of the top parts of L_{1} and L_{2} as mmHG score (π), defined above.
Results
Estimation of the mmHG pvalue – introducing first upper bound – B1
Since the number of permutations is huge, the process described above is very far from feasible. Therefore, we seek a computationally tractable upper bound, preferably tight.
A trivial upper bound is the Bonferroni corrected mmHG score defined by s · N^{2}. A more subtle upper bound was suggested by Steinfeld et al.[11] and is briefly described later as bound B3. In this work we introduce tighter bounds that are polynomially computable.
as we first choose b elements from the range [1,…,n_{1}] to appear at the first n_{2} entries of the permutation (there are $\left(\begin{array}{c}{\mathit{n}}_{1}\\ \mathit{b}\end{array}\right)$ possibilities). Then, we choose the positions amongst the first n_{2} entries that are occupied by these b elements, while considering all internal arrangements (for each choice of b elements there are $\left(\begin{array}{c}{\mathit{n}}_{2}\\ \mathit{b}\end{array}\right)\mathit{b}!$ possibilities). We next choose n_{2}b elements from the range [n_{1} + 1,…,N] to appear at the rest of the first n_{2} entries of the permutation (there are $\left(\begin{array}{c}\mathit{N}{\mathit{n}}_{1}\\ {\mathit{n}}_{2}\mathit{b}\end{array}\right)$ possibilities for that) and consider all possible (n_{2}  b) ! arrangements. Finally, we take into account all possible (N  n_{2}) ! arrangements of the remaining Nn_{2} entries of the permutation.
This upper bound is simple and requires O(N^{3}) HGT calculations. An HGT calculation takes O(N) time, assuming binomial coefficients can be calculated in constant time. Constant time computation can be achieved using Stirling’s approximation [27]: $\sqrt{2\mathrm{\pi n}}{\left(\frac{\mathrm{n}}{\mathrm{e}}\right)}^{\mathrm{n}}{\mathrm{e}}^{\frac{1}{12\mathrm{n}+1}}\le \mathrm{n}!\le \sqrt{2\mathrm{\pi n}}{\left(\frac{\mathrm{n}}{\mathrm{e}}\right)}^{\mathrm{n}}{\mathrm{e}}^{\frac{1}{12\mathrm{n}}}$, which is tight for large factorials.
A refined upper bound for the pvalue – B2
The upper bound introduced in the previous section counts the number of permutations for which the value HGT(N, n_{1}, n_{2}, b) is calculated when taking the nonpractical exhaustive approach that enumerates over all N! permutations. Ideally, we wish to count the number of permutations for which the value HGT(N, n_{1}, n_{2}, b) is also their mmHG score, as a permutation may correspond with many HGT values that are better than s, so it can be counted more than once. This explains why the formula introduced earlier is an upper bound and not an exact pvalue. A second observation that follows is that the smaller the mmHG score s, the tighter the bound, because a permutation will have fewer combinations (N, n_{1}, n_{2}, b) having HGT values better than s.
Therefore, if we can reduce the extent of multiple counting of the same permutation, we will get a tighter bound. We do this by looking one step backwards. If, for example, HGT(N, n_{1}, n_{2}, b) ≤ s, we can exclude from the counting permutations that contain b elements from the range [1,…,n_{1}1] at their first n_{2} entries because they are already taken into account in Λ(N, n_{1}  1, n_{2}, b) (because necessarily HGT(N, n_{1}  1, n_{2}, b) ≤ s, as we will later explain).
The properties of the hypergeometric distribution imply that given a tuple (N, n_{1}, n_{2}, b), the permutations in ψ_{1}, ψ_{2}, ψ_{4} can be disregarded from the current counting iteration. To explain why, we will demonstrate the argument on ψ_{1}. The permutations in ψ_{1} contain b elements from the range [1,…,n_{1}1] at their first n_{2} entries. Recall that we also assume that HGT(N, n_{1}, n_{2}, b) ≤ s. Therefore HGT(N, n_{1}  1, n_{2}, b) ≤ s also holds, as the same intersection is observed for even a smaller set. Thus, the permutations in ψ_{1} have already been counted as having HGT value better than s when handling the triplet n_{1}1, n_{2} and b, and can be disregarded for the combination of n_{1}, n_{2} and b. Similar arguments hold for ψ_{2} and ψ_{4}.
Since Λ* is recursive, we need to define a base case. Recall that given N, n_{1} and n_{2}, b can be any integer in the range [max(0, n_{2}  N + n_{1}), min(n_{1}, n_{2})], hence determining a base case for n_{1} and n_{2} is sufficient (N is known). The base case here is that when n_{1} ≤ 1 or n_{2} ≤ 1, Λ*(N, n_{1}, n_{2}, b) is defined the same as Λ(N, n_{1}, n_{2}, b).
This upper bound uses more delicate counting than the bound B1 introduced in the previous section. In the following sections we assess the tightness of this bound. In later sections we demonstrate an application for PWM motif search.
Comparison to a different mmHG variant – B3
Assessment of tightness
In order to assess the quality of our bound B2, we compared it to the pvalue, which can be calculated exactly for small values of N (that is, in cases where N! is not too large) and empirically for larger values of N (by randomly sampling permutations). Evidently, our bound B2 was significantly tighter than the Bonferroni bound for N = 10 (Figure 1A) and N = 20 (Figure 1B). We also observed that the smaller the mmHG scores – the tighter the bound, consistent with lesser overcounting for smaller scores as explained in previous sections. Furthermore, our refined bound B2 is tighter than the bound B1 (Figure 1C), and the latter is significantly better than the Bonferroni bound. Both bounds B1 and B2 are derived by enumerating HGT scores rather than enumerating permutations in S_{ N }. The refinement of this approach produced by reducing the extent of multiple counting of permutation further improves the upper bound. In addition, the bound B2 was almost always observed to be tighter than the bound B3 (Figure 1D).
An upper bound which balances between tightness and computational cost – B4
Performance of various bounds
Protein  N  mmHG score  Bound B2  Bound B3  Bound B4  Bonferroni bound 

TNWMNG  500  2.17e18  2.75e14  7.5e14  6.28e14  5.43e13 
0.274 min  0.0028 min  0.079 min  
CTNNNAT  500  2.86e27  1.32e28  3.66e28  2.37e28  2.86e27 
0.155 min  0.0029 min  0.059 min  
MMMMMMMM  500  1.08e43  1.07e39  3.47e39  1.69e39  2.71e38 
0.104 min  0.003 min  0.048 min  
REB1  4000  1.66e137  9.18e133  1.19e131  1.54e132  1.67e131 
17.25 min  0.04 min  2.753 min  
CBF1  4000  1.95e80  9.15e76  4.62e75  1.84e75  1.96e74 
26.05 min  0.03 min  3.409 min  
UME6  4000  5.42e88  2.62e83  3.04e82  5.11e83  5.43e82 
23.81 min  0.03 min  3.374 min  
TYE7  4000  1.62e43  5.63e39  2.83e38  1.39e38  1.62e37 
34.25 min  0.02 min  4.05 min  
GCN4  4000  2.04e50  7.66e46  4.62e45  1.80e45  2.04e44 
35.43 min  0.03 min  3.95 min  
Puf5  4795  7.91e85  3.38e80  5.60e79  6.95e80  7.93e79 
31.51 min  0.027 min  4.51 min  
Pub1  4251  1.49e84  6.86e80  1.33e78  1.37e79  1.5e78 
27.74 min  0.033 min  3.81 min  
Pab1  4142  2.46e11  3.57e7  5.17e7  1.37e6  2.46e5 
48.46 min  0.007 min  5.41 min  
Khd1  4773  2.74e20  5.09e16  1.46e14  1.73e15  2.74e14 
47.58 min  0.015 min  5.84 min  
Nab2  4101  2.09e11  3.08e7  1.48e5  1.18e6  2.09e5 
48.7 min  0.016 min  5.34 min  
Vts1  1787  1.44e10  4.74e6  1.33e5  1.4e5  1.45e4 
21.94 min  0.003 min  2.07 min  
Pin4  4261  8.16e14  1.32e9  8.08e9  4.83e9  8.18e8 
49.38 min  0.011 min  5.48 min  
Nrd1  3947  5.72e12  9.09e8  5.71e6  3.36e7  5.74e6 
47.67 min  0.014 min  5.11 min  
Yll032c  2286  1.06e9  2.62e5  1.61e4  8.3e5  0.001 
35.58 min  0.003 min  2.77 min 
Application in PWM motif search
In this section we discuss mmHG as a framework for assessing the significance of PWM motifs in ranked lists. Given a ranked list of sequences and a PWM motif, by using the mmHG statistics and the bounds introduced earlier, we can assign a pvalue to represent the significance of that PWM being enriched at the top of the list. To apply this approach for denovo motif search, one needs to theoretically consider all possible PWMs. However, the search space  when considering position weight matrix motifs – is huge. Assuming the probabilities in the matrix are multiples of 0.1 and the alphabet is of size 4, there are 286^{ k } possible candidate PWMs of length k (since each column must sum to 1, the number of combinations in each column of the matrix is equal to the number of integer solutions for the equation X_{1} + X_{2} + X_{3} + X_{4} = 10, which is $\left(\begin{array}{c}13\\ 10\end{array}\right)$). Our approach to navigating in this search space is to narrow the search using the IUPAC alphabet, which considers all possible combinations of letters in the alphabet, and then represent the motif as a PWM based on its actual occurrences in the list. This heuristic approach, called mmHGFinder, takes as input a ranked list of DNA or RNA sequences and returns significant motifs in PWM format. In cases where sequence ranking is not relevant or not available, it allows the use of positive and negative sets of sequences, searching for enriched motifs in the positive set using the negative set as the background.
We next describe the methodology implemented in mmHGFinder. The input consists of a ranked list of sequences (or, alternatively, two sets of sequences representing target and background), as well as the motif width, given as a range [k_{1}, k_{2}].
 1.
Build a generalized suffix tree for the input sequences.
 2.
For each k=k _{1},…,k _{2}

Traverse the tree to find all kmers

Sort the kmers according to their enrichment at the top of the list (this is done using the mHG statistics), as explained in Leibovich et al.[8].

Take the most significant fifty kmers, to be used as starting points for the next step. This set of candidates is chosen such that the members are quite different. Note that this is a heuristic approach and the number 50 is somewhat arbitrary, chosen to succeed in catching the best performing PWMs without heavily paying in complexity.

For each starting point, we iteratively replace one position in the kmer by considering all possible IUPAC replacements and taking the one that improves the enrichment the most. We repeat this process for all positions several times, and eventually we get a motif in the IUPAC alphabet. We note that given an IUPAC pattern P, the occurrences of P in the list are extracted efficiently by traversing the paths in the suffix tree that agree with P.

Each IUPAC word is then expanded through a heuristic approach which is based on the Hamming neighbors of that word. Hamming neighbors are added as long as the new addition improves the enrichment pvalue of the set of words, and as long as the overall similarity between the members in the set does not decrease below a similarity threshold. Since the neighbors are defined as exact words, they usually help in finetuning the correct weights of each letter in each position of the resulting PWM. Finally, the expanded motif is converted to a PWM.
 3.
The PWMs found in the previous step are assessed using the mmHG statistics and the best PWMs are returned as output, together with their pvalue. The score assigned by a PWM to a string S is the maximal score obtained for a substring of S. To obtain the likelihood of a substring of length k (where k is the PWM width), we simply multiply the scores assigned to each letter in each of the positions in that substring.
We provide an efficient implementation of the algorithm described above as publicly available software. Our application takes as input a ranked list of sequences and returns significant PWM motifs. It is compatible with all operating systems and can be freely downloaded from http://bioinfo.cs.technion.ac.il/people/zohar/mmHGFindercode/.
PWM motif search in longnoncoding RNA sequences
 1.
We searched for denovo PWM motifs in the promoter sequences of the tissuespecific lncRNAs using mmHGFinder (introduced in the previous section). Promoter sequences were defined as 1000 bp upstream the transcription start site.
 2.
We scanned the promoter sequences of the tissuespecific lncRNAs with PWMs corresponding to known transcription factors, downloaded from the JASPAR database [35].
 3.
Independently of sequence, we calculated the statistical enrichment of measured transcription factor binding events within our lists of loci. Transcription factor binding events within lncRNAs were downloaded from ChIPBase database [36], which aggregates highthroughput sequencing data taken from hundreds of ChIPSeq experiments.
CpG hypomethylation in tissuespecific lncRNA promoters
Tissue  Mutual enrichment between GC content and tissuespecificity  Mutual enrichment between hypomethylation and tissuespecificity  

Normal/subnormal cells  Cancer cells  
Thyroid  3.89e31  No methylation data  No methylation data 
Prostate  5.76e22  4.16e11 (PrEC)  0.002 (LNCaP) 
Adrenal  5.46e20  No methylation data  No methylation data 
Brain  1.57e14  1.21e8 (NHA)  5.55e5 (U87) 
Ovary  8.80e12  No methylation data  0.0085 (ovcar3) 
Lymph node  3.64e6  No methylation data  No methylation data 
Adipose  9.25e6  No methylation data  No methylation data 
Foreskin  2.25e5  0.72 (BJ)  No methylation data 
Breast  4.40e5  5.08e5 (HMEC)  8.45e5 (MCF7) 
2.0e10 (MCF10A)  0.0065 (T47D)  
Kidney  6.34e5  1.56e5 (HEK293)  No methylation data 
White blood cell  3.78e4  0.6 (GM12878)  0.21 (Jurkat) 
Placenta  0.011  No methylation data  No methylation data 
Colon  0.012  No methylation data  1.0 (Caco2) 
Skeletal muscle  0.04  0.34 (SKMC)  No methylation data 
Lung  0.33  No methylation data  No methylation data 
Heart  1.0  1.0 (HCM)  No methylation data 
1.0 (HCF)  
Liver  1.0  1.0 (Hepatocytes)  1.0 (HepG2) 
Testes  1.0  No methylation data  1.0 (NT2D1) 
Lung fibroblasts  1.0  1.0 (IMR90)  No methylation data 
1.0 (AGO4450) 
Furthermore, by intersecting the results of the second and the third tests together, we identified transcription factors that may regulate lncRNAs, mainly in thyroid and prostate. This set includes NRF1, E2F1, E2F3, E2F4, E2F6, EGR1, SP1, SP2 and ZBTB33. Moreover, the consensus recognition sites of EGR1, SP1 and E2F3 were found to be similar to the motifs returned by mmHGFinder in thyroid, prostate and other tissues (Figure 3; the comparison was done using the motif discovery tool TOMTOM [37]). The full output of the second and the third tests are summarized in Additional file 2.
As GCrich motifs may be associated with CpG methylation, and due to the possible binding of SP1 which has been suggested to protect CpG islands from de novo methylation [17, 38], we further tested the association between hypomethylation and tissue specificity. For that, we downloaded genome wide (450 K) CpG Methylation data from UCSC Table Browser [39] (ENCODE/HAIB). We intersected lncRNA promoter regions with CpG methylation data, and continued only with the 1099 loci that were covered by the methylation experiment. For them, we calculated the mutual enrichment between hypomethylation and tissuespecificity (the results are summarized in Table 2). Thyroid cells were not covered in this experiment, however two cell lines corresponding to prostate were tested (normal prostate epithelial cell line and cancerous prostate endodermal cell line). We observed that prostatespecific lncRNA promoters were less methylated than nonprostatespecific lncRNAs, and this was much stronger in normal cells than in cancer (mmHG pvalue ≤ 4.16e11 in normal prostate cells, and 0.002 in prostate adenocarcinoma cells). We observed strong mutual enrichment between CpG hypomethylation and tissuespecificity also in brain, ovary, breast and kidney. That is, significant mutual enrichment values were found for tissues where tissuespecific lncRNAs had GCrich promoter sequences, but these values were not significant for tissues that did not show such GCbias (heart, liver, testes, and lung fibroblasts). Additionally, in most cases the significance in normal cells was stronger than in cancer, which may be related to changes in methylation patterns acquired during carcinogenesis [40, 41].
Discussion
The assessment of mutual enrichment in ranked lists is often required to support the analysis of biological measurement data, such as in the case of identifying sequence motifs that are involved in regulation processes. Relative ranking can be represented by using permutations over the measured elements. Therefore – the statistical assessment of mutual enrichment can be modelled by characterizing properties of random permutations. Due to the size of the measure space, statistics over S_{ N }, the group of permutations over N elements, is difficult to perform and implement. Mutual enrichment is more informative from the point of view of practical biological science than simple correlation measures, as it focuses on the top of the lists and not on the overall agreement, which may be weak even in cases where extremities agree. In this work we derive polynomially computable bounds for the associated tail distribution of mutual enrichment in ranked lists. Namely – we provide methods for computing an upper bound on the pvalue of mutual enrichment at the top of two permutations uniformly and independently drawn over S_{ N }. Naïve approaches to computing such bounds include variants of the Bonferroni approach. These do not provide tight bounds and may lead to mislabeling results as nonsignificant. For several representative datasets, we note that our bound improves the Bonferroni derived pvalue estimates by a factor of almost 40, on average. Nevertheless, these improvements become relevant only for high pvalues  for which significant scores should be treated with care anyway. We therefore note that the Bonferroni correction is applicable in many cases, as demonstrated in Table 1. Using our bounds is highly beneficial in borderline cases but is also important in cases where an accurate estimate of the pvalue is desired, even if nuances do not affect the final biological research conclusions.
We use our statistical/algorithmic framework to support PWM motif searches and demonstrate the application to biological data. We identify motifs that characterize tissue specificity of lncRNA in thyroid and in prostate. Specifically, we find the EGR1 binding motif to be enriched in the promoter regions of lncRNAs which are thyroidspecific. EGR1 was observed to be highly expressed in thyroid (Additional file 1, taken from [36]), consistent with our stronger motif findings. Similarly, EGR1 is highly expressed in adipose tissue and its transcription factor binding sites are enriched in lncRNAs specific to this tissue. We do not have methylation data for the latter two tissues types. However – we do observe the promoters of lncRNAs that are specific to breast to have enriched occurrences of motifs that are similar to EGR1 transcription factor binding sites (pvalue of similarity according to TOMTOM = 3.52∙10^{5}). EGR1 is also highly expressed in breast. Finally, the promoters of lncRNAs that are specific to breast are less methylated in breast (MCF10A and HMEC cells) than other promoters. This suggests the role of EGR1 in driving tissue differentiation by transcribing tissuespecific lncRNAs and by protecting the associated promoters from methylation. EGR1 has been previously shown to recognize GCrich consensus sequences located in CpG island promoters of active genes [42]. Generally, we observed that tissuespecific lncRNA promoters tend to be less methylated than those of nontissuespecific lncRNAs in prostate, brain, ovary, breast and kidney, which may be associated with the GCrich patterns enriched among their tissuespecific lncRNA promoter sequences.
Thresholdfree alternatives to mmHG include the work of McLeay and Bailey, in which a linear regression method is applied [43]. It was shown to achieve high accuracy on a benchmark comprising 237 ChIPchip datasets, which was higher than all other data driven methods tested, and specifically higher than Spearman’s rank correlation. We note that applying linear regression or Spearman correlation to PWM motif search in ranked lists requires that for significant motifs we observe an overall agreement between the biological measurement and the PWM score. Nevertheless, the standard PWM formulation fails to predict binding affinity when the latter decreases to the point of nonspecific binding [44]. In other words, the overall agreement between the PWM score and the binding affinity may be relatively weak. High correlation between the PWM score and the binding affinity needs to hold, in effect, only for sequences demonstrating highbinding affinity with respect to the protein of interest (that is, for sequences that are located at the top of the list) [45]. This weaker relationship is naturally addressed by the mmHG statistics. A combination of mmHG and a linear model, such as suggested in [43], applied to strong binders (top of the list), may yield an even more faithful and informative model.
Future research directions include more extensive application to biological data and the development of tighter and more efficient bounds. Our results show promise in enabling efficient and userfriendly PWM motif search in ranked lists. The software is freely available at http://bioinfo.cs.technion.ac.il/people/zohar/mmHGFindercode/. Finally, the full characterization of the distribution of mmHG as a random variable over S_{ N } remains an open question.
Conclusions
In this work we developed tight bounds on the tail distribution of mutual enrichment in ranked lists. Our bounds are computable in polynomial time and potentially add to the accuracy of reported results. We demonstrated the utility of mutual enrichment in motif search – specifically, when searching position weight matrix motifs in ranked lists, where the ranking can be according to binding affinity or according to any other biological measurement. Additionally, we used mutual enrichment to study tissuespecific long noncoding RNA regulation, and suggest that tissuespecific lncRNAs are regulated through GCrich elements located on their promoters, in several tissue types. We hypothesize that these GCrich patterns are associated with DNA hypomethylation.
Declarations
Acknowledgements
We thank Israel Steinfeld for critical and inspiring discussions. We thank Roy Navon for technical help with the software download service. We also thank the WABI anonymous reviewers for their useful comments. LL was partially supported by Israel Ministry of Science and Technology and by ISEF Fellowship.
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 credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
 Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP: Gene set enrichment analysis: a knowledgebased approach for interpreting genomewide expression profiles. Proc Natl Acad Sci U S A. 2005, 102: 1554515550. 10.1073/pnas.0506580102.View ArticlePubMedPubMed CentralGoogle Scholar
 Eden E, Navon R, Steinfeld I, Lipson D, Yakhini Z: GOrilla: a tool for discovery and visualization of enriched GO terms in ranked gene lists. BMC Bioinformatics. 2009, 10: 4810.1186/147121051048.View ArticlePubMedPubMed CentralGoogle Scholar
 GOrilla Webserver. [http://cblgorilla.cs.technion.ac.il/]
 RagleAure M, Steinfeld I, Baumbusch LO, Liestøl K, Lipson D, Nyberg S, Naume B, Sahlberg KK, Kristensen VN, BørresenDale AL, Lingjærde OC, Yakhini Z: Identifying intrans process associated genes in breast cancer by integrated analysis of copy number and expression data. PLoS ONE. 2013, 8: e5301410.1371/journal.pone.0053014.View ArticleGoogle Scholar
 Akavia UD, Litvin O, Kim J, SanchezGarcia F, Kotliar D, Causton HC, Pochanard P, Mozes E, Garraway LA, Pe’er D: An integrated approach to uncover drivers of cancer. Cell. 2010, 143: 10051017. 10.1016/j.cell.2010.11.013.View ArticlePubMedPubMed CentralGoogle Scholar
 Dehan E, BenDor A, Liao W, Lipson D, Frimer H, Rienstein S, Simansky D, Krupsky M, Yaron P, Friedman E, Rechavi G, Perlman M, AviramGoldring A, Izraeli S, Bittner M, Yakhini Z, Kaminski N: Chromosomal aberrations and gene expression profiles in nonsmall cell lung cancer. Lung Cancer. 2007, 56: 175184. 10.1016/j.lungcan.2006.12.010.View ArticlePubMedGoogle Scholar
 AlShahrour F, DíazUriarte R, Dopazo J: FatiGO: a web tool for finding significant associations of gene ontology terms with groups of genes. Bioinformatics. 2004, 20: 578580. 10.1093/bioinformatics/btg455.View ArticlePubMedGoogle Scholar
 Leibovich L, Yakhini Z: Efficient motif search in ranked lists and applications to variable gap motifs. Nucleic Acids Res. 2012, 40: 58325847. 10.1093/nar/gks206.View ArticlePubMedPubMed CentralGoogle Scholar
 Leibovich L, Paz I, Yakhini Z, MandelGutfreund Y: DRIMust: a web server for discovering rank imbalanced motifs using suffix trees. Nucleic Acids Res. 2013, 41: W174W179. 10.1093/nar/gkt407.View ArticlePubMedPubMed CentralGoogle Scholar
 DRIMust Webserver. [http://drimust.technion.ac.il/]
 Steinfeld I, Navon R, Ach R, Yakhini Z: miRNA target enrichment analysis reveals directly active miRNAs in health and disease. Nucleic Acids Res. 2013, 41: e45e45. 10.1093/nar/gks1142.View ArticlePubMedPubMed CentralGoogle Scholar
 miTEA Webserver. [http://cblgorilla.cs.technion.ac.il/miTEA/]
 Enerly E, Steinfeld I, Kleivi K, Leivonen SK, RagleAure M, Russnes HG, Rønneberg JA, Johnsen H, Navon R, Rødland E, Mäkelä R, Naume B, Perälä M, Kallioniemi O, Kristensen VN, Yakhini Z, BørresenDale AL: miRNAmRNA integrated analysis reveals roles for miRNAs in primary breast tumors. PLoS ONE. 2011, 6: e1691510.1371/journal.pone.0016915.View ArticlePubMedPubMed CentralGoogle Scholar
 Plis SM, Weisend MP, Damaraju E, Eichele T, Mayer A, Clark VP, Lane T, Calhoun VD: Effective connectivity analysis of fMRI and MEG data collected under identical paradigms. Comput Biol Med. 2011, 41: 11561165. 10.1016/j.compbiomed.2011.04.011.View ArticlePubMedPubMed CentralGoogle Scholar
 Eden E, Lipson D, Yogev S, Yakhini Z: Discovering motifs in ranked lists of DNA sequences. PLoS Comput Biol. 2007, 3: e3910.1371/journal.pcbi.0030039.View ArticlePubMedPubMed CentralGoogle Scholar
 Steinfeld I, Navon R, Ardigò D, Zavaroni I, Yakhini Z: Clinically driven semisupervised class discovery in gene expression data. Bioinformatics. 2008, 24: i90i97. 10.1093/bioinformatics/btn279.View ArticlePubMedGoogle Scholar
 Straussman R, Nejman D, Roberts D, Steinfeld I, Blum B, Benvenisty N, Simon I, Yakhini Z, Cedar H: Developmental programming of CpG island methylation profiles in the human genome. Nat Struct Mol Biol. 2009, 16: 564571. 10.1038/nsmb.1594.View ArticlePubMedGoogle Scholar
 Lee BK, Bhinge AA, Iyer VR: Wideranging functions of E2F4 in transcriptional activation and repression revealed by genomewide analysis. Nucleic Acids Res. 2011, 39: 35583573. 10.1093/nar/gkq1313.View ArticlePubMedPubMed CentralGoogle Scholar
 Rhee Ho S, Pugh BF: Comprehensive genomewide proteinDNA interactions detected at singlenucleotide resolution. Cell. 2011, 147: 14081419. 10.1016/j.cell.2011.11.013.View ArticlePubMedPubMed CentralGoogle Scholar
 Lebedeva S, Jens M, Theil K, Schwanhäusser B, Selbach M, Landthaler M, Rajewsky N: Transcriptomewide analysis of regulatory interactions of the RNAbinding protein HuR. Molecular Cell. 2011, 43: 340352. 10.1016/j.molcel.2011.06.008.View ArticlePubMedGoogle Scholar
 Hafner M, Landthaler M, Burger L, Khorshid M, Hausser J, Berninger P, Rothballer A, Ascano M, Jungkamp AC, Munschauer M, Ulrich A, Wardle GS, Dewell S, Zavolan M, Tuschl T: Transcriptomewide identification of RNAbinding protein and MicroRNA target sites by PARCLIP. Cell. 2010, 141: 129141. 10.1016/j.cell.2010.03.009.View ArticlePubMedPubMed CentralGoogle Scholar
 Staden R: Computer methods to locate signals in nucleic acid sequences. Nucleic Acids Res. 1984, 12: 505519. 10.1093/nar/12.1Part2.505.View ArticlePubMedPubMed CentralGoogle Scholar
 Stormo GD, Schneider TD, Gold L: Quantitative analysis of the relationship between nucleotide sequence and functional activity. Nucleic Acids Res. 1986, 14: 66616679. 10.1093/nar/14.16.6661.View ArticlePubMedPubMed CentralGoogle Scholar
 Hertz GZ, Stormo GD: Identifying DNA and protein patterns with statistically significant alignments of multiple sequences. Bioinformatics. 1999, 15: 563577. 10.1093/bioinformatics/15.7.563.View ArticlePubMedGoogle Scholar
 Tompa M, Li N, Bailey TL, Church GM, De Moor B, Eskin E, Favorov AV, Frith MC, Fu Y, Kent WJ, Makeev VJ, Mironov AA, Noble WS, Pavesi G, Pesole G, Regnier M, Simonis N, Sinha S, Thijs G, van Helden J, Vandenbogaert M, Weng Z, Workman C, Ye C, Zhu Z: Assessing computational tools for the discovery of transcription factor binding sites. Nat Biotech. 2005, 23: 137144. 10.1038/nbt1053.View ArticleGoogle Scholar
 Sinha S: On counting position weight matrix matches in a sequence, with application to discriminative motif finding. Bioinformatics. 2006, 22: e454e463. 10.1093/bioinformatics/btl227.View ArticlePubMedGoogle Scholar
 Abramowitz M, Stegun IA: Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. 1964, New York: Dover Publications, Inc.Google Scholar
 Bailey TL, Elkan C: Unsupervised learning of multiple motifs in biopolymers using expectation maximization. Mach Learn. 1995, 21: 5180.Google Scholar
 Bailey TL: DREME: motif discovery in transcription factor ChIPseq data. Bioinformatics. 2011, 27: 16531659. 10.1093/bioinformatics/btr261.View ArticlePubMedPubMed CentralGoogle Scholar
 Luehr S, Hartmann H, Söding J: The XXmotif web server for eXhaustive, weight matriXbased motif discovery in nucleotide sequences. Nucleic Acids Res. 2012, 40: W104W109. 10.1093/nar/gks602.View ArticlePubMedPubMed CentralGoogle Scholar
 Smeenk L, van Heeringen SJ, Koeppel M, van Driel MA, Bartels SJJ, Akkers RC, Denissov S, Stunnenberg HG, Lohrum M: Characterization of genomewide p53binding sites upon stress response. Nucleic Acids Res. 2008, 36: 36393654. 10.1093/nar/gkn232.View ArticlePubMedPubMed CentralGoogle Scholar
 Harbison CT, Gordon DB, Lee TI, Rinaldi NJ, Macisaac KD, Danford TW, Hannett NM, Tagne JB, Reynolds DB, Yoo J, Jennings EG, Zeitlinger J, Pokholok DK, Kellis M, Rolfe PA, Takusagawa KT, Lander ES, Gifford DK, Fraenkel E, Young RA: Transcriptional regulatory code of a eukaryotic genome. Nature. 2004, 431: 99104. 10.1038/nature02800.View ArticlePubMedPubMed CentralGoogle Scholar
 Hogan DJ, Riordan DP, Gerber AP, Herschlag D, Brown PO: Diverse RNAbinding proteins interact with functionally related sets of RNAs. Suggesting an extensive regulatory system. PLoS Biol. 2008, 6: e25510.1371/journal.pbio.0060255.View ArticlePubMedPubMed CentralGoogle Scholar
 Cabili MN, Trapnell C, Goff L, Koziol M, TazonVega B, Regev A, Rinn JL: Integrative annotation of human large intergenic noncoding RNAs reveals global properties and specific subclasses. Genes Dev. 2011, 25: 19151927. 10.1101/gad.17446611.View ArticlePubMedPubMed CentralGoogle Scholar
 Mathelier A, Zhao X, Zhang AW, Parcy F, WorsleyHunt R, Arenillas DJ, Buchman S, Chen Cy, Chou A, Ienasescu H, Lim J, Shyr C, Tan G, Zhou M, Lenhard B, Sandelin A, Wasserman WW: JASPAR 2014: an extensively expanded and updated openaccess database of transcription factor binding profiles. Nucleic Acids Res. 2013, 42: D142D147.View ArticlePubMedPubMed CentralGoogle Scholar
 Yang JH, Li JH, Jiang S, Zhou H, Qu LH: ChIPBase: a database for decoding the transcriptional regulation of long noncoding RNA and microRNA genes from ChIPSeq data. Nucleic Acids Res. 2013, 41: D177D187. 10.1093/nar/gks1060.View ArticlePubMedPubMed CentralGoogle Scholar
 Gupta S, Stamatoyannopoulos J, Bailey T, Noble W: Quantifying similarity between motifs. Genome Biol. 2007, 8: R2410.1186/gb200782r24.View ArticlePubMedPubMed CentralGoogle Scholar
 Brandeis M, Frank D, Keshet I, Siegfried Z, Mendelsohn M, Names A, Temper V, Razin A, Cedar H: Sp1 elements protect a CpG island from de novo methylation. Nature. 1994, 371: 435438. 10.1038/371435a0.View ArticlePubMedGoogle Scholar
 UCSC Table Browser. [http://genome.ucsc.edu/cgibin/hgTables?command=start]
 Bert SA, Robinson MD, Strbenac D, Statham AL, Song JZ, Hulf T, Sutherland RL, Coolen MW, Stirzaker C, Clark SJ: Regional activation of the cancer genome by longrange epigenetic remodeling. Cancer Cell. 2013, 23: 922. 10.1016/j.ccr.2012.11.006.View ArticlePubMedGoogle Scholar
 Nejman D, Straussman R, Steinfeld I, Ruvolo M, Roberts D, Yakhini Z, Cedar H: Molecular rules governing de novo methylation in cancer. Cancer Res. 2014, 74: 14751483. 10.1158/00085472.CAN133042.View ArticlePubMedGoogle Scholar
 Kubosaki A, Tomaru Y, Tagami M, Arner E, Miura H, Suzuki T, Suzuki M, Suzuki H, Hayashizaki Y: Genomewide investigation of in vivo EGR1 binding sites in monocytic differentiation. Genome Biol. 2009, 10: R4110.1186/gb2009104r41.View ArticlePubMedPubMed CentralGoogle Scholar
 McLeay R, Bailey T: Motif enrichment analysis: a unified framework and an evaluation on ChIP data. BMC Bioinformatics. 2010, 11: 16510.1186/1471210511165.View ArticlePubMedPubMed CentralGoogle Scholar
 Frank DE, Saecker RM, Bond JP, Capp MW, Tsodikov OV, Melcher SE, Levandoski MM, Record MT: Thermodynamics of the interactions of lac repressor with variants of the symmetric lac operator: effects of converting a consensus site to a nonspecific site. J Mol Biol. 1997, 267: 11861206. 10.1006/jmbi.1997.0920.View ArticlePubMedGoogle Scholar
 Benos PV, Lapedes AS, Stormo GD: Is there a code for proteinDNA recognition? Probab(ilistical)ly. Bioessays. 2002, 24: 466475. 10.1002/bies.10073.View ArticlePubMedGoogle Scholar