A generic reference defined by consensus peaks for single-cell ATAC-seq data analysis

Domain: neuroscience
Published: 2026-02-09
Editorial outcome: major_revision

Abstract
The rapid advancement of transposase-accessible chromatin using sequencing (ATAC-seq) technology, particularly with the emergence of single-cell ATAC-seq (scATAC-seq), accelerates the studies of gene regulation. However, the absence of a generic feature reference for ATAC-seq data limits single-cell analyses and hinders the development of comprehensive cell atlases. To address this, we construct a generic chromatin accessibility reference by aggregating peaks from 624 high-quality bulk ATAC-seq datasets, defining about 1.4 million consensus peaks (cPeaks). Leveraging a deep neural network model, we expand cPeaks to include previously unobserved regions, enhancing their coverage across diverse tissues and cell types. cPeaks exhibit consistent shapes across tissue types, sequencing technologies, and peak-calling methods, indicating that they represent inherent genomic features. Compared to existing feature-defining methods and references, cPeaks show superior performance in scATAC-seq analyses, improving cell annotation and rare cell type identification. Additionally, cPeaks provide insights into chromatin dynamics during cellular differentiation and tumor progression. cPeaks can serve as a robust reference for chromatin accessibility studies to promote cross-dataset consistency and accelerate biological discoveries.

Introduction
Chromatin accessibility refers to the availability of chromatinized DNA to directly interact with other macromolecules1. The accessible regions2,3, also known as open regions1,2, often associate with regulatory elements such as promoters, enhancers, insulators, and transcription factor (TF) binding sites (TFBSs)4,5. Studying accessible regions is crucial for deciphering gene regulations2,6–8 across various biological processes, including cell differentiation9,10, organ development11,12, and disease manifestation13–16. The advancements in high-throughput sequencing technologies3, particularly Assay for Transposase-Accessible Chromatin using sequencing (ATAC-seq)17,18, have revolutionized our ability to measure chromatin accessibility on a genome-wide scale by offering more efficient and simplified approaches. Especially, the emergence of single-cell ATAC-seq (scATAC-seq)19–24 now allows for the exploration of chromatin accessibility at single-cell resolution, complementing single-cell RNA sequencing (scRNA-seq) for investigating the molecular properties of cells of different types or states and for building the biomolecular atlases of cells5,16,25,26.

Unlike RNA-seq, which benefits from established reference transcriptomes27, ATAC-seq lacks a standardized reference for localizing and quantifying chromatin accessibility across the genome. Traditional approaches in bulk ATAC-seq involve creating sample-specific feature sets, using peak-calling algorithms such as MACS228 and Cell Ranger ATAC29 to identify regions with significantly enriched read signals as accessible regions or peaks5,25,30. When integrating data from different samples, current practices usually first call peaks as features in each dataset, and then integrate these features with the peak merging31,32 or iterative removal5,16 approach. The merging approach combines overlapping peaks into broader regions, potentially merging adjacent peaks inaccurately. The iterative removal approach selects the most significant peak from overlapping areas, risking the loss of valuable information. Recently, an alternative approach avoids peak calling by tiling the genome into fixed-width bins (e.g., 500 bp) for accessibility comparison across cells and samples33. This strategy avoids the complexity of merging non-aligned peaks and enables easier comparison across datasets. While it has proven effective for analyzing gene regulatory programs in large cohorts, using uniform bins may overlook true biological peak boundaries and reduce resolution for detecting fine regulatory elements.

The advent of scATAC-seq provides a single-cell resolution of chromatin accessibility information and brings substantial increases in the scale of the data, often comprising tens to hundreds of thousands of cells5,26,30. However, peak-calling at the single-cell level is impractical, necessitating the merging of multiple cells into pseudo-bulk samples to define accessible regions5,25,26,30,34,35. This results in the loss of single-cell resolution and potentially obscures unique chromatin features of rare cell types. To improve this, iterative techniques for clustering5 and optimized latent semantic indexing (LSI)36 have been developed to improve the identification of accessible regions, especially in rare cell types. Nevertheless, these techniques still involve multiple rounds of feature integration, which increases with scATAC-seq data scales, leading to labor-intensive processes and aggravated information loss, especially in atlas-scale analyses.

These challenges highlight the critical need for a generic reference of accessible regions to standardize analyses across the increasing volume of sequencing data, and to enhance the consistency and accuracy of research. Valentina et al. have shown that using predefined sets of genomic regions can facilitate the mapping of scATAC-seq data37. This underscores the potential transformative impact of a comprehensive reference that encompasses all or at least the major potential accessible regions shared across datasets on scATAC-seq data analysis. Currently, several pre-defined sets of accessible regions have been proposed to identify and characterize possible gene regulatory elements5,38,39. For example, the consensus DNase I hypersensitive site (cDHS)38 offers a standard index of DNase I hypersensitive sites identified through DNase I hypersensitive seq (DNase-seq) data, and the cis-element ATLAS (CATLAS)5 was introduced as a repository of candidate cis-regulatory elements (cCREs) across the genome. cDHS is specifically based on DNase-seq technology, and CATLAS cCREs do not retain information on the length of accessible regions. Consequently, these resources are often employed for quality control in data analysis rather than as comprehensive references for data analysis3.

To construct a robust chromatin accessibility reference, we collect 624 high-quality public bulk ATAC-seq datasets, covering diverse tissue types. These datasets are used to define approximately 1.4 million observed consensus peaks (cPeaks) as a chromatin accessibility reference. To evaluate the generalizability of cPeaks as a universal feature set for scATAC-seq data, we validate them using an independent set of 231 scATAC-seq datasets. We find consistent patterns in cPeak shapes across tissue types, sequencing technologies, and peak-calling methods, suggesting that cPeaks represent inherent genomic features. To address the limitation of unseen tissue-specific accessible regions in existing data, we develop a deep neural network model to predict new cPeaks, adding ~0.28 million predicted regions to the reference. We use an additional independent set of 789 bulk ATAC-seq datasets to validate the reliability of predicted cPeaks. The combined observed and predicted cPeaks serve as a standardized set of chromatin accessibility references, enhancing cross-dataset consistency in ATAC-seq analyses. Compared to existing feature-defining methods and references, cPeaks significantly improve scATAC-seq cell annotation and rare cell type identification. Additionally, cPeaks facilitate the exploration of chromatin dynamics in complex biological states such as cellular differentiation and tumor progression. We provide comprehensive resources and bioinformatics tools to integrate cPeaks into multiple analysis pipelines to enable their effective use in all types of chromatin accessibility studies.

Results
Defining cPeaks by aggregating ATAC-seq data
We used bulk ATAC-seq data to establish a robust chromatin accessibility reference. We collected all available bulk ATAC-seq datasets derived from healthy adult samples from two widely used repositories, ENCODE40 and ATACdb41. After quality control, the collection resulted in a total of 624 high-quality datasets, covering over 40 human organs and all body systems (Methods, Supplementary Data 1, Supplementary Fig. 1a). Each dataset provided a set of pre-identified peaks, each marking an accessible region on the genome observed in respective dataset (Fig. 1a). Peaks were iteratively merged, one dataset at a time, with non-overlapping peaks directly added to the pool and overlapping peaks combined into broader regions (Methods). As more datasets were integrated, the total number of peaks and their genomic coverage initially grew rapidly, followed by a gradual deceleration in growth (Fig. 1b), suggesting that the total number of all possible accessible regions on the genome may be limited.

After processing all datasets, we identified approximately 1.2 million merged peaks, covering about 30% of the whole genome. Each merged peak represents a genomic region that could be accessible in one or more cell types. For each position within a merged peak, we calculated the proportion of datasets in which the position was covered by peaks, defining this proportion as the “accessible score” for the position. The series of accessible scores along the 5’-to-3’ direction of a merged peak was summarized into a profile referred to as the “shape” of the peak (Methods).

We then investigated whether the shapes of merged peaks remained consistent under varying conditions. As an internal evaluation, we first grouped the datasets used to construct merged peaks based on tissue origin and source databases. About 51% of the datasets came from blood samples, while the rest were from solid tissues such as lung, heart, and brain. We applied the same merging procedure to each group independently, and the majority of peaks overlapped across blood and solid tissues (87% and 76%, respectively) showed highly similar shapes, as illustrated in Fig. 1c. We quantified shape similarity using the Pearson correlation coefficient (PCC; Methods, Supplementary Fig. 1b), yielding a median PCC of 0.916 (Supplementary Fig. 1c). To assess the sensitivity and significance of this similarity, we compared the PCCs against two baselines: shuffled pairings and 250 bp–shifted pairs (Methods, Supplementary Fig. 1b). As expected, the paired peaks showed substantially higher PCCs than either baseline: 61.2% had empirical p < 0.05 compared to shuffled controls, and 77.7% compared to shifted controls (Methods, Supplementary Fig. 1d). These results support the robustness and biological relevance of peak shape consistency across tissue types.

We then evaluated the robustness of peak shapes across data sources and technologies. First, we generated merged peaks separately from ENCODE and ATACdb datasets, which had been processed using MACS228 with different peak-calling parameters (Methods). Once again, we observed a high level of shape similarity (median PCC = 0.868, Supplementary Fig. 1c), with 25.3% of pairs significant against shuffled controls, and 93.1% against shifted controls (Supplementary Fig. 1d). We then compared merged peaks from ATAC-seq (ENCODE) and DNase-seq (cDHS)38 datasets (Methods). Despite differences in assay principles and processing pipelines, we observed high shape similarity at corresponding genomic regions (median PCC = 0.902, Supplementary Fig. 1c, d). Since both ENCODE and ATACdb, datasets used for generating merged peaks in this study were analyzed by MACS2, we specifically compared peaks generated by all ENCODE datasets with those generated by CellRanger29 from a 10x Genomics scATAC-seq dataset (Methods). We clustered the single cells into 17 groups, derived pseudo-bulk peaks for each group, and then constructed merged peaks using the same procedure. Due to the limited number of cell types, the resulting merged peaks had relatively low resolution, which likely contributed to the lower observed PCC (median = 0.665). Nevertheless, the shape similarity remained significantly higher than both the shuffled and shifted controls (19.3% and 22.9% of pairs significant against shuffled controls and shifted controls, respectively), supporting the robustness of peak shapes (Supplementary Fig. 1c, d).

Together, these findings highlight the consistent shapes of merged peaks, regardless of tissue type, sequencing technology, or peak-calling methodology. This suggests that merged peaks may represent an inherent genomic feature and can serve as a basis for generating a chromatin accessibility reference. These findings also suggest that the generation of cPeaks is largely unaffected by batch effects, including those arising from differences in sample origin, sequencing technology, or computational processing. This robustness likely stems from the design of our approach: cPeaks are constructed from pre-identified peak sets rather than raw sequencing reads, allowing much of the technical noise to be filtered out during initial peak calling.

We found that merged peaks exhibited two distinct shapes: unimodal or multimodal (Fig. 1a, Methods). We found that multimodal merged peaks were composed of adjacent unimodal peaks and therefore we decomposed them into multiple unimodal peaks (Fig. 1a, Methods). These unimodal merged peaks were designated as the basic unit of the chromatin accessibility reference, termed observed consensus peaks (cPeaks). In this way, we obtained approximately 1.4 million observed cPeaks in total (an example provided in Supplementary Fig. 1e). The majority of observed cPeaks were shorter than 2000 bp, with a median length of 525 bp, closely aligning with previous findings that ATAC-seq peaks typically span around 500 bp42 (Supplementary Fig. 1f). To assess the robustness of the final output, we also evaluated shape consistency on the set of observed cPeaks (Methods). This analysis confirmed that cPeaks retained high shape similarity across different conditions (Fig. 1d and Supplementary Fig. 1g), consistent with the trends observed in the merged peaks.

We assessed the feasibility of cPeaks as a chromatin accessibility reference for processing scATAC-seq data. To ensure broad validation, we collected 231 scATAC-seq datasets from published studies (Supplementary Data 2, Supplementary Fig. 1a), each containing peaks pre-identified using different peak-calling methods28,29. For each dataset, we calculated the proportion of the pre-identified peaks overlapping with observed cPeaks, referred to as cPeak coverage. Most datasets exhibited coverage above 0.9, with an average of 0.98 and a median of 0.99 (Fig. 1e). Observed cPeaks consistently maintained high coverage across different tissue types, pre-identified peak numbers, and peak-calling parameters without clear bias (Supplementary Fig. 1h).

To further evaluate the enrichment of sequencing reads on observed cPeaks, we used the 27 scATAC-seq datasets from CATLAS5 to investigate cPeaks’ performance across different organs. We found that the sequencing reads were predominantly enriched near the centers of observed cPeaks (Fig. 1f, Methods). On average, 68.8% of sequencing reads in these scATAC-seq data were mapped to observed cPeaks, surpassing the typical mapping rate of ~65% reported in the technical reports of 10x Genomics (Methods). These findings establish cPeaks as a reliable reference for capturing the signals in scATAC-seq data across diverse organs.

We compared the performance of cPeaks against other pre-defined accessible region sets, including cDHS and CATLAS cCREs of human tissues. First, we assessed the differences in coverage between cPeaks and cDHS on accessible regions identified by ATAC-seq and DNase-seq. To do this, we constructed three unbiased pairs of ATAC-seq and DNase-seq datasets from the GM12878 cell line using data in ChIP-Atlas43 (Methods). We found that compared to cDHS, cPeaks provided higher coverage of ATAC-seq-specific accessible regions while maintaining similar coverage for DNase-seq-specific regions in these datasets (Fig. 1g, Methods). This demonstrates that cPeaks encompass a broader range of accessible regions detected by ATAC-seq compared to cDHS.

We next evaluated how well different pre-defined accessible region sets capture chromatin accessibility signals in real scATAC-seq data, to assess whether these region sets could serve as effective references for scATAC-seq analysis. To do so, we used pseudo-bulk peaks identified by the commonly used MACS228 pipeline as a practical baseline, reflecting common analysis practices in current workflows. Using a hematopoietic differentiation dataset9, we mapped sequencing reads to each region set and sorted accessible regions by the number of mapped reads. As shown in Fig. 1h, cPeaks captured a comparable proportion of sequencing reads to pseudo-bulk peaks, and substantially outperformed both cDHSs and CATLAS cCREs. We note that cDHSs were originally generated using Hotspot2, while cPeaks and cCREs were called with MACS2, which may introduce minor tool-related differences. However, as MACS2 is commonly adopted as the default peak-calling tool for ATAC-seq in community pipelines32, this result still suggests that cPeaks may serve as a more comprehensive and appropriate reference for scATAC-seq data analysis.

Expanding cPeaks through deep learning
cPeaks were derived from a finite set of ATAC-seq datasets. Despite the considerable sample size, concerns remain about their ability to handle unseen cell types with unique accessible regions. To address this limitation and enhance cPeaks as a universal chromatin accessibility reference, we employed a deep learning model to expand cPeaks.

We hypothesized that cPeak positions are related to their DNA sequences. This assumption is supported by the non-random distribution of cPeaks across the genome (Fig. 1b) and their highly consistent shapes regardless of tissue type, sequencing technology, or peak-calling methodology. These observations suggest a sequence-dependent mechanism underlying cPeak positioning.

We thus adapted the architecture of scBasset44, a one-dimensional convolutional neural network (CNN) originally designed to predict chromatin accessibility in a single-cell context using DNA sequence information. However, our objective differs from the original design: rather than predicting cell-type-specific accessibility profiles, we aim to identify genomic regions with an intrinsic potential for accessibility, independent of tissue or cell identity. Accordingly, we modified the output layer of the model to perform binary classification, assigning each input sequence to either cPeaks or non-cPeaks (Fig. 2a). Observed cPeaks served as positive samples, while randomly selected regions without cPeaks were used as negative controls (Methods).

The model achieved a high accuracy with an AUC of 94.72% on test data (Fig. 2b), highlighting a strong relationship between cPeak locations and their DNA sequence signatures. These findings reinforce the hypothesis that cPeaks are an inherent genomic feature and demonstrate the feasibility of predicting unobserved cPeaks using DNA sequence information.

We evaluated whether the CNN model could predict cPeaks specific to cell types or tissues that were not included during model training. Each cPeak in our reference is annotated with the cell types or tissues in which it is accessible, and we defined cell-type- or tissue-specific cPeaks as those observed in only one context. To assess generalization, we designed a leave-one-context-out strategy, where all cPeaks specific to a given cell type or tissue were removed from the training data and reserved as a held-out test set. The remaining cPeaks were used as positive training examples, and negative samples were randomly selected from genomic regions without cPeaks, following the same procedure illustrated in Fig. 2a. The CNN model was then trained from scratch on each revised dataset (Methods).

As shown in Fig. 2c, the model consistently achieved high prediction scores for distinguishing excluded cell-type/tissue-specific cPeaks from an equal number of negative samples, achieving a mean AUC of 91.67% across all evaluations. For instance, heart-specific cPeaks, which represent 0.97% of all cPeaks, were effectively identified with an AUC of 94.70% after training the model on the remaining 99.03% of cPeaks. These results indicate that the CNN model can effectively recover cPeaks that were not used during training, including those identified in tissues or cell types excluded from the training set.

We then assessed the robustness of the model using a leave-one-chromosome-out evaluation (Fig. 2d). In each round, the model was trained on all chromosomes except one, which was held out as an independent test set. The model consistently achieved high AUC scores (~0.95) on most held-out chromosomes, closely matching its performance on the original test set derived from the remaining chromosomes. When chromosome Y (chrY) was held out, the AUC dropped slightly to ~0.90, suggesting that the relationship between DNA sequence and chromatin accessibility on chrY may differ from that on other chromosomes. This observation is consistent with prior reports indicating that chrY harbors a complex repeat structure compared to other chromosomes45,46.

We further employed NeuronMotif47 to interpret the learned features of the CNN model (Methods). By analyzing the last convolutional layer, we identified 342 TF motifs that could robustly activate at least one neuron. Notably, motifs from the bZIP family and CTCF/CTCFL were found to be highly influential across many neurons (Supplementary Fig. 1i). These findings provide additional confidence that the model captures generalizable regulatory patterns from DNA sequence, rather than relying on dataset-specific artifacts.

To expand the cPeak set, we used the trained CNN model to scan the entire human genome, applying a sliding window of 500 bp across all background regions which are longer than 100 bp. This approach assigned an accessibility score to each genomic region, representing its probability of being a cPeak (Fig. 2e, Methods). Regions were ranked by accessibility scores, with those scoring highest designated as predicted cPeaks, resulting in a total of approximately 0.28 million predicted cPeaks (Methods). We provided the complete list of cPeaks and their related information (see Code Availability for details). In all downstream analyses, we consistently used the full set of cPeaks as the reference.

To evaluate the performance of predicted cPeaks, we collected an independent set of bulk ATAC-seq datasets from ChIP-Atlas43. After quality control, we obtained 789 datasets covering eight cell types or tissues: B cells, T cells, fibroblasts, brain, hippocampus, liver, lung, and retina (Supplementary Data 2, Supplementary Fig. 1a). These datasets were chosen to include both tissue types represented in the initial cPeak construction and a previously unseen tissue, retina. Importantly, none of these datasets were used in the initial construction process, ensuring an unbiased evaluation of predicted part. We also randomly generated a set of background regions with lengths and numbers matching the predicted cPeaks, referred to as random regions, to serve as random controls for comparison.

We evaluated the coverage of observed cPeaks alone, observed cPeaks plus random regions, and observed cPeaks plus predicted cPeaks across the 789 datasets. We found that observed cPeaks already covered more than 99.21% of the reported peaks for most datasets. For retina, the mean coverage of observed cPeaks was slightly lower at 92.23%. Adding random regions to observed cPeaks did not significantly improve coverage, whereas incorporating predicted cPeaks significantly increased coverage across most cell types and tissues (Fig. 2f). In previously observed tissues or cell types, datasets with initially lower observed coverage exhibited the greatest improvement, reducing coverage variance among samples after the inclusion of predicted cPeaks. For the unseen tissue, retina, the model demonstrated the most notable enhancement, with mean coverage increasing from 92.23% to 93.95% upon adding predicted cPeaks. These results indicate that the inclusion of predicted cPeaks enhances the ability of cPeaks to encompass accessible regions across diverse datasets, particularly benefiting cell types or tissues not initially represented in cPeak construction.

Association between cPeaks with known regulatory elements and cellular functions
We explored the associations between cPeaks and known regulatory elements and biological functions. We first validated the regulatory relevance of cPeaks using ChIP-seq data from ReMap 202248 (Methods). Compared to cDHSs and CATLAS cCREs, cPeaks showed the highest coverage of ReMap peaks (88.7%) and the strongest enrichment over matched background regions (13.5-fold) (Supplementary Fig. 2a). At the individual TF level, cPeaks also consistently outperformed CATLAS cCREs and slightly exceeded cDHSs, confirming their robustness in capturing regulatory elements (Supplementary Fig. 2b). Footprint analysis also confirmed that cPeaks cover abundant TF binding sites and can support footprint-based analyses when base-resolution data are available (Supplementary Fig. 2c, Methods).

