MOST+: A de novo motif finding approach combining genomic sequence and heterogeneous genome-wide signatures

Background Motifs are regulatory elements that will activate or inhibit the expression of related genes when proteins (such as transcription factors, TFs) bind to them. Therefore, motif finding is important to understand the mechanisms of gene regulation. De novo discovery of regulatory elements, like transcription factor binding sites (TFBSs), has long been a major challenge to gain insight on mechanisms of gene regulation. Recent advances in experimental profiling of genome-wide signals such as histone modifications and DNase I hypersensitivity sites allow scientists to develop better computational methods to enhance motif discovery. However, existing methods for motif finding suffer from high false positive rates and slow speed, and it's difficult to evaluate the performance of these methods systematically. Result Here we present MOST+, a motif finder integrating genomic sequences and genome-wide signals such as intensity and shape features from histone modification marks and DNase I hypersensitivity sites, to improve the prediction accuracy. MOST+ can detect motifs from a large input sequence of about 100 Mbs within a few minutes. Systematic comparison method has been established and MOST+ has been compared with existing methods. Conclusion MOST+ is a fast and accurate de novo method for motif finding by integrating genomic sequence and experimental signals as clues.


Background
Gene expression can be activated or inhibited when proteins (transcription factors, TFs) bind to a genomic sequence segment (i.e. motif) dispersed in promoter and enhancer regions. To gain insights into the mechanisms of gene regulation, the patterns and locations of functional TF binding sites should be identified. Nevertheless, they remain to be challenging problems until high-throughput approaches like ChIP-seq (chromatin immunoprecipitation followed by massively parallel sequencing) [1,2] and ChIP-chip (ChIP on chip) appeared. However, ChIP-seq methods have poor resolution ranging from 200 bps to 300 bps. Besides, their results depended on high quality antibodies for the TFs, which are hard to obtain in general. Moreover, despite the possibility of examining multiple TFs in a single experiment [3,4], it is still prohibitively expensive to test thousands of TFs in various cell types and conditions. Therefore, computational methods are still in a great demand as complementary means to analyze TF binding sites, like PWM (Position Weighted Matrix) guided TFBS identifying and motif discovery.
The major task of motif discovery (or "motif finding") could be viewed as deciphering hidden patterns significantly over-represented in a given genome. Alternative binding motifs, spacing in motifs, palindromes, tandems of single motifs and binding via cofactors make this task complicated. Besides, some cis-regulatory elements are highly degenerated or involve complicated dependency among positions, which is presumably derived from biochemical structure [1], rendering it difficult to discover them. As a result, a motif found by existing algorithms is an approximation instead of an exact reflection of functional binding patterns.
Two types of motif identifying approaches were developed to detect motifs: one uses orthologous sequence alignments extracted from multiple genomes, and the other uses genome regions enriched with particular motifs in a single genome [3]. The first type of methods, such as MEME [5], NestedMICA [6], and Gibbs sampler [7] utilize Bayesian model, EM algorithm, or Gibbs sampling techniques to obtain a Maximum A Posteriori (MAP) for motifs with a given length. Based on these methods, they iteratively update the motif and position with an expectation step followed by a maximization step until some termination criteria are reached [8]. These settings fit quite well with problems of extracting shared sequences from multiple genomes or a small subset of promoter regions (or enhancer regions believed to share some motifs). However, the majority of these methods are slow when processing big datasets. For example, it takes MEME [9] a few days to deal with input sequences of 0.5 Mbps. In addition, due to the probabilistic nature of the models, they are more or less sensitive to noises in input data. Efforts have been made to improve these classic methods over the last few years. STEME [10] speeds up MEME by indexing sequences with a suffix tree. ChIPMunk [11] searches a gapless multiple local alignment using a greedy algorithm to deal with a large set of inputs.
On the other hand, the other type of approaches, word-frequency-based approaches, have been developed, such as WEEDER [12], MDScan [13], Trawler [14], Amadeus [15], DREME [16], and CisFinder [17]. These tools demonstrate significant reduction in computational time and are often better devised for single genome analysis with a large input. Typically, they exhaustively enumerate all possible arrangements of all nucleotides up to a user-specified length and select those words (k-mers containing no degenerated nucleotide) that occur significantly more frequent than those in the Markov background. Some methods iterated by masking stronger motifs first and to find weaker motifs later, while others try to identify all significant motifs, which will then be refined (or clustered) elaborately by expectation-maximization, Self Organizing Map (SOM) [18] or hierarchical clustering. Suffix tree, a data structure to represent the organization of all suffixes of a string, has often been utilized to index input sequences and speed up sequence searching in some word-count based methods like Trawler [14], STEME [10], WEEDER (v1.4.2) and DRI-MUST [19]. With the Ukkonen algorithm, the whole process can be completed in linear time in terms of the size of input [10]. One major challenge for typical wordcount methods is that motifs with weak patterns may be discarded and many spurious motifs may be generated. As a result, a motif can barely be differentiated from the background if its occurrence is relatively low. Limitations like this can hardly be eliminated especially when the size of the input sequences is small [17].
We hope to enhance the word-count methods by utilizing genome-wide signatures from experimental data. Many in vivo events (iVEs) [20], such as open chromatin accessibility or epigenetic modification, strongly associated with the TF binding process, indicating their potential use as biological context to identify TFBSs. For instance, with genome-wide nucleosome occupancy profiles (identified by DNase-Seq or FAIRE-Seq [4,21] in a same cell line, the search for TFBSs on genomic sequence could be confined to only regions with open chromatin accessibility. Meanwhile, epigenetic modifications, especially several kinds of methylations on histones, have also been proved very informative for TFBS finding [4,21]. These modifications can either be associated with activation or repression of corresponding genes in different cell types [22]. Some recent algorithms have already taken advantage of histone mark information to identify individual occurrences of known motifs, such as CENTIPEDE [4], CHRO-MIA [23], FIMO [21] and HINT [24]. Considering a large collection of TFBSs and other CRMs (cis-regulatory modules) are still unrevealed, and sources like histone modification signals can be informative, there is some space to develop better motif finding approaches by taking advantage of epigenetic profiles around TFBS instances.
In this paper, we present MOST+, which can de novo predict motifs and their loci in a genome with the guidance of histone modification and DNase I hypersensitivity information. MOST+ first finds overrepresented words using a suffix tree, then incorporates external genome-wide signals to enhance the search for significant seeds and the aggregation of seeds into motifs. We show that our algorithm integrates successfully the intensity and shape features of external signals to promote the accuracy of motif finding.

Results
We have created fast motif finding methods MOST and MOST+. The system diagram of MOST/MOST+ was shown in Figure 1 and described in Methods section. MOST uses only sequences as the inputs while MOST+ can integrate additional experimental data.
Comparison of MOST+/MOST with existing algorithms MOST+ and MOST were compared with several prevalent motif finding systems, such as MEME, Trawler, WEEDER, CisFinder and HOMER, in terms of speed and accuracy. Different fragments of the mESC ChIP-seq dataset and the genome-wide promoter regions of the mouse genome were used as the gold standard datasets (see Methods for more details).

Comparison of speed
MOST outperformed all tested algorithms ( Figure 2, Table S1) in terms of speed for all datasets with various sizes. MEME [9] had not produced any result after 4 days running on a 16 Mbps dataset. WEEDER failed to produce a result when the input sequence length exceeds 20 Mbps. CisFinder was fast, but had a restriction on the input size (up to 25 Mbps). In short, MOST +/MOST and HOMER [25] (version 4.2) were able to produce results for large datasets in a reasonable time.
For all methods, default parameter settings were used when possible (details can be found in Additional file 1). Word lengths were chosen to a range from 4 bps-7 bps, with an exception of WEEDER (up to 8 bps). Figure 1 The pipeline of MOST+ system. A set of target genomic sequences are extracted from a genome then indexed by a suffix tree to count occurrence of each word (or K-mer). If under MOST+ mode, histone modification marks and/or DNase I hypersensitivity (referred as tag signals in this schema or mark distribution) of each word are used to yield mark distribution scores. Top ranked words are put into clustering and motifs are generated from the resulted clusters. The strategy for clustering is illustrated on the right panel of this figure.
MOST+ is insensitive to the input size. Counting and locating words takes little time. Memory consumption remains relatively low when data size grows up to 100 Mbps (linear spatial complexity, MOST+ requires 500 megabyte memory for an 8 Mbps input).

MOST/MOST+ can find complicate motif patterns
MOST/MOST+ was designed to bridge gaps and find long gap patterns. For instance, Sox2 (canonically represented by CATTGTT) and Oct4 (or known as POU5F1, canonically represented by ATGCAAAT) often occur as a heterodimer (characterized by motifs of OCT-N-SOX or OCT-3N-SOX). These motifs had been detected simultaneously ( Figure S1). MOST+ also successfully detected motif of SMAD1, which has relatively low occurrences in tested ChIP-seq peak regions. For OCT4, MOST+ successfully detected a tandem repeat of OCT4 core motif ( Figure S1).
Alternative TFBSs and palindromes were found for ESrrb (Estrogen-related receptor beta) when a larger word width (K = 11) was employed to capture more sophisticated co-occurrence and a more stringent clustering threshold was set to discriminate sub-patterns ( Figure S1).
Incorporating epigenetic marks helps filter out non-functional motifs MOST reported a similar result compared to the prevalent algorithms. After epigenetic signals were added (MOST+ mode), motifs were better aligned with known motifs in databases and the rank of essayed TF was lifted (Table 1).
We observed that spuriously over-represented k-mers (like some tandem repeat that may not be motif in our dataset) were more likely associated with a higher level of noise and asymmetry (Figure 3a, and Figure S2). Indicating these features may help eliminating false positives. Furthermore, similar motifs can have distinct mark   distribution patterns ( Figure 3b). If guided by epigenetic marks, it would be possible to distinguish similar motifs. In general, MOST+ found more motifs that exist in motif database than MOST, which means more essayed TF or co-factors are found. According to Bieda et al., [26] the motif of E2f1 is very hard to be identified maybe due to indirect binding. However, MOST+ did report a motif (featured by CGCCAT) that ranked second and resembled the motif of E2f-family member E2f3. Detections of n-Myc, c-Myc, zfx and SMAD1 would also benefit from epigenetic marks in terms of highlighting the assayed motifs.
Comparison of motif accuracy on ChIP-seq datasets for mouse and human MOST+ and exiting motif-finding tools were compared using the ChIP-seq dataset of 13 TFs of mouse ES cell.
General comparison was conducted based on alignments with known motifs in the reference database. Comparison on other 5 algorithms showed that MOST+ was among the best algorithms in both capturing major TFBSs (the binding sites of the assayed TFs) and detecting co-factors. MOST+ and HOMER identified 11 of 13 major motifs (with the major motif ranked the first in results) with significant e-value of alignment whereas DREME reported 10 motifs (Table 2). MOST+ can recover a comparable amount of co-factors with DREME, which was devised for finding co-factors (Table 3). Like CisFinder, MOST+ could automatically determine a self-adaptive length for each cluster.
Assessment on site-level accuracy of motif finding methods was conducted by using the pipeline described in Methods section ( Figure 4). Figure S3 showed a comparison of found motifs for different methods. Results show MOST+ has the highest AUROC sum over 13 TFs, though the situation varies from TF to TF ( Figure 5, Figure S4). With parameters learned by part of the dataset, MOST+ achieved better accuracy on recovering motif positions with validation data, suggesting that motifs could reflect actual binding sites better if external signals were available. The AUC (Area Under the Curve) of ROC increased when word counts and mark distribution scores were combined under the appropriate parameter setting. This supported the idea that epigenetic marks can be helpful to cluster words [4]. With this parameter setting learned from training and validation datasets, we compared MOST+ with other 5 algorithms on the remaining partition data. The learned parameters were given in Table S2 (see Additional file 1 for more detail). Each time, we included one feature in our model to test its contribution to the accuracy improvement. Results show that the mere use of asymmetry feature contributes the most to the overall improvement of AUROC, while the utilization of original signals of essayed TF ChIP-seq data failed to show any advantage.
To validate the power of our algorithm, we also applied the model to the human data. On average, MOST+ outperforms MOST by better AUROC, however, with an exception of JUND 1 and 2 ( Figure S5). When epigenetic signals were added to the model, MOST+ found some novel co-factors that MOST could not find. For instance, with the aid of DNase I hypersensitivity information, additional motifs that strongly resembled STAT1, GABPA, LM4 and LM130 (both were long motifs reported by Xie et al.) were found in VDR (Vitamin D Receptor) datasets ( Figure 6).

Genome-wide promoter analysis
To date, large portions of TFBSs are still undetected. However, poor performance has often been reported on capturing motifs in generic promoters [27]. With the power of epigenetic marks, we hoped our overall finding of motifs in all regions close to the transcription start sites (TSSs) might reveal novel motifs. We took an exhaustive search on 89 Mbps mouse promoter regions. When more hints from ChIP-seq data were available, we put on a more stringent threshold of word clustering and masked repeats (see Additional file 1 for MOST+ command lines).
Results showed that 117 motifs discovered by MOST+ (16% of total found motifs) could be aligned with one or more motifs in our mouse motif databases (JASPAR, UNIPROBE and Chen et al., 541 in total, Figure S6). Left panel of Figure 7 showed an example of motifs discovered by MOST+. MOST+ also output hundreds of novel motif candidates, which were not found in mouse motif databases. However, the mark distributions centred by some of these motifs show non-trivial shapes. This indicates either a putative novel motif or a  potential sequence mark of epigenetic modifications (right panel of Figure 7).

Discussion
Our results demonstrated that chromatin accessibility and epigenetic information should not be underestimated for their potential to identify motifs and their co-factors.
In the near future, it is possible to utilize these marks better by introducing different model frameworks, such as conditional random field (CRF [28]) framework, which can estimate the weight for each feature automatically and predict the binding loci more precisely. Figure 4 The diagram of pipeline for parameter optimization and method comparison. A motif-finding step is followed by a TFBS identification step (by CisFinder) using motifs and genomic sequences as input. Training data (8 of 10 folds) are fed into motif finding tools, and then accuracy is evaluated based on how well the motifs recovered can pinpoint TFBSs. AUROC is used to represent the accuracy of each method. This algorithm can be improved in many ways. A parameter advisory phase can be included before MOST+ is executed. More sophisticated analysis like wavelet analysis may be employed to extract oscillation features of mark distributions. The improvement on resolution of ChIP-seq will make the distribution more informative, which can potentially boost the accuracy of the algorithm. Some aspects of the algorithm improvement are discussed as follows.

Setting parameters for epigenetic signals
Parameters for different epigenetic signals and different genomes vary. To make the most out of external signals, the optimal parameters should be learned from the data. However, different genomes may share some common epigenetic features. Users are encouraged to use our pre-set parameters for their genomes of interest with signals like DNase I, H3k4me1 and H3K4me3. Although H3K4me2 is informative in human datasets, it was not informative in mouse dataset (result not shown). Therefore, it was not included as the default signals.

Impact of the window size K
Since the results and running time of MOST+ depend on the window size K, we examined its performances on detecting major motifs of CTCF with distinct K values. Often but not necessarily, bigger K values could capture sophisticated patterns like gaps or alternative binding motifs. However, a very large K could be neither spatially nor temporally efficient in terms of computing. Therefore, we set K = 9 as default for typical ChIP-seq data (with 100 bps peaks) to quickly depict a whole landscape  of a certain motif without losing much sensitivity (see Table 4 for details).

Selecting the background model
The last parameter is the background model for motif finding. Background of the 1 st -order Markov chain barely showed any difference with the 0 th -order (i.e. without considering the dependency of neighbouring bases) in terms of detection ratio. Either the whole genome sequence or merely the target regions can be selected as the background model. During the accuracy assessment, we set the length of flanking region to be 650 bps. However, when it was set as 500 bps and 800 bps, the results were quite similar and robust.

Constructing a motif from overrepresented words
Clustering on primitive motifs was necessary for our methods. In the current version of MOST/MOST+, word instances were extended by substituting bases one by one. Therefore, only a limited space instead of the whole possible space was explored. For instance, in searching motif characterized by ATGCAAAT, it enumerated 4*8 = 32 single nucleotide substitutions to simulate position frequency matrices (PFM) whereas the actual possible solution space might be larger (for example, word instances with more than one substitutions, like ACCCAAAT, were not considered).

Conclusion
Genome-wide epigenetic or DNase I signals can improve motif finding significantly. In this paper, we present MOST+, an accurate and fast motif discovering system through combining the genomic sequence and these genome-wide signals.
There are four highlights for MOST+: 1). It is one of the first de novo motif finding systems to combine the word frequency information and genome-wide signals from different sources such as histone modification information and DNase I hypersensitivity information. 2). It is a fast system for motif finding. MOST+ is 2 orders of magnitude faster than DREME. 3). It can deal with a large size of input data. For example, MOST+ can find motifs from an input data of 100 Mbps. 4). Its accuracy is at least comparable with the best existing systems.

