Improving rigor and reproducibility in chromatin immunoprecipitation assay data analysis workflows with Rocketchip

Implementation

Rocketchip is an automated bioinformatics workflow written in the Python-based (v3.10.12) workflow manager, Snakemake (v7.32.4)²¹ ( Figure 1). Rocketchip takes full advantage of Snakemake’s capabilities to ensure a highly flexible and adaptable workflow. By consolidating Snakemake rules from various software tools into a single file tailored to the user’s data, Rocketchip supports the seamless incorporation of new tools and updates. This design is crucial in the rapidly evolving field of bioinformatics, allowing Rocketchip to stay current with the latest advancements without requiring a complete overhaul of the existing pipeline. Additionally, Rocketchip efficiently handles large-scale genomic studies by leveraging Snakemake’s ability to process hundreds of files and utilize parallelization based on available computational resources. This makes it suitable for handling extensive datasets on high-performance computing (HPC) systems, ensuring that users can perform complex analyses with optimal resource allocation. Furthermore, Rocketchip utilizes Snakemake’s error-handling capabilities, including a log of standard error outputs per sample, enabling users to investigate and address specific data processing issues as well as utilizing the capability to resume workflows from the point of failure rather than re-running the entire pipeline. This error-handling capability is further enhanced by the inclusion of two FastQC steps in Rocketchip’s workflow, one at the beginning and one after deduplication, to provide additional checkpoints for assessing sequence quality and detecting potential issues.

Figure 1. Rocketchip pipeline for ChIP-seq data analysis.

Rocketchip begins by downloading raw sequencing data and the reference genome, followed by an initial quality control check. Next, the raw data is aligned to the reference genome using the user’s choice of alignment software, and PCR duplicates are optionally removed. A second quality control check is performed after deduplication to assess alignment quality and any remaining artifacts. Aligned data is then converted into bigwig files for visualization, and peak calling is performed using customizable options for peak type and control usage. The final outputs include processed data ready for downstream analysis or visualization in genome browsers.

Operation

Rocketchip downloads ChIP-seq data directly from the SRA, the largest publicly available sequence database, using the SRA Toolkit (sra-tools, 3.0.9).^2,22 The SRA Toolkit downloads data via the prefetch and fasterq-dump functions, then splits the raw DNA sequence read file into respective paired-end (PE) read files using the flag --split-files and converts them into FASTQ file formats, whereas for single-end (SE) read files, the SRA Toolkit is utilized solely for the download and conversion of raw read data into FASTQ format. Rocketchip also presents the option for users to use local (i.e. non-SRA sourced) ChIP-seq data by providing their own FASTQ files. In parallel, Rocketchip downloads and processes a reference genome of the user’s choice from the UCSC Genome Browser.²³ After downloading and converting files, files are stored in the local file system in FASTA and FASTQ formats. Users also have the option of storing their own custom genomes and sequencing reads that are not yet publicly available.

In Rocketchip, raw sequence data undergoes a quality control step using FastQC (v0.12.1) with default parameters to assess potential problems that could impact downstream analysis, including levels of data duplication, base pair quality, adapter contamination, and overall sequence quality.²⁴ This initial FastQC step provides users with the ability to investigate samples if downstream analysis fails, which can occur if there are problems involving sequence quality, contamination, or high levels of data duplication. Raw sequence data is then aligned to the reference genome using the user’s choice of alignment software from BWA-MEM (v0.7.17), Bowtie2 (v2.5.2), or STAR (v2.7.11a), each with the default parameters.^25–27 Intermediate data files are processed as necessary for deduplication. PCR duplicates are removed from sequence data using a user’s choice of deduplication software from Samtools (v1.18) with --mode s (i.e. standard PCR duplicate detection), Picard (v3.1.1) using MarkDuplicates, Sambamba (v1.0.0) with default parameters, or no deduplication.^28–30 Deduplicated data files are subject to another quality control step using FastQC to provide further information, such as alignment quality, validating processing steps (i.e. success of deduplication), and detecting any remaining biases or artifacts (e.g. GC content, adapter contamination, etc.) to ensure overall data integrity. Next, Deeptools’ (v3.5.4) bamCoverage function with default parameters is used to convert data from the BAM file format to the bigwig file format, which can be used for visualization of ChIP-seq data in the UCSC Genome Browser or other visualization tools.³¹ Finally, in Rocketchip, peaks are called with a user’s choice of software from MACS3 (v3.0.0b3) with the --bdg option, Genrich (v0.6.1) using default parameters, PePr (v1.1.24) using default parameters, or CisGenome (v2.0).^32–35 For CisGenome, if a control was used, the seqpeak command was employed with the default options. Without a control, CisGenome ran two rounds of peak calling using the hts_peakdetectorv2 command with options “-w 100 -s 25 -c 10 -br 1 -brl 30 -ssf 1.”