By comparing observed cPeaks with genomic annotations, we found that observed cPeaks are enriched near transcriptional start sites (TSSs) and transcription end sites (TESs) (Fig. 3a, Methods). Approximately 13% of observed cPeaks are localized at promoter regions (±2000 bp around TSS), while the remaining cPeaks are primarily distributed at exons or distal intergenic regions (Fig. 3b, Methods). We further collected multiple types of epigenetic annotations, including histone modification and CTCF ChIP-seq data (Supplementary Data 3). Using these annotations, we found that most observed cPeaks contain or overlap with one or more regulatory elements, such as enhancers, promoters (both known and predicted), untranslated regions (UTRs), or CTCF binding sites (CTCFBSs) (Fig. 3c, Methods). Further genomic annotations and regulatory role inferences for predicted cPeaks revealed properties similar to those of observed cPeaks (Fig. 3b, c).

For each cPeak, we defined the maximum value of its shape as its “accessible ratio” to reflect the proportion of datasets where the cPeak is accessible. We defined a cPeak with an accessible ratio > 90% as a “housekeeping” cPeak (Fig. 3d). Among the observed cPeaks, 24,538 (1.8%) were identified as housekeeping cPeaks. The median length of housekeeping cPeaks was 1441 bp, with most located at promoter regions (Fig. 3e). These promoter regions overlap with 77.4% of the promoters of housekeeping genes49 (Fisher’s exact test p-value < 1 × 10−5, Methods). Besides, the proportion of housekeeping cPeaks and the proportion of housekeeping genes across chromosomes were strongly correlated (Supplementary Fig. 2d). Gene Ontology (GO) enrichment analysis of genes linked to housekeeping cPeaks revealed significant enrichment in essential biological pathways (Methods), including ribonucleoprotein complex biogenesis, cell cycle checkpoint signaling, and methylation (Fig. 3f). These findings align with those of a recently published study50 which proposed a similar concept Housekeeping CREs (HK-CREs) and reported that most HK-CREs are located in core promoters and associated with housekeeping gene-related biological processes. We found that housekeeping cPeaks encompass 98.4% of HK-CREs, underscoring their functional relevance (Fig. 3g).

We further examined the relationship between cPeak shapes and gene regulation. We defined the genomic region corresponding to the maximum value of each cPeak shape as its core accessible region. The distance from the endpoints of cPeak to the nearest boundary of the core accessible region is defined as the edge length. This metric reflects the consistency of accessible region boundaries across different cells. Smaller edge lengths indicated more consistent boundaries, suggesting the presence of well-positioned flanking nucleosomes, while larger edge lengths implied weakly-positioned flanking nucleosomes1,51,52. To ensure reliability, we filtered cPeaks with an accessible ratio greater than 0.3. As shown in Fig. 3h, the distribution of edge lengths exhibited three clusters with boundaries at 170 bp and 340 bp, approximately corresponding to the lengths of DNA wrapped around one and two nucleosomes, respectively. Edges shorter than 170 bp were categorized as well-positioned, while longer edges were classified as weakly-positioned. Based on edge type, cPeaks were further divided into three patterns: “well-positioned,” “asymmetrically-positioned,” and “weakly-positioned” (Fig. 3i).

We focused on well-positioned edges as they represent more precisely regulated regions. Two potential biological mechanisms were hypothesized to underlie the occurrence of well-positioned edges: the presence of specific TFBSs at the inner boundary and well-positioned flanking nucleosomes outside the boundary. To investigate the first hypothesis, we conducted motif enrichment analysis around well-positioned edges to identify well-positioned-associated TFs. These identified TFs included CTCF, BORIS, Fosl2/Fra-2, Fli1, and ETS1 (Supplementary Fig. 2e, Methods). The top 20 well-positioned-associated TFs shared common DNA-binding domains, including the basic leucine zipper (bZIP), Zinc finger (Zf), and erythroblast transformation specific (ETS) domains (Fig. 3j). Beyond the top 20, all well-positioned-associated TFs were also enriched for additional domains, such as homeobox and basic helix-loop-helix (bHLH) (Supplementary Fig. 2f). The GO enrichment analysis result showed that these TFs were related to essential biological processes, including the regulation of transcription, cell development, and cell differentiation (Supplementary Fig. 2g, Methods). We assumed that these well-positioned-associated TFs with these domains may possess stronger nucleosome remodeling capabilities. Besides, we found pioneer transcription factors (PTFs)53–55 were enriched in these TFs with a p-value of 0.016 (Fisher’s exact test, Methods). For the second hypothesis, we performed de novo motif enrichment analysis for the sequence of well-positioned edges, identifying sequence motifs associated with well-positioned flanking nucleosome positioning around cPeak boundaries (Fig. 3k, Methods). Together, well-positioned edges likely arise from the interplay of specific TFBSs and well-positioned flanking nucleosomes, highlighting their role in precise genomic regulation.

Most cPeaks are accessible in only a limited number of tissues or cell types (Fig. 3d), leaving 92.62% of cPeaks with unlabeled shapes due to low accessible ratios. To address this, we developed a neural network, a modified version of the CNN model used for cPeak prediction, to capture the DNA sequencing patterns corresponding to specific cPeak shapes (Supplementary Fig. 2h, Methods). The network independently predicts whether each edge of a cPeak is well-positioned and achieved an AUC of 86.56% on the test set when trained on labeled cPeaks. We applied the model to cPeaks with unlabeled shapes and used the newly labeled, well-positioned cPeaks for TF enrichment analysis. The enriched TF results from this analysis strongly aligned with those obtained from previously labeled cPeaks (Spearman correlation coefficient = 0.76) (Supplementary Fig. 2i), demonstrating the network’s capability to classify cPeak shapes based on DNA sequence patterns. Incorporating these predictions, we classified 312,497 (18.86%) as “well-positioned”, 366,716 (22.13%) as “asymmetrically-positioned”, and 977,981 (59.01%) as “weakly-positioned”. These annotations were integrated into the cPeak resource to support further analyses and applications.

Enhancing cell type annotation of scATAC-seq data with cPeaks
We evaluated the effectiveness of using cPeaks as a reference for annotating scATAC-seq data, comparing it against three commonly used approaches for defining chromatin accessibility features and two pre-identified accessible region sets. The first approach divides the genome into 2000 bp bins as basic chromatin accessibility feature units56–58, referred to here as “genomic bins”. The second approach aggregates all sequenced single cells into pseudo-bulk samples and calls peaks from these aggregates29,30, referred to as “pseudo-bulk peaks”. The third approach clusters cells based on ground-truth cell type labels, call peaks by combining cells within each cell type, and then unites the peaks from all clusters together as the reference features5,59, referred to as “united peaks”. These approaches, along with two pre-defined feature sets—cDHS and CATLAS cCREs—were evaluated against cPeaks for their performance in cell annotation tasks. For each approach, we tested different levels of feature selection, including all features, the top 50,000 highly variable features (HVFs), and the top 5000 HVFs (Methods, Fig. 4a).

We used a scATAC-seq dataset of hematopoietic differentiation9, comprising 2,034 cells and 10 cell types with FACS-sorting labels as ground truth (Supplementary Fig. 3a). For each approach, we generated cell-by-feature matrices, applied multiple embedding methods for dimensionality reduction19,29,42,56,60,61, and used the Louvain method62 for clustering and annotation (Fig. 4a, Methods). Annotation results were assessed using Adjusted Rand Index (ARI), Adjusted Mutual Information (AMI), and homogeneity (H)63 (Fig. 4a, Methods). Across all evaluation metrics, cPeaks consistently delivered the best or among the best performances compared to the other feature-defining approaches (Fig. 4b and Supplementary Fig. 3b, c). The superiority of the cPeaks became more pronounced when using the top 5000 or top 50,000 HVFs. Considering that the definition of cPeaks was independent of this particular dataset, this performance demonstrates that cPeaks can effectively serve as a robust general reference for scATAC-seq data.

To further validate the performance of cPeaks, we employed the iterative method of ArchR36, which employs an iterative reduction approach to select appropriate features and refine cell embedding results. We conducted 1, 2, and 4 rounds of iteration, comparing cPeaks to the five other feature sets (Fig. 4a, Methods). The results showed that cPeaks outperformed other feature sets after just one iteration and continued to improve with additional iterations (Fig. 4c). The findings demonstrate that cPeaks not only streamline the ArchR analysis workflow but also provide superior results compared to other feature sets, making cPeaks an optimal choice for cell type annotation of scATAC-seq datasets.

We further tested the effectiveness of cPeaks on a peripheral blood mononuclear cell (PBMC) dataset profiled using DOGMA-seq24, referred to as PBMC-DOGMA. This dataset contains 7,623 cells annotated with 8 cell types, as provided by the original study (Supplementary Fig. 3d). Using the same evaluation framework, we applied the well-performing embedding methods from the hematopoietic dataset analysis (Supplementary Fig. 3e). The results showed that cPeaks again outperformed the other feature sets in most conditions (Supplementary Fig. 3f, g), reinforcing the generalizability across datasets and platforms, and demonstrating their advantage as a robust, pre-defined feature set for scATAC-seq data analysis.

Facilitating the discovery and analysis of rare cell types with cPeaks
An important expected advantage of using a dataset-independent reference like cPeaks, rather than relying on dataset-specific peak calls, is the enhanced ability to resolve low-abundance cell types in the dataset (i.e., rare cell types). In conventional scATAC-seq workflows, peaks are typically called from pseudo-bulk profiles aggregated across all cells. In this process, peaks specific to rare populations are often missed due to their low abundance, leading to information loss early in the analysis. In contrast, cPeaks are constructed across a wider range of tissues and cell types, not dominated by their prevalence in any single dataset. This avoids frequency-driven filtering and preserves many informative regions that can distinguish rare populations. In addition, our CNN model can predict cPeaks for unseen or masked cell types and tissues (Fig. 2c), thereby recovering accessible regions that were not directly observed in the training data. Together, these properties ensure that rare or missing cell types are not overlooked in downstream analysis.

To evaluate this advantage, we analyzed a PBMC scATAC-seq dataset (PBMC-8k) from 10x Genomics, which contains 8728 cells (Methods). We constructed cell-by-feature matrices using five different references: cPeaks, cDHSs, CATLAS cCREs, 2000-bp genomic bins, and pseudo-bulk peaks (called from the same dataset). Each matrix underwent normalization, feature selection, and dimensionality reduction using standard procedures (Methods). We transferred cell type labels from a corresponding scRNA-seq PBMC dataset through cross-modality integration using the anchor method64 in Seurat65. For each cell in the scATAC-seq data, we obtained an inferred cell type with a prediction score reflecting the annotation confidence (Supplementary Fig. 4a).

Uniform Manifold Approximation and Projection (UMAP) visualizations showed that most feature sets successfully resolved major cell types, including three rare subtypes: plasmacytoid dendritic cells (pDCs), dendritic cells, and NKbright cells, each representing less than 1% of total cells (Fig. 5a, Supplementary Fig. 4b). Genomic bins, however, failed to separate these rare subtypes. To quantify clustering accuracy, we calculated five metrics: the Calinski–Harabasz index, mean cell-type Local Inverse Simpson’s Index (cLISI), Davies–Bouldin index, balanced k-nearest neighbor (KNN) accuracy, and mean silhouette score (Methods). cPeaks consistently achieved the best or second-best performance across all metrics (Supplementary Fig. 4c). In addition, when transferring scRNA-seq labels to scATAC-seq data, cPeaks yielded higher prediction scores than all other references (Fig. 5b), further supporting their effectiveness in improving annotation quality.

To assess the utility of cPeaks in downstream analyses, we examined differentially accessible regions identified using each reference set. Taking pDCs as an example, we compared their accessibility profiles to all other cells (Methods). We found that cPeaks enabled the identification of slightly more differentially accessible regions than cDHSs, and substantially more than CATLAS cCREs and pseudo-bulk peaks. We compared the differential accessible regions identified by cPeaks against pseudo-bulk peaks, and found that many of them were specifically accessible in pDCs (Fig. 5d and Supplementary Fig. 4d), highlighting the limitations of pseudo-bulk peaks in detecting accessible regions for rare cell types. These differential accessible regions were important for investigating the functions of rare cell types. For example, we identified the promoter region of the PTCRA gene as a pDC-specific accessible region (Fig. 5e, f). Consistently, scRNA-seq data showed that PTCRA expression was restricted to pDCs (Fig. 5g), aligning with previous studies that identified PTCRA as a marker of human pDCs66,67.

We repeated this analysis on two additional PBMC datasets: PBMC-10k from 10x Genomics (10,412 cells) and PBMC-48k68 (48,487 cells; Supplementary Fig. 5, Methods). The results were consistent: cPeaks again ranked among the top in clustering metrics (Supplementary Fig. 5c and Supplementary Fig. 5g), and enabled the identification of comparable or more DARs compared to other references (Supplementary Fig. 5d and Supplementary Fig. 5h). These findings further support the robustness and generalizability of cPeaks for facilitating the analysis of rare cell types across datasets.

We extended this analysis to rare subtypes within cell groups, focusing on gastrointestinal (GI) epithelial cells in CATLAS5. They included nine subtypes, four of which represented less than 2% of the cells (Supplementary Fig. 6a). Using the same analysis procedure (Methods), we found that CATLAS cCREs, cDHSs, and 2,000-bp genomic bins were insufficient to clearly distinguish the rare subtypes from other epithelial cells (Supplementary Fig. 6b). In contrast, both cPeaks and pseudo-bulk peaks enabled clearer separation of these subtypes (Fig. 5h and Supplementary Fig. 6b), with cPeaks offering the additional advantage of requiring no dataset-specific customization or additional peak calling. We quantified all the five clustering accuracy metrics and found that cPeaks outperformed the other feature sets (Supplementary Fig. 6c). In particular, cPeaks achieved the highest KNN accuracy in identifying each of the four rare subtypes (Supplementary Fig. 6d), demonstrating its enhanced ability to resolve fine-grained cell identity.

We further employed cPeaks to reconstruct the whole CATLAS, which encompasses 30 adult tissue types5. CATLAS includes a pre-defined unified feature set, CATLAS cCREs, which was derived from pseudo-bulk peak-calling and peak integration. To evaluate cPeaks, we performed one-round cell clustering using both cPeaks and CATLAS cCREs as features, following the strategies in the original paper5 (Methods). We assessed the results using 30 major cell groups and 111 cell types annotated by the original paper as ground truth. Both cPeaks and CATLAS cCREs exhibited comparable performance for the 30 major groups (Supplementary Fig. 6e).

When clustering was refined to 111 cell types, cPeaks showed slightly superior to CATLAS cCREs on all evaluations. For example, ARI increased by 0.03, and silhouette score increased by 0.02 (Supplementary Fig. 6f). Depending on the proportion of each cell type in the whole dataset, we categorized cell types into three parts, each with 37 cell types: extremely rare cell types (average proportion 0.1%), rare cell types (average proportion 0.3%), and common cell types (average proportion 1.2%) (Fig. 5i, Methods). We randomly sampled 50% of the dataset 10 times and calculated clustering metrics for each feature set on each subset. cPeaks consistently achieved slightly better performance than CATLAS cCREs across most metrics in all three abundance groups, including those with low cell numbers. Given that CATLAS cCREs were derived from the same dataset, these results suggest that cPeaks offer a robust and transferable feature set for analyzing diverse cell populations, including rare subtypes.

Overall, the results demonstrated that cPeaks can simplify scATAC-seq data processing and facilitate the discovery of rare cell types and the analysis of their specific features, which may otherwise be difficult to detect using other approaches. This capability reinforces its value as a universal reference for scATAC-seq analyses, enabling robust exploration of cellular diversity.

Investigating cell state transitions during development with cPeaks
cPeaks were defined using the data of ATAC-seq datasets of adult healthy samples. We investigated their usability and performance in analyzing developmental data. We utilized a human fetal retinal scATAC-seq dataset encompassing early (8 weeks) and late (13 weeks) stages of retina organogenesis34 (Supplementary Fig. 7a). The late stage included two anatomically distinct regions: the center and the periphery of the retina. Using cPeaks as the reference, we generated cell-by-feature matrices, applied UMAP for dimensionality reduction, and labeled cells according to the original study (Methods).

Using cPeaks as features, we were able to reproduce the developmental trajectories that had been reported in the original study, supporting their utility for analyzing developmental processes despite being constructed solely from adult samples. For early-stage cells, three developmental trajectories were observed, leading from multipotent progenitor cells to horizontal cells, cones, and retinal ganglion cells (Fig. 6a). The pseudo-time values of these cells displayed a clear progression order in the UMAP (Fig. 6b), similar to the trend reported in the original study. For later stages, we separately analyzed samples collected from the center and the periphery of the retina, and the results showed clear clusters of cell types and their transition states (Supplementary Fig. 7b, c).

While multiple reference sets may also capture the overall trajectories (Supplementary Fig. 7d), only cPeaks provide the additional property of being “well-positioned.” Building on our previous finding that TFs associated with well-positioned cPeaks are significantly enriched for GO terms related to development (Supplementary Fig. 2g) and show substantial overlap with known pioneer TFs, we hypothesized that well-positioned cPeaks play critical roles in cellular differentiation and development. To test this, we analyzed the accessibility dynamics of well-positioned cPeaks during development by calculating the well-positioned-cPeak ratio in all accessible cPeaks (Methods). The results revealed increased well-positioned-cPeak ratio during transition periods, followed by a decrease in accessibility at later developmental stages compared to earlier ones (Fig. 6c and Supplementary Fig. 7a–c). To further support this observation, we examined the ratio along pseudo-time in the early-stage data and found a marked increase during cell state transition (Fig. 6d, Methods), reinforcing the association of well-positioned cPeaks with dynamic regulatory activity during development.

To further investigate differences in well-positioned-cPeak ratio between fetal and adult cells, we examined two large single-cell chromatin accessibility datasets: a human fetal cell atlas26 and CATLAS5 (Methods). We first visualized the distribution of well-positioned-cPeak ratios across cells in both datasets together using UMAP projections. The results showed that fetal cells generally exhibited higher well-positioned-cPeak ratios compared to adult cells (Supplementary Fig. 7e, f). To quantify this difference, we calculated the ratio for each cell type. We found that for most cell types shared by fetal and adult samples, the well-positioned-cPeak ratio consistently decreased in adult cells relative to their fetal counterparts (Fig. 6e). These findings underscore the association of well-positioned cPeaks with cellular differentiation and development.

Characterizing tumor cell subtypes and tumor progression with cPeaks
We then investigated the use of cPeaks in analyzing tumor cells. We re-analyzed a human gynecologic malignancy scATAC-seq dataset of 74,621 cells in the endometrium or ovary of 11 treatment-naive cancer patients35, using cPeaks as the reference (Methods). The analysis demonstrated clear clustering of cells by type, with tumor cells from gastrointestinal stromal tumor (GIST), endometrial cancer, and ovarian cancer distinctly separated from normal cells (Fig. 7a and Supplementary Fig. 8a). Labeling the cells with copy number variations (CNVs) (Fig. 7b) and the expression of cancer marker genes such as KIT, WFDC2, and MUC16 (Supplementary Fig. 8b) further validated these clusters.

Beyond replicating the findings from the original study, we observed that ovarian cancer cells in sample 8 formed two distinct clusters (Fig. 7a and Supplementary Fig. 8a). Both clusters expressed the cancer marker WFDC2 (Supplementary Fig. 8b), but only one cluster exhibited significant CNV enrichment (Fig. 7b). This observation suggested the presence of two subclones with distinct genomic properties within the tumor, highlighting the potential of cPeaks to identify subclusters in tumor cells. The generality of cPeaks as a reference may reveal subtle cellular differences that could be missed using dataset-derived references.

We further examined the accessibility dynamics of well-positioned cPeaks in various tumor cells. Figure 7c highlights the variation in the proportion of well-positioned cPeaks across tumor stages, showing a non-monotonic trend: an initial increase followed by a slight decrease at later stages. To validate this pattern, we analyzed the malignant cells in an independent hepatocellular carcinoma dataset69 and observed a similar trend (Fig. 7d and Supplementary Fig. 8c, Methods). We thus hypothesized that well-positioned cPeaks may reflect dynamic tumor progression.

cPeaks as a generic reference for ATAC-seq data analysis
We developed a comprehensive resource for cPeaks, including basic information and annotations (see Methods for details). The basic information includes cPeak IDs and genomic locations, while the annotations describe the housekeeping status and shape patterns.

To enhance usability, we linked cPeaks with the other widely used accessible region sets and regulatory databases, namely cDHS38, CATLAS cCREs5, and ReMap48,70, offering relevant information about regulatory vocabularies, cell-type specificity, and binding TFs (Methods).

To facilitate the practical application of cPeaks, we developed multiple computational interfaces that integrate with existing scATAC-seq analysis platforms. Specifically, cPeaks have been incorporated into SnapATAC271 and ArchR36, two widely used packages for single-cell epigenomics analysis. SnapATAC2, implemented in Python, supports efficient construction of cell-by-feature matrices, where cPeaks serve as a standardized reference, enabling consistent feature selection across datasets. In ArchR, an R-based platform, cPeaks can replace default peak or bin sets to generate a cell-by-feature matrix. Both platforms allow for the direct integration of cPeaks into their workflows, streamlining the analysis of chromatin accessibility in single-cell datasets.

In addition to these platform-specific integrations, we provide a standalone Python script for researchers requiring additional flexibility in their workflows. This script enables the transformation of fragment files into cPeaks-based data matrices and the transformation of cell-by-peak matrices into cell-by-cPeak matrices, offering a straightforward method to incorporate cPeaks into custom pipelines.