Mouse dataset
A benchmark ChIP-seq experiment of 13 mESC TFs associated with maintenance of pluripotency in mouse embryonic stem cell (cell type EF0_000462) [29] were collected. Nine epigenetic marks [22] in the same cell line were tested: DNase I, H3K27me3, H3K36me3, H3K4me1, H3K4me2, H3K4me3, H3K9me3, H3 and H4K20me3. Three of them (H3K4me1, H3K4me3, DNase I) were selected for subsequent analysis in which each TF against these three marks showed strong and distinct patterns of either kurtosis or single spike ( Figure S7 and S8). We also included a raw ChIP-seq signal track. ChIP-seq peaks of each TF (controlled by removing peaks also found in anti-GFP ChIP-seq under the same cell line) are called using MACS [30]. All tracks were lifted and aligned to reference genome mm9.
For mESC date set, we included 146 motifs from JAS-PAR, 386 from UNIPROBE and 13 by Chen et al., 541 motifs in total, to form our mouse motifs database. We identified motifs using TOMTOM, a motif comparing and visualization tool [31]. Motif alignments with E-values less than 0.05 were considered statistically significant.

Human dataset
We also evaluated our method on human lymphoblastic cell line GM12878. 4 TFs were tested: CTCF, JUND, MAX (From UCSC genome browser) and VDR ChIP-seq (From Ramagopalan et al. [32]). DNase I hypersensitivity signals and histone modification marks H3k4me1, H3k4me2 and H3k4me3 (all from ENCODE project [33]) were included. ChIP-seq peaks were called by MACS. All tracks were mapped to human reference genome hg18. For human data we used 1,268 known motifs selected from JASPAR and UNIPROBE as the motif database (See  table S3).

