Mutual enrichment in ranked lists and the statistical assessment of position weight matrix motifs
© Leibovich and Yakhini; licensee BioMed Central Ltd. 2014
- Received: 1 December 2013
- Accepted: 30 March 2014
- Published: 5 April 2014
Statistics in ranked lists is useful in analysing molecular biology measurement data, such as differential expression, resulting in ranked lists of genes, or ChIP-Seq, 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.
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, mmHG-Finder, which is publicly available, to study PWM motifs in several datasets. In addition to validating known motifs, we found GC-rich strings to be enriched amongst the promoter sequences of long non-coding RNAs that are specifically expressed in thyroid and prostate tissue samples and observed a statistical association with tissue specific CpG hypo-methylation.
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 prostate-specific long non-coding RNAs are regulated by transcription factors that bind GC-rich sequences, such as EGR1, SP1 and E2F3. We further suggest that this regulation is associated with DNA hypo-methylation.
- Statistical enrichment
- Position weight matrices
- High-throughput sequencing data analysis
- Tissue specific methylation patterns
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 , 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  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 ChIP-Seq.
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 ChIP-Seq 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 ). To address mutual enrichment in two ranked lists over the same set of N elements, miTEA  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,…,n1] is enriched in the first n2 elements of π, that is in the set π(1), …, π(n2). The statistics introduced by miTEA is called mmHG (minimum-minimum-Hyper-Geometric). 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 , 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  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 non-parametric 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 p-value). In this section we describe the mmHG statistics and in later sections we suggest tight bounds for the p-value. Our definition of the mmHG statistics varies slightly from that of Steinfeld et al., 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 ≤ n1, n2 ≤ N such that the observed intersection between the numbers 1,…,n1 and the first n2 elements of π – namely, π(1), …, π(n2) - is the most surprising in terms of the hypergeometric p-value.
In the Results section we present significantly tighter bounds.
Position weight matrix motifs
Data produced by techniques such as ChIP-Seq , ChIP-exo , CLIP , PAR-CLIP  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 . 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 = S1…SK is defined as , 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 position-specific scores for each symbol in the substring. This definition can be generalized to yield a score for a sequence S = S1…SM longer than the PWM by calculating . 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 .
mmHG statistics for PWM motifs
Given a set of sequences that were tested in a high throughput experiment such as ChIP-Seq , CLIP  and others, they can be ranked according to the measured binding affinities, yielding a ranked list L1. 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 de-novo 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 L2, different from L1. 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 L1 and L2. Given two ranked lists L1 and L2 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 L1 and L2 as mmHG score (π), defined above.
Estimation of the mmHG p-value – 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 · N2. A more subtle upper bound was suggested by Steinfeld et al. 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,…,n1] to appear at the first n2 entries of the permutation (there are possibilities). Then, we choose the positions amongst the first n2 entries that are occupied by these b elements, while considering all internal arrangements (for each choice of b elements there are possibilities). We next choose n2-b elements from the range [n1 + 1,…,N] to appear at the rest of the first n2 entries of the permutation (there are possibilities for that) and consider all possible (n2 - b) ! arrangements. Finally, we take into account all possible (N - n2) ! arrangements of the remaining N-n2 entries of the permutation.
This upper bound is simple and requires O(N3) 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 : , which is tight for large factorials.
A refined upper bound for the p-value – B2
The upper bound introduced in the previous section counts the number of permutations for which the value HGT(N, n1, n2, b) is calculated when taking the non-practical exhaustive approach that enumerates over all N! permutations. Ideally, we wish to count the number of permutations for which the value HGT(N, n1, n2, 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 p-value. 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, n1, n2, 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, n1, n2, b) ≤ s, we can exclude from the counting permutations that contain b elements from the range [1,…,n1-1] at their first n2 entries because they are already taken into account in Λ(N, n1 - 1, n2, b) (because necessarily HGT(N, n1 - 1, n2, b) ≤ s, as we will later explain).
The properties of the hypergeometric distribution imply that given a tuple (N, n1, n2, 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,…,n1-1] at their first n2 entries. Recall that we also assume that HGT(N, n1, n2, b) ≤ s. Therefore HGT(N, n1 - 1, n2, 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 n1-1, n2 and b, and can be disregarded for the combination of n1, n2 and b. Similar arguments hold for ψ2 and ψ4.
Since Λ* is recursive, we need to define a base case. Recall that given N, n1 and n2, b can be any integer in the range [max(0, n2 - N + n1), min(n1, n2)], hence determining a base case for n1 and n2 is sufficient (N is known). The base case here is that when n1 ≤ 1 or n2 ≤ 1, Λ*(N, n1, n2, b) is defined the same as Λ(N, n1, n2, 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 p-value, 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 over-counting 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
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 p-value to represent the significance of that PWM being enriched at the top of the list. To apply this approach for de-novo 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 X1 + X2 + X3 + X4 = 10, which is ). 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 mmHG-Finder, 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 mmHG-Finder. 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 [k1, k2].
Build a generalized suffix tree for the input sequences.
For each k=k 1,…,k 2
Traverse the tree to find all k-mers
Sort the k-mers according to their enrichment at the top of the list (this is done using the mHG statistics), as explained in Leibovich et al..
Take the most significant fifty k-mers, 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 k-mer 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 p-value 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 fine-tuning the correct weights of each letter in each position of the resulting PWM. Finally, the expanded motif is converted to a PWM.
The PWMs found in the previous step are assessed using the mmHG statistics and the best PWMs are returned as output, together with their p-value. 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/mmHG-Finder-code/.
PWM motif search in long-non-coding RNA sequences
We searched for de-novo PWM motifs in the promoter sequences of the tissue-specific lncRNAs using mmHG-Finder (introduced in the previous section). Promoter sequences were defined as 1000 bp upstream the transcription start site.
We scanned the promoter sequences of the tissue-specific lncRNAs with PWMs corresponding to known transcription factors, downloaded from the JASPAR database .
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 ChIP-Base database , which aggregates high-throughput sequencing data taken from hundreds of ChIP-Seq experiments.
CpG hypo-methylation in tissue-specific lncRNA promoters
Mutual enrichment between GC content and tissue-specificity
Mutual enrichment between hypo-methylation and tissue-specificity
No methylation data
No methylation data
No methylation data
No methylation data
No methylation data
No methylation data
No methylation data
No methylation data
No methylation data
No methylation data
No methylation data
White blood cell
No methylation data
No methylation data
No methylation data
No methylation data
No methylation data
No methylation data
No methylation data
No methylation data
No methylation data
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 mmHG-Finder in thyroid, prostate and other tissues (Figure 3; the comparison was done using the motif discovery tool TOMTOM ). The full output of the second and the third tests are summarized in Additional file 2.
As GC-rich 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 hypo-methylation and tissue specificity. For that, we downloaded genome wide (450 K) CpG Methylation data from UCSC Table Browser  (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 hypo-methylation and tissue-specificity (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 prostate-specific lncRNA promoters were less methylated than non-prostate-specific lncRNAs, and this was much stronger in normal cells than in cancer (mmHG p-value ≤ 4.16e-11 in normal prostate cells, and 0.002 in prostate adenocarcinoma cells). We observed strong mutual enrichment between CpG hypo-methylation and tissue-specificity also in brain, ovary, breast and kidney. That is, significant mutual enrichment values were found for tissues where tissue-specific lncRNAs had GC-rich promoter sequences, but these values were not significant for tissues that did not show such GC-bias (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].
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 p-value 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 mis-labeling results as non-significant. For several representative datasets, we note that our bound improves the Bonferroni derived p-value estimates by a factor of almost 40, on average. Nevertheless, these improvements become relevant only for high p-values - 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 p-value 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 thyroid-specific. EGR1 was observed to be highly expressed in thyroid (Additional file 1, taken from ), 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 (p-value 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 tissue-specific lncRNAs and by protecting the associated promoters from methylation. EGR1 has been previously shown to recognize GC-rich consensus sequences located in CpG island promoters of active genes . Generally, we observed that tissue-specific lncRNA promoters tend to be less methylated than those of non-tissue-specific lncRNAs in prostate, brain, ovary, breast and kidney, which may be associated with the GC-rich patterns enriched among their tissue-specific lncRNA promoter sequences.
Threshold-free alternatives to mmHG include the work of McLeay and Bailey, in which a linear regression method is applied . It was shown to achieve high accuracy on a benchmark comprising 237 ChIP-chip 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 non-specific binding . 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 high-binding affinity with respect to the protein of interest (that is, for sequences that are located at the top of the list) . This weaker relationship is naturally addressed by the mmHG statistics. A combination of mmHG and a linear model, such as suggested in , 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 user-friendly PWM motif search in ranked lists. The software is freely available at http://bioinfo.cs.technion.ac.il/people/zohar/mmHG-Finder-code/. Finally, the full characterization of the distribution of mmHG as a random variable over S N remains an open question.
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 tissue-specific long non-coding RNA regulation, and suggest that tissue-specific lncRNAs are regulated through GC-rich elements located on their promoters, in several tissue types. We hypothesize that these GC-rich patterns are associated with DNA hypo-methylation.
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.
- 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 knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005, 102: 15545-15550. 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: 48-10.1186/1471-2105-10-48.View ArticlePubMedPubMed CentralGoogle Scholar
- GOrilla Webserver. [http://cbl-gorilla.cs.technion.ac.il/]
- Ragle-Aure M, Steinfeld I, Baumbusch LO, Liestøl K, Lipson D, Nyberg S, Naume B, Sahlberg KK, Kristensen VN, Børresen-Dale A-L, Lingjærde OC, Yakhini Z: Identifying in-trans process associated genes in breast cancer by integrated analysis of copy number and expression data. PLoS ONE. 2013, 8: e53014-10.1371/journal.pone.0053014.View ArticleGoogle Scholar
- Akavia UD, Litvin O, Kim J, Sanchez-Garcia 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: 1005-1017. 10.1016/j.cell.2010.11.013.View ArticlePubMedPubMed CentralGoogle Scholar
- Dehan E, Ben-Dor A, Liao W, Lipson D, Frimer H, Rienstein S, Simansky D, Krupsky M, Yaron P, Friedman E, Rechavi G, Perlman M, Aviram-Goldring A, Izraeli S, Bittner M, Yakhini Z, Kaminski N: Chromosomal aberrations and gene expression profiles in non-small cell lung cancer. Lung Cancer. 2007, 56: 175-184. 10.1016/j.lungcan.2006.12.010.View ArticlePubMedGoogle Scholar
- Al-Shahrour F, Díaz-Uriarte R, Dopazo J: FatiGO: a web tool for finding significant associations of gene ontology terms with groups of genes. Bioinformatics. 2004, 20: 578-580. 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: 5832-5847. 10.1093/nar/gks206.View ArticlePubMedPubMed CentralGoogle Scholar
- Leibovich L, Paz I, Yakhini Z, Mandel-Gutfreund Y: DRIMust: a web server for discovering rank imbalanced motifs using suffix trees. Nucleic Acids Res. 2013, 41: W174-W179. 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: e45-e45. 10.1093/nar/gks1142.View ArticlePubMedPubMed CentralGoogle Scholar
- miTEA Webserver. [http://cbl-gorilla.cs.technion.ac.il/miTEA/]
- Enerly E, Steinfeld I, Kleivi K, Leivonen S-K, Ragle-Aure 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ørresen-Dale A-L: miRNA-mRNA integrated analysis reveals roles for miRNAs in primary breast tumors. PLoS ONE. 2011, 6: e16915-10.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: 1156-1165. 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: e39-10.1371/journal.pcbi.0030039.View ArticlePubMedPubMed CentralGoogle Scholar
- Steinfeld I, Navon R, Ardigò D, Zavaroni I, Yakhini Z: Clinically driven semi-supervised class discovery in gene expression data. Bioinformatics. 2008, 24: i90-i97. 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: 564-571. 10.1038/nsmb.1594.View ArticlePubMedGoogle Scholar
- Lee B-K, Bhinge AA, Iyer VR: Wide-ranging functions of E2F4 in transcriptional activation and repression revealed by genome-wide analysis. Nucleic Acids Res. 2011, 39: 3558-3573. 10.1093/nar/gkq1313.View ArticlePubMedPubMed CentralGoogle Scholar
- Rhee Ho S, Pugh BF: Comprehensive genome-wide protein-DNA interactions detected at single-nucleotide resolution. Cell. 2011, 147: 1408-1419. 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: Transcriptome-wide analysis of regulatory interactions of the RNA-binding protein HuR. Molecular Cell. 2011, 43: 340-352. 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 A-C, Munschauer M, Ulrich A, Wardle GS, Dewell S, Zavolan M, Tuschl T: Transcriptome-wide identification of RNA-binding protein and MicroRNA target sites by PAR-CLIP. Cell. 2010, 141: 129-141. 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: 505-519. 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: 6661-6679. 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: 563-577. 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: 137-144. 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: e454-e463. 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: 51-80.Google Scholar
- Bailey TL: DREME: motif discovery in transcription factor ChIP-seq data. Bioinformatics. 2011, 27: 1653-1659. 10.1093/bioinformatics/btr261.View ArticlePubMedPubMed CentralGoogle Scholar
- Luehr S, Hartmann H, Söding J: The XXmotif web server for eXhaustive, weight matriX-based motif discovery in nucleotide sequences. Nucleic Acids Res. 2012, 40: W104-W109. 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 genome-wide p53-binding sites upon stress response. Nucleic Acids Res. 2008, 36: 3639-3654. 10.1093/nar/gkn232.View ArticlePubMedPubMed CentralGoogle Scholar
- Harbison CT, Gordon DB, Lee TI, Rinaldi NJ, Macisaac KD, Danford TW, Hannett NM, Tagne J-B, 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: 99-104. 10.1038/nature02800.View ArticlePubMedPubMed CentralGoogle Scholar
- Hogan DJ, Riordan DP, Gerber AP, Herschlag D, Brown PO: Diverse RNA-binding proteins interact with functionally related sets of RNAs. Suggesting an extensive regulatory system. PLoS Biol. 2008, 6: e255-10.1371/journal.pbio.0060255.View ArticlePubMedPubMed CentralGoogle Scholar
- Cabili MN, Trapnell C, Goff L, Koziol M, Tazon-Vega B, Regev A, Rinn JL: Integrative annotation of human large intergenic noncoding RNAs reveals global properties and specific subclasses. Genes Dev. 2011, 25: 1915-1927. 10.1101/gad.17446611.View ArticlePubMedPubMed CentralGoogle Scholar
- Mathelier A, Zhao X, Zhang AW, Parcy F, Worsley-Hunt R, Arenillas DJ, Buchman S, Chen C-y, 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 open-access database of transcription factor binding profiles. Nucleic Acids Res. 2013, 42: D142-D147.View ArticlePubMedPubMed CentralGoogle Scholar
- Yang J-H, Li J-H, Jiang S, Zhou H, Qu L-H: ChIPBase: a database for decoding the transcriptional regulation of long non-coding RNA and microRNA genes from ChIP-Seq data. Nucleic Acids Res. 2013, 41: D177-D187. 10.1093/nar/gks1060.View ArticlePubMedPubMed CentralGoogle Scholar
- Gupta S, Stamatoyannopoulos J, Bailey T, Noble W: Quantifying similarity between motifs. Genome Biol. 2007, 8: R24-10.1186/gb-2007-8-2-r24.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: 435-438. 10.1038/371435a0.View ArticlePubMedGoogle Scholar
- UCSC Table Browser. [http://genome.ucsc.edu/cgi-bin/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 long-range epigenetic remodeling. Cancer Cell. 2013, 23: 9-22. 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: 1475-1483. 10.1158/0008-5472.CAN-13-3042.View ArticlePubMedGoogle Scholar
- Kubosaki A, Tomaru Y, Tagami M, Arner E, Miura H, Suzuki T, Suzuki M, Suzuki H, Hayashizaki Y: Genome-wide investigation of in vivo EGR-1 binding sites in monocytic differentiation. Genome Biol. 2009, 10: R41-10.1186/gb-2009-10-4-r41.View ArticlePubMedPubMed CentralGoogle Scholar
- McLeay R, Bailey T: Motif enrichment analysis: a unified framework and an evaluation on ChIP data. BMC Bioinformatics. 2010, 11: 165-10.1186/1471-2105-11-165.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 non-specific site. J Mol Biol. 1997, 267: 1186-1206. 10.1006/jmbi.1997.0920.View ArticlePubMedGoogle Scholar
- Benos PV, Lapedes AS, Stormo GD: Is there a code for protein-DNA recognition? Probab(ilistical)ly. Bioessays. 2002, 24: 466-475. 10.1002/bies.10073.View ArticlePubMedGoogle Scholar