By integrating cPeaks into established analysis platforms and offering flexible standalone solutions, we provide a robust framework for analyzing scATAC-seq data. Comprehensive documentation, tutorials, and code examples for all tools are available at https://mengqiuchen.github.io/cPeaks/Tutorials.

Methods
cPeaks generation
Collecting ATAC-seq data for cPeak generation. We collected bulk ATAC-seq data from ENCODE40 and ATACdb41, focusing on datasets from healthy adult samples with more than 10,000 pre-identified peaks. For each dataset, we excluded peaks on chromosome Y (ChrY) if fewer than 20 ChrY peaks were identified, retaining the remaining peaks for further analysis. For datasets from ATACdb originally generated by the hg19 reference, we used CrossMap73 (version 0.6.5) to convert peak positions to the hg38 reference. Peaks that underwent length changes during the conversion were excluded for further analysis.

Defining merged peaks. Merged peaks were generated by joining overlapping peaks across datasets. To leverage the high-quality datasets from ENCODE and the broad cell type coverage of ATACdb, we generated merged peaks for ENCODE and ATACdb separately. The final merged peaks consisted of two parts: all ENCODE merged peaks, and all ATACdb merged peaks that were not overlapped with any ENCODE merged peak.

Defining shapes of merged peaks. For each position within a merged peak, we calculated an accessible score, defined as the normalized proportion of datasets in which the position was covered by peaks, reflecting its accessibility across datasets. These accessible scores were arranged sequentially along the 5’-to-3’ direction of the merged peak, forming a continuous profile that we referred to as the “cPeak shape”. This shape captures the distribution of accessibility across the peak, providing a detailed representation of both consistent and variable accessibility patterns across datasets. When defining shapes of merged peaks on ChrY, we accounted only for the number of datasets retaining ChrY as the denominator.