Genome-wide promoter regions of the mouse genome
An exhaustive search is conducted on the promoter regions of the whole mouse genome (upstream/downstream 1000bps of transcription start sites). The total size of the whole mouse promoter region is 89 Mbps. Transcription start site information (of mm9) was retrieved from UCSC RefSeq track. For time complexity comparison, different sizes of datasets are generated from this dataset by cutting this promoter region dataset. The test was run on mouse CTCF dataset (8 Mbp). Columns 3-5: 3) Times: running time under each word width K; 4) Motifs predicted: number of unique motifs found by each schema; and 5) E-value: the lowest e-value of alignments between a major motif found by each schema and motifs in mouse motif database.

MOST+ pipelines
MOST+ contains two modes: one for word frequency counting only (MOST) and the other for combining different genome-wide signatures (MOST+). The general design of MOST+ is shown in Figure 1. MOST finds overrepresented words only from genomic sequences (such as a set of ChIP-seq peak regions of 200 bps), while MOST+ incorporates signals like histone modification marks and DNase I hypersensitivity sites to enhance the search of significant seeds. Finally, the program exports results including a list of found motifs and their corresponding position frequency matrices (PFM) (See Additional file 1 for more details).

Exploiting epigenetic signals
Mark distributions nearby one particular site can indicate whether this site tends to be bound by TFs. However, it is weak and noisy. Thus we aggregated the aligned signals for all occurrences of a particular k-mer seed. We denote this overall distribution as mark distribution for this kmer seed. The mark distribution is discretized over a set of fixed-size bins (illustrated by Figure S9).

