TrawlerWeb: an online de novo motif discovery tool for next-generation sequencing datasets

Web implementation

Input file

The web interface for TrawlerWeb offers users the option of uploading input sequences as BED indexed format, in addition to FASTA files, which is the most common form of input offered by most motif discovery tools. BED files are lists of genomic intervals and are the direct output from ChIP-seq experiments. Users can therefore directly input the data of their ChIP-seq experiments in TrawlerWeb, without having to retrieve corresponding FASTA sequences. Note that PSCAN [21] also offers BED file input but performs motif discovery on known motifs rather than de novo. DeepSEA [22] accepts input files as BED format, however, currently only one genome (hg19) is supported and input sequence length is limited to 1000 bp. Since peaks obtained from ChIP-seq on histone modifications, for example, could exceed 1000 bp, this option might be restrictive to users. At present, TrawlerWeb remains the only online motif discovery tool accepting input files in BED format for a large range of species.

Background file

For each set of input BED regions, a randomised background specific to the dataset is automatically generated to match the distribution of the genomic locations of the input regions. This input-matched background allows TrawlerWeb to account for sequence biases present in certain genomic regions (for example promoter regions are known to be enriched in CpG islands [29]).

To generate a customised background dataset, first, TrawlerWeb calculates the distribution of the distances of the input regions with respect to the nearest transcription start site (TSS). Each input region is associated to a ‘nearest gene’ and the distance of this region to the gene’s TSS is calculated using gene coordinates downloaded from EnsEMBL using BioMART [26, 30] for the given organism. Distances are then plotted across discrete ranges (e.g. -5000 bp to 0, 0 to 5000 bp, etc.) to produce a frequency table representing the input regions. Next, genes are randomly selected from the entire genome and genomic regions are extracted upstream or downstream of the TSS, so that the distribution of the selected regions match the frequency table generated for input regions. This frequency table is displayed on the results page as TrawlerWeb is running. The amount of randomly selected regions to generate this background dataset is eight times the amount of regions in the input, which we have previously demonstrated to robustly provide adequate background dataset [5, 6].

As a new background will be generated for each new input submission, users have the option to download the background FASTA sequences generated for a specific FASTA input, should they need to re-run the exact same analysis using the same background.

Input options

TrawlerWeb comes with an array of options to optimise the user’s search results. By default, Trawler will search for motifs that are at least 8 bp in length and at most 20 bp. However, the user can reduce the minimum motif length in order to allow motifs of shorter length to be identified.

The wildcard option allows for mismatches in the identified motif. For a minimum motif length of 8, two wildcards are used by default. Should the user choose to reduce the minimum motif length (e.g. 6), one wildcard is recommended to maintain sensitivity.

If the final list of clustered motifs (named families) in the results page retain some similarity and should be clustered together, the “percentage overlap” option can be reduced so that the amount of similarity required between instances to be clustered is reduced.

The frequency at which a motif is expected to occur in a sample can vary depending on the type of data. Generally, a minimum of 10 to 20 occurrences is suitable for most ChIP experiments. If this parameter is set too low, Trawler will be unable to identify any significant results whereas setting the parameter too high when the desired motif is not present in all sample sequences will produce a motif with low information content. This can be configured using the “occurrence” option.

The number of motifs used for clustering is determined by the “number of motifs” option. By default, Trawler takes the top 200 ranked motifs for clustering. To allow for identification of secondary and possibly tertiary motifs, the number of motifs can be increased, however this will also increase computation time.

After running Trawler initially with default settings once, the final number of motifs identified can be fine-tuned by indicating a value for “number of clusters” option. By default, Trawler will cluster by strongly connected component (SCC) when “number of clusters” is set to zero. However, k-means clustering can be used to cluster the motifs into n clusters by providing an integer value (n) to this option.

Output files and data download

The final results are summarised as a web page displaying all discovered motifs in a table along with z-scores for statistical over-representation [5]. Putative matches against known TFBSs and maximum sequence conservation of the motif’s instances will be displayed in this summary web page. By default, motifs are ordered in decreasing order of conservation score, however any of the columns can be used for sorting the output. For each identified motif, further information is available upon clicking on the motif name of logo. First, the distribution of the motif locations within the input sequences can be visualised as a histogram. Consensus sequence, length and identity of the TFBS match are displayed in the next table where a mouse-over on the column titles will provide a description of the properties of the match. A link to the original TFBS will allow to visualise the putative hit from either UniPROBE, Jaspar and HOCOMOCO databases [31–34]. The final table lists the location of every instance of the predicted motif in the input sequences, along with an average and maximum conservation score within the instance. Each location is linked to a dynamic view of the region in the UCSC genome browser [35]. Similarly, a mouse-over on the column’s title will provide further information. All tables in this detailed page are sortable by column, searchable and are dynamically filtered for the searched term. The entire results webpage and files generated by TrawlerWeb are also available for download in a single zip file.

TrawlerWeb runs the fastest amongst popular web-based motif discovery tools

For all five species tested (Table 3), TrawlerWeb ran the fastest (Fig. 1a) in accordance with its standalone version [5, 6]. TrawlerWeb generally returned the discovered motifs in less than 2 min (min) for the four smaller datasets Danio rerio (Dr), Arabidopsis thaliana (At), Mus musculus (Mm) and Drosophila melanogaster (Dm). Only for the human dataset (Hs), being the largest input file tested, running time ranged from 4 min to over 21 min, averaging 10 min overall. RSAT peak-motifs identified motifs after 2-10 min with only a few outliers. Motif discovery for the human dataset ran for 17 to 22 min, however, one run was completed in under 7 min. Nonetheless, we experienced larger variations in processing times with RSAT peak-motifs compared to TrawlerWeb. MEME-ChIP had the longest running time among the three tested tools. Typical motif discovery was completed on average after about 30 min, with the exception of Hs which finished after up to almost 1 h. Of note, MEME-ChIP provides by default two different algorithms, MEME and DREME. DREME [37] discovers short and ungapped motifs. Since MEME did not find any motifs for the zebrafish dataset (Dr) we used the results provided by DREME, which can explain the shorter time compared to the performance of MEME.

