
==== Front
DNA Res
DNA Res
dnares
DNA Research: An International Journal for Rapid Publication of Reports on Genes and Genomes
1340-2838
1756-1663
Oxford University Press UK

38102723
10.1093/dnares/dsad026
dsad026
Research Article
AcademicSubjects/MED00774
AcademicSubjects/SCI01140
AcademicSubjects/SCI01140
Churros: a Docker-based pipeline for large-scale epigenomic analysis
https://orcid.org/0000-0003-3110-0605
Wang Jiankang School of Biomedical Sciences, Hunan University, Changsha, Hunan, China
Institute for Quantitative Biosciences, The University of Tokyo, Bunkyo-ku, Tokyo, Japan

https://orcid.org/0000-0003-3019-5817
Nakato Ryuichiro Institute for Quantitative Biosciences, The University of Tokyo, Bunkyo-ku, Tokyo, Japan

To whom correspondence should be addressed. Tel. +81 3 5841 1471. Fax: +81 3 5841 7308. Email: rnakato@iqb.u-tokyo.ac.jp
2 2024
15 12 2023
15 12 2023
31 1 dsad02628 9 2023
23 11 2023
13 12 2023
05 1 2024
© The Author(s) 2023. Published by Oxford University Press on behalf of Kazusa DNA Research Institute.
2023
https://creativecommons.org/licenses/by-nc/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution-NonCommercial License (https://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com

Abstract

The epigenome, which reflects the modifications on chromatin or DNA sequences, provides crucial insight into gene expression regulation and cellular activity. With the continuous accumulation of epigenomic datasets such as chromatin immunoprecipitation followed by sequencing (ChIP-seq) data, there is a great demand for a streamlined pipeline to consistently process them, especially for large-dataset comparisons involving hundreds of samples. Here, we present Churros, an end-to-end epigenomic analysis pipeline that is environmentally independent and optimized for handling large-scale data. We successfully demonstrated the effectiveness of Churros by analyzing large-scale ChIP-seq datasets with the hg38 or Telomere-to-Telomere (T2T) human reference genome. We found that applying T2T to the typical analysis workflow has important impacts on read mapping, quality checks, and peak calling. We also introduced a useful feature to study context-specific epigenomic landscapes. Churros will contribute a comprehensive and unified resource for analyzing large-scale epigenomic data.

epigenomics analysis
bioinformatics pipeline
Docker
T2T genome
large-scale ChIP-seq
==== Body
pmc1. Introduction

Epigenomics using next-generation sequencing can depict various modifications to DNA sequences in a genome-wide manner.1 Epigenomic data are crucial for elucidating biological mechanisms and disease-associated pathologies.2,3 The most common epigenomic technologies include: chromatin immunoprecipitation followed by sequencing (ChIP-seq),4 which is used to study histone modifications, transcription factors and chromatin regulators; assay for transposase accessible chromatin using sequencing (ATAC-seq),5 which assesses chromatin accessibility; and DNA methylation sequencing,5 such as whole-genome bisulfite sequencing and reduced representation bisulfite sequencings. A more recent technology, namely cleavage under targets & tagmentation (CUT&TAG),6 enables the detection of specific epigenomic features with relatively greater precision and sensitivity. With the rapidly decreasing cost of sequencing, epigenomic datasets are now widely used in diverse areas of biomedical research. For example, international projects including ENCODE7 and IHEC8 have generated and analyzed extensive epigenomic data to understand the functional elements within the human genome.

Despite major efforts to develop methods for systematically exploration of the ever-increasing scale and diversity of epigenomic data, computational analysis for such large-scale comparisons involving hundreds of samples is still not straightforward. Epigenomic analysis comprises multiple steps, each using different tools and requiring laborious computational configurations.4 Users also need to experiment with various combinations of tools, which requires specific knowledge and experience in conjunction with each step.9 Format-conversion problems often arise, as the output file of one step may not be compatible with the input of the next step. In some cases, data filtering between steps brings additional complications. Moreover, researchers often encounter practical problems, as handling large datasets is not trivial but rather highly error-prone. For instance, the writing and maintaining an abundance of custom scripts for the multistep processing of large datasets; naming output files and scripts in a consistent way; keeping all parameters fixed and generating comparable results across samples during multiple updates of the entire workflow; and automating all streamlined steps with customized intermediate steps while maintaining reproducibility of published results. Consequently, an epigenomic analysis pipeline equipped with pre-compiled computational environments, automated processes with customizable details, and optimizations for large-scale datasets would be particularly helpful for addressing these challenges.

