Comparative analysis of commonly used peak calling programs for ChIP-Seq analysis

Abstract

Chromatin immunoprecipitation coupled with high-throughput DNA sequencing (ChIP-Seq) is a powerful technology to profile the location of proteins of interest on a whole-genome scale. To identify the enrichment location of proteins, many programs and algorithms have been proposed. However, none of the commonly used peak calling programs could accurately explain the binding features of target proteins detected by ChIP-Seq. Here, publicly available data on 12 histone modifications, including H3K4ac/me1/me2/me3, H3K9ac/me3, H3K27ac/me3, H3K36me3, H3K56ac, and H3K79me1/me2, generated from a human embryonic stem cell line (H1), were profiled with five peak callers (CisGenome, MACS1, MACS2, PeakSeq, and SISSRs). The performance of the peak calling programs was compared in terms of reproducibility between replicates, examination of enriched regions to variable sequencing depths, the specificity-to-noise signal, and sensitivity of peak prediction. There were no major differences among peak callers when analyzing point source histone modifications. The peak calling results from histone modifications with low fidelity, such as H3K4ac, H3K56ac, and H3K79me1/me2, showed low performance in all parameters, which indicates that their peak positions might not be located accurately. Our comparative results could provide a helpful guide to choose a suitable peak calling program for specific histone modifications.

Keywords: ChIP-Seq; histone modification; human embryonic stem cell; peak calling program.

Fig. 1.

Overview of the analysis. (A) Histone modifications are classified as narrow, broad, and mixed types. (B) Five programs were used for data preparation and peak calling. MACS2 was executed with the default or broad option. (C) The called peaks were compared in terms of enrichment, consistency, and specificity.

Fig. 2.

Pairwise comparison of shared regions. (A) Percentage of peaks recaptured by programs shown pairwise. Each panel shows the percentage of total peaks from one method (column) that was recaptured by another peak caller (row) after filtering blacklist peaks. (B) The concordance rate of peak regions derived from two peak callers. The ranked coincidence weas calculated and the values of percentage and correlation coefficients were denoted after filtering blacklist peaks.

Fig. 3.

Peak consistency between replicates (A) Jaccard correlation coefficient between biological duplicates for each histone ChIP-Seq data. (B) Reproducible peak numbers passing the IDR threshold of 0.01%.

Fig. 4.

Peak coverage with different sequencing depths. The genomic coverage of the regions was shown by sampling with different sequence reads (0.5, 0.75, 1.0, 2.5, 5.0, 7.5, 10, 15, 20, and 30 million).

Fig. 5.

Specificity of peak calling against the noise signal. The specificity of each program was calculated by sampling with different noise levels. Fifty percent (A), 100% (B), and 150% (C) of control reads added.