In conclusion, in all tested cases, TrawlerWeb outperformed popular web-based de novo motif discovery tools in terms of speed.

TrawlerWeb accurately identifies the expected binding site

To measure TrawlerWeb’s accuracy against other online motif discovery tools, we tested TrawlerWeb against a yeast ChIP benchmark dataset [38]. This same dataset was previously used to test Trawler_standalone algorithm against other software [5, 6] but lacked RSAT peak-motifs [20], which we have included in this run. On this yeast dataset, TrawlerWeb outperformed all other methods in terms of accuracy, identifying 54 out of 65 motifs previously found (85%; Fig. 1b and Additional file 1: Table S1). These results were expected since TrawlerWeb and Trawler_standalone share the same underlying algorithm.

In order to further characterise TrawlerWeb’s accuracy to identify the expected binding sites, we compared the motifs predicted by each tool, across eight users, to the known binding site of the transcription factor of interest, for the five ChIP-seq datasets (Table 3). To identify which predicted motif best matched the known binding site, all predicted motifs with the two highest z-scores (TrawlerWeb) or lowest e-values (RSAT peak-motifs, best of 6 nucleotide or 7 nucleotide length, and MEME/DREME) were aligned and clustered with the expected binding sites (Table 3) using STAMP (default options) [39]. For each ChIP-seq dataset, and for each program, motifs with the shortest distance to the expected binding site were identified as the closest motif. Across all users, the closest primary and secondary motifs (Fig. 2a-e), and corresponding distance to the expected binding site (Fig. 2f), were downloaded from STAMP in Newick format and visualised with iTOL using the expected motif as a reference. In two cases out of five, TrawlerWeb and RSAT predicted the known binding sites (for Zic3.2 (Fig. 2a), TOC1 (Fig. 2b)). TrawlerWeb is the only tested tool that identified the expected binding site for MEF2A (Fig. 2c), albeit as a secondary motif. In the case of Su(H) (Fig. 2d), TrawlerWeb and RSAT peak-motifs identified the same primary motif, which interestingly, is different to the expected binding site. Finally for Sox15.1 (Fig. 2e), RSAT peak-motifs discovered a similar motif to the known binding site, whereas TrawlerWeb found motifs that are quite different. However, these motifs resemble the PWM of Sox15.2 [32] suggesting that the reported binding site of Sox15 [40] is Sox15.2 rather than Sox15.1. In two out of five cases (Fig. 1b, d) the primary motifs discovered by MEME-ChIP agreed with the other motif discovery tools and the expected binding site. For MEF2A and Sox15.1, the motifs identified by MEME-ChIP have low similarity to the known binding sites. For Zic3.2 MEME-ChIP did not find any motifs, hence we used DREME, which found two motifs that are quite distant from the expected binding site.

Of note, in three out of five datasets, TrawlerWeb identifies the predicted motif from the highest number of input sequences. When comparing over-represented motifs (based on the percentage of motif occurrences of the highest score) discovered using DREME, MEME, RSAT peak-motifs, and TrawlerWeb (Fig. 1c), TrawlerWeb returns the highest percentage occurrence in dm3, hg19 and tair9 test datasets (65.5%, 37.6% and 36.6% respectively), and delivers occurrences comparable to both MEME and RSAT peak-motifs in the mm9 and zv9 datasets (1.6% and 0.3%) (Additional file 2: Table S2). Variations in motif occurrence could be due to the presence of a secondary motif corresponding to a co-factor of the TF of interest [11]. Indeed, for TOC1, MEF2A and Su(H) datasets, TrawlerWeb identified a motif different to the known binding site (Fig. 2b-d). Overall, TrawlerWeb robustly identifies the primary motif with high similarity to the expected binding site (Fig. 2f).

TrawlerWeb offers the unique option of displaying motif conservation scores

Motif discovery tools deliver a list of over-represented putative DNA binding sites, usually ranked by over-representation score. Motifs are then often selected for experimental validation, for instance by verifying whether the transcription factor of interest effectively binds to the predicted motif. Amongst all of the instances of the predicted binding site in the submitted sequences, selecting the ones for experimental validation is not trivial. In order to prioritise identified motifs for downstream analysis, evolutionary conservation has been used as a proxy to select for the binding sites which are likely to be functional. Indeed, TFBSs harbouring an essential function are under strong evolutionary constraint compared to neutrally evolving non-coding sequences, and will therefore display higher sequence conservation [41]. We have implemented this feature in TrawlerWeb whereby when input files are provided in BED format, corresponding genomic coordinates are matched against the reference genome selected. This provides the unique advantage for TrawlerWeb over other tools to display conservation scores for every instance of motifs discovered. To display conservation scores, PhastCons scores were downloaded from UCSC [35] and stored along with the chromosome sequences for each organism. For each instance of the identified motif in the input sequences, the average and maximum conservation score is calculated using the “bigWigOverAverage” tool provided in kentUtils [42] package from UCSC. Version 305 of kentUtils source is used as it offers “minMax” option to calculate the maximum value for conservation. The average and maximum conservation scores calculated from PhastCons scores are automatically displayed in the output (Fig. 3). Although this option is only available with the BED input option, it provides scientists the opportunity to filter for the evolutionarily conserved predicted binding sites for downstream biological validations.