Evaluating the consistency of peak shapes. We assessed the consistency of merged peak shapes across various biological and technical conditions. For evaluating different tissue types, we divided the ENCODE bulk ATAC-seq datasets into two groups based on tissue origin: 48 datasets from blood and 101 from solid tissues (e.g., lung, heart, brain). We applied the same merging procedure independently to each group. For evaluating different peak-calling parameters, we generated merged peaks separately from ENCODE and ATACdb datasets, These datasets were processed using MACS228 with different peak-calling parameters. For evaluating different sequencing technologies, we analyzed merged peaks generated by all ENCODE bulk ATAC-seq samples with merged peaks generated by 733 DNase-seq samples used for generating for cDHS38 (downloaded from https://www.encodeproject.org/annotations/ENCSR857UZV/). For evaluating different peak-calling methods, we compared merged peaks generated by all ENCODE bulk ATAC-seq samples with merged peaks generated by a 10x Genomics scATAC-seq sample, which used MACS228 and CellRanger29 for peak calling, respectively. The samples were downloaded from https://www.10xgenomics.com/datasets/10k-human-pbmcs-atac-v2-chromium-controller-2-standard. We first clustered the single cells into 17 groups, constructed pseudo-bulk samples for each group, and then performed peak calling followed by merging to obtain the merged peak set.

We selected merged peaks derived from each approach with a shape maximum greater than 0.3 for further analysis. In each comparison, we designated a baseline group and a comparison group, and focused on the shared merged peaks between them. For tissue-type-based evaluation, peaks generated from blood samples served as the baseline group, and those from solid tissues as the comparison group. For all other evaluations, ENCODE-derived peaks were used as the baseline group.

To assess shape similarity, we used the genomic window defined by merged peaks generated from all ENCODE bulk ATAC-seq samples as a reference. The shape profile of each peak was aligned to this reference window, with excess regions trimmed and shorter peaks zero-padded to ensure comparability. We then treated the base-pair–level signal values across the aligned window as observations and calculated the PCC to quantify similarity (Methods, Supplementary Fig. 1b).

To evaluate the significance of the observed PCC values, we introduced two negative control groups (Supplementary Fig. 1b). In the shuffled control, peaks from one group were randomly paired with peaks from the other group, aligned by center, and their lengths were adjusted as above before calculating PCC. This control captures the background similarity expected from unrelated genomic regions. In the shifted control, each peak was displaced 250 bp downstream relative to its partner before PCC calculation, and their lengths were adjusted as above before calculating PCC. These two complementary baselines allowed us to evaluate whether observed PCC values for paired peaks were significantly higher than expected by chance or misalignment.

Empirical p-values were calculated by comparing each observed PCC to the corresponding null distribution. Specifically, we defined the empirical p-value as the proportion of shuffled or shifted control PCC values greater than or equal to the observed PCC.

For evaluating shape consistency on the set of observed cPeaks, we used the observed cPeaks to define the reference genomic windows for length unification. All other procedures remained the same as described above.

Generation of cPeaks from merged peaks. We used the Hartigan’s dip test74 (the “dip.test” function in the diptest package (version 0.76-0) in R) to measure the multimodality of merged peak shapes. We defined merged peaks with p-values less than 0.1 as multimodal ones. These multimodal peaks were then split into multiple unimodal peaks based on their shapes. For each position within a multimodal peak, we calculated the number of datasets where it was accessible and created an equivalent number of tuples. The first dimension of each tuple corresponds to the position, while the second dimension is a value randomly sampled from a standard normal distribution. In this way, we transformed the shape of each peak into a point cloud in a two-dimensional space. Then, we performed Hierarchical Density-Based Spatial Clustering of Applications with Noise (HDBSCAN)75 to cluster the point cloud using the “hdbscan” function in the dbscan package (version v1.1-11) in R with variable minimum size of clusters (minPts) set as the number of points divided by the merged peak length, multiplied by 50. We used the cluster ID calculated by HDBSCAN to split the merged peak into several separate regions. For each separated region, we decreased the value of minPts by 5 and repeatedly performed clustering and splitting until all separated regions were unimodal or achieved the minimum length of 800 bp. After, all the separated regions were regarded as cPeaks generated from the multimodal merged peak.

Assessing the availability of cPeaks as a chromatin accessibility reference
To analyze the distribution of sequencing reads around cPeak centers, we aligned the centers of all cPeaks and collected reads mapped within ±1000 bp of these centers to generate an aggregated distribution. The frequency of read counts at each position was calculated and normalized by dividing it by the total number of aggregated reads.

For comparing cPeaks with cDHS in terms of coverage over ATAC-seq accessible regions, we downloaded DNase-seq and ATAC-seq data of GM12878 from the ChIP-Atlas website43. To ensure objectivity and avoid selection bias, we investigated all available DNase-seq and ATAC-seq datasets for GM12878 in H. sapiens (hg38) from the ChIP-Atlas database. Data files were filtered by setting the “significant” parameter to 50 and restricting peak numbers to between 30,000 and 100,000. A total of 65 datasets were obtained, comprising 62 ATAC-seq datasets and 3 DNase-seq datasets.

Due to the limited number of high-quality, comparable DNase-seq/ATAC-seq datasets for the same cell type, constructing well-matched pairs was not straightforward. To ensure objectivity, we paired each DNase-seq dataset with the ATAC-seq dataset that had the closest number of peaks, without any manual selection. This resulted in three unbiased DNase–ATAC pairs: SRX5766093 & SRX069089, SRX7124440 & SRX069213, and SRX5766134 & SRX069213. Within each pair, regions uniquely identified by ATAC-seq or DNase-seq were defined as ATAC-seq- or DNase-seq-specific accessible regions, respectively.

The typical mapping rate reported in the technical reports of 10x Genomics can be found at https://cf.10xgenomics.com/samples/cell-atac/2.1.0/10k_pbmc_ATACv2_nextgem_Chromium_X/10k_pbmc_ATACv2_nextgem_Chromium_X_web_summary.html and https://cf.10xgenomics.com/samples/cell-atac/2.1.0/10k_pbmc_ATACv2_nextgem_Chromium_Controller/10k_pbmc_ATACv2_nextgem_Chromium_Controller_web_summary.html.

Expanding cPeaks through deep learning
Model architecture. The architecture of the deep CNN model was modified based on scBasset44, a one-dimensional (1D) CNN that included six 1D convolutional blocks followed by two linear blocks. Each convolutional block includes a 1D convolution layer, batch normalization, dropout, max pooling, and GELU activation function, while each linear block consists of a linear layer, batch normalization, dropout, and GELU activation function. The convolutional layers have 288, 323, 362, 406, 456, and 512 kernels, with kernel sizes of 17, 5, 5, 5, 5, and 5, respectively. The final convolutional output of 5,120 dimensions is reduced to 32 through a linear layer, and then to 2 dimensions by another linear layer, enabling binary classification. A softmax activation function is applied at the final layer to generate class probabilities.

Generation of samples. Observed cPeaks shorter than 100 bp were filtered out, and longer cPeaks were chunked into segments of 2000 bp, defining the positive samples. Using the “bedtools shuffle” function from BEDTools version 2.31.176, the positions of these positive samples were randomized across the human genome, avoiding blacklist regions, to create the negative samples. The “bedtools nuc” function76 was applied to compute the nucleotide composition of the genomic positions corresponding to both sample types, excluding any samples containing ‘N’. The remaining samples were converted into their respective base sequences for further analysis.

Model training. The sequences of cPeaks (positive samples) and negative samples were encoded using one-hot encoding. Their reverse-complement sequences were also included in the training process with the corresponding labels to enhance model robustness. The data was randomly split into training, validation, and test datasets in an 8:1:1 ratio. The model was trained to classify the samples, assigning positive samples a label of 1 and negative samples a label of 0. The Adam optimizer was used with a learning rate of 0.001 and a binary cross-entropy loss (BCELoss) function (1):

L = 1/N * sum(y_i * log(p_i) + (1 - y_i) * log(1 - p_i)) (1)

Where N is the number of samples, y_i in {0,1} is the ground truth label of the i-th sample, and p_i is its predicted probability of class 1 (i.e., is a cPeak). This formulation enables a consistent label structure between training and testing, which makes the model easily extensible to new, unseen tissues or cell types. Training was performed in batches of size 64. Early stopping was applied based on an increase in validation loss, and performance was evaluated on the test dataset.

Evaluating model’s performance in predicting missing cPeaks. For each cell type or tissue, we masked its specific cPeaks, treating the remaining cPeaks as positive samples. Negative samples were generated from background regions, including those corresponding to the masked cell type or tissue-specific cPeak regions, to simulate their absence. After preparing the positive and negative samples, we trained the CNN model and used it to predict the accessibility scores of the “lost” cell type or tissue-specific cPeaks.

Determination of the threshold for predicting additional cPeaks. Building on the masking experiments described above, where cell type– or tissue-specific cPeaks were withheld and then predicted by the CNN, we used the resulting score distributions to derive a decision threshold. Because the masked regions are genuine positives, their predicted scores were typically concentrated near 1. For each masking experiment, we examined the histogram of predicted scores and identified an elbow point (the value at which the frequency of scores began to rise sharply). This procedure was repeated for every masked cell type or tissue, and the resulting elbow points were averaged to yield a unified threshold of 0.872. This threshold was then applied to classify additional predicted cPeaks.

Evaluating the model’s performance by leaving one chromosome out. We first generated the full set of positive and negative samples across the genome. For each held-out chromosome, we removed all samples (both cPeaks and background regions) located on that chromosome and used the remaining samples for model training. The samples from the held-out chromosome were reserved as an additional independent test set to evaluate generalizability.

Interpreting the model with NeuronMotif. To interpret sequence patterns learned by the CNN model, we employed NeuronMotif47 to analyze the last convolutional layer of the trained CNN model. Using default parameters, we identified sequence motifs that activated each neuron and aligned them to known TF motifs in the HOCOMOCO v11 database77 using TOMTOM78. A motif was considered associated with a neuron if it matched with a q-value < 0.05. We then calculated, for each TF motif, the proportion of neurons it activated and its average frequency across neurons.

Predicting additional cPeaks. The trained model was subsequently used to predict additional cPeaks. Candidate sequences were generated by sliding a 500 bp window with a 250 bp stride across background regions, excluding regions on chrX, chrY, or shorter than 100 bp. These candidate sequences were fed into the trained model to compute accessibility scores, and sequences with scores exceeding the threshold defined by the mean elbow point were selected as predicted cPeaks. Importantly, the CNN model was trained only once using the observed cPeaks, and the predicted cPeaks were not fed back into the training data. This one-pass design was intentionally adopted to prevent potential error accumulation from iterative retraining, which could amplify boundary noise due to the fixed-window predictions scheme.

Predicting cPeak shapes. The base sequence corresponding to the genomic positions of cPeaks was used as input, with shape categories at each endpoint serving as labels. To ensure label balance, we included all well-positioned samples and randomly selected an equal number of weakly-positioned ones. The CNN model used the same architecture as before, excluding the final activation function. Input sequences were one-hot encoded, and the data were split into training, validation, and test sets in an 8:1:1 ratio. The network was trained using the AdamW optimizer and BCELoss criterion, with a batch size of 128 and a learning rate of 0.001. Early stopping was employed when validation loss increased, and performance was assessed on the test set.

Using the trained CNN, we predicted the shape categories of all unlabeled cPeaks. Outputs for each endpoint were plotted across the genome, and a kernel density estimate (KDE) plot was generated. A threshold was determined based on the trough between two peaks in the KDE plot. This threshold was applied to classify the shape categories of all cPeaks, including predicted cPeaks.

Annotating cPeaks
Comparison between cPeaks and ReMap ChIP-seq peaks. We downloaded the non-redundant peak set of ReMap 202248 from https://remap.univ-amu.fr/storage/remap2022/hg38/MACS2/remap2022_nr_macs2_hg38_v1_0.bed.gz. We compared three chromatin accessibility references (cPeaks, cDHSs, and CATLAS cCREs) by calculating the proportion of ReMap peaks they overlapped. To assess enrichment beyond random expectation, we generated background peak sets for each reference: for each reference, we randomly sampled the same number of peaks, with the same length distribution, from the genome after excluding blacklist regions and the original reference peaks. We then computed overlaps between ReMap peaks and both the reference and its background. Fold enrichment was calculated as the ratio of overlap with the true reference versus its background. Additionally, to examine overlap specificity at the TF level, we calculated the percentage of each individual TF’s ReMap peaks that overlapped with each reference.

Footprinting analysis. To illustrate the presence of transcription factor (TF) binding sites within observed cPeaks, we used TOBIAS79 to analyze bulk ATAC-seq BAM files from a classical study80. We used the GM12878 sample as an example, downloaded from the Gene Expression Omnibus (GEO) database under accession code GSE47753. We selected the canonical motifs CTCF and BATF as representative examples, and TF-footprinting analysis was performed within the accessible chromatin regions defined by cPeaks on chromosome 4.

Genomic annotations of cPeaks. We used the ChIPSeeker package81 (version 1.35.1) to annotate cPeaks with genomic annotations. We used the “plotPeakProf2” function to compare the observed cPeaks with genomic annotations of transcriptional start sites (TSSs) and transcriptional end sites (TESs). We discretized the genomic region into 800 bins and using transcript annotations (TxDb) to map the cPeaks to TSSs and TESs. Then, we defined the promoter region as ±2000 bp around TSS and annotated cPeaks using the “annotatePeak” function. The annotation database was EnsDb.hg38.v86 provided by the ensembldb package27. ChIPSeeker categorized cPeaks into seven types: promoter, 5’ UTR, 3’ UTR, exon, intron, downstream (≤300 bp downstream of the gene end), and distal intergenic regions.

Inferring regulatory roles of cPeaks. We categorized cPeaks based on their overlap with genomic and epigenetic features. cPeaks overlapping any gene promoter region were defined as known promoters, while those overlapping 3’ UTRs were annotated as UTRs. To identify cPeaks associated with CTCF binding sites, we analyzed 210 CTCF ChIP-seq data from ENCODE40 and classified cPeaks overlapping any CTCF ChIP-seq signals as CTCF binding sites. Histone modification data from ENCODE40 were also incorporated, using their called peaks as signals. cPeaks overlapping H3K4Me1 signals were predicted as enhancers, while cPeaks overlapping H3K4Me3 signals but not overlapping known promoter regions were categorized as predicted promoters.

Annotating housekeeping cPeaks. We used a pre-defined housekeeping gene set49, downloaded from https://www.tau.ac.il/~elieis/HKG/HK_genes.txt. A contingency table was constructed to quantify the number of housekeeping and non-housekeeping genes associated with housekeeping and non-housekeeping cPeaks, considering only protein-coding genes. We used Fisher’s exact test to calculate the significance of the contingency table. For housekeeping cPeaks annotated as known promoters, the associated genes were defined as those driven by these promoters. Functional enrichment analysis of GO biological process terms was performed on these genes using the “enrichGO” function in the ChIPSeeker package81. Enriched pathways were further summarized using the “simplify” function in the clusterProfiler package82.

Motif enrichment analysis. Only the endpoints of unimodal peaks were selected for this analysis. We extracted the 200 bp sequences from the inner sides of the cPeak endpoints. We used “findMotifsGenome.pl” in Homer to perform motif enrichment analysis, with sequences derived from well-positioned edges as the experimental group and those from weakly-positioned edges as the background. We used “annotation.pl” in Homer to annotate the enriched motifs83. A KDE plot was generated using seaborn.kdeplot84 to visualize the distribution of motifs around the endpoints. We filtered TFs associated with the enriched motifs with an adjusted p-value of less than 0.1 for domain and GO enrichment analysis. We used the enrichR library85 for enrichment analysis. For GO enrichment, we used the GO_Biological_Process_2023 database86 and applied an adjusted p-value threshold of 1e–5. Subsequently, we clustered the enriched GO terms using the “GO_similarity” function in the “simplifyEnrichment” library87 to group similar terms.

Evaluating the cell type annotation performance of cPeaks
Data preprocessing. We obtained BAM files for the FACS-sorting hematological differentiation scATAC-seq data from https://github.com/pinellolab/scATAC-benchmarking/tree/master/Real_Data/Buenrostro_2018/input/sc-bams_nodup, which was provided by a benchmark study63. The dataset included 2034 cells and 10 cell types with the hg19 reference genome. We obtained fragment files for the PBMC-DOGMA dataset from the GEO database under accession code GSE156478. The dataset included 7623 cells and 8 cell types with the hg38 reference genome. The cell type label was obtained from the original study24.

We generated cell-by-feature matrices using six feature sets: genomic bins, pseudo-bulk peaks, united peaks, cPeaks, cDHSs, and CATLAS cCREs. Each feature set provided genomic regions defined as accessible regions. For genomic bins, the genome was divided into fixed-length 2000 bp bins, each treated as a feature.

For pseudo-bulk peaks of the hematopoietic dataset, we used the pre-defined peaks of the original study9 from the GEO database under accession code GSE96769, where peaks were called on all aggregated cells using MACS228. United peaks were derived from the peak set provided by the benchmark study63 from https://github.com/pinellolab/scATAC-benchmarking/blob/master/Real_Data/Buenrostro_2018/input/combined.sorted.merged.bed. For pseudo-bulk peaks of the PBMC-DOGMA dataset, we obtained the pseudo-bulk peaks from the original study. We generated the united peaks of this dataset by calling the pseudo-bulk peaks for each cell type, and merged all the peaks as the united peaks.

For cPeaks, we excluded housekeeping cPeaks and extra-long cPeaks (> 2,000 bp), and converted the cPeaks from hg38 to hg19 using CrossMap73 for the hematopoietic dataset. For cDHSs, we downloaded the hg19 version from https://zenodo.org/records/3838751/files/DHS_Index_and_Vocabulary_hg19_WM20190703.txt.gz, and the hg38 version from https://zenodo.org/records/3838751/files/DHS_Index_and_Vocabulary_hg38_WM20190703.txt.gz. For CATLAS cCREs, we retrieved the hg38 version from http://catlas.org/catlas_downloads/humantissues/cCRE_hg38.tsv.gz and converted it to hg19 using CrossMap for the hematopoietic dataset.

For each cell in the hematopoietic dataset, sequencing reads in the BAM file were mapped to the corresponding feature set, and the number of reads overlapping each feature was counted as chromatin accessibility. The chromatin accessibility values across all features and cells were compiled into a cell-by-feature matrix for downstream analyses. For each cell in the PBMC-DOGMA dataset, we mapped each fragment to the reference to obtain a cell-by-feature matrix.

Selecting HVFs. We first applied TF-IDF transformation to normalize the cell-by-feature matrix. We then calculated the variance of each feature across all cells after normalization. We sorted the features by their variance and selected the top features as HVFs. We selected 5000 and 50,000 HVFs for the hematopoietic dataset, respectively, and selected 5000 and 10,000 HVFs for the PBMC-DOGMA dataset, respectively, as the number of pseudo-bulk peaks of the PBMC-DOGMA dataset is less than that of the hematopoietic dataset.

Embedding methods for cell representations. We utilized various embedding methods to process cell-by-feature matrices.

For ChromVAR-motif42, we used the ChromVAR package (version 1.16.0) in R, adjusted for GC bias, resized features to 500 bp, and matched motifs from the JASPAR database88 using the “matchMotifs” function. The resulting cell-by-motif matrix was transformed into a cell-by-deviation Z-score matrix via “computeDeviations” in ChromVAR. ChromVAR-motif-PCA extended this method by performing principal component analysis (PCA) on the deviation matrix and selecting the top 10 principal components (PCs) for the final embedding. For ChromVAR-kmer, instead of motifs, we used 6-mers as features with the “matchKmers” function in ChromVAR, generating a cell-by-kmer deviation Z-score matrix. ChromVAR-kmer-PCA followed the same steps but added PCA to select the top 10 PCs. cisTopic60 (version 0.3.0) in R applied latent Dirichlet allocation (LDA) on the cell-by-feature matrix using the “runCGSModels” function, generating models with varying topic numbers as 10, 20, 25, 30, 35, and 40. The best model was selected using “selectModel”, and the final cell-by-topic probability matrix was normalized via “modelMatSelection” with “method=probability”. Scasat61 computed cell-cell Jaccard distances from the cell-by-feature matrix and reduced dimensions to 15 via multidimensional scaling (MDS). Signac89 used Seurat (v4.3.0)65 to create a cell-by-feature matrix, applied TF-IDF transformation, and performed latent semantic indexing (LSI) via SVD, retaining PCs 2~50 for reduced dimensions. SnapATAC271 generated a cell-by-cPeaks matrix using “snapatac2.pp.make_peak_matrix” and reduced dimensions with “snapatac2.tl.spectral”. Neighbor graphs and cell clusters were computed using “snapatac2.pp.knn” and “snapatac2.tl.leiden”, respectively, yielding 10 cell types. TF-IDF-PCA19,56,90 filtered sex chromosome features, applied TF-IDF transformation, and performed SVD to select the top 50 PCs as embedding features. ArchR36 used cPeaks as features, reduced dimensions via “addIterativeLSI”, and identified 10 cell types using “addClusters”, with iterations set to 1, 2, and 4 for comparison.

Evaluation. After embedding, we performed Louvain clustering62 with the number of clusters set equal to the number of FACS-sorted labels. The clustering results were evaluated using three metrics: Adjusted Mutual Information (AMI), Adjusted Rand Index (ARI), and Homogeneity (H). All these indexes were calculated using the “sklearn.metrics.cluster” library in Python.

Analysis of the PBMC datasets for rare cell type discovery and analysis
We download the fragment file of the PBMC-8k dataset from the 10x Genomics website (https://cf.10xgenomics.com/samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_fragments.tsv.gz). Cells with extremely high total number of fragments may represent doublets and were removed. We mapped all fragments to 2000-bp genomic bins, cPeaks, cDHSs, and CATLAS cCREs, respectively, and counted the number of fragments that overlapped with each reference peak as its accessibility. This process generated a cell-by-feature matrix for each reference. For comparison, we downloaded the matrix using pseudo-bulk peaks as features from the same website (https://cf.10xgenomics.com/samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5). Then, we followed the Signac89 PBMC vignette (https://stuartlab.org/signac/articles/pbmc_vignette.html) to preprocess these cell-by-feature matrices separately. Low-dimensional clustering and UMAP visualization were performed using dimensions 2~30.

We annotated cells in the scATAC-seq data by integrating them with a labeled PBMC scRNA-seq dataset provided by 10x Genomics. The raw scRNA-seq data was downloaded from https://support.10xgenomics.com/single-cell-gene-expression/datasets/3.0.0/pbmc_10k_v3, and the pre-processed scRNA-seq Seurat object was downloaded from https://signac-objects.s3.amazonaws.com/pbmc_10k_v3.rds. We utilized the “FindTransferAnchors” function in Seurat65 for cross-modality integration and used the “TransferData” function in Seurat to infer cell type labels in the scATAC-seq dataset along with their corresponding prediction scores. After label transformation, cluster 18 was identified as having low prediction scores, with an average below 0.5. This indicated that the cells in this cluster were likely of low quality. Consequently, this cluster was excluded from subsequent differential accessibility analyses.

To quantitatively assess the clustering performance under different reference peak sets, we used five metrics: Calinski–Harabasz index, Davies–Bouldin index, mean silhouette score, balanced k-nearest neighbor (KNN) accuracy, and mean cLISI (cell-type Local Inverse Simpson’s Index). The Calinski–Harabasz index evaluates the ratio of between-cluster dispersion to within-cluster dispersion; a higher value indicates more distinct clusters. The Davies–Bouldin index measures the average similarity between each cluster and its most similar one; lower values imply better cluster separation. The silhouette score reflects how well each cell fits within its assigned cluster relative to other clusters; higher scores represent more coherent clustering. Balanced KNN accuracy was computed using k=10. Each cell’s predicted label was inferred from a majority vote among its k-nearest neighbors in the PCA space (KNN accuracy). Balanced accuracy averages accuracies across all cell types, thus adjusting for label imbalance. Mean cLISI was used to quantify local label purity. For each cell, cLISI measures the diversity of neighboring cell types, with lower values indicating less label mixing and purer clusters. All metrics were calculated in the first 30 dimensions of the LSI space generated from the TF-IDF-normalized matrix using singular value decomposition (SVD).

We performed differential accessibility analysis to find differentially accessible regions in pDCs compared to other cell types. We used “FindMarkers” in the Signac package89 with test.use = “LR” and latent.vars=“peak_region_fragments” to perform differential analysis.

The PBMC-10k dataset was downloaded from the 10x Genomics website (https://support.10xgenomics.com/single-cell-multiome-atac-gex/datasets/1.0.0/pbmc_granulocyte_sorted_10k). We followed the processing steps outlined in the Seurat v5 vignette (https://satijalab.org/seurat/articles/seurat5_atacseq_integration_vignette). As this is a multi-omics dataset, we did not perform label transfer from scRNA-seq to scATAC-seq. Cell type labels were directly taken from the original dataset metadata. The pseudo-bulk peaks were also derived from the original dataset. We followed the same procedures used for the PBMC-8k dataset to conduct all downstream analyses.

The PBMC-48k dataset68 was downloaded from GEO (accession code GSE190992), containing paired scRNA-seq and scATAC-seq data. We processed the scATAC-seq data using the ArchR pipeline36, identifying 9.08% of cells as doublets using addDoubletScores and removing them with filterDoublets. For the scRNA-seq data, cells with predicted.celltype.l2.score > 0.9 were retained as high-confidence annotations. These cells were then used to transfer labels to scATAC-seq data via the addGeneIntegrationMatrix function in ArchR, which also provides a predictedScore for each assigned label. We selected 48,487 scATAC-seq cells with a predictedScore > 0.9 for downstream analyses, following the same procedures as for the PBMC-8k dataset. We found that cPeaks did not outperform other methods when doublets were not removed (Supplementary Fig. 5i), This observation emphasizes that rigorous data quality control is an essential step in ensuring the reliability of analyses.

Analysis of the CATLAS dataset for rare cell type discovery and analysis
The fragment file and metadata of CATLAS data were downloaded from http://catlas.org/catlas_downloads/humantissues/. We first analyzed GI epithelial cells from CATLAS5, constructing cell-by-feature matrices with 2000-bp genomic bins, cPeaks, cDHSs, CATLAS cCREs, and pseudo-bulk peaks, separately. For cPeaks, we used non-housekeeping cPeaks accessible in more than 50 cells as features. Pseudo-bulk peaks for GI epithelial cells were called using MACS2 with the parameters --shift 0 --nomodel --call-summits --nolambda --keep-dup all. Each matrix was preprocessed using the same steps applied to PBMC-8k, followed by graph-based clustering and UMAP visualization across dimensions 2 to 30. The evaluation of clustering performance was the same as PBMC-8k.

For the analysis of the whole CATLAS, we imported 92 fragment files from CATLAS into snapATAC271. Using the barcodes in the metadata download from http://catlas.org/catlas_downloads/humantissues/fragment/, we got 615,998 cells. Batch correction was performed using the “snapatac2.pp.mnc_correct” function. We used the “snapatac2.pp.select_features” function to select 500,000 features. We reduced the dimensions to 30 through the “snapatac2.tl.spectral” function and generated a neighbor graph through “snapatac2.pp.knn”. UMAP visualization was carried out using the “snapatac2.tl.umap” function. Adjusting the resolution parameter in “snapatac2.tl.leiden”, we identified either 30 or 111 cell clusters.

We calculated the proportion of each of the 111 cell types annotated in the CATLAS paper5. Depending on the ratio, we categorized cell types into three parts, each with 37 cell types. Based on the ranking of cell proportions, we evenly divided these cell types into three groups, with 37 cell types in each group. We sampled 50% from the whole dataset for 10 times, and calculated AMI, ARI, and H for each group, using the 111 cell type labels as ground truth, to evaluate clustering performance.

Analysis of the development datasets
We retrieved the fragment files and metadata for the retinal fetal dataset34 from GEO (accession code GSE184386) and processed the data following the filtering steps outlined in the original study. The generation of cell-by-feature matrices, clustering, and visualizations were conducted using the same methodology as applied to the PBMC dataset analysis. Cell labeling was based on the provided metadata, and all analyses were performed independently for each sample.

To quantify the well-positioned-cPeak ratio in all accessible cPeaks, we first binarized the cell-by-feature matrices, and then computed the normalized accessibility of all well-positioned cPeaks which is defined as the number of accessible well-positioned cPeaks divided by the total number of well-positioned cPeaks. To visualize the relationship between pseudotime and the well-positioned-cPeak ratio, we performed outlier filtering based on pseudotime values within each cell type cluster. Specifically, cells with pseudotime values below Q1−1.5×IQR or above Q3+1.5×IQR were identified as outliers and removed, where Q1 and Q3 denote the first and third quartiles of pseudotime, respectively, and IQR=Q3−Q1. This IQR-based filtering ensured a robust fit by reducing the influence of extreme pseudotime values. We then plotted the well-positioned-cPeak ratio along pseudotime using a smoothed curve fitted by a generalized additive model, using the default settings of the geom_smooth function in the ggplot2 R package.

We downloaded the human fetal cell atlas data26 from http://catlas.org/catlas_downloads/humantissues/fragment/. The analysis of these data followed the procedure as that of CATLAS data. For visualization, we used the UMAP coordinates provided in the supplementary files of the CATLAS paper5, allowing us to directly map and compare the well-positioned-cPeak ratios across fetal and adult cell types in a consistent embedding space.

Analysis of the tumor datasets
We downloaded the fragment files and metadata for the gynecologic malignancy dataset35 from GEO (accession code GSE173682). Cells from all 11 patients were combined to generate a unified cell-by-feature matrix. We followed the preprocessing steps in the original study to filter cells. We downloaded the fragment files and metadata for the hepatocellular carcinoma dataset69 from GEO (accession code GSE227265). We selected all samples diagnosed with hepatocellular carcinoma that contained more than 50 malignant cells. Malignant cells from these samples were then combined to construct a unified cell-by-feature matrix.

The generation of the cell-by-feature matrix, clustering, and visualizations were performed using the same methodology as applied to the PBMC dataset analysis. We labeled each cell using labels from the metadata file. The calculation of the well-positioned-cPeak ratio was performed using the same methodology as applied to the retinal fetal dataset.

The resources for cPeaks
The resource for cPeaks includes basic information and cPeak annotation, organized into 14 columns:
- ID: A unique 13-character string representing each cPeak, starting with “CP”, followed by “HS” (human), and ending with a nine-digit number. For example, the first cPeak is encoded as “CPHS000000001”.
- source: Indicates the origin of the cPeak, either “observed” or “predicted”.
- chr_hg38: The chromosome where this cPeak is located in the hg38 reference genome.
- start_hg38: The start position of this cPeak in the hg38 reference genome.
- end_hg38: The end position of this cPeak in the hg38 reference genome.
- housekeeping: Specifies whether the cPeak is accessible across nearly all datasets (“TRUE” or “FALSE”).
- shape_pattern: The shape pattern of the cPeak, categorized as “well-positioned”, “asymmetrically-positioned”, or “weakly-positioned”.
- inferredElements: The inferred regulatory elements associated with the cPeak, such as “CTCF”, “TES”, “TSS”, “Enhancer” or “Promoter”.
- chr_hg19: The chromosome where this cPeak locates in the hg19 reference genome.
- start_hg19: The start position of this cPeak in the hg19 reference genome.
- end_hg19: The end position of this cPeak in the hg19 reference genome.
- cDHS_ID: The ID of the overlapped cDHS region in the hg38 reference, formatted as “chr_start_end”.
- CATLAS_ID: The ID of the overlapped CATLAS cCREs region in the hg38 reference, formatted as “chr_start_end”.
- ReMap_ID: The ID of the overlapped ReMap region in the hg38 reference, formatted as “chr_start_end”.

We provide the complete list of cPeaks in https://mengqiuchen.github.io/cPeaks/Tutorials and https://zenodo.org/records/1779936791. In addition, we offer a dataset-level matrix recording the presence or absence of each cPeak across all datasets analyzed in this study, available at https://zenodo.org/records/17799367/files/cpeaks-by-dataset.tsv.gz. To further enhance usability for researchers interested in organ- or tissue-level analyses, we also aggregated this matrix by tissue type to generate a tissue-by-cPeak matrix, which is available at https://zenodo.org/records/17799367/files/cpeaks-by-tissue.tsv.gz. While they do not represent a comprehensive atlas of chromatin accessibility for all human tissues, they provide a practical and scalable reference for both biological interpretation and computational modeling.

All data analysis and visualization were performed using R (version 4.1.3) and Python (version 3.10). Schematic illustrations and final figure assembly and layout were completed using Adobe Illustrator 2024.

Statistics & Reproducibility
The sample size was determined by the number of samples included in the analyzed datasets. No selection or manual intervention was applied to samples. Statistical analyses were conducted using standard methods implemented in R version 4.1.3. We used two-sided statistical tests for all related analysis. The specific methods for statistic tests used are described in the corresponding sections of the Methods and figure legends.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Meng et al developed a generic chromatin accessibility reference by aggregating peaks from 624 high-quality bulk ATAC-seq datasets, defining more than 1 million consensus peaks (cPeaks). They demonstrated cPeaks outperforms existing feature defining methods and references by improving cell annotation and rare cell type identification, and provide insights into chromatin dynamics during cellular differentiation and tumor progression. The collected large consortium of ATAC-seq would be a good resource in studying gene regulation. I have some comments as follows,\n1. How the batch effect like ATAC-seq generated from different labs, protocols, etc are corrected before the integration?\n2. For the developed CNN model is claimed to be modified based on scBasset. However, the difference and similarity of CNN and scBasset is unclear, and many technical details of CNN is missing. scBasset aims to predict the peaks for all cells in the single-cell setting, while the CNN is predict peaks for multiple tissues/organs. Please clarify this point and how the loss function is defined. Moreover, when the CNN performs zero-shot prediction, where training set does not include the tissue/organ/cell types in the testing set, please clarify the dimension of output labels in the training set and testing set as supervised learning need the same dimension of training and testing sets. Furthermore, are the training, validation and testing set split by different chromosomes? It is highly encouraged the detailed mathematical formula need to be present in the method section.\n3. Further comments have for using cPeaks for facilitating the discovery and analysis of rare cell types with cPeaks or Investigating cell state transitions during development with cPeaks. It is unclear which peak set of CPeaks for the data analysis, tissue-matched for the single-cell data or all peaks? What is the peak overlap between scATAC peak-calling and the peak set provided by CPeaks in the data analysis? It is also unclear additional peak information provided by CPeaks with the same read information can achieve both tasks. In another word, what is the intuition to using the mapped reads from studied scATAC-seq data in a broad collection of peaks across more tissues and cell types, provided cPeaks, would improve the identification of rare cell types?\n4. The github seems only provide the pipeline and software code. Would the CPeak data resources (peaks across all tissues/organs) can be provided in a format of online database as well as standard-alone download package? This will facilitate the users from both biology background and computational modeling.", "reviewer_2": "This manuscript presents cPeaks, a reference list of chromatin-accessible regions created by merging peaks from 624 bulk ATAC-seq datasets. The authors aim to provide a shared peak set for scATAC-seq analysis and employ a deep learning model trained on DNA sequence to predict additional peaks in unobserved tissues or cell types. While the study addresses an important challenge in the field, it has three major limitations. First, the cPeak list is not tissue- or cell-type-specific, which limits its applicability for studies requiring context-specific regulatory information. Second, using only DNA sequence to predict chromatin accessibility overlooks critical biological factors such as transcription factor binding and chromatin context, raising doubts about the validity of the predicted peaks. Third, the CNN model is trained on data from known tissues, and it remains unclear how it can reliably predict peaks in tissues it has not encountered during training. The detailed comments are as below.\nMajor Concerns\n1. The cPeaks are derived by merging peaks identified using various peak-calling tools, each with different characteristics and sensitivities. This can introduce substantial technical variation. Although the authors report shape similarity across methods, the Pearson correlation coefficients range from 0.4 to 1.0, indicating inconsistency. The authors should clarify how they account for this heterogeneity and whether peak-calling parameters were harmonized prior to merging.\n2. While the merged cPeaks serve as a universal reference, they may obscure important tissue- or cell-type-specific regulatory elements. Including annotations that flag cPeaks unique to specific contexts would enhance the utility of the resource. Alternatively, providing separate cPeak lists for each tissue or cell type would support more targeted and biologically meaningful analyses.\n3. The deep learning model used to predict additional cPeaks is trained solely on DNA sequence, without considering key biological factors that influence chromatin accessibility, such as transcription factor binding, histone modifications, and chromatin context. As a result, it is unclear whether the predicted peaks truly capture cell-type-specific regulatory activity. Moreover, the model is trained on cPeaks derived from known tissues, yet the manuscript does not provide sufficient justification for its ability to generalize to unseen tissues. Finally, the prediction framework presents a potential circular logic issue: newly predicted peaks could be added to the shared peak list and used for further model training, allowing the process to expand indefinitely without clear stopping criteria or external validation. The authors should clarify how they prevent such iterative propagation of prediction error.\n4. In benchmarking the coverage of scATAC-seq reads, the authors use MACS2-called peaks as ground truth. Since cPeaks may be partially derived from MACS2-called datasets, this could introduce bias in favor of cPeaks. Moreover, if MACS2 output is considered the gold standard, it raises the question of why not simply recommend MACS2-based peak calling per dataset. The authors should report the proportion of cPeak-contributing datasets that used MACS2. If MACS2 is dominant, it should be treated as a baseline method—like cDHS and CATLAS cCREs—rather than as the default reference. A side-by-side comparison of performance across all methods would be more appropriate.\n5. The reliability of the curated cPeaks would be strengthened through validation with independent chromatin accessibility assays, such as ChIP-seq, CUT&RUN, or CUT&Tag. Such validation would help confirm that cPeaks represent biologically relevant regulatory elements.", "reviewer_3": "Summary\nIn this manuscript, the authors assembled 1.4 million observed consensus peaks (cPeaks) on chromatin openness from bulk ATAC-seq data in ENCODE and ATACdb and developed a deep CNN model to predict an additional 0.28 million cPeaks across the human genome. They further annotated the cPeaks and demonstrated utilities of the cPeaks in labeling scATAC-seq data, identifying “rare” cell types, investigating cell state transitions during development, and characterizing tumor cell subtypes. On the positive side, the topic is interesting, the assembled cPeaks can be a valuable resource to the field, and most Results sections are well written. On the negative side, crucial details are missing in Methods sections, some benchmarking analyses are neither thorough nor convincing, and some statements are over stretched. The manuscript needs a significant improvement.\nMajor concerns\n1. The authors did not provide enough details on how Pearson correlation coefficient (PCC) is calculated. It can be easily misused and inflates similarity. What are the start and end points for PCC calculations when comparing two peaks? For example, comparing the peak in Fig. 1c with one shifting to the right by 250bp, how will the PCC be calculated? What is the PCC value?\n2. Details on permutation tests need to be provided and the PCC distribution from 1000 permutations should be provided. Otherwise, Fig. 1d is not informative. The authors may want to convert the PCC values in Fig 1d to empirical p values, using the distribution from the permutation tests as the null distribution, and show the histogram of the empirical p values.\n3. The similarity analysis should be performed on consensus peaks (cPeaks), not on multimodal peaks. Users care more about the end product than the intermediate results. Basically, the order between lines 120-126 and lines 104-119 should be switched.\n4. Why is there no evidence for footprinting from transcription factor binding in Fig. 1c or Fig. 1f?\n5. Fig. 4: This benchmarking analysis is very weak: a) The dataset of 2034 cells is very small and may miss complexities from more complicated datasets. The authors should include additional datasets generated from DOGMA-seq (Mimitou et al. 2021) and/or TEA-seq (Swanson et al. 2021) for benchmarking. b) The authors should provide more details on “clustering and annotation” (such as how many clusters and how many cell populations were used for benchmarking) and show a few examples of UMAPs.\n6. The benchmarking on PBMC dataset in Fig. 5 is very problematic: a) The authors benchmarked against five different feature (peak) sets in Fig. 4, why did they limit it to pseudo-bulk peaks here? It is not obvious other feature sets won’t be able to identify the three “rare” cell types. b) There are much larger datasets in the literature with tens to hundreds of thousands of PBMCs. Why did the authors benchmarked on a small dataset of 8728 cells? c) Did they consider the existence of doublets? d) The “rare” cell types (pDC, DC, NK bright) are not that rare and are thus poor examples for “discovering” rare cell types.\n7. The benchmarking on gastrointestinal (GI) epithelial cells in Fig. 5 is also problematic: a) Again, why did the authors limit the benchmarking to one feature set instead of five? b) The claim that “(p)seudo-bulk peaks in the original study failed to separate these subtypes” (lines 384-385) is not supported by the evidence in Fig. s5b. It is a purity issue, not a failure. c) It appears that there are impurities in the four rare cell types in Fig. 5h as well. The authors should quantify and compare impurity between methods. d) The improvements in Fig. 5i by cPeaks were small. Were they statistically meaningful?\n8. The presentation on Fig. 6 can be improved: a) Similar analysis with respect to Fig. 6ab can be done with other feature sets. What is the novelty here? b) Lines 451-453: A scatter plot of pseudo-time versus well-positioned-cPeak ratio or a plot like Fig. 6d can be informative here. c) Lines 454-458: The presentation on the two datasets should be separated and mirror the analysis on the retinal data. As is, it feels rushed and incomplete.\n9. Lines 483-488: Provide p values to assess significance.\nMinor concerns\n10. Lines 45-50: The literature review is incomplete. It ignores a different approach by a recent publication (Rachid Zaim et al. 2024).\n11. In addition to Table S1-2, summary pie charts showing tissue and cell type compositions can be informative.\n12. Fig. 1b: The curves grow linearly without reaching a plateau from 200 datasets. Descriptions should be accurate.\n13. How representative is Fig. 1g or 1h? Can the authors provide similar results based on more datasets? As is, it is hard to know whether it is “cherry picking”.\n14. The CNN model predicts cPeaks based on genome sequence regardless of tissue types. In other words, predicted cPeaks may not be open in specific tissues. The authors should distinguish testing prediction accuracy on known tissue-specific cPeaks from assigning predicted cPeaks as open to specific tissues. The latter is beyond the CNN model. The present description is ambiguous when claiming “that the CNN model can effectively predict missing cPeaks for previously unseen tissues or cell types”.\n15. Fig. 3h: Remove the extra “the distribution of” from the legend.\n16. Fig. S2b: Provide a legend.\n17. Lines 369-370: The anchor method in Seurat has limitations in labeling scATAC-seq data. The authors are encouraged to explore newer methods for labeling scATAC-seq data in the literature.\n18. It is unclear what is the difference between Fig. S5b and Fig. S5c. More details can be helpful.\n19. Fig. 6 and Fig. S6: a) Lines 447-450: Results on enriched GO terms should be provided. b) A UMAP by life stage can be informative. c) All UMAPs should use the same projection and thus look similar and comparable.\n20. Line 464: “Each box represents a single cell type …”. Did you mean each “line” here? How many total cell types were presented?\n21. Line 479: It should be “Fig. 7b” not “Fig. 6b”.\n22. Lines 489-494: These statements are stronger than what the evidence showed. The authors should weaken them to align with the evidence. Such statements also belong to Discussion, not Results.\n23. Many readers may be interested in Table S4. It would be helpful to provide a brief description in Results (lines 508-512) and list it explicitly in Data Availability.\n24. Line 846: There are 14 columns in Table S4, not 15 columns. This is also a good place to list Table S4.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Responses to the Reviewer Comments\nReviewer #1\nComment 1-0\nMeng et al developed a generic chromatin accessibility reference by aggregating peaks from 624 high-quality bulk ATAC-seq datasets, defining more than 1 million consensus peaks (cPeaks). They demonstrated cPeaks outperforms existing feature defining methods and references by improving cell annotation and rare cell type identification, and provide insights into chromatin dynamics during cellular differentiation and tumor progression. The collected large consortium of ATAC-seq would be a good resource in studying gene regulation. I have some comments as follows,\nReply:\nThank you so much for providing the excellent summary of our manuscript. We appreciate your feedback and have addressed all the points raised, as detailed in the responses below.\nComment 1-1\nHow the batch effect like ATAC-seq generated from different labs, protocols, etc are corrected before the integration?\nReply:\nThank you for the question. We have considered the potential batch effects in defining cPeaks. We apologize for not providing the relevant details in the previous version of the manuscript.\nFirst, cPeaks are not derived from raw sequencing reads, which are indeed subject to batch-to-batch variability, but from pre-identified peaks that have already undergone peak calling within each dataset. This peak calling process effectively filters out substantial noise and retains biologically meaningful accessibility signals, thereby minimizing the influence of batch effects at the source.\nWe also evaluated the robustness of our approach by generating cPeaks separately from different batches. We have demonstrated that the shape profiles of merged peaks remain highly consistent across diverse experimental conditions. We observed strong concordance in peak shape across tissues of different origins, across sequencing technologies such as ATAC-seq and DNase-seq, and across different peak calling methods and parameters. These results suggest that the merged peaks represent intrinsic genomic features rather than artifacts driven by experimental variation.\nIn this revision, we reorganized the presentation of these results and added additional analyses, including the distributions of shuffled peak pairs and 250 bp–shifted peak pairs, to help readers better interpret the consistency and robustness of the peak shapes (Figure R1). We added the results in Line 119-141 and Figure 1d. We also expanded the discussion on this point in the revised manuscript to clarify the robustness of our method to batch effects (Line 144-148).\nFigure R1 Correlations of peak shapes between different conditions.\nComment 1-2\nFor the developed CNN model is claimed to be modified based on scBasset. However, the difference and similarity of CNN and scBasset is unclear, and many technical details of CNN is missing. scBasset aims to predict the peaks for all cells in the single-cell setting, while the CNN is predict peaks for multiple tissues/organs. Please clarify this point and how the loss function is defined. Moreover, when the CNN performs zero-shot prediction, where training set does not include the tissue/organ/cell types in the testing set, please clarify the dimension of output labels in the training set and testing set as supervised learning need the same dimension of training and testing sets. Furthermore, are the training, validation and testing set split by different chromosomes? It is highly encouraged the detailed mathematical formula need to be present in the method section.\nReply:\nThank you for the question. We apologize for not providing a detailed description of the CNN model structure in previous manuscript.\nscBasset aims to predict cell type-specific chromatin accessibility. Each output dimension in scBasset corresponds to number of cell types. In contrast, our goal is to predict the intrinsic chromatin accessibility potential of a genomic region. Rather than focusing on cell type- or tissue-specific accessibility patterns, we aim to identify a universal set of consensus peaks (cPeaks) shared across diverse biological contexts.\nTo achieve this, we modified the architecture of scBasset accordingly. We retained the core one-dimensional convolutional structure to capture DNA sequence patterns but simplified the output to a binary classification task. The model predicts whether a given genomic region corresponds to a cPeak or not. The output layer contains two dimensions representing probabilities of the two classes (cPeak / non-cPeak). We use the standard Binary Cross-Entropy Loss (BCELoss) as the training objective:\n\"\n1\n!=− %[' log( +(1−')log(1−()]\n$ ! ! ! !\n!#1\nWhere N is the number of samples, ' ∈{0,1} is the ground truth label of the i-th sample, and ( is\n! !\nits predicted probability of class 1 (i.e., is a cPeak).\nThis whole design enables a consistent label structure between training and testing, which makes the model easily extensible to new, unseen tissues or cell types. When simulating a zero-shot scenario by holding out the tissue- or cell-type-specific cPeaks from training and using them as the test set, as the label space remains binary, the model does not require any structural adjustment for prediction on unseen biological contexts.\nIn the revised version, we added a leave-one-chromosome-out approach to further evaluate the robustness of the model (Figure R2). Specifically, for each held-out chromosome, we trained the model on all other chromosomes by splitting peaks into training, validation, and test sets with an 8:1:1 ratio. The held-out chromosome was then used as an additional test set. We found that the model consistently achieved high AUC scores (~0.95) on most left-out chromosomes except for chrY, closely matching performance on the original test set. When holding out chrY, the model achieved a slightly lower AUC (~0.90), suggesting that the relationship between DNA sequence and chromatin accessibility on chrY may differ from that on other chromosomes. This observation is consistent with prior reports indicating that chrY harbors a complex repeat structure compared to other chromosomes1,2.\nWe added a detailed description of the model architecture, loss function, and training/test design in Line 218-223, 766-769. We added the chromosome-based evaluation in Line 243-249 and 779-782 and Figure 2d.\nFigure R2. The AUC for distinguishing peaks on the original test set derived from the remaining chromosomes and for peaks on the held-out chromosome.\nNature, 2003, 423(6942): 825-837.\nComment 1-3\nFurther comments have for using cPeaks for facilitating the discovery and analysis of rare cell types with cPeaks or Investigating cell state transitions during development with cPeaks. It is unclear which peak set of CPeaks for the data analysis, tissue-matched for the single-cell data or all peaks? What is the peak overlap between scATAC peak-calling and the peak set provided by CPeaks in the data analysis? It is also unclear additional peak information provided by CPeaks with the same read information can achieve both tasks. In another word, what is the intuition to using the mapped reads from studied scATAC-seq data in a broad collection of peaks across more tissues and cell types, provided cPeaks, would improve the identification of rare cell types?\nReply:\nThank you for the comment. We are sorry that the writing in the previous version of the manuscript caused misunderstandings regarding the usage of cPeaks. We would like to clarify this point as follows.\n(1) Peak set used for analysis\nWe employed the complete set of cPeaks (observed cPeaks + predicted cPeaks), without restricting to tissue-matched or dataset-specific subsets in all analyses using cPeaks throughout this study, including the identification of rare cell types (e.g., Fig. 5) and the investigation of developmental transitions (e.g., Fig. 6). This is a key motivation for defining cPeaks in the first place: to establish a unified, high-confidence reference derived from diverse adult tissues, which enables researchers to map newly generated scATAC-seq reads directly to a common peak set. We added more detailed clarification in the revised manuscript (Line 260-261).\n(2) Overlap with standard scATAC peak calling\nWe agree that the overlap between cPeaks and dataset-specific peaks is an important factor. We quantified this in two representative analyses. In the PBMC-8k analysis, MACS2-called pseudo-bulk peaks from all PBMCs totaled 133,648 peaks, of which 98.74% were covered by cPeaks. In the GI epithelial cell analysis, the fetal retina pseudo-bulk peak set contained 132,172 peaks, with 98.95% covered by cPeaks. These results indicate that cPeaks effectively encompass most accessible regions identified through conventional peak calling.\n(3) Rationale for improved performance using cPeaks\nThe advantage of using cPeaks, particularly for analyzing rare cell types, lies in avoiding the limitations of pseudo-bulk peak calling. Traditional scATAC-seq analysis pipelines typically require aggregating cells into pseudo-bulk profiles before calling peaks. This step tends to overlook signals from rare cell types, as their accessibility patterns are often drowned out by more abundant populations and treated as noise. Accurately recovering such rare subtypes often requires iterative clustering and peak refinement.\nIn contrast, using cPeaks provides a comprehensive set of high-confidence potential accessible regions. We can quantify reads over a consistent reference from all cells, including rare subpopulations. This design avoids the information loss and bias introduced during dataset-specific peak calling, leading to improved resolution of rare cell types. In this revision, we revised the beginning of the relevant Results section to better highlight this motivation (Line 426-432).\nAlso, we added an additional benchmark analysis included in the revision, to support our explanation. Specifically, we constructed cell-by-feature matrices using five different references: cPeaks, cDHSs, CATLAS cCREs, 2,000-bp genomic bins, and pseudo-bulk peaks (called from the same dataset). We quantified clustering accuracy using five metrics, and cPeaks consistently achieved the best or second-best performance across all metrics (Figure R3). In addition, when transferring scRNA-seq labels to scATAC-seq data, cPeaks yielded higher prediction scores than all other references (Figure R4), further supporting their effectiveness in improving annotation quality. Furthermore, to assess the utility of cPeaks in downstream analyses, we examined differentially accessible regions identified using each reference set. Taking pDCs as an example, we compared their accessibility profiles to all other cells. We found that cPeaks enabled the identification of slightly more differentially accessible regions than cDHSs, and substantially more than CATLAS cCREs and pseudo-bulk peaks (Figure R5). All the results support the use of cPeaks for facilitating the analysis of rare cell types. These results are included in the revised manuscript (Line 440-462, Fig. 5b-c, Fig. S4c).\nThank you again for raising these points.\nFigure R3 Performance of clustering accuracy using different reference peak sets in the PBMC-8k dataset.\nFigure R4 Distributions of the prediction scores of the cross-modality integration between the scRNA-seq data and scATAC-seq data analyzed with different reference sets. The prediction score reflects the confidence of predicted labels.\nFigure R5 The Venn plot of the number of differential accessible regions calculated by cPeaks and other reference sets for pDCs.\nComment 1-4\nThe github seems only provide the pipeline and software code. Would the CPeak data resources (peaks across all tissues/organs) can be provided in a format of online database as well as standard-alone download package? This will facilitate the users from both biology background and computational modeling.\nReply:\nThank you for the question. In our previous manuscript, besides pipeline and software code, we have provided the detailed information of cPeaks, including genome location, annotation, etc. Our goal was to provide a consistent reference peak set rather than tissue-specific accessibility data. Due to the limited dataset and technical noise, the data may not fully reflect the true chromatin accessibility patterns within each tissue. So, we did not include cPeak-by-tissue resources you mentioned. However, we recognize that such resources are indeed valuable, and Reviewer 2 also raised a similar point. Therefore, we have now included the cPeak-by-tissue data as a supplementary resource. We also emphasize in the manuscript that this dataset reflects only the observed signals, which may contain noise and may not fully represent the true accessibility landscape within each tissue.\nSpecifically, we have uploaded a matrix that records the presence of each cPeak across all datasets used in our study, available at:\nhttps://zenodo.org/records/15666169/files/cpeaks-by-dataset.tsv.gz\nTo further improve usability for users interested in organ- or tissue-level summaries, we aggregated the dataset-level matrix by tissue type to generate a tissue-by-cPeak matrix, which is available at:\nhttps://zenodo.org/records/15666169/files/cpeaks-by-tissue.tsv.gz\nThese matrices reflect the distribution of observed cPeaks in the datasets analyzed in this study. While they do not constitute a complete atlas of chromatin accessibility across all human tissues, they provide a practical reference for analyses or modeling.\nWe added descriptions of these datasets in the revised manuscript (Line 1043-1049). We will continue working to make the cPeaks data more broadly accessible through additional database interfaces in future updates.\nReviewer #2\nComment 2-0\nThis manuscript presents cPeaks, a reference list of chromatin-accessible regions created by merging peaks from 624 bulk ATAC-seq datasets. The authors aim to provide a shared peak set for scATAC-seq analysis and employ a deep learning model trained on DNA sequence to predict additional peaks in unobserved tissues or cell types. While the study addresses an important challenge in the field, it has three major limitations. First, the cPeak list is not tissue- or cell-type-specific, which limits its applicability for studies requiring context-specific regulatory information. Second, using only DNA sequence to predict chromatin accessibility overlooks critical biological factors such as transcription factor binding and chromatin context, raising doubts about the validity of the predicted peaks. Third, the CNN model is trained on data from known tissues, and it remains unclear how it can reliably predict peaks in tissues it has not encountered during training. The detailed comments are as below.\nReply:\nThank you very much for the comments. We would like to clarify the following:\n(1) On context specificity: The main purpose of cPeaks is to provide a universal, robust reference set to support general scATAC-seq analysis. While not tissue-specific, cPeaks demonstrate strong performance across diverse datasets, including developmental and cancer contexts. Nevertheless, we now additionally provide tissue-level cPeak matrices to support context-aware analyses.\n(2) On using only DNA sequence: Our CNN model is intentionally designed to predict the potential for chromatin accessibility based on intrinsic sequence features, not to capture context-dependent activity. As shown in our model interpretation (e.g., NeuronMotif analysis), the model captures biologically meaningful motifs.\n(3) On prediction in unseen tissues: We evaluated the model’s generalization using held-out tissue-specific peaks as test sets and achieved high prediction performance. The model predicts whether a region could be accessible, but not in which tissue it will be so. We have clarified this in the revised manuscript. Also, the CNN predictions were not fed back into training, avoiding circular logic.\nWe appreciate your suggestions and have revised the manuscript to better clarify the design and intended use of cPeaks and the CNN model.\nFor detail, please see our point-by-point responses below.\nComment 2-1\nThe cPeaks are derived by merging peaks identified using various peak-calling tools, each with different characteristics and sensitivities. This can introduce substantial technical variation. Although the authors report shape similarity across methods, the Pearson correlation coefficients range from 0.4 to 1.0, indicating inconsistency. The authors should clarify how they account for this heterogeneity and whether peak-calling parameters were harmonized prior to merging.\nReply:\nThank for the comment. All cPeaks were derived from peaks identified using the same peak-calling tool, MACS2. Sorry that we didn’t write it clearly in previous manuscript, and we added the clarification in cPeak generation part.\nFor the report of shape similarity across methods, referring to Fig. 1d, we compared peak shapes generated by MACS2 (used in our pipeline) with those produced by CellRanger, to examine whether cPeak shapes are sensitive to the choice of peak-calling method. While some PCCs were relatively low, the overall median PCC remained around 0.665, indicating that cPeak shapes are generally stable across different algorithms. As we only evaluated a single sample produced by CellRanger, due to the limited number of cell types, the resulting merged peaks had relatively low resolution, which likely contributed to the lower observed PCC. Nevertheless, the shape similarity remained significantly higher than controls, supporting the robustness of peak shapes (Figure R6).\nIn addition, the ENCODE and ATACdb datasets were processed using different parameter settings of MACS2. We compared the cPeak shapes generated from these two sources and found high concordance, with a median PCC of 0.868 (Figure R6). Results show that cPeak shapes are largely robust to parameter variations within MACS2.\nTo improve clarity and reduce potential confusion, we revised the relevant descriptions in the manuscript (Line 110-141). Besides, we updated some words to explicitly state that all peaks used for generating cPeaks were called using MACS2 (Line 133-134). Moreover, we have split Fig. 1d into two parts (Figure R6): one showing comparisons internal to our pipeline (e.g., ENCODE vs. ATACdb with different MACS2 parameters), and the other showing external comparisons (e.g., MACS2 vs. CellRanger), to help readers better understand the robustness across different peak-calling methods.\nThank you again for helping us improve the presentation of our work.\nFigure R6 Correlations of peak shapes between different conditions (new Fig. 1d).\nComment 2-2\nWhile the merged cPeaks serve as a universal reference, they may obscure important tissue- or cell-type-specific regulatory elements. Including annotations that flag cPeaks unique to specific contexts would enhance the utility of the resource. Alternatively, providing separate cPeak lists for each tissue or cell type would support more targeted and biologically meaningful analyses.\nReply:\nIn previous manuscript, our goal was to provide a consistent reference peak set rather than tissue-specific accessibility data. Due to the limited dataset and technical noise, the data may not fully reflect the true chromatin accessibility patterns within each tissue. So, we did not include cPeak-by-tissue resources you mentioned. However, we recognize that such resources are indeed valuable, and Reviewer 1 also raised a similar point. Therefore, we have now included the cPeak-by-tissue data as a supplementary resource. We also emphasize in the manuscript that this dataset reflects only the observed signals, which may contain noise and may not fully represent the true accessibility landscape within each tissue.\nWe have now provided two additional matrices:\nl A dataset-level matrix showing the presence or absence of each cPeak across the 624 bulk ATAC-seq datasets used in our construction: https://zenodo.org/records/15666169/files/cpeaks-by-dataset.tsv.gz\nl A tissue-level matrix obtained by aggregating the dataset matrix by organ or tissue type: https://zenodo.org/records/15666169/files/cpeaks-by-tissue.tsv.gz\nThese resources allow users to identify cPeaks that are broadly shared or enriched in specific tissues, supporting both general and context-specific analyses. We note that these annotations are based on the datasets included in this study and may not cover the full diversity of human tissues. We have updated the manuscript to describe these additions (Line 1043-1049), and we plan to expand more datasets and annotations in future versions of the resource.\nComment 2-3\nThe deep learning model used to predict additional cPeaks is trained solely on DNA sequence, without considering key biological factors that influence chromatin accessibility, such as transcription factor binding, histone modifications, and chromatin context. As a result, it is unclear whether the predicted peaks truly capture cell-type-specific regulatory activity. Moreover, the model is trained on cPeaks derived from known tissues, yet the manuscript does not provide sufficient justification for its ability to generalize to unseen tissues. Finally, the prediction framework presents a potential circular logic issue: newly predicted peaks could be added to the shared peak list and used for further model training, allowing the process to expand indefinitely without clear stopping criteria or external validation. The authors should clarify how they prevent such iterative propagation of prediction error.\nReply:\nThank you for the comment. We apologize that our original description may have led to misunderstandings regarding the design and scope of our prediction framework. We clarify the following points as below:\n(1) Considering key biological factors that influence chromatin accessibility\nOur prediction task aims to identify genomic regions with an intrinsic potential for chromatin accessibility, rather than predicting accessibility status in specific cell types or tissues. While factors such as transcription factor binding, histone modifications, and chromatin context influence accessibility in a context-dependent manner, many accessibility patterns are encoded in the DNA sequence itself. Therefore, our deep learning model is intentionally based only on sequence features, which remain constant across conditions. Recent studies, such as DeepSEA and Enformer, have shown that some of these regulatory factors can, to a certain extent, be inferred from sequence. We thus believe that a well-trained CNN model can capture these patterns and accurately predict regions with the potential to become accessible.\nIn this revision, we further explored the interpretability of the CNN model using the NeuronMotif tool3 to identify sequence patterns that activate specific neurons in the network (Figure R7). By analyzing the last convolutional layer, we identified 342 TF motifs that could robustly activate at least one neuron. Notably, motifs from the bZIP family and CTCF/CTCFL were found to be highly influential across many neurons. These results suggest that the model has captured biologically meaningful features based on DNA sequence alone.\nTo improve clarity, we added a more explicit explanation of the model’s objective and underlying assumptions (Line 214-217). We also incorporated the model interpretation results into the revised manuscript (Line 250-254, Fig. S1h).\nPNAS, 2023, 120(15): e2216698120. https://doi.org/10.1073/pnas.2216698120\nFigure R7 TF motifs discovered by NeuronMotif in the last convolutional layer of the CNN model for predicting cPeaks. The x-axis represents the proportion of neurons that can be activated by the TF motif. The y-axis represents the average frequency of the TF motif in the activated neurons.\n(2) Generalizing the model to unseen tissues\nIn the original manuscript, to evaluate the model’s ability to generalize to unobserved tissues or cell types, we conducted leave-one-context-out tests. We excluded cPeaks specific to a given tissue or cell type from the training set and used them as an independent test set. As shown in Figure 2c, the model was able to recover these held-out cPeaks with high performance, demonstrating its ability to learn general sequence features associated with chromatin accessibility from limited tissues.\nWe added further explanation of this evaluation strategy in the revised manuscript (line 229-236) to give a better explanation of the design of our generalization experiment.\n(3) potential circular logic issue\nWe apologize that our previous description was unclear and caused misunderstanding. Our original framework did not involve any iterative expansion using CNN-predicted cPeaks. Specifically, we did not incorporate predicted peaks back into the training set at any stage. The CNN model was trained once using the observed cPeaks, and the prediction was performed in a single forward pass. To prevent misunderstanding, we have revised the manuscript to explicitly clarify this point (Line 793-795). Thank you again for helping us better present the design of deep-learning model.\nComment 2-4\nIn benchmarking the coverage of scATAC-seq reads, the authors use MACS2-called peaks as ground truth. Since cPeaks may be partially derived from MACS2-called datasets, this could introduce bias in favor of cPeaks. Moreover, if MACS2 output is considered the gold standard, it raises the question of why not simply recommend MACS2-based peak calling per dataset. The authors should report the proportion of cPeak-contributing datasets that used MACS2. If MACS2 is dominant, it should be treated as a baseline method—like cDHS and CATLAS cCREs—rather than as the default reference. A side-by-side comparison of performance across all methods would be more appropriate.\nReply:\nThank you for the comment. We apologize for the lack of clarity in the manuscript, which led to a misunderstanding. All cPeaks were generated using MACS2. In the benchmarking experiment shown in Fig. 1h, we used peaks called by MACS2 on pseudo-bulk scATAC-seq data as a practical baseline. This strategy is the common practice in current scATAC-seq analysis pipelines.\nOur goal in this experiment was to evaluate how well different predefined accessible region sets (cPeaks, cDHS, and CATLAS cCREs) recover accessibility signals in real scATAC-seq data, and thereby assess their utility as reference sets. To address your concern regarding potential bias, we examined the processing steps of each reference. All datasets used to generate cPeaks were indeed originally processed using MACS2, as were those used in the CATLAS cCREs. In contrast, the cDHS reference was derived using Hotspot2. We therefore believe the comparison is fair, as MACS2 was used across both cPeaks and CATLAS, and thus any advantage is not unique to cPeaks.\nFurthermore, the MACS2-called pseudo-bulk peaks used as the comparison baseline in Fig. 1h were generated from an independent scATAC-seq dataset not used in cPeak construction. We found that cPeaks captured a similar proportion of reads as the MACS2 pseudo-bulk peaks, while substantially outperforming CATLAS cCREs and cDHSs. This suggests that cPeaks generalize well beyond the datasets used in their construction and provide a robust and comprehensive representation of accessible regions. We have clarified these points in the revised manuscript (Line 179-186, 681-683).\nLastly, we agree that for analyzing individual homogeneous datasets, MACS2 (or other conventional peak calling methods) is generally sufficient. However, our motivation for constructing cPeaks to facilitate large-scale integrative analysis across heterogeneous samples and improving performance in downstream analysis, such as rare cell detection. These advantages are demonstrated in the later benchmarking and case study sections of the manuscript.\nComment 2-5\nThe reliability of the curated cPeaks would be strengthened through validation with independent chromatin accessibility assays, such as ChIP-seq, CUT&RUN, or CUT&Tag. Such validation would help confirm that cPeaks represent biologically relevant regulatory elements.\nReply:\nThank you for this valuable suggestion. To assess whether cPeaks represent biologically meaningful regulatory elements, we validated their overlap with an independent ChIP-seq dataset from ReMap2022. Specifically, we downloaded the non-redundant ReMap peak set, and analyzed three pre-defined accessibility references: cPeaks, cDHS, and CATLAS cCREs, and calculated the proportion of ReMap peaks each covered. For each reference, we also generated a background peak set with matching peak length distribution and total peak number, excluding the original reference and ENCODE blacklist regions. Enrichment was assessed by comparing observed versus background overlap ratios.\nOur results (Figure R8) showed that cPeaks covered 88.7% of ReMap peaks, slightly higher than cDHS (88.2%) and substantially more than CATLAS cCREs (72.9%). Relative to their respective background peak sets, the enrichment fold-changes were 13.5 (cPeaks), 7.2 (cDHS), and 8.4 (CATLAS), supporting the strong regulatory relevance of cPeaks.\nTo further assess specificity, we evaluated overlap between each reference and the individual TF ChIP-seq peak sets in ReMap. Among all the TFs evaluated, cPeaks showed greater coverage than cDHS for 54.8% of TFs, and greater coverage than CATLAS for 99.8% of TFs (Figure R9). On average, cPeaks covered 87.5% of TF peaks (median = 92.6%), slightly outperforming cDHS (mean = 87.0%, median = 92.0%) and clearly exceeding CATLAS (mean = 73.0%, median = 77.9%).\nThese results collectively support that cPeaks are highly enriched for biologically relevant regulatory regions and provide a more comprehensive representation than existing references.\nWe have revised the manuscript accordingly (Line 297-302, Figure S2a-b). Thank you again for your suggestion.\nFigure R8 The proportion of each reference set covering ReMap ChIP-seq peaks.\nFigure R9 The proportion of each reference set covering peaks of each individual TF in ReMap ChIP-seq data. Each dot represented a TF.\nReviewer #3\nComment 3-0\nIn this manuscript, the authors assembled 1.4 million observed consensus peaks (cPeaks) on chromatin openness from bulk ATAC-seq data in ENCODE and ATACdb and developed a deep CNN model to predict an additional 0.28 million cPeaks across the human genome. They further annotated the cPeaks and demonstrated utilities of the cPeaks in labeling scATAC-seq data, identifying “rare” cell types, investigating cell state transitions during development, and characterizing tumor cell subtypes. On the positive side, the topic is interesting, the assembled cPeaks can be a valuable resource to the field, and most Results sections are well written. On the negative side, crucial details are missing in Methods sections, some benchmarking analyses are neither thorough nor convincing, and some statements are over stretched. The manuscript needs a significant improvement.\nReply:\nThanks for the summary and comments. We appreciate the recognition on the positive sides of our work. We carefully studied all the pointed negative issues and revised the work accordingly both in the writing and in some experiments: We added missing methodological details, refined several benchmarking analyses to make them more rigorous and convincing, and removed or rephrased overstatements in the previous writing. We are grateful for your questions that really helped to significantly improve the manuscript.\nComment 3-1\nThe authors did not provide enough details on how Pearson correlation coefficient (PCC) is calculated. It can be easily misused and inflates similarity. What are the start and end points for PCC calculations when comparing two peaks? For example, comparing the peak in Fig. 1c with one shifting to the right by 250bp, how will the PCC be calculated? What is the PCC value?\nReply:\nThank you for the suggestion. We apologize for the lack of clarity in our original description of the Pearson correlation coefficient (PCC) calculation. We provide a detailed clarification below and in revised manuscript regarding how PCC was calculated in our analysis, including the genomic window used and how shifts in peak position affect the correlation.\nIn our analysis, we calculated the PCC between two peak shapes using the genomic region defined by the merged peak boundaries from ENCODE as a fixed reference window. For any other peak set being compared, we aligned the signal within the same window. If the compared peak was shorter than the reference window, we padded the missing positions with zeros. If it extended beyond the reference window, we truncated the excess. This ensured that all comparisons were made over a consistent length and position.\nFollowing your suggestion, we conducted an additional test to examine the impact of positional shifts. We introduced a 250 bp shift to peaks and recalculated the PCC values. The results, shown as green density plots in Figure R10, demonstrate that the PCC decreases under such shifts, as expected (We changed the violin plot in the original version to the density plot to better illustrate these distributions in limited space). This supports the sensitivity of our method to spatial misalignment and clarifies how PCC behaves under positional variation.\nWe have added this clarification and the new analysis to the revised manuscript (Line 110-148, 678-712, new Fig. 1d). Thank you again for your helpful suggestion to help us improve the evaluation of cPeak shapes.\nFigure R10 Correlations of peak shapes between different conditions (new Fig. 1d).\nComment 3-2\nDetails on permutation tests need to be provided and the PCC distribution from 1000 permutations should be provided. Otherwise, Fig. 1d is not informative. The authors may want to convert the PCC values in Fig 1d to empirical p values, using the distribution from the permutation tests as the null distribution, and show the histogram of the empirical p values.\nReply:\nThank you for the suggestion. Following your recommendation, we revised and improved our permutation testing procedure to better assess the significance of peak shape similarity. This updated procedure differs from the original version and is now described in detail in the revised Methods section (Line 703-710).\nWe re-designed the generation of shuffled controls. We now randomly pair each peak in the baseline group with a different peak from the comparison group, align them by center position, and harmonize their lengths by trimming or zero-padding to match the baseline window before calculating the Pearson correlation coefficient (PCC).\nTo improve interpretability, we have added the distribution of the permuted PCC values in Figure R10 (the revised Fig. 1d). We also converted the observed PCC values into empirical p-values using the permutation-derived null distribution (Figure R11). In addition, we calculated empirical p-values for the shifted-250bp peaks mentioned in your Comment 3-1 to further illustrate how peak misalignment affects correlation significance. The results supported the high consistency of shapes across different conditions.\nThese updates are included in the revised Fig. 1d, Fig. S1b and described in the corresponding text (Line 110-129, 703-710). Thank you for helping us enhance the statistical rigor of this analysis.\nFigure R11 Empirical p-values of PCCs for true pairs of merged peaks, using shuffled controls (green) and shifted controls (blue) as null distributions. Colored numbers indicate the proportion of true peak pairs with p-value < 0.05 under each control.\nComment 3-3\nThe similarity analysis should be performed on consensus peaks (cPeaks), not on multimodal peaks. Users care more about the end product than the intermediate results. Basically, the order between lines 120-126 and lines 104-119 should be switched.\nReply:\nThank you for the suggestion. Our previous logic was to first show that the intermediate peaks used in the merging process were stable under different conditions so that the feasibility of constructing a consensus peak set is justified, and then to introduce the way of constructing the cPeaks. We agree that users are more interested in the final cPeaks and that evaluating their consistency is important. So, we kept the original logic flow but revised the paragraph according to this comment and the next one, and we added more results on the final cPeaks after this paragraph to better serve the readers’ interests.\nWe agree with your suggestion that the similarity analysis should include the final cPeaks. In the revised manuscript, we performed an additional shape consistency analysis on the final cPeaks and included the results (Figure R12). The results were highly consistent with those obtained using the intermediate peaks. We also revised the presentation of the relevant sections (Line 155-158, Fig. S1e-f) to make the logic of our workflow clearer and to help readers better understand both the process and outcome. Thank you again for helping us improve the logic and clarity of this section.\nFigure R12 Left: Correlations of shapes of observed cPeaks between different conditions. Right: Empirical p-values of PCCs for true pairs of observed cPeaks, using shuffled controls (green) and shifted controls (blue) as null distributions. Colored numbers indicate the proportion of true peak pairs with p-value < 0.05 under each control.\nComment 3-4\nWhy is there no evidence for footprinting from transcription factor binding in Fig. 1c or Fig. 1f?\nReply:\nThe purpose of Fig. 1c and Fig. 1f was to analyze the shape and consistency of peak regions, rather than to examine raw read distributions. As these analyses were conducted at the peak level, without access to the underlying read-depth profiles, we were not able to perform transcription factor footprinting, which requires base-resolution read information. We agree that footprinting could provide additional insights, and we will consider integrating such analyses in future work when enough read-level data are available.\nComment 3-5\nFig. 4: This benchmarking analysis is very weak: a) The dataset of 2034 cells is very small and may miss complexities from more complicated datasets. The authors should include additional datasets generated from DOGMA-seq (Mimitou et al. 2021) and/or TEA-seq (Swanson et al. 2021) for benchmarking. b) The authors should provide more details on “clustering and annotation” (such as how many clusters and how many cell populations were used for benchmarking) and show a few examples of UMAPs.\nReply:\nThank you for the suggestions. We have addressed both points as follows:\n(a) To strengthen the benchmarking analysis, we included an additional dataset, PBMC-DOGMA, generated using the DOGMA-seq protocol (Mimitou et al., 2021). This dataset comprises 7,623 cells annotated into 8 immune cell types, providing a larger-scale system for evaluation. We applied the same feature extraction and evaluation framework to this dataset and found that cPeaks consistently ranked among the top-performing feature sets across multiple embedding methods and evaluation metrics (Figure R13 and Figure R14), supporting the generalizability and robustness of our approach across datasets and platforms. We added the results in the revised manuscript (Fig. S3f-g, Line 409-414).\nFigure R13 Performance of different reference sets in the PBMC-DOGMA dataset.\nFigure R14 UMAP plots of the PBMC-DOGMA dataset analyzed by different reference sets with Scaset, using 10,000 HVFs.\n(b) To improve clarity, we added more details about the clustering and annotation procedures in the Methods section (Line 868-871). We added figures to describe the distributions of cell types in both datasets (Figure R15, or Fig. S3a and S3d in the revised manuscript). We also included representative UMAP visualizations for both datasets (Figure R14 and Figure R16, or Fig. S3c and S3g in the revised manuscript) to illustrate the clustering performance.\nWe thank you again for these constructive suggestions, which helped improve the depth of our benchmarking analysis.\nFigure R15 The proportion of each subtype in the hematopoietic and PBMC-DOGMA datasets.\nFigure R16 UMAP plots of the hematopoietic dataset analyzed by different reference sets with Scaset, using 50,000 HVFs.\nComment 3-6\nThe benchmarking on PBMC dataset in Fig. 5 is very problematic: a) The authors benchmarked against five different feature (peak) sets in Fig. 4, why did they limit it to pseudo-bulk peaks here? It is not obvious other feature sets won’t be able to identify the three “rare” cell types. b) There are much larger datasets in the literature with tens to hundreds of thousands of PBMCs. Why did the authors benchmarked on a small dataset of 8728 cells? c) Did they consider the existence of doublets? d) The “rare” cell types (pDC, DC, NK bright) are not that rare and are thus poor examples for “discovering” rare cell types.\nReply:\nThank you for the comments. We apologize for not including a more comprehensive benchmarking analysis and larger-scale PBMC data in the previous version of the manuscript. In the revised version, we have added the relevant experiments to address these concerns.\n(a) Benchmarking with other feature sets\nIn the revision, we have extended the analysis to include additional reference sets (genomic bins, CATLAS cCREs, and cDHSs). We did not include the united peaks approach in this section, as it requires pseudo-bulk peak calling for each cell type using known labels. Such high-confidence ground truth labels are not available for these datasets, unlike the hematopoietic dataset used in our earlier benchmarking.\nUMAP visualizations showed that most feature sets successfully resolved major cell types, including three rare subtypes: plasmacytoid dendritic cells (pDCs), dendritic cells, and NK bright cells, each representing less than 1% of total cells (Figure R17). Genomic bins, however, failed to separate these rare subtypes. We then evaluated them using multiple clustering accuracy metrics (Calinski-Harabasz Index, Silhouette Score, cLISI, balanced KNN accuracy, and Davies-Bouldin Index), inspired by your Comment 3-7 (Figure R18). We also investigated the prediction scores when transferring scRNA-seq labels to scATAC-seq data (Figure R19). We calculated the number of differentially accessible regions of pDC cell type identified using each reference set and draw Venn plots (Figure R20). Compared to other reference sets, cPeaks enabled the identification of more differentially accessible regions (genomic bins were excluded from this comparison because they failed to cluster pDCs as a distinct population). In all these analyses, cPeaks outperformed other reference sets, supporting its utility in analyzing rare cell types. We added the results in Fig. 5b, 5c, S4b, S4c, and Line 433-453.\nFigure R17 The UMAP plots of the dataset analyzed by each reference set.\nFigure R18 Performance of clustering accuracy using different reference peak sets in the PBMC-8k dataset.\nFigure R19 Distributions of the prediction scores of the cross-modality integration between the scRNA-seq data and scATAC-seq data analyzed with different reference sets. The prediction score reflects the confidence of predicted labels.\nFigure R20 The Venn plot of the number of differential accessible regions calculated by cPeaks and other reference sets for pDCs.\n(b) Benchmarking on larger datasets\nIn this study, we initially chose the 10X PBMC-8k dataset (8,728 cells) as it is widely used in benchmarking. To address concerns about dataset size and generalizability, we repeated the same analysis using a larger PBMC dataset from 10X Genomics (PBMC-10k, 10,412 cells). The results were consistent: cPeaks again ranked among the top in clustering accuracy metrics and enabled the identification of more differentially accessible regions compared to other references (Figure R21). These results support the robustness of cPeaks across datasets of different sizes and confirm that our findings are not limited to a specific small dataset. We added the results in Line 459-462 and the new Fig. S5.\nFigure R21 Rare cell type discovery using cPeaks on the PBMC-10k dataset. (a) The proportion of each cell type in the dataset. (b) The UMAP plots of the dataset analyzed by each reference set. Each dot represented one cell. Different colors represented different clusters. (c) Performance of clustering accuracy using different reference peak sets. (d) The Venn plot of the number of differential accessible regions calculated by cPeaks and other reference sets for pDCs.\n(c) Considering the existence of doublets\nIn our study, we followed quality control procedures to minimize their impact. Specifically, we excluded cells with an abnormally high total number of fragments, which are commonly indicative of doublets. This filtering step was applied prior to downstream analysis to reduce the presence of technical artifacts. Although we did not use a dedicated doublet detection algorithm, the resulting clusters showed clear separation and consistent label transfer from scRNA-seq, suggesting that any remaining doublets did not substantially affect our results. We have now added a clarification of this filtering criterion in the revised Methods section (Line 928).\n(d) Issue of “rare” cell types\nWe apologize for not clearly stating the definition of “rare” in the previous version of the manuscript. We consider “rare” to refer to cell types that appear at low proportions within a given sequencing dataset. This definition may not correspond to biological rarity. The sampling bias or technical limitations in sequencing may also result in lower capture rates for certain cell types. The limited representation of these “rare” cells poses challenges for downstream analysis. For example, in the PBMC-8k dataset used in our analysis, pDC cells accounts for less than 2% of the total cells, and they cannot be well separated by genomic bins as reference, as shown in the new Fig. S4b. According to our definition, these cells fall within the category of rare cell types that are the focus of our analysis.\nIn the revised version, we retain the use of “rare” for consistency with common usage in the field but added explanation to more precisely define our intent (Line 427 and 487-489).\nWe thank you again for helping us to better clarify the definition of “rare” cells in our manuscript.\nComment 3-7\nThe benchmarking on gastrointestinal (GI) epithelial cells in Fig. 5 is also problematic: a) Again, why did the authors limit the benchmarking to one feature set instead of five? b) The claim that “(p)seudo-bulk peaks in the original study failed to separate these subtypes” (lines 384-385) is not supported by the evidence in Fig. s5b. It is a purity issue, not a failure. c) It appears that there are impurities in the four rare cell types in Fig. 5h as well. The authors should quantify and compare impurity between methods. d) The improvements in Fig. 5i by cPeaks were small. Were they statistically meaningful?\nReply:\nThank you for the comments. We apologize for not performing a comprehensive analysis and quantification on the GI cell experiments. We have revised this part of the manuscript to address your concerns.\n(a) In the revision, we expanded the benchmarking to include five representative feature sets: CATLAS cCREs, cDHSs, pseudo-bulk peaks, 2,000-bp fixed-size genomic bins, and cPeaks. This comparison now provides a more comprehensive evaluation of performance. We have updated the corresponding figures (revised Fig. S6b-c, see Figure R22 and Figure R23) and the corresponding text (Lines 463-471). The results showed that cPeaks outperformed the other feature sets, demonstrating improved resolution of rare cell subtypes.\nFigure R22 The UMAP plot of subtypes in GI epithelial cells using different reference peak sets. Each dot represented one cell. Each color represented a cell type provided by the original study. The black circles highlight the four rare subtypes.\nFigure R23 Performance of clustering accuracy using different reference peak sets.\n(b-c) As you pointed out, the issue is one of purity rather than complete failure. In the revision, we refined the language and provided quantitative assessments of clustering accuracy for all reference sets using five metrics (Calinski–Harabasz index, Davies–Bouldin index, mean silhouette score, balanced KNN accuracy, and mean cLISI). These results are now shown in revised Fig. S6c (see Figure R23). We calculated the KNN accuracy of rare cell subtypes and found that cPeaks consistently achieved the highest purity across most subtypes, supporting their utility in rare cell type detection. These results are now shown in revised Fig. S6d (see Figure R24).\nFigure R24 Clustering KNN accuracy of four rare subtypes using different reference peak sets.\n(d) To determine whether the improvements observed with cPeaks in Fig. 5i are statistically significant, we performed repeated subsampling experiments on the CATALS dataset. Specifically, we randomly sampled 50% of the dataset 10 times and calculated clustering metrics for each feature set on each subset. cPeaks consistently achieved slightly better performance than CATLAS cCREs across most metrics in all three abundance groups, including those with low cell numbers (Figure R25). Though the improvement is slight, given that CATLAS cCREs were derived from the same dataset, these results suggest that cPeaks offer a robust and transferable feature set for analyzing diverse cell populations, including rare subtypes. We have added this analysis and its summary to the revision (Fig. 5i, Line 482-486), and we thank you for prompting us to perform this important validation.\nFigure R25 The performance of cell clustering on different groups of cell types, analyzed with different references. The significant levels are calculated by the paired t-test (p-values denoted as ****, < 0.0001; **, < 0.01; *, < 0.05; ·, < 0.1; ns, ≥ 0.1).\nComment 3-8\nThe presentation on Fig. 6 can be improved: a) Similar analysis with respect to Fig. 6ab can be done with other feature sets. What is the novelty here? b) Lines 451-453: A scatter plot of pseudo-time versus well-positioned-cPeak ratio or a plot like Fig. 6d can be informative here. c) Lines 454-458: The presentation on the two datasets should be separated and mirror the analysis on the retinal data. As is, it feels rushed and incomplete.\nReply:\nThank you for the suggestions. We have revised the corresponding section and addressed each of your points below:\n(a) We apologize for not explaining this clearly in the original manuscript, which led to a misunderstanding. Our intention was not to emphasize the novelty of cPeaks compared to other feature sets in terms of dimensionality reduction or clustering performance. Instead, our focus was on the proportion of well-positioned peaks.\nMore specifically, our primary objective in this section is to demonstrate that cPeaks, although constructed from adult chromatin accessibility data, can be effectively applied to fetal developmental datasets to capture dynamic regulatory changes across cell states. This cross-context applicability is non-trivial, as many existing references are typically context-specific. While we acknowledge that other reference sets such as cDHSs or CATLAS cCREs may also allow for similar downstream trajectory visualization, in the revised version we shifted the focus to a more distinctive aspect: the analysis of chromatin accessibility dynamics using the ratio of well-positioned cPeaks along pseudo-time. This approach leverages cPeaks to quantify regulatory changes in a way that other features cannot support. We clarified this motivation in the revised text (Line 522-524), and adjusted the surrounding analysis and figures, according to your suggestion (b) and (c), to better highlight this strength of cPeaks beyond standard clustering or UMAP-based assessments.\n(b) We generated scatter plots and lineage-wise visualizations of pseudo-time versus well-positioned-cPeak ratio, which confirmed our previous conclusion that this ratio increases during lineage transitions and decreases afterward (Figure R26). These new results have been included in Fig. 6d and detailed in Line 535-537. Thank you for this helpful suggestion to improve the presentation of results.\nFigure R26 The well-positioned-cPeak ratio changed along with pseudotime. Three lineages are shown separately. The black curve represents a smoothed trend estimated using a generalized additive model, and the shaded area indicates the 95% confidence interval of the fitted curve.\n(c) Following your suggestion, we have included more analysis results and figures to make this part more comprehensive. We added UMAP-based visualizations of both atlases annotated with well-positioned-cPeak ratios, and provided detailed examples for representative cell types (Figure R27). We found that fetal cells consistently exhibit higher well-positioned-cPeak ratios than adult cells, supporting our hypothesis that these peaks are associated with differentiation potential. We have rewritten the corresponding text (Line 538-542) and included a new panel (Fig. S7e-g) to clarify and strengthen this analysis.\nWe thank you again for your suggestions, which has significantly improved the clarity and completeness of this section.\nFigure R27 (e–f) UMAP plots of scATAC-seq data from the fetal and adult cell atlases. Cells are colored by dataset origin in (e) and by well-positioned cPeak ratios in (f). (g) UMAP visualization of selected clusters from (e). Cells are colored by dataset origin, cell type, and well-positioned cPeak ratio, with corresponding distribution plots shown on the right. The top three most abundant cell types from fetal and adult samples are labeled, respectively.\nComment 3-9\nLines 483-488: Provide p values to assess significance.\nReply:\nThank you for your suggestion. In this revision, we performed statistical significance analysis on all data shown here. However, on the TCGA data, possibly due to the very limited sample size and potential batch effects in bulk sequencing data, only a few cancer types showed significant trends. Given that our observation of changes in the well-positioned cPeak ratio was originally made at the single-cell level, we conducted the same analysis using an independent single-cell dataset of hepatocellular carcinoma4. We observed a similar trend as shown in Fig. 7c: the well-positioned cPeak ratio increases at intermediate stages and decreases at later stages of tumor progression (Figure R28). We have replaced the original TCGA analysis with this new result in the manuscript (Line 573-577, 1014-1017, Fig. 7c-d and S8c), as we believe it more reliably reflects our observation. We also discuss potential reasons for this phenomenon in the Discussion section (Line 632-639). We appreciate your comment, which prompted us to further refine and strengthen our findings.\nFigure R28 Well-positioned-cPeak ratios in the two cancer scATAC-seq datasets across different tumor stages.\nMinor concerns\nComment 3-10\nLines 45-50: The literature review is incomplete. It ignores a different approach by a recent publication (Rachid Zaim et al. 2024).\nReply:\nThank you for pointing this out. We added a citation and brief discussion of the study by Rachid Zaim et al. (2024) in the introduction (Line 52-56).\nComment 3-11\nIn addition to Table S1-2, summary pie charts showing tissue and cell type compositions can be informative.\nReply:\nThank you for this suggestion. We added summary pie charts showing the tissue and cell type compositions and included them in the revised manuscript (Figure S1a or Figure R29).\nFigure R29 Tissue and cell type compositions of the datasets for generating and evaluating cPeaks.\nComment 3-12\nFig. 1b: The curves grow linearly without reaching a plateau from 200 datasets. Descriptions should be accurate.\nReply:\nThank you for pointing this out. We have changed the description in the revised manuscript (Line 103) to clarify that as more datasets were integrated, the total number of peaks and their genomic coverage initially grew rapidly, followed by a gradual deceleration in growth.\nComment 3-13\nHow representative is Fig. 1g or 1h? Can the authors provide similar results based on more datasets? As is, it is hard to know whether it is “cherry picking”.\nReply:\nThank you for the question. We apologize for not clearly explaining the data in the main text, which led to a misunderstanding. Due to the limited number of available datasets, we were unable to provide more results.\nIn Figure 1g and 1h, we systematically downloaded all available DNase-seq and ATAC-seq datasets for GM12878 in H. sapiens (hg38) from the ChIP-Atlas database. After filtering for quality and limiting peak counts to a comparable range (30,000–100,000), we identified 65 usable datasets, including 62 ATAC-seq and only 3 DNase-seq datasets.\nDue to the limited number of high-quality, comparable DNase-seq/ATAC-seq datasets for the same cell type, constructing well-matched pairs was not straightforward. To ensure objectivity, we paired each DNase-seq dataset with the ATAC-seq dataset that had the closest number of peaks, without any manual selection. This resulted in three unbiased DNase–ATAC pairs.\nWe emphasize that the examples shown in Fig. 1g and 1h were drawn from these three randomly matched pairs and not selected based on performance. Therefore, we believe they are representative and not the result of cherry-picking. We have clarified this procedure in the revised manuscript (Line 174-175, 733-741) to make our analysis more transparent.\nComment 3-14\nThe CNN model predicts cPeaks based on genome sequence regardless of tissue types. In other words, predicted cPeaks may not be open in specific tissues. The authors should distinguish testing prediction accuracy on known tissue-specific cPeaks from assigning predicted cPeaks as open to specific tissues. The latter is beyond the CNN model. The present description is ambiguous when claiming “that the CNN model can effectively predict missing cPeaks for previously unseen tissues or cell types”.\nReply:\nThank you for the comment. We agree that our original wording may have caused confusion regarding the scope of the CNN model’s predictions. Our model is designed to identify genomic regions with an inherent potential for chromatin accessibility based solely on DNA sequence, independent of any specific tissue or cell type. It does not incorporate tissue-specific chromatin features or predict accessibility in a context-dependent manner.\nThe purpose of evaluating the model on held-out, tissue- or cell-type-specific cPeaks was not to assign predicted peaks to specific tissues, but rather to test the model’s generalizability: whether it can recover potentially accessible regions that were not present in the training data. These regions, though observed in particular tissues, are treated in our analysis as sequence-defined elements whose accessibility may depend on cellular context.\nIn response to your comment, we have revised the relevant sentences in the manuscript to clarify this point. The previous statement has been revised to: “the CNN model can effectively recover cPeaks that were not used during training, including those identified in tissues or cell types excluded from the training set” (Line 240-242). We hope this better reflects the actual scope of the model.\nComment 3-15\nFig. 3h: Remove the extra “the distribution of” from the legend.\nReply:\nThank you for pointing this out. We have revised the figure legend accordingly.\nComment 3-16\nFig. S2b: Provide a legend.\nReply:\nWe apologize for the oversight. We have added the legend to Fig. S2b: “Examples of known motifs enriched around the well-positioned edges of cPeaks.”\nComment 3-17\nLines 369-370: The anchor method in Seurat has limitations in labeling scATAC-seq data. The authors are encouraged to explore newer methods for labeling scATAC-seq data in the literature.\nReply:\nIn this study, our primary objective for cell type annotation was to establish a consistent and reasonable reference framework for downstream analysis, rather than to comprehensively benchmark annotation strategies. To this end, we adopted the anchor-based label transfer method implemented in Signac5, following the widely used Seurat pipeline6. While we acknowledge that the anchor method has certain limitations when applied to scATAC-seq data, it remains a practical and broadly adopted approach, particularly in settings where ground truth labels are unavailable.\nImportantly, our downstream analysis centers on plasmacytoid dendritic cells (pDCs), a well-defined population with distinct clustering and marker profiles. As such, we believe the annotation method used does not materially affect our conclusions.\nNonetheless, we fully agree that the field of scATAC-seq cell annotation is rapidly evolving, and we will actively explore and adopt more advanced annotation strategies in future analyses as appropriate. Thank you again for your suggestion.\nComment 3-18\nIt is unclear what is the difference between Fig. S5b and Fig. S5c. More details can be helpful.\nReply:\nThank you for the comment. Sorry that the difference was not clearly describe in the previous version of the manuscript.\nThe key difference between Fig. S5b (new Fig. S6b, CATLAS cCREs part) and Fig. S5c (new Fig. S6b, pseudo-bulk peaks part) lies in the choice of feature sets used for clustering. Fig. S5b uses CATLAS cCREs as input features. These cCREs were defined globally by merging and trimming peaks across all cell types. In contrast, Fig. S5c is based on peaks called from a pseudo-bulk of all cells in the GI dataset. This dataset-specific reference better captures context-specific signals and partially improves the resolution of rare cell types. As shown in Fig. 5h in the original manuscript, cPeaks achieve comparable performance in identifying rare subtypes without requiring dataset-specific peak calling.\nWe added a more detailed description (Line 463-471) and change the style of figure (Fig. S6b) to help better understanding the results.\nComment 3-19\nFig. 6 and Fig. S6: a) Lines 447-450: Results on enriched GO terms should be provided. b) A UMAP by life stage can be informative. c) All UMAPs should use the same projection and thus look similar and comparable.\nReply:\nThank you for the suggestions.\n(a) The GO enrichment results mentioned here are based on our analysis of TFs associated with well-positioned cPeaks in the original manuscript. To improve clarity, we have revised the text to explicitly reference the corresponding results in Fig. S2f (Line 529-530), which shows that these TFs are enriched in development-related terms.\n(b) As suggested, we have added a UMAP visualization colored by developmental stage to better illustrate the life-stage-specific distribution of cell states (Fig. S7e-f, Line 538-542).\n(c) We apologize that our previous wording caused some misunderstanding. The current figures include several distinct datasets: the retina dataset contains three separate groups (early, late-central, and late-peripheral), while the fetal and adult atlases are analyzed together. For each dataset, we performed UMAP separately to preserve the internal developmental or anatomical structure of the cells. We have revised the figure legends and text (Line 518-519) to clarify this point and guide the reader through the multiple analyses.\nComment 3-20\nLine 464: “Each box represents a single cell type …”. Did you mean each “line” here? How many total cell types were presented?\nReply:\nThank you for pointing this out. Yes, we meant to say, “each line represents a single cell type.” We have corrected this wording in the revised manuscript. A total of 111 cell types were presented for each of the fetal and adult stages, respectively. We added the information in the figure legend (Line 553-554).\nComment 3-21\nLine 479: It should be “Fig. 7b” not “Fig. 6b”.\nReply:\nThank you for pointing out this error. We apologize for the incorrect figure reference and have corrected “Fig. 6b” to “Fig. 7b” in the revised manuscript (Line 565 and 569). We have also reviewed all figure citations to ensure their accuracy throughout the text.\nComment 3-22\nLines 489-494: These statements are stronger than what the evidence showed. The authors should weaken them to align with the evidence. Such statements also belong to Discussion, not Results.\nReply:\nThank you for the suggestion. We have revised the statements to better reflect the strength of the evidence and have moved them from the Results section to the Discussion (Line 632-639), where they are more appropriately placed. In addition, based on your suggestion in Comment 3-9, we have further elaborated on the observed trend in tumors, which helped us refine and expand our hypothesis.\nComment 3-23\nMany readers may be interested in Table S4. It would be helpful to provide a brief description in Results (lines 508-512) and list it explicitly in Data Availability.\nReply:\nThank you for the suggestion. We have added a brief description of Table S4 in the Results section (Line 259-260) to highlight its relevance, and we now explicitly include it in the Data Availability statement to ensure it is easily accessible to readers (Line 1056-1057).\nComment 3-24\nLine 846: There are 14 columns in Table S4, not 15 columns. This is also a good place to list Table S4.\nReply:\nWe are sorry for this mistake. We have corrected the description to reflect that Table S4 contains 14 columns (Line 1024). We have also added an explicit reference to Table S4 here (Line 1043).\nReviewer #4\nComment 4-0\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\nThe authors have provided a comprehensive tutorial on how to use the cPeak database.\nReply:\nThank you for your comment and for participating in the review process.\nResponses to Reviewer #2’s comments\nGeneral response to Reviewer #2\nWe sincerely thank you for your careful and constructive feedback, and especially for taking the additional effort to review our responses to Reviewer #1. We greatly appreciate the time and thought you invested in helping us strengthen the manuscript. Following your suggestions, we have substantially revised the text, clarified methodological details, and added new analyses where appropriate. Our detailed, point-by-point responses are provided below.\nComment 2-1. Clarification of Figure 1d and Peak Shape Similarity Analysis\nThe revised Figure 1d helps clarify the use of MACS2 and the overall strategy for generating cPeaks. I recommend that the authors provide additional detail on how the Pearson correlation was calculated. It would also be helpful to define the two negative control groups and the true peak group more clearly. Furthermore, I suggest avoiding the term “true peaks” when referring to the experimental group, since these peaks are still based on estimates from observed ATAC-seq or DNase-seq data and not experimentally validated ground truth.\nResponse:\nThank you for the helpful suggestions. To improve clarity, we have substantially revised the description of the peak shape similarity analysis in the Methods (lines 700–711). In the revised text, we describe how the genomic window was defined, how peaks were trimmed or zero-padded for alignment, and how base-pair-level values were used to calculate the PCC. We also provide clearer definitions of the two control groups. To further aid understanding, we added a new schematic (Figure R 1, or Fig. S1b in the revised manuscript) that illustrates the calculation procedure and the construction of the two control groups.\nIn addition, we replaced the term “true peak pairs” with “paired peaks” throughout the manuscript to avoid implying experimentally validated ground truth. We believe these revisions address your concerns and make the analysis more transparent and rigorous.\nFigure R 1 Diagram of PCC computation and control group definitions.\nComment 2-2. Concern Regarding Potential Model Circularity\nThe authors appear to have misunderstood the original question regarding circularity. The concern is whether the current training process, which uses DNA sequence and cPeaks derived from observed ATAC-seq data, might predict novel peaks that are not included in the existing cPeak list. If those novel peaks were then added to the cPeak list and used to retrain the model, it is unclear whether this would result in further predictions of additional peaks. This raises the possibility of an iterative feedback process without a clear stopping criterion. The authors are encouraged to clarify whether such retraining was considered and whether this approach would be methodologically appropriate.\nResponse:\nThank you for further clarifying this question. Indeed, iteratively expanding the training set with model-predicted positives is a potential strategy explored in other contexts, such as positive-unlabeled (PU) learning. Following your suggestion, we evaluated whether such an iterative retraining approach could improve our model. Specifically, we added all predicted cPeaks (~0.28 million) to the training set and retrained the CNN. The updated model subsequently predicted an additional ~0.22 million peaks with scores above the threshold. Repeating this procedure led to the identification of progressively more peaks until nearly all candidate regions were classified as cPeaks, indicating uncontrolled expansion of predictions (Figure R 2). The results suggested a clear behavioral drift in the retrained model: with each iteration, the CNN increasingly classified nearly all candidate regions as positive. We hypothesize that this is likely due to the fixed-window prediction scheme used by the CNN, which does not provide precise peak boundary information. These results indicate that iterative retraining may propagate boundary noise and accumulate errors rather than improve model performance. Therefore, we consider the one-time training strategy adopted in this study to be the most appropriate and stable approach.\nTo help readers better understand this point, we have added a brief explanation in the revised manuscript (Lines 656–657, 803–805), noting that improving the prediction procedure to more accurately capture peak boundaries will be an important direction for future work.\nFigure R 2 Number of newly/accumulated predicted cPeaks identified after each iteration of model retraining.\nComment 2-3. Benchmarking Comparison in Figure 1h\nThe comparison in Figure 1h may be affected by differences in the peak-calling tools used. Both cPeaks and CATLAS cCREs were generated using MACS2, but cDHSs were generated using Hotspot2. This makes it hard to tell whether the performance differences are due to biological signals or the choice of tool. Please clarify this point in the manuscript and consider discussing how tool differences may affect the results.\nResponse:\nThank you for the suggestion. We agree that the comparison in Figure 1h could potentially be influenced by differences in peak-calling algorithms, as cDHSs were originally generated with Hotspot2, whereas both cPeaks and CATLAS cCREs were produced using MACS2.\nTo directly address this concern, we attempted to reanalyse the data using Hotspot2 and perform comparison. However, despite our extensive technical efforts (including reimplementation of preprocessing scripts, setting computational environments on different servers, optimization of input data, and multiple recovery strategies), the analysis could not be completed successfully due to a fatal data-dependent bug in the core hotspot2_part1 module. This issue, which causes silent crashes even on validated datasets, reflects the fact that Hotspot2 has not been actively maintained for nearly a decade and is no longer compatible with nowadays computing environments.\nWe therefore approached the question conceptually. MACS2 is now commonly adopted as the default peak-calling tool for ATAC-seq in community pipelines (e.g., Currently, MACS2 is the default peak caller of the ENCODE ATAC-seq pipeline. from Yan et al., Genome Biology, 2020, https://doi.org/10.1186/s13059-020-1929-3). Since the dataset analysed here is ATAC-based, using MACS2 is both appropriate and consistent with community standards.\nWe acknowledge your point regarding the potential methodological bias and have added a clarification in the revised manuscript (Lines 179–183) to explicitly discuss this limitation and the rationale behind our tool selection.\nReviewer #2’s additional comments to cover Reviewer #1\nComment 2-4. Response to comment 2 by Reviewer #1\nThe authors have improved the text, but their response still lacks clarity on the rule of thumb or decision criterion for classifying DNA sequence patterns into cPeak vs. non-cPeak. It would strengthen the manuscript if they could provide more details on the exact procedure.\nResponse:\nThank you for the suggestion. We have now clarified the decision criterion for distinguishing cPeaks from non-cPeaks. Specifically, we added a new subsection titled “Determination of threshold for predicting additional cPeaks” (Lines 781–787), where we describe how the threshold was derived from the masking experiments. In brief, we examined the predicted score distributions of masked cell type- or tissue-specific cPeaks, identified an elbow point (the value at which the frequency of scores began to rise sharply) in each case, and averaged these values across all datasets to obtain a unified threshold of 0.872. This empirically derived threshold was then used to classify additional predicted cPeaks.\nComment 2-5. Response to comment 3 by Reviewer #1\nThe reviewer 1 raised a key concern: if the training dataset used to generate cPeaks does not include one cell type (e.g., a rare cell type), what is the rationale for expecting cPeaks to aid in the discovery of that cell type? The authors’ current response does not sufficiently address this. Intuitively, if the training data lacks signals from a particular cell type, then simply using the aggregated cPeak list would not necessarily enable its detection. The authors should provide a stronger explanation or empirical demonstration of why and how cPeaks could facilitate rare cell type identification under such conditions.\nResponse:\nThank you for the opportunity to clarify our rationale regarding rare cell types. We would like to clarity that the advantage of cPeaks lies in two aspects.\nFirst, cPeaks address the limitation of conventional pseudo-bulk peak calling, which is strongly biased toward prevalent signals. Accessibility sites from rare populations are often treated as noise and filtered out during peak calling, meaning that the information is lost before downstream analysis. In contrast, cPeaks are constructed across a wide range of tissues and cell types, independent of their abundance in any single dataset. This avoids prevalence-driven filtering and ensures that many informative regions relevant to rare populations remain in the reference feature space.\nSecond, to further address the scenario where a cell type or tissue is entirely absent from the training data, we use a CNN to predict additional cPeaks. As shown in Fig. 2c, the model achieved a mean AUC of 91.67% when predicting peaks specific to cell types or tissues that were masked during training. This demonstrates that even for unseen populations, most of their accessible regions can still be recovered as predicted cPeaks. Together, the broad construction of cPeaks and the predictive capacity of the CNN ensure that rare or missing cell types are not excluded from analysis. We have added further clarification on this point in Lines 427–432.\nResponses to Reviewer #4’s comments\nRemarks on code availability\nThe code has a well-written tutorial and introduction.\nResponse:\nThank you. We are glad that you found the tutorial and introduction clear and well written.\nResponses to Reviewer #5’s comments\nGeneral response to Reviewer #5\nWe sincerely thank you for carefully reviewing our manuscript and for also taking the extra step of examining how we addressed the comments from the previous reviewer. Your feedback helped us substantially improve both the clarity and the rigor of the manuscript. Following your suggestions, we revised the text, clarified details, and added new analyses. Our detailed, point-by-point responses are provided below.\nComment 3-3 (presentation order of similarity analysis)\nGeneral comment: You added cPeak-level shape-consistency results (Fig. R12), which is excellent, but the paragraph order was not reversed as originally suggested. While this is a stylistic matter, switching the intermediate- vs final-peak sections would make the logic flow more intuitive for readers.\nDetailed comment: Although the authors kept the original paragraph order (justifying the merging strategy with intermediate-peak results before introducing the new cPeak analysis), the added Fig. R12 and Suppl. S1e-f now demonstrate robust similarity for the final cPeaks, the scientific concern is resolved even if the stylistic request to swap those sections is unmet.\nResponse:\nThank you for the suggestion. We fully agree that readers are primarily interested in the final cPeaks. At the same time, we believe that showing the stability of the intermediate peaks is essential to justify the feasibility of constructing consensus peaks, which is a central part of our argument.\nTo better balance these two aspects, we have made two changes in the revision. First, we shortened the intermediate-peak section so that the key message is easier to grasp (Lines 110–135). Second, we moved the intermediate-peak results to the Supplementary Information (now Fig. S1c–d) and highlighted the final cPeak results in the main text (now Fig. 1d). We appreciate that this is largely a stylistic matter, and we hope that the current restructuring strikes a good balance between logical justification and readability.\nComment 3-4 (TF-binding footprints in Fig. 1c/f)\nGeneral comment: No footprinting analysis was added; the manuscript now only states that read-level data were unavailable. Even a limited demonstration (e.g., aggregate cut-site metaplots for one representative ENCODE BAM) would provide direct evidence of TF occupancy and fully satisfy the request.\nDetailed comment: The authors states that transcription-factor footprinting requires base-resolution read coverage, which they did not retain when summarising peaks; they therefore offer no new analyses or figure updates and only suggest footprinting for “future work.” A review of the revised manuscript confirms that no footprinting results (or TF-binding evidence) were added anywhere, so the original request remains unmet.\nAt minimum, the authors could have downloaded one high-depth ENCODE bulk ATAC-seq BAM that matches your example peaks, extracted single-base Tn5 cut sites, scanned those peaks for frequent JASPAR motifs (e.g., CTCF or AP-1), and produced an aggregate ±50 bp metaplot showing the expected central dip flanked by higher insertion density; a simple footprint-depth statistic and one or two IGV snapshots would have provided direct, base-resolution evidence of TF occupancy without re-processing all 624 samples.\nResponse:\nThank you. Following your advice, we performed an additional demonstration using publicly available read-level ATAC-seq data to directly visualize transcription factor (TF) footprints. Specifically, we downloaded the original bulk ATAC-seq BAM files from a classical study (Buenrostro et al., Nature Methods, 2013, https://doi.org/10.1038/nmeth.2688). To ensure computational efficiency, we focused our analysis on chromosome 4. Using the TOBIAS framework, we conducted TF-footprinting analysis within the open chromatin regions defined by cPeaks.\nWe obtained clear footprinting patterns and visualized the aggregated ATAC-seq cut-site profiles and heatmaps. Representative examples for the canonical motifs CTCF and BATF are shown, both displaying the expected central depletion of insertions flanked by regions of increased cut-site density, consistent with characteristic TF occupancy signatures (Figure R 3).\nThese results demonstrate that regions defined by cPeaks indeed cover TF-binding sites and that the cPeak framework can support footprint-based analyses when base-resolution data are available. Thank you again for your valuable suggestion.\nFigure R 3 Footprinting of CTCF (left) and BATF (right) motifs within the open chromatin regions defined by cPeaks.\nComment 3-6 (PBMC benchmark)\nGeneral comment: The additional PBMC-10k dataset and QC filters are welcome, but two gaps remain:\n(i) Scalability to a truly large PBMC collection (≥ 50 k cells) is still untested.\n(ii) A dedicated doublet-detection tool (e.g., Scrublet, DoubletFinder) has not been applied. Including one run and reporting the percentage of flagged cells would strengthen confidence in the benchmark.\nDetailed comment: The benchmarking section is substantially improved—additional feature sets, clearer clustering metrics, UMAP examples, and a second PBMC dataset now appear. However, the new PBMC-10k still falls short of the very large cell counts suggested, and doublet handling remains minimal. Thus Comment 3-6 is largely but not completely resolved.\nTo finish closing the gap you should still (i) benchmark on a genuinely large PBMC dataset—one of the public 50 k–100 k cell collections—so that scalability to tens-of-thousands of cells is demonstrated, (ii) apply a dedicated doublet-detection tool such as Scrublet, DoubletFinder or AMULET and report how many cells are removed and whether results change, and (iii optional) include at least one TEA-seq or other multimodal dataset (or an unsupervised variant of the united-peaks reference) to round out the cross-platform comparison.\nResponse:\n(i) Benchmarking on a larger-scale dataset\nThank you for highlighting the importance of testing scalability. In response, we downloaded the paired scRNA-seq and scATAC-seq dataset GSE190992. After quality control and doublet removal using the ArchR doublet-detection module, we obtained 48,487 high-confidence scATAC-seq cells for downstream analysis, which we refer to as PBMC-48k. Benchmarking results (Figure R 4) showed that cPeaks perform comparably or better in detecting rare cell types at this scale, supporting their robustness in large datasets. These results have been added to Fig. S5e–h and described in Lines 459–463 and 977–983 of the revised manuscript.\nFigure R 4 Rare cell type discovery using cPeaks on the PBMC-48k dataset. (e) The proportion of each cell type in the PBMC-48k dataset. Three rare cell types were less than 1% (the dashed black line): cDC2, Eryth, and B intermediate. (f) The UMAP plots of the PBMC-48k dataset analyzed by each reference set. Each dot represented one cell. Different colors represented different clusters. (g) Performance of clustering accuracy using different reference peak sets in the PBMC-48k dataset. (h) The Venn plot of the number of differential accessible regions calculated by cPeaks and other reference sets for pDCs in the PBMC-48k dataset.\n(ii) Doublet-detection tool\nSpecifically, we processed the scATAC-seq data using the ArchR pipeline, identified 9.08% of cells as doublets with the addDoubletScores function, and removed them using filterDoublets. To further evaluate the impact of doublet removal, we also performed the same analysis without filtering doublets and found that cPeaks no longer achieved the best clustering performance (Figure R 5). This indicates that cPeaks perform particularly well on high-quality datasets with proper doublet removal. Details of doublet removal have been added in Lines 978–979.\nIn addition, we confirmed that all datasets analyzed in the original manuscript had already undergone doublet filtering by their original authors, so the results presented previously were not affected by doublet contamination.\nFigure R 5 Performance of clustering accuracy using different reference peak sets on the PBMC-48k dataset without doublet removal.\n(iii) Benchmarking on TEA-seq\nRegarding TEA-seq, we explored available datasets (e.g., https://doi.org/10.7554/eLife.63632) but were unable to obtain reliable cell type labels from the original data to support a meaningful benchmark. We also contacted the corresponding authors for clarification but did not receive a response. In addition, our attempts to reproduce their published analysis pipeline did not yield results consistent with those reported in the article. For these reasons, we were unable to include a TEA-seq benchmark in the current revision.\nComment 3-8 (pseudotime analysis with alternative feature sets)\nGeneral Comment: The new lineage-wise plots using cPeaks clarify the trend, yet the analysis is still shown for cPeaks alone. Please add the same pseudotime/UMAP plots for at least one other reference (e.g., CATLAS cCREs) or provide quantitative evidence that those references cannot support an equivalent well-positioned-peak metric.\nDetailed Comment: The presentation is clearer and the new pseudotime plots meet requests (b) and (c). However, point (a) remains only partially resolved: the manuscript still lacks a side-by-side pseudotime/UMAP demonstration with at least one other reference set, so the claimed uniqueness of the well-positioned-cPeak metric is asserted rather than shown.\n(a) is still incomplete: to substantiate the claimed novelty, please include the same pseudotime/UMAP analysis with at least one alternative reference (e.g., CATLAS cCREs or cDHSs) or provide quantitative evidence that those feature sets cannot support an equivalent well-positioned-peak metric.\nResponse:\nThank you for this suggestion. We would like to emphasize that the well-positioned-cPeak metric is a property unique to cPeaks and cannot be directly derived from other reference sets. Our goal in this section was not to argue that only cPeaks can recover developmental trajectories, but to demonstrate the distinctive biological insights enabled by this metric.\nIn line with your comment, we repeated the pseudotime/UMAP analysis using CATLAS cCREs and cDHSs (Figure R 6). The overall trajectories were similar across reference sets; however, only cPeaks and cCREs revealed the two sub-branches within the RGC lineage that were noted in the original study (Figure R 7), whereas cDHSs failed to capture this structure. Since the original work did not provide detailed analysis of these sub-branches, we cannot ascertain their biological significance at this stage.\nTo avoid overstating the claims, we have revised the text to make clear that this analysis is intended to highlight the unique role of the well-positioned-cPeak metric in developmental studies, rather than to suggest exclusive superiority of cPeaks in trajectory reconstruction (Lines 523-524, 530-531).\nFigure R 6 UMAP visualization of the 8-week retina scATAC-seq data analyzed by different reference sets, colored by cell types (top) and pseudotime (bottom).\nFigure R 7 UMAP visualization of the 8-week retina scATAC-seq data from the original paper (https://doi.org/10.1016/j.celrep.2021.110294). Cells are colored by cell types (left), pseudotime (right).\nComment 3-17 (minor) (annotation method)\nGeneral Comment: The revision retains the Seurat/Signac anchor strategy without testing any newer scATAC annotation tools. Running a contemporary method such as ArchR label transfer, scGLUE or CellTypist—and showing that conclusions remain unchanged—would address the original concern.\nDetailed Comment: I appreciate the justification for retaining the Seurat/Signac anchor approach, but because no alternative annotation method was tested the limitation highlighted in Comment 3-17 persists; incorporating at least one contemporary scATAC-labelling tool (e.g., ArchR label transfer or scGLUE) would strengthen confidence in the annotation robustness.\nResponse:\nThank you. As part of the new PBMC-48k analysis described in our response to Comment 3-6, we applied ArchR’s label transfer approach to annotate scATAC-seq cells. Specifically, for the scRNA-seq data, cells with predicted.celltype.l2.score > 0.9 were retained as high-confidence annotations. These cells were then used to transfer labels to scATAC-seq data via the addGeneIntegrationMatrix function in ArchR, which also provides a predictedScore for each assigned label. We selected 45,761 scATAC-seq cells with a predictedScore > 0.9 for downstream analyses.\nAs shown in Figure R 4, the results were consistent with PBMC-8k which transferred labels with Seurat, thereby confirming that our conclusions are robust to the choice of annotation method. We have added the details of label transferring in Lines 977-983.\nComment 3-19c (minor) (comparability of UMAP projections)\nGeneral Comment: UMAPs are still calculated separately for each dataset, limiting direct visual comparison across life stages. A joint embedding (or at least side-by-side UMAPs generated with the same global coordinates) would make cross-dataset patterns easier to interpret.\nDetailed Comment: Points (a) and (b) are fully satisfied. Point (c) is only partially resolved because the figures still rely on distinct projections, limiting direct visual comparison across datasets.\nThe added GO-term heatmap reference and the new life-stage-coloured UMAPs clarify the biology nicely; however, because each dataset is still visualised with its own UMAP embedding, cross-dataset comparisons remain difficult. A joint projection (or at least side-by-side UMAPs created with identical global coordinates) would fully resolve the comparability issue noted in Comment 3-19.\nResponse:\nThank you for the suggestion. We have now performed a joint embedding by visualizing the different developmental stages of the retina together (Figure R 8), and the results have been added to Fig. S7a, S7d, and Lines 519, 537. Consistent with our previous observations, the well-positioned-cPeak ratio shows an overall decrease with developmental progression, supporting the biological trend described in the main text.\nFigure R 8 (a) The UMAP plot of all scATAC-seq data collected from retina. Cells are colored by sample labels (left) and cell types (right). (d) The UMAP plot of all scATAC-seq data collected from retina. Cells are colored by well-positioned-cPeak ratios.\nRemarks on code availability\nI did not run it, but looked at whether the repo was usable.\nThe repository is usable and largely reproducible for downstream analyses (mapping fragments, rerunning benchmarks). Re-generating the entire resource (merging 624 datasets or replicating the CNN training) would still require substantial external data and compute, but that is typical for large genomics references.\nSuggested very minor improvements:\n(i) Add a requirements.txt or environment.yml with tested versions.\n(ii) Upload the pre-trained CNN weights (or provide an S3/Zenodo link) and a small test fragment file for a quick end-to-end smoke test.\n(iii) Include a minimal CI job that runs map2cpeak on the toy dataset to ensure code health over time.\n(iv) With these tweaks, the repo would be an even more robust community resource.\nResponse:\nThank you for these practical suggestions. We have made the following updates:\n(i) added a requirements.txt file specifying tested versions at https://github.com/MengQiuchen/cPeaks/blob/main/deepCNN/requirements.txt and https://github.com/MengQiuchen/cPeaks/blob/main/requirements.txt;\n(ii) uploaded the pre-trained CNN weights at:\nhttps://github.com/MengQiuchen/cPeaks/tree/main/deepCNN/save_model;\n(iii) provided a small test fragment file for running a quick map2cpeak end-to-end test at:\nhttps://github.com/MengQiuchen/cPeaks/tree/main/map2cpeak/demo, and implemented a minimal CI job at https://github.com/MengQiuchen/cPeaks/blob/main/.github/workflows/smoke-test.yml.\nWe hope that these improvements can make the repository more robust and user-friendly for the community.\nResponses to Reviewer #2’s comments\nThe authors have addressed most of our comments. One minor issue remains: Reviewer 1’s question concerns unobserved cell types—i.e., cases with zero abundance in the training data. In that setting, how can cPeaks serve as a valid peak list for that cell type? The hold-out (masked) evaluation in Fig. 2c helps address this and should be emphasized. However, the newly added sentence “cPeaks are constructed across a wide range of tissues and cell types, independent of their prevalence in any single dataset…” is misleading. “Independent” is ambiguous and overstates the claim. We understand the intent is to say that aggregating across datasets reduces prevalence-driven filtering (so rare-accessible loci are less likely to be lost), not that cPeaks are statistically independent of training composition.\nReponses:\nThank you for the helpful feedback. We agree with your point regarding the ambiguity of the term “independent.” We have revised the sentence in Line 429 to read: “cPeaks are constructed across a wide range of tissues and cell types, and their composition is not dominated by the prevalence of any single dataset.” This change more accurately reflects our intended meaning.\nResponses to Reviewer #5’s comments\nAll six points are addressed to a standard I would consider acceptable for a second-round review. I’d like to thank the authors for addressing these comments.\nThe last tweak I would still ask for before acceptance, is to move some figures in the rebuttal into the Supplement. It would greatly benefit the manuscript in my opinion.\nFor instance you could migrate rebuttal Figs. R3, R5, and R6–R7 into the Supplement (either as new Supplementary Figs, or as additional panels) so that the footprinting evidence, the no-doublet control, and the alternative-feature pseudotime comparisons are part of the archived record.\nReponses:\nThank you very much for your suggestions. We have incorporated your recommendations as follows:\nl We have added Fig. R3 to the manuscript as Fig. S2c, cited it in Lines 299–230, and included an explanation of the footprinting experiment in Lines 830–834.\nl We have added Fig. R5 as Fig. S5i, cited it, and discussed the results in Lines 989–991.\nl We have added Fig. R6 as Fig. S7d, and cited it in Line 531.\nl We did not include Fig. R7, as it corresponds to figures from the original manuscript (Ref. [34] in our paper). However, we added a note in Line 528 indicating that the trajectories are “similar to the trend reported in the original study,” to help readers better interpret the results."}, {"round": 2, "reviewer_comments": {"reviewer_2": "The authors have addressed most of our comments. One minor issue remains: Reviewer 1’s question concerns unobserved cell types—i.e., cases with zero abundance in the training data. In that setting, how can cPeaks serve as a valid peak list for that cell type? The hold-out (masked) evaluation in Fig. 2c helps address this and should be emphasized. However, the newly added sentence “cPeaks are constructed across a wide range of tissues and cell types, independent of their prevalence in any single dataset…” is misleading. “Independent” is ambiguous and overstates the claim. We understand the intent is to say that aggregating across datasets reduces prevalence-driven filtering (so rare-accessible loci are less likely to be lost), not that cPeaks are statistically independent of training composition.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_5": "Remarks to the author.\nAll six points are addressed to a standard I would consider acceptable for a second-round review. I’d like to thank the authors for addressing these comments.\nThe last tweak I would still ask for before acceptance, is to move some figures in the rebuttal into the Supplement. It would greatly benefit the manuscript in my opinion.\nFor instance you could migrate rebuttal Figs. R3, R5, and R6–R7 into the Supplement (either as new Supplementary Figs, or as additional panels) so that the footprinting evidence, the no-doublet control, and the alternative-feature pseudotime comparisons are part of the archived record.\nComment 3-3 (presentation order of similarity analysis)\nIt was asked to swap paragraph order so final cPeaks precede intermediate-peak results.\n=>Acceptable.\nThe stylistic swap has suggested wasn’t done, but the scientific concern (robust similarity for final cPeaks) is now clear and appropriately placed (main vs. supplement).\nComment 3-4 (TF-binding footprints in Fig. 1c/f)\nIt was asked to add footprinting evidence (e.g., aggregate cut-site metaplots).\n=>Acceptable, *with a note*.\nThe analysis meets the request. However, the footprint panels appear only in the rebuttal (I don’t see the figure in the supplied Supplementary PDF). I’d recommend the authors port Figure R3 into the Supplement so the evidence is captured in the manuscript files.\nComment 3-6 (PBMC benchmark)\nIt was asked to show scalability (≥ 50k cells PBMC)\n=> Acceptable.\nThey met both scale and doublet-handling asks. TEA-seq wasn’t added (they explain label issues), which is reasonable for this revision.\nComment 3-8 (pseudotime analysis with alternative feature sets)\nIt was asked to show side-by-side pseudotime/UMAP using another reference (e.g., cCREs), or justify why others can’t support the metric.\n=> Acceptable\nThe side-by-side demonstration is provided (in rebuttal figures) and the main text is aligned with a balanced claim.\nComment 3-17 (minor) (annotation method)\nIt was asked to test a modern alternative to Seurat/Signac (e.g., ArchR label transfer).\n=>Acceptable.\nThis directly addresses robustness to the annotation framework.\nComment 3-19c (minor) (comparability of UMAP projections)\nIt was asked to provide a joint embedding (or at least shared coordinates) across developmental stages.\n=>Acceptable.\nThe comparability issue is resolved with a unified projection."}, "author_response": "Reponses:\nThank you very much for your suggestions. We have incorporated your recommendations as follows:\nl We have added Fig. R3 to the manuscript as Fig. S2c, cited it in Lines 299–230, and included an explanation of the footprinting experiment in Lines 830–834.\nl We have added Fig. R5 as Fig. S5i, cited it, and discussed the results in Lines 989–991.\nl We have added Fig. R6 as Fig. S7d, and cited it in Line 531.\nl We did not include Fig. R7, as it corresponds to figures from the original manuscript (Ref. [34] in our paper). However, we added a note in Line 528 indicating that the trajectories are “similar to the trend reported in the original study,” to help readers better interpret the results."}, {"round": 3, "reviewer_comments": {"reviewer_2": "The authors have addressed all my comments.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": ""}]

Decision letter