For this purpose, we developed Churros (https://churros.readthedocs.io/), a portable, downloadable, reproducible, and customizable pipeline for large-scale epigenomic analysis. Churros utilizes a Docker system (https://www.docker.com/), a platform that allows users to package, distribute, and run programs in a containerized environment.10 The docker environment of Churros is equipped with comprehensive pre-installed tools and custom scripts. Churros also provides pre-made genome references and all required databases for over 14 genome builds. For users new to computational analysis, Churros can automate all the steps required for a typical analysis and seamlessly produce standard outputs. For advanced users, Churros offers sophisticated parameters that allow customized analysis during each intermediate step. Regarding large-scale datasets, Churros implements several optimizations such as metadata input, parallel computing, unified output filenames, standard directory structure, aggregated analysis reports, and specialized visualizations.

To demonstrate the effectiveness of Churros, we applied it to large-scale ChIP-seq data. Especially, Churros supports epigenomic analysis with T2T (telomere-to-telomere), a recently introduced, complete human genome assembly.11 Using Churros, we systematically compared the results between T2T and the previous human reference genome hg38. We found that applying T2T to the typical analysis workflow slightly increased the mapping rate and affected the correlation analysis among samples. Churros visualizations suggested that T2T does provide additional epigenomic information compared with hg38. T2T reduced the ChIP-seq quality metric normalized strand coefficient (NSC), likely owing to the increased mapping rate. Unexpectedly, T2T could either decrease or increase the number of peaks called by MACS2, which is associated with changes in estimated fragment length. Finally, we introduced a function for clustering and visualizing large-scale epigenomic profiles. Churros and the results presented in this study should accelerate the data-driven epigenomic analysis that lead to the discover of new biological insights.

2. Materials and methods

2.1. Processing epigenomic data with Churros

Churros provides a one-stop solution for epigenomic analysis. Churros supports many popular programs, as well as original scripts and functions. In Churros, public data are fetched by SRAtoolkit.12 Mapping of sequencing reads is accomplished using Bowtie,13 Bowtie2,14 BWA15 and Chromap.16 Quality control of the sequencing data is performed by FastQC (https://github.com/s-andrews/FastQC) and fastp.17 Adapter trimming is conducted by Cutdapt18 and TrimGalore (https://github.com/FelixKrueger/TrimGalore). Peak calling is performed by MACS219 and DROMPAplus.4 Quality assessment of ChIP-seq data is performed by SSP.20 Motif identification is performed by Homer.21 Peak annotation is done by ChIPseeker.22 Other downstream analysis and visualization of epigenomic data are done by DROMPAplus4 and deepTools.23 Churros also includes other useful tools such as super-enhancer analysis by ROSE,24 chromatin state segmentation by ChromHMM11 and epilogos (https://github.com/meuleman/epilogos), data imputation by ChromImpute,25 and gene regulation prediction by STITCHIT.2 Bisulfite sequencing data is analyzed by Bismark.26 File processing and data format conversion are performed by SAMtools27 and BEDtools.28 In addition, Churros utilizes MultiQC29 to summarize the main analysis results for multiple tools and samples in a single report. The main code for the analysis in this paper is provided as Supplementary Data 4.

2.2. Installation and execution of Churros

The initial obstacle that researchers must confront is setting up specific software environments by installing the required tools. Installing dozens of tools with controlling versions across multiple computing environments is often problematic. Fortunately, Docker container, which is used by Churros, allows for an isolated and self-contained computational environment. Users can easily install Churros by downloading a single pre-built, ready-to-run image that contains all required programs and configurations, which eliminates the need to install numerous tools with complex dependencies.10

The Churros Docker image can be installed with the one-line command ‘docker pull rnakato/churros’. After that, users can use the full functionality of Churros with another one-line command ‘docker run rnakato/churros’. An alternative way to install Churros is Singularity,30 another containerization technology specifically designed for scientific research using high-performance computing systems. Singularity has two main advantages: it does not require root privileges, and the entire Singularity container can be packaged into a single file (e.g. ‘Churros.sif’). Users can create the Churros image with the one-line command ‘singularity build churros.sif docker://rnakato/churros’ and then access the full functionality with the command ‘singularity exec churros.sif’.

2.3. Preparing input data and reference databases

The raw sequencing data used in Churros can be in FASTQ or CSFASTQ format. Both single-end and paired-end data are supported. It is also possible to assign multiple separate sequencing files to one sample. In addition, Churros takes a metadata file containing information for all samples, such as FASTQ file locations, labels, and corresponding input samples, to eliminate manual errors and maintain reproducibility when processing large datasets.

Before starting an epigenome analysis, all dependent datasets of the reference genome (e.g. genome, genes, and mapping indices) must also be carefully constructed. These command-line-based tasks are often complex and time-consuming for novices. In Churros, the command ‘download_genomedata.sh’ can automatically download and generate reference genomes, gene annotations and relevant databases, while the command ‘build-index.sh’ can construct the mapping index files. Currently, Churros supports 14 types of genome builds: human (hg38, hg19, and T2T), mouse (mm39 and mm10), rat (rn7), fly (dm6), zebrafish (danRer11), chicken (galGal6), African clawed frog (xenLae2), Caenorhabditis elegans (ce11), Saccharomyces cerevisiae (sacCer3), Schizosaccharomyces pombe (SPombe), and Hydra vulgaris (HVAEP). Specifically, the T2T genome data is downloaded from https://genomeinformatics.github.io/CHM13v2; SPombe genome data is downloaded from https://www.pombase.org; Hydra genome data is downloaded from https://research.nhgri.nih.gov/HydraAEP/; and the data for other genome build are downloaded from UCSC and Ensembl.

2.4. Parallel computing

To facilitate the processing of large datasets, Churros employs two strategies for parallel computing: (i) for tools that support multithreading (e.g. Bowtie2, SAMtools), Churros passes the -p parameter directly to the corresponding option; (ii) for other tools that do not support multithreading (e.g. MACS2, bedtools), Churros executes multiple jobs for many samples in parallel.

2.5. Comparison of results from the T2T and hg38 analyses

To compare the analysis results between hg38 and T2T, we used a total of 196 ChIP-seq samples, including 156 samples for MCF-7 cells and 40 samples for other cell types. All data were obtained from the GEO and ENCODE databases, accession numbers for which are listed in Supplementary Data 1. We applied Churros to the FASTQ files with either the T2T or hg38 reference genome. All results were obtained by Churros without installing any additional tools. For MACS2 peak calling, we used the default parameters: sharp model, q-value cutoff of 0.05, and building the shifting model.

To reflect the change in values from hg38 to T2T, we defined a diff_ratio as:

diff_ratio=NT2T−Nhg38max(NT2T,Nhg38)×100

where N is a specific value of interest. Whereas a typical percentage of change (NT2T−Nhg38)/Nhg38×100 generates a value range from −100% to +∞%, diff_ratio ranges from −100% to 100%, providing a balanced representation of the degree of change in terms of decreases or increases.

2.6. Clustering of epigenomic profiles

A notable feature of Churros is the ‘classheat’ function (classification and heatmap) that clusters and visualizes epigenomic profiles of large datasets. This function takes genomic regions of interest (e.g. specific protein binding sites) as input 1 and epigenomic signal files as input 2. Classheat has two modes: binary mode and continuous mode. In the binary mode, classheat outputs a binary matrix (output 1) representing the overlap of epigenomic markers (usually ChIP-seq peaks) at given genomic regions. The binary matrix is then formatted and sorted by the user-defined column (i.e. the filename of the selected marker) to generate the processed matrix (output 2) and plot the sorted heatmap (output 3). Next, classheat uses principal component analysis followed by mini-batch k-means (a variation of k-means that is faster and more efficient for large datasets) clustering to produce the clustered matrix (output 4) and the clustered heatmap (output 5). Other clustering methods such as k-means, DBSCAN, spectral clustering, mean shift, and affinity propagation are also available. In the continuous mode, classheat calculates the averaged read density of each sample at given genomic regions (output 1). After logarithmic transformation, z-score normalization (optional method is 0-to-1 scaling), and sorting, classheat generates the remaining outputs in the same manner as the binary mode. All these steps can be done automatically with the single-line command ‘classheat.sh binary input1 input2’.

2.7. Code availability and reference documentation

The source code of Churros is available at https://github.com/rnakato/Churros. The docker image of Churros can be found at https://hub.docker.com/r/rnakato/churros. Reference documentation can be found at https://churros.readthedocs.io/. Churros will be continuously updated to incorporate the latest advances in epigenomic analysis.

3. Results

3.1. Structure and implementation of Churros

The current version of Churros supports epigenomic data from ChIP-seq, CUT&TAG, ATAC-seq, and DNA methylation sequencing. This paper primarily focuses on large-scale ChIP-seq analysis. To simplify the tedious task of preparing many input files and to improve reproducibility, Churros uses two formatted tables as input (Fig. 1A). The ‘samplelist.tsv’ table provides locations and corresponding labels for all sequencing files, whereas the ‘samplepairlist.csv’ table provides information about matches between ChIP and input samples, as well as other computational parameters. Another preparatory step is the construction of dependency databases.31 Notably, Churros can automatically download and appropriately construct all required databases, including reference genome sequences, alignment indexes, gene annotations, and other program-specific databases (Fig. 1B) for 11 species involving 14 genome builds (see Section 2).

Figure 1. Construction of Churros. (A) The two types of metadata inputs required by Churros. (B) Churros automatically prepares the reference data (e.g. genome assemblies, annotations, and indexes) for various species and genome builds. (C) The main data processing workflow in Churros. The pipeline can be executed automatically or incrementally. (D) Churros supports a variety of downstream analyses to further explore the epigenomic data.

After preparing the inputs and dependency databases for a specific genome build, users can execute all analysis steps in a streamlined manner using only a single-line command (Fig. 1C). Automated processes within Churros generate a series of outputs, including sequencing quality, mapping statistics, alignment results (in bam or cram format), read coverage (in bigwig format), ChIP-seq quality metrics, peak calling results, pairwise correlation analysis, reads distribution visualizations (genome-wide, 5-kb resolution and 5-Mb resolution), a formatted analysis report, and running logs. These analyses are valuable for both bench biologists and bioinformaticians who want to preliminarily check their datasets. For example, the formatted report based on MultiQC29 provides rich analysis results in a single html file (Supplementary Data 2). Alternatively, users could obtain the same output by executing the four core steps individually: churros_mapping, churros_callpeak, churros_visualize, and churros_compare (Fig. 1C). This is especially useful if users need to modify certain parameters and re-run the analysis starting with intermediate outputs.

Another feature of Churros is the capability to perform many customized analyses (Fig. 1D), which are provided either by original code or by optimized scripts based on existing tools. For example, users can cluster epigenomic profiles, visualize multi-omics data, investigate the aggregation of epigenomic markers, examine correlations and overlaps between samples, segment the genome by chromatin states, and so on.

3.2. Comparison of epigenomic analysis based on hg38 and T2T

In 2022, the Telomere-to-Telomere (T2T) consortium published the gapless assembly for the complete human reference genome.32 Despite the valuable insights provided by the newly introduced sequences,11 it remains unclear whether a current epigenomic analysis pipeline using the T2T reference genome would lead to significantly different results. Because Churros supports both T2T and hg38, we systematically compared the results obtained using these two genome builds. Churros successfully processed all 196 public ChIP-seq datasets obtained from the GEO and ENCODE databases, demonstrating its compatibility with large-scale epigenomic datasets. It should be noted that the highly repetitive sequences in T2T require specialized mapping algorithms and/or long-read sequencing data.32 Here we attempted to simply map single-end short-read data to both hg38 and T2T using Bowtie2, as our primary goal was to investigate whether a typical analytical workflow using T2T could outperform the results obtained with hg38. The scope of our study was not to investigate the biological insights from a T2T analysis because this aspect has been explored in previous studies.11,32,33

We first compared mapping rate considering only uniquely mapped reads (aligned exactly one time). Fig. 2A illustrates the difference in mapping rate after subtracting T2T from hg38 (Supplementary Data 3). We observed that all ChIP-seq samples had an increased mapping rate (~1%) regardless of data source or cell type. Notably, large increases in mapping rate were observed for several epigenetic markers, including H3K9me3, H3K4me3, BRD3, CBP, and P300. Since T2T has many heterochromatin sequences not present in hg38,11 the increase in mapping rate for the heterochromatin marker H3K9me3 was considered reasonable. BRD3 is an essential actor of pericentric heterochromatin activation,34 whereas CBP and P300 have been reported to be located at (peri-)centromeres and involved in heterochromatin remodeling.35,36 Increased mapping rate was also apparent for the promoter marker H3K4me3, which is consistent with the fact that 1,956 new genes were introduced into the T2T reference genome.32 To further investigate the influence of T2T on read mapping, we analyzed the pairwise correlation of ChIP-seq read coverage between samples, and subtracted the correlation matrix of T2T from that of hg38. As illustrated in Fig. 2B and C, the Pearson correlation coefficients were minimally affected by T2T (less than ± 0.05), which is consistent with the slight (~1%) increases in mapping rate. Coincident with the results in Fig. 2A, several markers (H3K9me3, H3K4me3, BRD3, CBP, and P300) showed notably larger changes in the correlation analysis, indicating that these markers are more enriched in T2T-specific sequences.

Figure 2. Comparison of read mapping between T2T and hg38 assemblies. (A) Differences in the rate of uniquely mapped reads comparing T2T to hg38. Publicly available ChIP-seq datasets for MCF-7 and other cell types were obtained from Encode and GEO databases. (B, C) Differences in the Pearson correlation matrix between T2T and hg38 for the Encode MCF-7 datasets (B) or GEO MCF-7 datasets (C).

In summary, even the simple application of T2T to a typical short-read analysis workflow could provide additional epigenomic reads, especially for specific epigenomic markers.

3.3. T2T introduces additional epigenomic information

To intuitively investigate the additional reads mapped to T2T, we used Churros to visualize the read intensity of many epigenomic markers on chromosome 15, which is one of the most affected chromosomes, where T2T introduced nearly 20 million additional reference bases.32Fig. 3A shows the read distribution with hg38, where large gaps (purple dashed-line rectangle) are apparent due to the absence of reference sequences. Fig. 3B presents the read distribution with T2T. Importantly, the signal distribution was similar in the regions available in both T2T and hg38, suggesting that the use of T2T did not significantly affect the observations made with hg38. On the other hand, a clear enrichment of various epigenomic markers was observed in the newly introduced regions of T2T. For example, we observed the enrichment of H3K4me3 and H3K9me3, as well as the absence of the enhancer markers H3K27ac and H3K4me1, reflecting the existence of active gene promoters at centromeres11,37,38 (Fig. 3B, green dashed-line rectangle). The acetylated SMC3 (SMC3ac) at centromeres is required to establish the cohesiveness of chromatin-loaded cohesin and is maintained until the subsequent deacetylation.39–41 The green rectangle (bottom) in Fig. 3B shows the centromeric enrichment of SMC3ac in different cell types. In addition, on the short arms of the acrocentric chromosomes (SAAC), we observed many active histone markers, RNA polymerase II and various transcriptional factors (Fig. 3B, blue rectangle), suggesting the presence of ‘euchromatin-like’ sequences within SAAC.42,43Fig. 3C shows an enlarged region from Fig. 3B (indicated by the gray bar at the bottom), where ‘peak-like’ enrichment of ChIP-seq signals could be observed. This suggests that the additional reads mapped in T2T do not simply represent ‘noise’ but rather specific epigenomic features.

Figure 3. Visualization of read intensities for T2T and hg38. (A) Visualization of reads for various epigenomic markers on chromosome 15 using reference genome hg38. All reads were normalized relative to the whole genome. ChIP/Input enrichment is shown on a scale of 0–2. The cell type was MCF-7 unless otherwise indicated. (B) Visualization of reads for various markers on chromosome 15 using reference genome T2T. (C) Enlarged representation of the region in panel B denoted by the gray bar at the bottom. Genomic locations are also indicated at the bottom.

Collectively, the use of T2T with the current analytical workflow in Churros indeed yields new epigenomic information without compromising the original observations of hg38.

3.4. T2T reduces the ChIP-seq quality metric NSC

Cross-correlation analysis is commonly used for quality assessment of ChIP-seq data.20 We can compute a ‘strand-shift profile’ as the correlation coefficient between the read densities mapped onto the forward and reverse DNA strands, after shifting one strand by d base pairs. Owing to the fact that ChIP-seq sequence reads accumulate on the forward and reverse strands centered around the protein binding site, a strand-shift profile typically peaks at the shift distance corresponding to the DNA fragment length dflen. Based on this observation, the normalized strand coefficient (NSC) is defined as the ratio of the strand-shift profile at the fragment-length shift distance dflen to the background shift distance dbg. NSC is a widely used ChIP-seq quality metric. We therefore compared NSC values between T2T and hg38 to determine any impact of T2T on ChIP-seq quality assessments. As shown in Fig. 4A, NSC values were lower with T2T than with hg38 for almost all samples. Given that T2T introduced only a few additional reads in our analysis (~1% compared with hg38), the extent of the changes in NSC is notably large (~30%).

Figure 4. Influence of T2T on the quality metric NSC. (A) Diff_ratio represents the difference in NSC values when comparing T2T to hg38. The order of samples is the same as that shown in Fig. 2A. (B) Strand-shift profiles of three representative samples. The horizontal axis represents the shifted distance (bp), and the vertical axis represents the correlation coefficient. Estimated fragment lengths are indicated. (C) Scatter plot showing the relationship between the changes in mapping rate and the changes in NSC, comparing hg38 to T2T.

According to the definition described above,20NSC=J(dflen)/J(dbg), where J is the correlation coefficient (Jaccard index) for a given shift distance. We then asked which part of the NSC equation was most affected by T2T. First, the background shift distance dbg was not considered a possible factor because it has a fixed range of 500 kb to 1 Mb, with steps of 5 kb. Second, dflen represents the fragment length estimated from the strand-shift profile. Although several samples experienced a >50% change in NSC (Supplementary Fig. S1A), their corresponding dflen values were minimally affected by T2T (Supplementary Fig. S1B). Thus, dflen was unlikely to be the primary cause of the decreased NSC. Third, J(dflen) and J(dbg) are the correlation coefficients at the estimated fragment length and at the background distance, respectively. We plotted the strand-shift profiles, which represent the correlation coefficients for each strand-shift distance d (Fig. 4B). Despite a slight decrease in J(dbg), J(dflen) at the estimated fragment length decreased more substantially, suggesting that J(dflen) was the main factor contributing to the markedly reduced NSC. Considering that T2T introduced new reads mainly at specific regions, the decreased NSC could be attributed to the low enrichment of signals for the newly mapped reads at T2T-specific regions (e.g. heterochromatin bindings). Incorporating Figs. 2A and 4A, we noticed that samples that had a greater increase in mapping rate also exhibited a greater decrease in NSC. Indeed, the data in Fig. 4C illustrate the strong correlation between changes in mapping rate and changes in NSC (R2 = 0.3547, P < 0.0001).

These results suggest that using T2T in a standard ChIP-seq analysis will largely decrease the NSC metrics. Users need to consider adjusting the NSC threshold, or improving the algorithm to better suit T2T analysis.

3.5. T2T substantially affects MACS2 peak-calling results

Peak calling is one of the most critical steps in ChIP-seq analysis, and MACS2 is the most popular tool for it.19 Considering that T2T always results in a slight increase in the mapped reads without affecting regions that are also present in hg38, we hypothesized that the inclusion of the additional reads by T2T would result in slightly more identified peaks. Fig. 5A presents the changes in MACS2 peak numbers from hg38 to T2T. Unexpectedly, we observed both increases and decreases in peak numbers. Moreover, the magnitude of these changes was remarkable, with some samples showing fluctuations of over ±50%. To confirm whether such large percentage changes were limited to samples with few peaks, we examined Fig. 5B and observed that many samples with a large number of peaks also exhibited substantial changes in percentages. Given that the use of T2T did not affect the reads in hg38 regions (Fig. 3), the dramatic changes in peak number were likely due to the specific algorithm of MACS2.

Figure 5. Influence of T2T on MACS2 peak calling. (A) Diff_ratio represents the difference in the number of MACS2 peaks comparing T2T to hg38. (B) Scatter plot illustrating the relationship between diff_ratio and the numerical change in peak number. Each dot represents a ChIP-seq sample. (C) Association between changes in the number of peaks and changes in mapping rate. (D) Diff_ratio shows the differences in the estimated fragment length between T2T and hg38. Upper panel: estimated by MACS2; lower panel: estimated by SSP. (E) Correlation analysis between the diff_ratio of MACS2 peak numbers and the diff_ratio of estimated fragment lengths; ‘abs’ represents absolute value. (F) Estimated fragment lengths for representative samples comparing T2T to hg38. (G) MACS2 peak numbers for representative samples comparing T2T to hg38.

We then analyzed which factors contributed to the large difference in MACS2 peak number between T2T and hg38. Considering the variables in the MACS2 workflow, we first focused on the most important step, i.e. peak detection by Poisson distribution. The number of reads at each locus (k) can be modeled by the one-parameter Poisson distribution, i.e. Pλ(X=k)=λk/k!∗e−λ, where parameter λ is the expected number of reads: λ=(read length ×total read number)/effective genome length. Whereas read length and effective genome length are fixed, the total number of reads only increased by ~1% as shown above (Fig. 2). Such a small change was inconsistent with the large variations in MACS2 peak numbers. A scatter plot was generated to examine the relationship between changes in mapped reads and MACS2 peaks (Fig. 5C), but no correlation was evident as indicated by the values R2 = 0.002 and P = 0.51. Consequently, it was unlikely that the peak-detection step was the primary cause of the dramatic changes in the MACS2 peak calling results.

Another key step in MACS2 is modeling the fragment length before peak detection. Fig. 5D (upper panel) shows that T2T indeed moderately affected the estimated fragment length determined by MACS2. Because DNA fragmentation was performed during the preparation of the ChIP-seq library,3 the distribution of fragment lengths should be independent of genome assembly. In this regard, MACS2 is not adept at estimating fragment length in T2T, whereas other tools such as SSP achieved more robust results across assemblies (Fig. 5D, lower panel). After estimating the fragment length d, MACS2 shifts all read by d/2 toward the 3ʹ ends and then slides across the genome using a window size of 2 d to find candidate peaks. Hence, the impact of changes in d on peak calling will be highly intricate. For example, we simulated peak calling with different fragment lengths d (Supplementary Fig. S2A). Varying d from 100 to 200 bp resulted in both increases and decreases in the number of called peaks. Fig. 5E illustrates the relationship between changes in MACS2 peak numbers (absolute value) and changes in estimated fragment length (absolute value), revealing a positive correlation (R2 = 0.60, P < 0.0001). Example data showed that samples with large variations in MACS2 peak numbers also exhibited considerable differences in the estimated fragment length (Figs. 5F and G). These results indicated that the abnormality in the estimated fragment length d is the main contributing factor to the unexpectedly large changes in the MACS2 peak calling analysis with T2T. In addition, peak calling using other tools was also affected by T2T in different ways (Supplementary Fig. S2B). For example, from hg38 to T2T, the number of H3K9me3 peaks called by MACS2 increased greatly, whereas those called by Drompa34 decreased substantially. This result suggests that the impact of using T2T on peak calling depends on the specific algorithm.

Overall, the peak-calling results of MACS2 were substantially affected by T2T owing to variability in the estimated fragment length, highlighting the need to consider the peak-calling method in T2T analysis and underscoring the importance of developing a specialized algorithm for T2T.

3.6. Clustering of large-scale epigenomic profiles

Churros also includes many custom functions for epigenomic analysis. Here, we present a function that allows clustering and visualization of large-scale epigenomic profiles. Despite some tools, such as DeepTools,23 offering features to visualize several ChIP-seq samples, these features are not straightforward for large-scale epigenomics. Inspired by previous studies,44 Churros provides a ‘classheat’ function (classification and heatmap) for analyzing many epigenomic samples at given genomic regions. As inputs, classheat uses a file representing regions of interest and a directory containing multiple epigenomic signals (binary or continuous) (Fig. 6A); it then generates five output files, including the integrated raw matrix, processed matrix, sorted heatmap, clustered matrix, and clustered heatmap (also see Materials and methods). The classheat function is particularly useful for classifying specific regions into several different types and determining the epigenomic features enriched for each type. Moreover, because the regions of interest (input 1) can be arbitrary, such as specific chromatin structures (3D genome) or differentially expressed genes (transcriptome), the classheat function holds promise for multiomics integration.

Figure 6. A function to cluster and visualize large-scale epigenomic profiles. (A) Workflow of the classheat function. (B) Clustered heatmap with the continuous mode. The results show the enrichment of epigenomic features on different clusters of cohesin binding sites. Labels on the left represent clusters of genomic locations. Labels on the top represent user-defined groups. (C) Clustered heatmap with the binary mode. (D) Read distributions on the example regions show the promoter-like and enhancer-like cohesin sites.

To illustrate the effectiveness of classheat, we applied it to study the function of cohesin, a key protein complex for transcriptional regulation and chromatin structure.45,46 In the continuous model, the inputs included: (i) A bed file for all cohesin binding sites; and (ii) a directory of bigwig files for transcription factors, histone markers, RNA polymerase II, and CTCF. Supplementary Fig. S3A shows the sorted heatmap (output 3) based on the processed matrix (output 2). The noisy signals of heterochromatin markers served as the negative control. Consistent with previous research describing CTCF-cohesin and non-CTCF-cohesin,47,48 we observed that some cohesin sites colocalized with intense signals for CTCF and RNA polymerase II, whereas other cohesin sites lacking CTCF signals exhibited high enrichment of many cis-regulatory elements. However, these observations were vague and contained much noise. Fig. 6B shows the heatmap based on the clustered matrix (see Materials and methods), allowing us to observe the same biological meaning but with clear clusters and distinct feature enrichment. On the other hand, the binary model is sometimes more useful because continuous ChIP-seq signals can suffer from biases related to signal-to-noise, studies, or sequencing libraries.49 In the binary model, we used similar inputs as with the continuous model but utilized a folder of peak files instead of read coverages. From the sorted heatmap (Supplementary Fig. S3B), we could also observe the distinction between non-CTCF cohesin and CTCF cohesin. This sorted heatmap can only yield two clusters because the values of the sorted epigenomic features only have 0 or 1. Alternatively, more than two clusters could be obtained (Fig. 6C) because the clustered heatmap also considered colocalization with other factors. In addition to the CTCF cohesin sites, we further observed promoter-like cohesin sites, which co-bind with H3K4me3 and RNA polymerase II, and enhancer-like cohesin sites, which co-bind with H3K4me1, H3K27ac, and many transcription factors. Fig. 6D shows example regions of promoter-like and enhancer-like cohesin sites. These results demonstrated that classheat is effective for exploring biological significance from large-scale epigenomic profiles.

4. Discussion

While advances in sequencing technologies have made it possible to generate large-scale datasets, analyzing them poses substantial practical challenges. The rapidly increasing scale and diversity of epigenomic data have arisen in the need to develop portable, reproducible, and scalable analytic workflows. Although several computational frameworks such as Nextflow50 and snakemake51 also allow flexibility in data processing, these frameworks are still challenging for novice users, as they require training in specific syntax and programming languages; snakePipes simplifies the epigenomic analysis workflow by using a repository of modular rules,52 but users still need to tinker with the specific configurations of each tool; ChIP-AP supports multiple analysis steps from raw sequencing files, but its primary focus is on peak calling9; Most existing pipelines, including nf-core/ChIP-seq,53 do not include visualization and detailed sample comparison; AutoChIP and others54,55 provide automated ChIP-seq pipelines but lack downstream analysis and are not suitable for large datasets. Chilin56 can process many datasets, but it does not have sufficient optimization for hundreds of samples (e.g. only the alignment step supports parallel computation) and does not allow for much custom analysis. CoBRA57 uses Docker container to solve software dependencies, but it lacks support for comprehensive steps (e.g., alignment and quality check) and large-scale datasets.

In this regard, Churros constitutes an all-in-one epigenomic analysis platform with a complete analysis workflow, i.e. from read mapping to function annotations. The main advantages of Churros include pre-compiled computational environments, automated processes with customizable details, and optimizations for large-scale datasets. Churros can be distributed across computers or cloud servers as long as the machinery supports Docker. Churros is amenable to both biologists with no programming experience and bioinformatics experts who want to quickly interpret their sequencing data. Churros is a long-term maintenance project, and new functions will continue to be added over time.

Churros was constructed based on a combination of many published tools and our custom scripts. Although the main Churros workflow follows best practices, Churros also includes other tools that are pre-installed without modification, as these tools are already easy to use. For example, in Churros, users can use Homer21 to extract DNA motifs, and use MAnorm58 for differential peak analysis. In this sense, Churros serves as a platform that allows users to seamlessly implement many tools without the need to struggle with installation and configuration. Nevertheless, it is nearly impossible to encompass all the hundreds of tools, especially for complex downstream epigenomic analysis depending upon biological questions.52 Moreover, despite the development of next-generation sequencing-based epigenomics over the past 20 years, many challenges remain concerning data analysis, such as integrating epigenomic data with other omics data (3D genomics, transcriptomics, etc.) and exploring epigenetic heterogeneity using single-cell technologies. While developing new algorithms is beyond the scope of Churros, we will continue to follow cutting-edge tools and incorporate the useful ones into Churros.

We demonstrated the effectiveness of Churros by analyzing large-scale epigenomic datasets. In the comparison between T2T and hg38, we applied canonical analysis strategies to short-read data. Not surprisingly, we only observed approximately 1% more mapped reads in T2T, which differs from the fact that T2T introduces almost 8% new sequences.32 Despite this shortcoming, we still obtained additional epigenomic information with T2T (Figs. 2 and 3). In particular, our results underscore that users must heed caution regarding the influence of T2T on ChIP-seq quality metrics and peak calling. Given the massive amount of existing short-read epigenomic datasets, conducting exploratory analysis for reference genome T2T with Churros is highly valuable and practicable. Undoubtedly, the full power of the T2T assembly will be unleashed with the increasing availability of long-read datasets and specialized analysis methods. On the other hand, we introduced classheat for clustering and visualizing large-scale epigenomic profiles. Although classheat also supports other clustering methods, we only utilize principal component analysis followed by mini-batch k-means, which generate clear clusters. However, as epigenomic studies become increasingly complex, alternative clustering algorithms should be compared and considered to effectively address evolving challenges.

In summary, we developed a comprehensive epigenomic analysis pipeline, Churros. We demonstrated its validity via the analysis of real-world datasets. Churros will be of great benefit to all researchers exploring the vast ocean of epigenomics.

Supplementary Material

dsad026_suppl_Supplementary_Datas_4

dsad026_suppl_Supplementary_Datas_1

dsad026_suppl_Supplementary_Datas_2

dsad026_suppl_Supplementary_Datas_3

dsad026_suppl_Supplementary_Figures_S1-3

Acknowledgements

This work was supported by the Japan Agency for Medical Research and Development under grant number JP23gm6310012h0004.

Conflict of Interest

The authors declare no competing interests.

Author Contributions

R.N. and J.W. conceived the project. R.N. and J.W. wrote the code for Churros. J.W. performed the bioinformatic analysis. J.W. and R.N. wrote the manuscript. R.N. supervised the project. All authors read and approved the final manuscript.
==== Refs
References

1. Luo, Y., Hitz, B.C., Gabdank, I., et al. 2020, New developments on the encyclopedia of DNA elements (ENCODE) data portal, Nucleic Acids Res., 48 , D882–9.31713622
2. Schmidt, F., Marx, A., Baumgarten, N., et al. 2021, Integrative analysis of epigenetics data identifies gene-specific regulatory elements, Nucleic Acids Res., 49 , 10397–418.34508352
3. Nakato, R. and Shirahige, K. 2017, Recent advances in ChIP-seq analysis: from quality management to whole-genome annotation, Brief. Bioinform., 18 , 279–90.26979602
4. Nakato, R. and Sakata, T. 2021, Methods for ChIP-seq analysis: a practical workflow and advanced applications, Methods, 187 , 44–53.32240773
5. Buenrostro, J.D., Wu, B., Chang, H.Y., and Greenleaf, W.J. 2015, ATAC-seq: a method for assaying chromatin accessibility genome-wide, Curr. Protoc. Mol. Biol., 109 , 21 29 21–9.
6. Kaya-Okur, H.S., Wu, S.J., Codomo, C.A., et al. 2019, CUT&Tag for efficient epigenomic profiling of small samples and single cells, Nat. Commun., 10 , 1930.31036827
7. Consortium, E.P., Moore, J.E., Purcaro, M.J., et al. 2020, Expanded encyclopaedias of DNA elements in the human and mouse genomes, Nature, 583 , 699–710.32728249
8. Cell editorial, t. 2016, A cornucopia of advances in human epigenomics, Cell, 167 , 1139.27863229
9. Suryatenggara, J., Yong, K.J., Tenen, D.E., Tenen, D.G., and Bassal, M.A. 2022, ChIP-AP: an integrated analysis pipeline for unbiased ChIP-seq analysis, Brief. Bioinform., 23 , bbab537.34965583
10. Di Tommaso, P., Palumbo, E., Chatzou, M., Prieto, P., Heuer, M.L., and Notredame, C. 2015, The impact of Docker containers on the performance of genomic pipelines, PeerJ, 3 , e1273.26421241
11. Gershman, A., Sauria, M.E.G., Guitart, X., et al. 2022, Epigenetic patterns in a complete human genome, Science, 376 , eabj5089.35357915
12. Leinonen, R., Sugawara, H., and Shumway, M.; International Nucleotide Sequence Database Collaboration. 2011, The sequence read archive, Nucleic Acids Res., 39 , D19–21.21062823
13. Langmead, B. 2010, Aligning short sequencing reads with Bowtie, Curr. Protoc. Bioinformat., Chapter 11 , Unit 11 17.
14. Langmead, B. and Salzberg, S.L. 2012, Fast gapped-read alignment with Bowtie 2, Nat. Methods, 9 , 357–9.22388286
15. Li, H. and Durbin, R. 2009, Fast and accurate short read alignment with Burrows-Wheeler transform, Bioinformatics, 25 , 1754–60.19451168
16. Zhang, H., Song, L., Wang, X., et al. 2021, Fast alignment and preprocessing of chromatin profiles with Chromap, Nat. Commun., 12 , 6566.34772935
17. Chen, S., Zhou, Y., Chen, Y., and Gu, J. 2018, fastp: an ultra-fast all-in-one FASTQ preprocessor, Bioinformatics, 34 , i884–90.30423086
18. Martin, M. 2011, Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet, J., 17 , 10–2.
19. Feng, J., Liu, T., Qin, B., Zhang, Y., and Liu, X.S. 2012, Identifying ChIP-seq enrichment using MACS, Nat Protoc, 7 , 1728–40.22936215
20. Nakato, R. and Shirahige, K. 2018, Sensitive and robust assessment of ChIP-seq read distribution using a strand-shift profile, Bioinformatics, 34 , 2356–63.29528371
21. Heinz, S., Benner, C., Spann, N., et al. 2010, Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities, Mol. Cell, 38 , 576–89.20513432
22. Yu, G., Wang, L.G., and He, Q.Y. 2015, ChIPseeker: an R/Bioconductor package for ChIP peak annotation, comparison and visualization, Bioinformatics, 31 , 2382–3.25765347
23. Ramirez, F., Dundar, F., Diehl, S., Gruning, B.A., and Manke, T. 2014, deepTools: a flexible platform for exploring deep-sequencing data, Nucleic Acids Res., 42 , W187–191.24799436
24. Whyte, W.A., Orlando, D.A., Hnisz, D., et al. 2013, Master transcription factors and mediator establish super-enhancers at key cell identity genes, Cell, 153 , 307–19.23582322
25. Ernst, J. and Kellis, M. 2015, Large-scale imputation of epigenomic datasets for systematic annotation of diverse human tissues, Nat. Biotechnol., 33 , 364–76.25690853
26. Krueger, F. and Andrews, S.R. 2011, Bismark: a flexible aligner and methylation caller for bisulfite-Seq applications, Bioinformatics, 27 , 1571–2.21493656
27. Li, H., Handsaker, B., Wysoker, A., et al. ; 1000 Genome Project Data Processing Subgroup. 2009, The sequence alignment/map format and SAMtools, Bioinformatics, 25 , 2078–9.19505943
28. Quinlan, A.R. and Hall, I.M. 2010, BEDTools: a flexible suite of utilities for comparing genomic features, Bioinformatics, 26 , 841–2.20110278
29. Ewels, P., Magnusson, M., Lundin, S., and Kaller, M. 2016, MultiQC: summarize analysis results for multiple tools and samples in a single report, Bioinformatics, 32 , 3047–8.27312411
30. Kurtzer, G.M., Sochat, V., and Bauer, M.W. 2017, Singularity: scientific containers for mobility of compute, PLoS One, 12 , e0177459.28494014
31. Pabinger, S., Dander, A., Fischer, M., et al. 2014, A survey of tools for variant analysis of next-generation genome sequencing data, Brief. Bioinform., 15 , 256–78.23341494
32. Nurk, S., Koren, S., Rhie, A., et al. 2022, The complete sequence of a human genome, Science, 376 , 44–53.35357919
33. Hoyt, S.J., Storer, J.M., Hartley, G.A., et al. 2022, From telomere to telomere: the transcriptional and epigenetic state of human repeat elements, Science, 376 , eabk3112.35357925
34. Col, E., Hoghoughi, N., Dufour, S., et al. 2017, Bromodomain factors of BET family are new essential actors of pericentric heterochromatin transcriptional activation in response to heat shock, Sci. Rep., 7 , 5418.28710461
35. Piacentini, L., Marchetti, M., Bucciarelli, E., et al. 2019, A role of the Trx-G complex in Cid/CENP-A deposition at Drosophila melanogaster centromeres, Chromosoma, 128 , 503–20.31203392
36. Naughton, C., Huidobro, C., Catacchio, C.R., et al. 2022, Human centromere repositioning activates transcription and opens chromatin fibre structure, Nat. Commun., 13 , 5609.36153345
37. Mellor, J., Dudek, P., and Clynes, D. 2008, A glimpse into the epigenetic landscape of gene regulation, Curr. Opin Genet. Dev., 18 , 116–22.18295475
38. Morrison, O. and Thakur, J. 2021, Molecular complexes at euchromatin, heterochromatin and centromeric chromatin, Int. J. Mol. Sci. , 22 , 6922.34203193
39. Hoencamp, C. and Rowland, B.D. 2023, Genome control by SMC complexes, Nat. Rev. Mol. Cell Biol., 24 , 633–50.37231112
40. Dauban, L., Montagne, R., Thierry, A., et al. 2020, Regulation of cohesin-mediated chromosome folding by eco1 and other partners, Mol. Cell, 77 , 1279–1293 e1274.32032532
41. Moronta-Gines, M., van Staveren, T.R.H., and Wendt, K.S. 2019, One ring to bind them—cohesin’s interaction with chromatin fibers, Essays Biochem., 63 , 167–76.31015387
42. Lyle, R., Prandini, P., Osoegawa, K., et al. 2007, Islands of euchromatin-like sequence and expressed polymorphic sequences within the short arm of human chromosome 21, Genome Res., 17 , 1690–6.17895424
43. Antonarakis, S.E. 2022, Short arms of human acrocentric chromosomes and the completion of the human genome sequence, Genome Res., 32 , 599–607.35361624
44. Faure, A.J., Schmidt, D., Watt, S., et al. 2012, Cohesin regulates tissue-specific expression by stabilizing highly occupied cis-regulatory modules, Genome Res., 22 , 2163–75.22780989
45. Wang, J. and Nakato, R. 2023, CohesinDB: a comprehensive database for decoding cohesin-related epigenomes, 3D genomes and transcriptomes in human cells, Nucleic Acids Res., 51 , D70–9.36162821
46. Wang, J. and Nakato, R. 2023, Comprehensive multiomics analyses reveal pervasive involvement of aberrant cohesin binding in transcriptional and chromosomal disorder of cancer cells, iScience, 26 , 106908.37283809
47. Schmidt, D., Schwalie, P.C., Ross-Innes, C.S., et al. 2010, A CTCF-independent role for cohesin in tissue-specific transcription, Genome Res., 20 , 578–88.20219941
48. Wang, J., Bando, M., Shirahige, K., and Nakato, R. 2022, Large-scale multi-omics analysis suggests specific roles for intragenic cohesin in transcriptional regulation, Nat. Commun., 13 , 3218.35680859
49. Zang, C., Schones, D.E., Zeng, C., Cui, K., Zhao, K., and Peng, W. 2009, A clustering approach for identification of enriched domains from histone modification ChIP-Seq data, Bioinformatics, 25 , 1952–8.19505939
50. Di Tommaso, P., Chatzou, M., Floden, E.W., Barja, P.P., Palumbo, E., and Notredame, C. 2017, Nextflow enables reproducible computational workflows, Nat. Biotechnol., 35 , 316–9.28398311
51. Koster, J. and Rahmann, S. 2018, Snakemake—a scalable bioinformatics workflow engine, Bioinformatics, 34 , 3600.29788404
52. Bhardwaj, V., Heyne, S., Sikora, K., et al. 2019, SnakePipes: facilitating flexible, scalable and integrative epigenomic analysis, Bioinformatics, 35 , 4757–9.31134269
53. Ewels, P.A., Peltzer, A., Fillinger, S., et al. 2020, The nf-core framework for community-curated bioinformatics pipelines, Nat. Biotechnol., 38 , 276–8.32055031
54. Kim, T., Lee, W., Han, K., and Kang, K. 2015, An automated analysis pipeline for a large set of ChIP-seq data: AutoChIP, Genes & Genomics, 37 , 305–11.
55. Park, S.J., Kim, J.H., Yoon, B.H., and Kim, S.Y. 2017, A ChIP-seq data analysis pipeline based on bioconductor packages, Genomics Inform, 15 , 11–8.28416945
56. Qin, Q., Mei, S., Wu, Q., et al. 2016, ChiLin: a comprehensive ChIP-seq and DNase-seq quality control and analysis pipeline, BMC Bioinf., 17 , 404.
57. Qiu, X., Feit, A.S., Feiglin, A., et al. 2021, CoBRA: containerized bioinformatics workflow for reproducible ChIP/ATAC-seq analysis, Genomics Proteomics Bioinformatics, 19 , 652–61.34284136
58. Shao, Z., Zhang, Y., Yuan, G.C., Orkin, S.H., and Waxman, D.J. 2012, MAnorm: a robust model for quantitative comparison of ChIP-Seq data sets, Genome Biol., 13 , R16.22424423