Finding over-represented seeds
MOST and MOST+ find each word's occurrences by utilizing a suffix tree. In order to reduce the space complexity, a hash table is utilized to represent the suffix tree. Ukkonen algorithm is employed to construct the suffix tree. Counting and locating each occurrence can be done by traversing the whole tree in asymptotically linear time (see Additional file 1 for more details). The null model we adopted to find over-represented kmer seeds was generated from an n th -order Markov chain [34] derived from the target regions. We assume the counts of each k-mer are in Poisson distribution. The pvalues can be estimated using Gaussian approximation.
Where wc and λ denote the observed and simulated average counts of certain k-mer respectively. φ is the cumulative distribution function of the standard normal distribution. For simplicity the calculation does not correct for the effects of self-overlapping words. We directly use word counts averaging on 5 simulated sequences with the same length of the query sequence, which was reported to have no substantial effect on result [14]. Nor do we correct for multiple testing, which is intractable since null odds ratio score distributions vary from word to word. Instead, we choose a small p-value threshold.
Ranking words MOST+ chooses a list of candidate motif seeds from all possible k-mers by applying a ranking step. If MOST+ mode is selected, aside from over-representation score, external experimental signals are decomposed into 3 subfeature scores for each candidate seed: intensity scores (characterizing average signal intensity over all occurrence), kurtosis score (characterizing the peak pattern of a mark distribution) and asymmetry score (characterizing the irregularity of a mark distribution). To calculate kurtosis and asymmetry score the mark distribution was first normalized, then smoothed (by taking Discrete Fourier Transform and removing high frequency components). Finally we measured the kurtosis and asymmetry level. An overall score integrating all 4 scores above to rank candidate seeds is calculated as follows.
Where S I , S K ,S A and S w represent score of intensity, kurtosis, asymmetry and over-representing odds ratio. ϕ(X) is the aggregated distribution of the genome-wide signals. f(X) represents normalized and smoothed distribution. The β values are weights for each score that need to be found. Symmetrised Kullback-Leibler divergence (KL-divergence) is adopted to measure the asymmetry of distributions. lhs and rhs denote left and right half of a distribution. wo and λ are defined as in equation (1).

