ChIPdig: a comprehensive user-friendly tool for mining multi-sample ChIP-seq data

Implementation

ChIPdig¹⁵ is developed in R using package Shiny and relies on multiple R/Bioconductor packages, including QuasR and BSgenome¹⁶, for mapping reads to a reference genome assembly¹⁷, BayesPeak for peak calling¹⁸, csaw¹² and edgeR¹³ for differential enrichment analysis, ChIPseeker for annotation of genomic regions of interest¹⁹, and EnrichedHeatmap for generating coverage heatmaps²⁰. Packages used in several analysis modules implemented in ChIPdig are GenomicRanges²¹, valr²², shinyFiles, GenomicFeatures, ggplot2, ggsignif, reshape2 and circlize. The GUI allows integrative and interactive usage of these powerful libraries without requiring programming or statistical experience from the user. The required packages may be installed manually by the user or automatically upon launching the tool for the first time. A step-by-step tutorial on how to use ChIPdig is provided at https://github.com/rmesse/ChIPdig.

Operation

R and RStudio versions 3.4.1 and 1.0.44 or higher, respectively, are required, and both memory usage and execution times are contingent on the specificities of the intended analysis (e.g. size of input files, sequencing depths and size of the reference genome assembly).

ChIPdig has the following capabilities: (1) alignment of reads to a reference genome; (2) normalization and comparison of distinct ChIP-seq data sets; (3) annotation of genomic regions; (4) generation of comparative heatmaps and metaplots for visualization of normalized coverage in user-defined specific regions. Each capability corresponds to a specific analysis module which can be loaded separate from other modules and using different input files. For the alignment and annotation modules, the tool automatically fetches the desired reference assembly and gene model, respectively, from public repositories. Analysis of mammalian genome data is supported.

