The aim of peak finding in ChIP-seq analysis is identification of genomic regions with a high density of mapped sequence tags relative to a measured or estimated background. A simple approach to achieving this goal consists of two steps. Firstly, a sequence of mapped tags along the genome is extracted. Secondly, every contiguous sequence of base pairs with more than a predefined threshold number of tags is selected as an enriched region or binding site. However, the experimental noise and inherent complexities of the tags require more sophisticated algorithms. Numerous solutions have been designed following different statistical models and enrichment measures, including window-based models such as Erange [1], Hidden Markov Model-based methods such as HPeak [2] and others, for example, FindPeaks [3]. Differences in the characteristics of the algorithms result in identification of different sets of peaks for the same ChIP-seq data set. These algorithmic diversities provide an opportunity to analyze ChIP-seq datasets under different conditions, but the problem of deciding which method is optimal for a given data set remains unsolved [4]. Here, in the spirit of protein structure prediction meta-servers, e.g. [5], we present the first meta-solution as a method that combines results from peak-finders chosen by the user.

The implemented meta-server collects results from several peak finders and from these extracts the final result. The peak finders currently included in PFMS are: MACS v1.3.7 [6], CisGenome v2.0 [7], SISSRs v1.4 [8], Erange v2.1 [1], SeqSite v1.0 [9], FindPeaks v3.1.9.2 [3], and HPeak v1.1 [2]. The system can be configured to include any combination of these peak finders. PFMS is implemented in a multi-threading manner. The peak finders can be run in parallel, or sequentially depending on user’s hardware.

Data processing

Seven peak identification methods are currently integrated in PFMS. They are run separately to identify lists of putative peaks. The peaks are collected, normalized and converted to either BED or WIG, as selected by the user. The BED option is more convenient for additional downstream analysis, while WIG is better for visualization purposes since it retains information about the peak shapes.

BED option

PFMS selects an integrated list of significant peaks from the combined results by any of the following methods: voting mechanism, minFP or minFN. The first method does not use any normalization, while the other two use normalized scores of peaks obtained from the individual peak finders. The normalization is carried out using one of five different normalization methods described below.

WIG option

WIG files store information about the shape of the peaks. Firstly, peaks are divided into steps and any operation on the file is in fact done at the level of steps, not the whole peaks. Therefore PFMS normalizes the peak scores using the Average or Naive Quantile method as described below. Ranking the steps or normalizing them using Normal method does not make sense from the statistical point of view, thus these two methods have not been implemented for WIG option. Finally, PFMS integrates the overlapping peaks and sums up their scores. The results are easy to visualize using genome browsers, such as UCSC Genome Browser either in WIG or BED format. An overview of the implemented procedure is shown in Figure 1.

Result generation

BED option

Prior to peak-selection, overlapping regions of the peaks obtained from each peak finder are aggregated and weighted by summing up the normalized scores of the regions. Apart from the aggregated score, the number of votes for each region is kept. The number of votes is defined as the number of peak finders that called the region. Finally, the consensus peaks are selected based on one of the following methods.

In-degree Centrality Voting Mechanism

Is our simplest mechanism of peak selection. All the aggregated regions are considered as candidate sites. A putative binding site is reported if the number of votes for a region is higher than or equal to a predefined threshold (min_rank). Note, that in this method the scores generated by the peak finders are not utilized. For example, given a set of regions 1A: [0 - 7], 1B: [9 - 15], 2: [5 - 25], and 3: [10 - 20] generated by PF 1, PF1, PF2 and PF3, respectively, with the threshold value 2, the regions [5 - 7] and [9 - 20] are reported as binding sites. If the threshold was set to 3, the output region would be [10 - 15] (Figure 2).

minFP peak selection method

Attempts to minimize the number of false positive peaks, at the cost of losing some true positive peaks. In this method the user defines min_rank, which is the minimum number of votes that a region has to receive in order to be identified as a putative peak. We calculate MaxP, the highest score of a peak having fewer than min_rank votes. Each peak that has a score higher than MaxP is reported as a putative binding site.

minFN peak selection method

Attempts to minimize the number of false negative peaks, at the cost of including some false positive peaks. min_rank is defined as in minFP. MinP is the lowest score of the peak having min_rank or more votes. Each peak that has a score higher than MinP is reported as a putative binding site.

WIG option

With the WIG option, peaks are divided into steps, each characterized by a score. The number of scores and the step size determine the width of the peak. All step sizes are unified to the smallest value among all WIG files returned by individual peak finders. All identified peaks are collected and the overlapping regions integrated. The integrated regions are weighted with the normalized score of the overlapping peaks. Optionally, the highest weighted peaks may be selected from the list of the integrated regions by setting a cut-off value. The goal is to select the regions that are called by the majority of the peak finders and that have high scores. This leads to a set of enriched regions that are easy to visualize.