Peak-calling options in Rocketchip include narrow versus broad peak calling and use of a control where appropriate. Software version control is handled using Conda to ensure reproducibility of results and compatibility between software versions.³⁶ Software options were chosen based on options for command line use, ease of installation, and their standard use in the field.

Users interested in using Rocketchip can install the source code via PyPi (i.e. `pip install Rocketchip`). Rocketchip requires a minimum Python version >=3.6. Associated dependencies are Conda installable and are stored in an environment file (https://github.com/vhaghani26/rocketchip/blob/main/environment.yml) with recommended software versions compatible with Rocketchip. This allows user flexibility in installing only dependencies of interest (e.g. only MACS3 instead of MACS3, Genrich, and PePr) and changing dependency versions if desired for algorithm benchmarking purposes or to ensure compatibility with the user’s environment during installation. After installation, the user must create a project YAML file containing relevant project details, including which genome to use for alignment, sample names and paths if using local data, and software selections for the workflow to execute. Rocketchip will use this input to generate a Snakefile tailored specifically to the user’s data, project details, and requested software combination with appropriate parameters. Then, the user will run the Snakefile using Snakemake’s command-line interface arguments, which are briefly detailed in the Rocketchip GitHub repository (https://github.com/vhaghani26/rocketchip). Overall, this approach offers users flexibility to efficiently test various software combinations for their data. It simplifies the execution of an entire workflow into just a few steps: installation, creating a project file, running the Rocketchip command, and executing the workflow with Snakemake. Additionally, it enhances the reproducibility of results by documenting all inputs, outputs, and software versions, while leveraging the automated nature of Snakemake’s workflow execution.

Conducting ChIP-seq data analysis via all software combinations

In order to ensure that all combinations of software can be integrated seamlessly and used for ChIP-seq data analysis, we selected a deeply sequenced experimental ChIP-seq study targeting a transcription factor, MeCP2, in mouse main olfactory epithelium conducted by Rube et al. 2016, hereafter referred to as “Rube”,³⁷ and a ChIP-seq study targeting NRF2 in human non-small lung cells conducted by Namani et al. 2019, hereafter referred to as “Namani”.³⁸ To facilitate this comprehensive assessment, these data sets were run through Rocketchip using all four peak callers, namely MACS3, CisGenome, Genrich, and PePr. This was used in combination with each of the three aligners, BWA-MEM, Bowtie2, and STAR and four deduplication techniques, Samtools, Sambamba, Picard, and no deduplication. Additionally, to validate that Rocketchip can be run with or without a control (i.e. input or IgG control), this analysis was conducted with both the usage and omission of the corresponding control for each data set, with the exceptions of CisGenome and PePr, which can only be run if a control is used in the analysis. Each test was run three times to assess Rocketchip’s ability to replicate experimental results.

Each algorithm demonstrated varying peak-calling efficiency, influenced by several key factors ( Figure 2, Supplementary Table 1). Most notably, the source of the data (i.e, Namani vs. Rube) revealed dramatic differences in the performance of peak-calling algorithms. For the Rube data, MACS3 consistently yielded the highest number of called peaks. All methods of deduplication produced comparable peak counts. Bowtie2 and BWA-MEM performed similarly, identifying slightly more peaks compared to STAR. CisGenome called the fewest peaks, significantly influenced by the deduplicator and aligner used. Without deduplication, CisGenome identified the fewest peaks across any software combination, while peak counts increased with any deduplication method. When STAR was used as the aligner with CisGenome, peak counts were the lowest. BWA-MEM yielded a higher peak count than Bowtie2 but exhibited more non-deterministic behavior compared to Bowtie2 and STAR. For Genrich, peak counts increased when the control was omitted during peak calling. In contrast to CisGenome, Genrich had the highest peak counts when no deduplication was used and when STAR was used for alignment. The other deduplicators and aligners showed negligible differences in peak counts. When using PePr for peak calling on the Rube data, no deduplication yielded the highest peak count, with other deduplication methods yielding similar results. Unlike the other algorithms, each aligner performed differently with PePr, with STAR identifying the most peaks, followed by BWA-MEM and then Bowtie2. Both the narrow- and broad-peak-calling algorithms yielded negligible differences, except for CisGenome, which had less deterministic results when peaks were defined as narrow and run through Rocketchip.

Figure 2. Peak counts for all software combinations.

ChIP-seq data from Namani et al. 2019 and Rube et al. 2016 were run through all software combinations in Rocketchip three times each. The data source (i.e. Rube vs. Namani) is yellow, with the light yellow corresponding to the Namani data and darker yellow corresponding to the Rube data. A subset of analyses were run with a control vs. no control. This is represented by light pink for a control being used during peak-calling and dark pink if a control was omitted during peak-calling. The data was run through the software combinations as both a broad peak (light orange) and narrow peak (dark orange) to evaluate algorithm performances, though the Namani data uses NRF2, which is characterized best by narrow peaks, and the Rube data uses MeCP2, which is characterized best by broad peaks. The aligners are represented in red, with Bowtie2 being a light red, BWA-MEM being a medium shade of red, and STAR being a dark red. The peak-caller is highlighted in green with shades of green getting progressively darker for Cisgenome, Genrich, MACS3, and PePr. Similarly, the deduplicator is represented by shades of purple, with no deduplication (“None”) being the lightest, followed by Picard, Sambamba, and then Samtools, which is the darkest. Because each combination was run three times, each of the three trials is plotted by shades of blue. The first trial is a light shade of blue, the second trial is a medium shade, and the third trial is the darkest shade. Raw peak counts were log2 transformed and plotted in the heatmap. Darker red corresponds to higher peak counts while darker blue corresponds to lower peak counts. Gray corresponds to “NA” values, as PePr and CisGenome cannot be run without a control. Raw peak counts can be found in Supplementary Table 1. The heatmap was created using R (v4.2.3) with a kernel (r-irkernel v1.3.2) in Jupyter Notebook (v1.0.0) using the following packages: ComplexHeatmap (v2.14.0), dplyr (v1.1.4), tidyr (v1.3.1), reshape2 (v1.4.4), stringr (v1.5.1), and MASS (v7.3.60.0.1).

For the Namani data, notable contrasts in peak-calling outcomes were observed across various algorithms. Genrich consistently yielded the highest number of peaks, contrasting with the Rube data where MACS3 showed higher peak counts. Interestingly, PePr consistently produced the lowest peak counts for the Namani dataset despite performing second best for the Rube data. The choice of aligner did not significantly impact peak counts overall; however, omitting deduplication generally resulted in slightly higher peak counts across all peak-callers. Unlike the Rube data, where the distinction between narrow- and broad-peak calling showed minimal differences across algorithms, the Namani data exhibited variability based on this distinction. Specifically, MACS3 identified higher peak counts under the assumption of broad peaks compared to narrow peaks.

In addition to assessing the results of the different algorithms, we also assessed Rocketchip’s ability to replicate experimental findings. Both the Namani and Rube data sets were run through every software combination in Rocketchip three times, totaling 288 software combinations with three trials per combination. Among the 288 software combinations tested, 274 combinations (95.14%) demonstrated perfect replication of peak counts across all three trials (i.e. every trial yielded identical peak counts). Surprisingly, despite identical software versions, computational resources, and inputs, certain combinations of data and software were non-deterministic, with 4.86% (14 out of 288) exhibiting variability in peak counts ( Table 1). Given that non-deterministic results were observed exclusively for the Rube data, this data was downsampled to 5,000 reads per sample using Seqtk (v1.4) and rerun through every combination 100 times to further assess the source of variation. While variation in peak counts persisted, the maximum range between the minimum and maximum observed peaks was limited to one across all combinations (see Supplementary Table 2). Notably, the software combinations that produced non-deterministic results differed from those in the initial analysis. Therefore, we hypothesize that the low level of variation observed in our results may stem from runtime errors.

Overall, these findings underscore four critical points. First, the choice of software combination, including the algorithm, deduplication method, and aligner, has a significant impact on peak-calling outcomes. Second, even with strict control of factors impacting the analysis, we still observe some variability in peak calling, albeit in a small subset of cases. Third, dataset-specific nuances significantly impact the performance of different software combinations, resulting in a differing consensus on what the “best” software combination is. Finally, the results validate that all available software combinations can be successfully executed within Rocketchip, ensuring flexibility and robustness in ChIP-seq analyses.

Benchmarking peak detection on simulated data and evaluating algorithm performance using F1 scores, precision, and sensitivity

Benchmarking algorithm performance on real biological data is challenging, as true binding sites are not known. To illustrate Rocketchip’s ability to benchmark algorithm performance, six simulated test data sets were designed for each combination of PE and SE reads with narrow and broad peaks. Each data set has two replicates and a synthetic genome composed of uniformly distributed base pairs. The dataset is modeled to have precisely 1000 peaks with varying numbers of 150 base pair long reads per peak to represent various peak coverage depths for different data sets: 4X, 8X, 16X, 32X, 64X, and 128X. Peaks were evenly distributed 2500 base pairs apart for narrow peaks and 5000 base pairs apart for broad peaks. Read locations within peaks were simulated with a Gaussian distribution centered on the peak summit, spanning 400 base pairs for narrow peaks and 1500 base pairs for broad peaks. Reads were dispersed with a standard deviation of 0.1 to ensure variability around the peak summit. A study assessing PCR duplication in ChIP-seq found that PCR duplication ranged from 2.94-55.38% among their tested ChIP-seq data sets, ultimately concluding that a rate of 20% duplication or less is most effective for read mapping.³⁹ Thus, PCR duplicates were simulated at a 20% rate. Input controls, which represent sonicated DNA fragments, were modeled as randomly distributed genomic reads.

This dataset was run through all combinations in Rocketchip, including the inclusion and exclusion of the input control, which accounts for genomic noise and is believed to improve accuracy. CisGenome and PePr require a control, so they were only tested with the input control. We varied read type (PE and SE) and peak type (narrow and broad) to systematically evaluate software performance across all combinations. The objective of this experiment was to assess accuracy, precision, and sensitivity of each combination, where accuracy, which is measured using F1 scores, measures how well peaks match true binding sites, precision indicates the proportion of true peaks among identified peaks, and sensitivity measures the proportion of true peaks detected.

F1, precision, and sensitivity scores were calculated based on true positive (TP), false positive (FP), or false negative (FN) base pairs within peaks. TP represents overlapping base pairs between expected and observed peaks, FP represents base pairs in observed peaks not matching expected peaks, and FN represents base pairs in expected peaks not observed. F1 scores were calculated as [2TP/(2TP + FP + FN)], precision as [TP/(TP + FP)], and sensitivity as [TP/(TP + FN)]. For Cisgenome and PePr, TP, FP, and FN were set to zero, as they cannot be run without a control, indicating an inability to be accurate, precise, or sensitive without using a control.

The results of this experiment revealed that, generally, narrow peak-calling had higher accuracy, precision, and sensitivity than broad peak-calling ( Figure 3). Interestingly, the peak callers with the highest F1, precision, and sensitivity scores, Genrich and MACS3, performed better when run without an input control compared to with a control. MACS3 and Genrich had the highest accuracy, precision, and sensitivity, followed by PePr and Cisgenome, which is consistent with the analyses using the Rube and Namani data. No major differences were observed based on the aligner or deduplicator used in any of the analysis. Overall, this analysis provides further insights into the algorithm performances, allowing for more robust and informed decisions in ChIP-seq data analysis.

Figure 3. Accuracy, precision, and sensitivity measures for rocketchip software combinations.

Simulated ChIP-seq data was run through all software combinations in Rocketchip. For each heatmap, the endedness of the data (PE vs. SE) is pink, with the light pink corresponding to PE data and dark pink corresponding to SE data. The data was run using broad peak (light orange) and narrow peak (dark orange) algorithms depending on whether broad or narrow peaks were modeled, respectively. The aligners are represented in red, with Bowtie2 being a light red, BWA-MEM being a medium shade of red, and STAR being a dark red. The peak-caller is highlighted in green with shades of green getting progressively darker for Cisgenome, Genrich, MACS3, and PePr. Similarly, the deduplicator is represented by shades of purple, with no deduplication (“None”) being the lightest, followed by Picard, Sambamba, and then Samtools, which is the darkest. The peak coverage, ranging from 4X-128X is depicted using shades of blue, with the shade getting darker as peak coverage increases for the given data set. Darker red shading in the plotted heatmap values corresponds to higher accuracy, precision, or sensitivity, while darker blue corresponds to lower accuracy, precision, or sensitivity. (A) The plotted values represent the calculated F1 scores when incorporating a control during peak-calling. (B) The plotted values represent the calculated F1 scores when omitting a control during peak-calling. (C) The plotted values represent the calculated precision scores when incorporating a control during peak-calling. (D) The plotted values represent the calculated precision scores when omitting a control during peak-calling. (E) The plotted values represent the calculated sensitivity scores when incorporating a control during peak-calling. (F) The plotted values represent the calculated sensitivity scores when omitting a control during peak-calling. The heatmaps were created using R (v4.2.3) with a kernel (r-irkernel v1.3.2) in Jupyter Notebook (v1.0.0) using the following packages: ComplexHeatmap (v2.14.0), dplyr (v1.1.4), tidyr (v1.3.1), reshape2 (v1.4.4), stringr (v1.5.1), and MASS (v7.3.60.0.1).

Demonstrating rocketchip’s compatibility with experimental ChIP-seq data from varying genomes

Rocketchip was tested for compatibility with a range of genomes from the UCSC Genome Browser, which includes the seven commonly used genomes displayed by default in the “Genomes” tab: human (hg38), mouse (mm10), rat (rn6), zebrafish (danRer11), fruitfly (dm6), worm (ce11), and yeast (sacCer3). Three published ChIP-seq data sets on the SRA with varying read coverage were selected per genome and run through Rocketchip ( Figure 4).^40–51 All experiments used PE data and were conducted using BWA-MEM for alignment, Samtools for deduplication, and MACS3 for peak-calling. These experiments were run on an HPC with 64 CPUs and 250 GB of memory available, running one job at a time without parallelization. To ensure consistency, genome copies were deleted between runs using the same genome to ensure that the run time accounts for the full workflow. The results of this experiment validated Rocketchip’s compatibility with the seven primary genomes from the UCSC Genome Browser and use of sequence data hosted on the SRA. Additionally, it provides users with estimates of how long Rocketchip should run for different genomes to allow for appropriate computational resource requests. Run times ranged from 0.55 hours to 20.54 hours depending on the genome size and read data used.

Figure 4. Rocketchip execution times per genome.

This is a horizontal bar plot representing how long (in hours) each sample per genome took to run. The label on the bars represents the time in strictly hours (i.e. not hours and minutes) that each sample took to run through the pipeline from beginning to end assuming no prior storage of any data required of the analysis. The color of the bars corresponds to which organism the sample comes from and which genome it was run with: red for dm6, orange for hg38, green for mm10, turquoise for rn6, blue for ce11, purple for sacCer3, and pink for danRer11. The bar plot was created using R (v4.2.3) with a kernel (r-irkernel v1.3.2) in Jupyter Notebook (v1.0.0) using the following packages: ggplot2 (v3.4.4) and dplyr (v1.1.4).

Conducting CUT&RUN and CUT&Tag data analysis with rocketchip

As of June 18, 2024, there are 363,213 ChIP-seq, 494,128 CUT&RUN, and 19,492 CUT&Tag data sets available on the SRA. Due to the increasing use of CUT&RUN and CUT&Tag, we wanted to assess Rocketchip’s ability to effectively process data generated by these mapping techniques. Therefore, we applied Rocketchip to CUT&RUN and CUT&Tag data generated by Akdogan-Ozdilek et al. that sought to characterize the zebrafish epigenome during embryogenesis.⁵² This data set was chosen due to the thorough documentation of the results and high sequence data quality.

In evaluating the performance of Rocketchip for CUT&RUN and CUT&Tag data analysis, we assessed several key metrics, including the percentage of reads aligned, sequence duplication percentages, and peak counts ( Table 2). The alignment percentage was chosen as a metric due to its significance in assessing the overall data processing efficiency and alignment accuracy. Alignment percentage serves as a key indicator of how effectively Rocketchip handles the unique characteristics of CUT&RUN and CUT&Tag data sets, ensuring that a substantial proportion of reads are appropriately mapped to the reference genome and available for further analysis. Sequence duplication percentages were also evaluated, as they indicate Rocketchip’s ability to reduce data redundancy while maintaining data quality throughout the pipeline. Additionally, peak counts provide insight into Rocketchip’s ability identify and quantify expected protein binding sites. Comparisons between processed read counts and peak counts demonstrate that Rocketchip scales appropriately with the number of reads remaining after data processing, reflecting its capacity to accurately assess and interpret the genomic data.

The CUT&RUN data set consisted of nine samples. Six samples were SE reads and corresponded to two replicates each for detection of H3K4me3, H3K27me3, and H3K9me3 (SRR14850825 and SRR14850826, SRR14850827 and SRR14850828, and SRR14850829 and SRR14850830, respectively). The remaining three sets were PE read samples that corresponded to two replicates for the detection of RNA polymerase II and mock IP control using the IgG antibody (SRR14850831 and SRR14850832, and SRR14850833, respectively). The study that originally produced these data sets employed Bowtie2 to align sequences, Samtools to filter aligned sequences, and HOMER for peak-calling. Rocketchip was run using Bowtie2, Samtools, and MACS3 for broad peak-calling. As the data was partially SE and partially PE reads, analyses were conducted separately, with the control being used for the PE analysis due to compatibility. Sequence duplication percentages varied across samples, with some samples exhibiting high duplication rates, which is indicative of data redundancy that Rocketchip managed effectively. Comparisons between processed read counts and peak counts demonstrated that Rocketchip scales appropriately with the number of reads remaining after data processing, reflecting its capacity to accurately assess and interpret genomic data. Original alignment percentages were compared to those obtained via Rocketchip. A paired t-test was conducted using each SRA input as an observation, yielding a p-value of 0.00302, with Rocketchip consistently producing significantly better alignment percentages for the CUT&RUN data compared to the original analysis. The mean difference in the alignment percentages was 32.66%. We also noted that Rocketchip detected a substantial number of peaks in relation to the processed reads, indicating its ability to effectively parse and process signal (peaks) versus biological noise. We hypothesize that the alignment accuracy differences are related to the use of different Bowtie2 versions, as the original study used Bowtie 2.4.1 to align their CUT&RUN data, whereas Rocketchip used version 2.5.2. The GitHub change log for Bowtie2 version 2.5.1, one version earlier than the one Rocketchip uses, notes: “fixed an issue affecting bowtie2 alignment accuracy” and “fixed a segmentation fault that would occur while aligning SRA data.” This highlights the need to revisit and utilize publicly available data, as software updates can improve data processing and thus yield more accurate results. This also highlights Rocketchip’s ability to accurately and effectively process CUT&RUN data.

The CUT&Tag data set used for testing in Rocketchip was comprised of six PE read samples. There are three replicates each of H2A.Z at 6 hours and 24 hours post fertilization (SRR14870792, SRR14870793, and SRR14870794 and SRR14870795, SRR14870796, and SRR14870797, respectively). The original study used Bowtie2 for alignment, Samtools for filtering, Picard for deduplication, and MACS2 for peak-calling. We ran Rocketchip using Bowtie2, Samtools, and MACS3. Similar to the CUT&RUN data, sequence duplication percentages varied, and Rocketchip was able to efficiently process the data without encountering any errors. Original alignment percentages were compared to those obtained via Rocketchip. A paired t-test was conducted using each SRA input as an observation, yielding a p-value of 0.00015, with Rocketchip resulting in lower alignment percentages for the CUT&Tag data compared to the original analysis. The magnitude of difference in alignment percentages, however, was less than that of the CUT&RUN data, as the mean difference by Rocketchip was -8.41% for CUT&Tag as opposed to +32.66% for CUT&RUN. However, the alignment percentages achieved in our Rocketchip CUT&Tag analysis still surpassed every alignment percentage reported in the original study for their CUT&RUN data, suggesting that Rocketchip is suitable for analyzing CUT&Tag data. We hypothesize that the alignment difference may be due to the usage of both Samtools and Picard in the original CUT&Tag analysis as opposed to just using Samtools in Rocketchip.