An earlier version of this article can be found on bioRxiv (DOI: https://doi.org/10.1101/220079).

Alignment of ChIP-seq reads

The aligning module of ChIPdig relies primarily on the QuasR package, which supports the analysis of single-read and paired-end deep sequencing experiments¹⁷. If desired, read preprocessing can be performed to prepare the input sequence files prior to alignment, e.g. removal of sequence segments corresponding to adapters and low-quality reads. Sequence files in FASTQ format for each input DNA and ChIP sample corresponding to the H3K4me3 and H3K36me3 ChIP-seq experiments were aligned to the WS220/ce10 reference genome (available through the BSgenome package¹⁶) (Figure 2), producing an average of 10 million mapped reads per sample. The execution time was approximately 2 h 30 min.

Figure 2. Alignment module.

The user supplies a text file listing the input files corresponding to unmapped reads in either uncompressed (e.g.: '.fq', '.fastq') or compressed (e.g.: '.gz', '.bz2', '.xz') format. Both single-end and paired-end data are supported. A scroll-down menu displays the reference genome assemblies available for alignment.

ChIP-seq data normalization, peak calling and differential enrichment analysis

A recurrent problem in ChIP-seq data analysis lies at the comparison of multiple coverage profiles generated from different experiments or corresponding to different conditions. This problem is addressed in the second analysis module of ChIPdig. Upon its selection, the user is prompted to provide a tab-delimited text file with the names of each ChIP sample (treatment) mapped reads file in BAM format, the matching input DNA (control) file, a sample ID, the condition or target name, and a color designation for the output peak and coverage files (Figure 3a). Any of the 657 R built-in colors can be chosen. A bin size parameter set to 50 bp (base pairs) by default is used in both peak calling and genome-wide library normalization for differential enrichment analysis. If desired, duplicate reads can be removed prior to normalization and sequences can be extended to a median fragment length that can be estimated either computationally or experimentally (e.g. average size of fragments generated by chromatin shearing) (Figure 3b). The initial mapped read processing outputs a table with library sizes and median fragment sizes, a chromosome size list, and a multidimensional scaling plot representing the similarity of samples in the data set (Figure 3c). Data are normalized based on library sizes and on the trimmed means of M values (TMM) approach²⁴, which is implemented in the edgeR package¹³. Initial mapped read processing for the H3K4me3 and H3K36me3 ChIP-seq data took approximately 10 min.

Figure 3. Initial processing of mapped reads.

The user supplies a tab-delimited text file listing the input DNA and ChIP sample mapped read files (a). If desired, duplicate read removal and fragment size extension can be performed (b). The application outputs a table with library sizes and median fragment sizes, a chromosome size list, and a multidimensional scaling plot representing the similarity of samples in the data set (c).

Completion of the initial processing of mapped reads releases options for downstream processing which are posted to the sidebar panel, namely export of files representing sequencing coverage, peak calling and differential enrichment analysis (Figure 4). Coverage is expressed in log2-transformed counts per million (cpm) for each genomic bin and, for each sample ID indicated by the user, three files with bedGraph format are generated: treatment (ChIP sample), control (input DNA sample) and the treatment-to-control ratio. Each file can be loaded onto a genome browser for visualization (Figure 5). Generation of coverage files for H3K4me3 and H3K36me3 data took approximately 8 min.

Figure 4. Downstream processing of mapped reads.

The user may generate coverage files in bedGraph format, perform peak calling for each sample, and assess differential enrichment for a set of genomic regions of interest.

Figure 5. H3K4me3 and H3K36me3 ChIP-seq coverage.

Values represent log2-transformed counts per million (cpm) for the ChIP sample subtracted by log2-transformed cpm for the input sample. Files in bedGraph format were exported from ChIPdig and loaded onto the UCSC Genome Browser. These profiles exemplify the well-known scenario whereby H3K4me3 accumulates around the transcription start site and H3K36me3 is enriched at the gene body.

Peak calling is performed via the BayesPeak package by using a hidden Markov model and Bayesian statistical methodology¹⁸. The user specifies a posterior probability (PP) threshold and genomic regions with PP above such threshold are identified as peaks. If replicates are indicated, commonly enriched genomic regions representing replicated peaks can be derived, and a track definition line can be added, allowing each file to be loaded onto the UCSC genome browser for visualization. In addition, a consensus peak set representing all candidate enriched regions across the full data set may be obtained. The user may be interested in choosing this option for posterior differential enrichment analysis. Peak calling was performed for the H3K4me3 and H3K36me3 ChIP-seq data sets, yielding 4994 replicated peaks for H3K4me3 and 25855 replicated peaks for H3K36me3. Peak calling execution time was approximately 12 h.

Differential enrichment analysis resorts to functions implemented in the csaw¹² and edgeR¹³ packages. Following library size and TMM normalization, bin-level coverage is computed and, if desired, bins in which coverage for input DNA samples exceeds that of the corresponding ChIP samples are filtered off. If the user is interested in differential enrichment analysis in a specific set of regions, the corresponding file in BED format has to be provided. To illustrate this feature of ChIPdig, differential enrichment analysis was performed for comparing H3K4me3 coverage with that of H3K36me3 using either H3K4me3 replicated peak coordinates (Figure 6a) or those corresponding to H3K36me3 peaks (Figure 6b), with a false discovery rate (FDR) threshold of 0.1. As expected, at H3K4me3 genomic peak coordinates, H3K4me3 coverage is greater than that of H3K36me3, and the opposite is observed for H3K36me3 peaks. Differential enrichment analysis for both comparisons took approximately 5 min.

Figure 6. Differential enrichment analysis results.

H3K4me3 ChIP-seq coverage was compared with that of H3K36me3 using either (a) H3K4me3 peak genomic coordinates or those of (b) H3K36me3 peaks with a false discovery rate threshold of 0.1. For each analysis, a mean-different plot representing the library size-adjusted log2-transformed fold change (the difference) against the average log2-transformed coverage (the mean), as well as a box-and-whisker plot showing the global change in normalized coverage between the two conditions, were originated.

Annotation of genomic regions

Annotation of genomic regions of interest is performed via the ChIPseeker package¹⁹. The user supplies the file with regions in BED format and chooses the reference genome assembly, as well as the distance upstream and downstream of the annotated transcription start site to be considered for assignment of promoter regions (Figure 7a). Replicated H3K4me3 and H3K36me3 peaks were annotated and, characteristically of these marks²³, most H3K4me3 peaks were assigned to promoters (Figure 7b), whereas H3K36me3 peaks lie predominantly at gene bodies (Figure 7c). The execution time for both annotation operations was approximately 3 min.

Figure 7. Annotation of genomic regions.

(a) The region file in BED format is uploaded and reference genome assembly, as well as distances upstream and downstream of transcription start site, are chosen. (b) Replicated H3K4me3 peaks are mostly promoter-proximal, whereas (c) H3K36me3 peaks are found predominantly at gene bodies.

Comparative heatmaps and metaplots

The user can supply multiple coverage files in bedGraph format and generate heatmaps and metaplots to visualize coverage in a specific region set in a comparative manner. This module of ChIPdig relies on a custom algorithm which builds a coverage matrix based on bin size and reference coordinates (start, end or both) selected by the user. Such matrix is then plotted in the form of a comparative metaplot and a set of heatmaps. H3K4me3 and H3K36me3 coverage files were loaded onto ChIPdig, along with a BED file with C. elegans transcription units in the WS220/ce10 reference assembly. Coverage was computed in 25 bp bins and gene bodies were either expanded or compressed to 500 bp. Windows of 250 bp upstream and downstream of each region were selected (Figure 8a). Typical of H3K4me3, coverage is higher in the vicinity of the transcription start site, whereas H3K36me3 is enriched at transcription unit bodies (Figure 8b, c). Heatmap and metaplot generation took approximately 4 min.

Figure 8. Generation of heatmaps and comparative metaplot for visualizing coverage.

Coverage files for H3K4me3 and H3K36me3 were supplied to ChIPdig, as well as a file with coordinates of C. elegans transcription units. (a) The heatmaps and metaplot were generated by considering a 25-bp bin size, and gene bodies were resized to 500 bp. Windows of 250 bp upstream and downstream of each region were selected. As observed in both (b) the heatmaps and (c) the comparative metaplot, coverage for H3K4me3 is overall higher in the vicinity of transcription start sites, while that of H3K36me3 is enriched at transcription unit bodies.

Comparison with publicly available tools available through the Galaxy framework

In the past decade, several tools have been made available through Galaxy, an interactive system that provides a simple Web portal enabling users to analyze genomic data derived from high-throughput sequencing techniques²⁵. We have compared results obtained by ChIPdig with those obtained within Galaxy using Bowtie³ for read alignment and MACS2^5,6 for peak calling. The percentage of mapped reads was higher for ChIPdig, but it should be noted that no read pre-processing was performed for read mapping in Galaxy (Figure 9a). For H3K4me3 peak calling, the number of called peaks using ChIPdig was similar to that in Galaxy, and, for H3K36me3, more peaks were called in ChIPdig than in Galaxy (Figure 9a). This can be attributed to the fact that the peak calling algorithm used by ChIPdig, and which is implemented in the BayesPeak package¹⁸, generates smaller peaks than those generated by Bowtie. More importantly, the spatial localization of ChIPdig-generated peaks overlapped very significantly with that of Galaxy-generated peaks (Figure 9b). We have also compared coverage tracks generated by ChIPdig (background-subtracted) with those generated by the Galaxy version of DeepTools and observed that the corresponding profiles are very similar (Figure 9c). These observations indicate that ChIP-seq data analysis using ChIPdig compares favorably with that in Galaxy.

Figure 9. Comparison of ChIPdig analysis results with those obtained through the Galaxy framework.

(a) The percentage of mapped reads using the alignment module of ChIPdig is higher than that using Bowtie within Galaxy, but it should be noted that read pre-processing was not included in the latter analysis. Peak calling in ChIPdig and in Galaxy resulted in similar numbers of peaks called for H3K4me3, and more H3K36me3 peaks were called using ChIPdig compared with Galaxy-generated peaks. (b) Importantly, ChIPdig-generated peaks overlapped very significantly with Galaxy-generated peaks. (c) Coverage tracks generated by ChIPdig have a profile similar to those generated by the Galaxy version of DeepTools.