Our method was evaluated against the benchmark datasets published by Rye et al. [12]. They analyzed ChIP-seq reads of three transcription factors (MAX, NRSF and SRF) using five different peak finders (MACS, SISSRs, PeakSeq, FindPeaks and QuEST). Furthermore, the authors visually inspected a number of detected peaks, classifying them as True Positive, False Positive or Ambiguous, providing an excellent resource for evaluation of peak finder performance. In our analysis we ran all seven peak finders included in PFMS on the three datasets incorporating control data with those peak finders that support experimental background measurements. PFMS with various combinations of parameters and peak finders was evaluated. The peaks obtained by the peak finders and by PFMS were intersected with the Rye’s results using BEDTools [10], keeping track of counts of True Positives and False Positives for each data set and parameter setting. The performance results of all runs are shown in Figure 3.

In order to clearly compare the performance of single peak finders and various configurations of PFMS three different measures were applied. Firstly, we calculated Euclidian distance of the points on the Figure 3 to point (0,1) for each of the investigated TFs. Secondly, for each configuration, we calculated Average Euclidean Distance for the three TFs results (see Additional file 1: Table S1). However, in the first measure the largest impact on the overall output have the results for the most 'difficult’ TF. A ranked-based approach was applied to compensate for this factor. We ordered increasingly the Euclidean distances separately for each TF. For each configuration we summed up the ranking position from all three TFs obtaining Combined Ranking Score. The third measure implemented in order to compare the performance of the peak finders and various configurations of PFMS was the Average Normalized Distance. To calculate this measure we applied studentization on the distances for each of the TF datasets. Then, for each peak finder and PFMS configuration we calculated mean value of the measure for the three datasets. For all three measures, the lower the score, the better configuration is.

In our validation analysis, the peak finders were ran with the default settings. However, tuning settings of the peak finders for the given datasets may improve the results. Configuration of the individual peak finders may be easily modified within PFMS (for detailed information see http://bioinf.icm.uu.se/~pfms/). CisGenome, MACS and SISSR were found to perform much better than the remaining four peak finders. Not surprisingly, PFMS performed best with only the three peak finders included. The top three choices for transcription factor dataset analysis utilize only these three peak finders with the min_rank parameter set to 2. We recommend either voting or minFP with Rank Normalization or Top Rank Normalization as these three configurations outperform each of the single peak finders as well as any other configuration of PFMS. Depending on the quality measure, one of the three configurations was the optimal one. However using any of the three quality measures the above mentioned configurations make the top three choices. As predicted, the minFP option returns results with only few False Positives but quite often misses True Positives. Apparently, this leads to a very high specificity and can be useful for applications such as selecting the strongest putative TF target genes for biological validation [13]. We recommend using the minFP option with CisGenome, MACS, SISSR and optionally with SeqSite.

In contrast, minFN has very high sensitivity, i.e. it returns most of the True Positives at the cost of including a considerable number of False Positives. This option also performs best with the four peak finders. The choice of normalization type did not seem to be crucial for the quality of the results when using minFN or minFP. Nevertheless, we do not recommend using the Normal Normalization unless one is certain that the peak scores obtained from all peak finders have a bell-shaped distribution. Otherwise, it is much safer to use Naive Quantile Normalization (e.g. quantile 0.75), Average or Rank Normalization. These methods do not require any assumptions about the distribution of the scores.

The configurations described in this section proved to be the best for transcription factor ChIP-seq datasets. Other types of ChIP-seq data such as histone modifications were not tested. Further investigation needs to be carried out and more validated datasets have to be provided to reveal the optimal settings for both PFMS itself as well as the individual peak finders. As an example, histone modifications ChIP-seq datasets are likely to have varying peak widths and shapes depending on the pattern of the modification (e.g. single or consecutive nucleosomes) and the density of the chromatin. Therefore, different approaches and options might be better for different cases.

PFMS can be used as a single interface for analyzing ChIP-seq datasets employing several peak finders simultaneously since users may choose a set of peak finders amongst the ones currently integrated in this application. In addition to the list of putative peaks identified by PFMS, the results of each peak finder may be stored in the output directory of PFMS.

We present Peak Finder Metaserver - a novel tool for finding peaks in ChIP-seq data. The tool combines the results from various widely used methods and generates consensus results. We investigated seven peak finders and identified three that perform best for transcription factor ChIP-seq datasets, i.e. CisGenome, MACS and SISSR. Applying only these three peak finders and setting the voting peak selection method and the minrank parameter to 2. To the best of our knowledge this is the best method of finding peaks in transcription factor ChIP-seq datasets. Meta-Server approach proved to be successful and PFMS with the above mentioned settings generates results of a better quality than any of the individual peak finders. Different configurations of our tool can be optimal for different types of analyses, but identification of optimal settings requires other validated datasets.