Generalizing seeds to motifs and Clustering
For each high-ranked seed, we generalized them into motifs of the same length by substituting each position of a word with other 3 oligonucleotides, and calculated each variant's occurrence, which produced a probability matrix of this word.
In order to remove redundancy and form full-length motifs, these primary motifs were further clustered. Here we use a similar approach as CisFinder [17] which clusters words by PFMs since PFMs contain more sitespecific details about each motif.
Our clustering procedure initially calculated a dissimilarity matrix according to primary motifs and corresponding mark distributions, which have a time complexity of O(M 2 L 2 ) for M motifs with length L. All possible offsets and orientations of two motifs were checked to decide how they overlapped, and the one with maximum Pearson correlation was chosen as the distance of two motifs. If under MOST+ mode, the KL-divergence between mark distributions of these two motifs could be also incorporated to calculate this distance ( Figure S10). Pseudocounts could be added to avoid zero probability.
After that, an agglomerative hierarchical clustering was applied using average-linkage approaches ( Figure  S11). Gap statistics [35] was employed to determine the optimal cluster number. At last, clusters were trimmed to remove ambiguous or low-count ends.

Assessment of site level accuracy
It is not straightforward to evaluate the accuracy of a motif discovery algorithm [2,8]. We evaluated the accuracy of our motif finding algorithm implicitly according to how well the found motifs can recall real sites.
With motifs found by MOST and other tested algorithms as inputs, we used CisFinder to predict bound sites on target regions. These predicted bound sites were compared with ChIP-seq signals as implicit evaluation of motif accuracy. However, as gold standard annotations, ChIP-seq has limited resolution. Therefore, we considered predicted sites that fell nearby ChIP-seq summits as true positives. Specifically, we flanked each ChIP-seq summit with 650bps at both sides (total 1,300bps for a single TFBS annotation, discretized into 13 100bp bins). The central bin of 100 bps was denoted as a TFBS bin (centred with the annotated TFBS summit), while other 12 bins at each end were denoted as flanking bins ( Figure S12). We evaluated the performance with a nine-fold cross validation (Figure 4). We randomly divided all data into ten partitions with equal sizes. Eight training partitions were used for MOST+/MOST and other methods to de novo find motifs. Then these found motifs were used to pinpoint TFBS (using CisFinder) on the other two partitions -One validation partition for determining appropriate parameters, and the other partition for making comparison over algorithms. In case the essayed motif ranks lower in motif finding result compared with its cofactors, we consider the top 3 motifs for each TF dataset. Finally we compared the predictions with ChIP-seq annotation, and further drew ROC (Receiver Operator Characteristic) curve using above defined positives and negatives. AUROC (Area Under ROC) was derived and used as the measurement of accuracy. The pipeline was repeated on different assignment of training and validation partitions to calculate the mean and standard variation of AUROC.

MOST+ parameter determination
During the training procedure, parameters of MOST+ need to be decided. MOST+ includes 4 parameters (intensity, kurtosis, asymmetry and clustering) concerning external tag information, 3 (intensity, kurtosis and asymmetry) during ranking phase and the remaining one in the clustering phase. To find sets of parameters for various TFs, the parameter spaces are searched systematically for a collection of discrete points, and the best parameter set is selected from thousands of combinations (See Additional file 1 for details).

Comparison with existing methods
Existing methods, such as DREME, Trawler, WEEDER, HOMER2 and CisFinder, have been compared with MOST/MOST+. The command lines used to run these tools are listed at the end of the Additional file 1.

Accessibility and requirement
MOST+ is implemented in C++. The source code is compiled with G++ 4.1.2. All tests were performed in a Linux server with CentOS 5.7, Intel Xeon E5520 @ 2.27GHz processors and a total memory of 16 Gigabytes. MOST+ and all supporting data are freely available at: http://cbb.sjtu.edu.cn/~ccwei/pub/software/ MOST/MOST.php.

Additional material
Additional file 1: This file contains supporting information for MOST+, including result figures S1 to S12, Tables S1 to S3 and some other detailed information like the parameter settings for the programs compared in this paper.