Peak Calling in Single-Cell ATAC-Seq: A Comparison of Methods and Best Practices

By Dr. Zubair Khalid, DVM, MS, PhD · 2026-08-20

Key Takeaways

• Data sparsity is a fundamental challenge in scATAC-seq, with only 1-10% of peaks detected per cell, necessitating aggregation strategies like pseudobulking (e.g., with MACS2) or cluster-specific calling to identify robust genomic features.
• Traditional peak calling methods (e.g., MACS2) are not inherently designed for scATAC-seq's sparsity and fragment structure; newer approaches like SnapATAC and ArchR offer integrated workflows and address these limitations through binning or pseudo-bulk replicates.
• Alternative frameworks like CREscendo and SCARlink are emerging, moving beyond de novo peak calling by leveraging Tn5 cleavage frequencies with regulatory annotations or employing regression on tile-level data, respectively, to improve portability and gene-accessibility linkage.
• Quality control metrics such as fragment size distribution (revealing nucleosome positioning), reads in peaks (FRiP), and peak width distribution are critical for assessing data quality and informing parameter selection for peak callers.
• Reproducibility in scATAC-seq peak calling requires meticulous documentation of software versions, parameters, and quality metrics, alongside consideration of batch effect correction (e.g., with epiConv) when integrating multiple datasets.

Single-cell ATAC-seq (scATAC-seq) maps chromatin accessibility across thousands of individual cells, but the data are inherently sparse because each cell carries only two copies of the genome. Peak calling, the process of identifying regions of open chromatin, is a foundational step that shapes all downstream analyses including clustering, cell-type annotation, differential accessibility, and transcription factor activity inference. This article compares the most widely used peak calling approaches for scATAC-seq data, explains their underlying assumptions, and provides practical guidance for selecting methods and parameters based on your experimental goals and data characteristics.

The choice of peak calling strategy is not a trivial technical detail. It determines which genomic regions are used as features for quantifying accessibility, and different methods can produce substantially different results from the same input data. A benchmarking study of 10 computational methods for scATAC-seq analysis found that the choice of feature definition and processing approach significantly affected the ability to discriminate cell types, with methods varying in their sensitivity to sequencing coverage and noise levels [<a href="#ref-1">1</a>]. More recent work has questioned whether traditional peak-based quantification is even the most appropriate framework for single-cell data, given the sparsity and sampling limitations inherent to the assay [<a href="#ref-2">2</a>].

This article is written for biology students, researchers, laboratory professionals, and life-science practitioners who need to make informed decisions about their scATAC-seq analysis pipelines. We cover the major peak calling tools including MACS2, SnapATAC, and ArchR, discuss their strengths and limitations in the single-cell context, and provide concrete recommendations for parameter selection and quality assessment.

Understanding scATAC-seq Data Characteristics

Data Sparsity and Its Implications

The fundamental challenge in scATAC-seq analysis is data sparsity. Unlike scRNA-seq where each cell can express thousands of genes, scATAC-seq samples a tiny fraction of the accessible chromatin in each cell. A benchmarking study reported that only 1 to 10 percent of peaks are detected per cell in scATAC-seq data, compared to 10 to 45 percent of expressed genes detected per cell in scRNA-seq data [<a href="#ref-1">1</a>]. This sparsity arises because each cell contains only two copies of the genome, and the transposase enzyme can only tag a limited number of accessible regions before the cell is lysed and sequenced.

This sparsity has direct consequences for peak calling. When you aggregate reads across all cells to call peaks, you are relying on the assumption that accessible regions are consistently accessible across many cells in your population. Regions that are accessible only in rare cell types may fall below the detection threshold and be missed entirely. Conversely, regions that are accessible in many cells will produce strong peak signals, but the per-cell quantification of those peaks will still be sparse and noisy.

A 2025 review of computational challenges in scATAC-seq emphasized that low genomic coverage per cell results in intrinsic data sparsity and missing-data issues that present unique methodological challenges [<a href="#ref-3">3</a>]. These challenges affect every stage of analysis, but they are most acute at the feature definition stage because errors made here propagate through all downstream steps.

Fragment Structure and Tn5 Bias

ATAC-seq relies on the Tn5 transposase to cleave and tag accessible chromatin. The enzyme has sequence preferences that introduce bias into the data, and understanding this bias is important for interpreting peak calling results. The Tn5 cleavage frequency itself carries information about chromatin structure and transcription factor binding, and some newer methods leverage this information directly instead of relying on peak boundaries.

A 2025 study introduced CREscendo, a method that uses Tn5 cleavage frequencies and regulatory annotations to identify differential usage of candidate regulatory elements across cell types [<a href="#ref-2">2</a>]. The authors argued that conventional peak-based methods can fail to capture precise cell-type-specific regulatory signals, producing results that are difficult to interpret and lack portability across datasets. This represents a shift away from traditional peak-based quantification toward a more robust framework that relies on standardized references of annotated regulatory elements.

For most researchers, however, peak calling remains the standard approach, and understanding the fragment structure of your data is essential for choosing appropriate parameters. Fragment size distributions, insertion site patterns, and mitochondrial contamination levels all provide information about data quality that should inform your peak calling strategy.

Library Complexity and Sequencing Depth

The number of unique fragments sequenced per cell, often called library complexity, varies substantially across scATAC-seq protocols. A 2024 benchmarking study of eight scATAC-seq methods across 47 experiments using human peripheral blood mononuclear cells found significant differences in sequencing library complexity and tagmentation specificity across protocols [<a href="#ref-4">4</a>]. These differences impacted cell-type annotation, genotype demultiplexing, peak calling, differential region accessibility, and transcription factor motif enrichment.

This finding has practical implications for peak calling. If your library complexity is low, you may need to aggregate more cells to call peaks reliably, or you may need to use a reference-based approach instead of calling peaks de novo from your own data. If your tagmentation specificity is poor, you may see elevated background signal that makes peak calling more difficult and requires more stringent thresholds.

Peak Calling Methods for scATAC-seq

MACS2: The Bulk ATAC-seq Standard Applied to Single-Cell Data

MACS2 is the most widely used peak caller for bulk ATAC-seq and ChIP-seq data, and it is frequently applied to scATAC-seq data after aggregating reads across cells. The method models the distribution of tag counts along the genome, identifies regions where the signal exceeds the local background, and estimates false discovery rates using a Poisson or negative binomial model.

In the single-cell context, MACS2 is typically applied to a pseudobulk profile created by merging all fragments from all cells or from cells within a cluster. This approach leverages the depth of the aggregated data to identify peaks that would be undetectable in individual cells. A 2026 pipeline for single-cell chromatin accessibility analysis describes a workflow that uses MACS2 for peak calling after preprocessing with scATAC-pro or Cell Ranger ATAC, followed by differential accessibility analysis to detect open chromatin regions and highlight regulatory differences among cell populations [<a href="#ref-5">5</a>].

The main advantage of MACS2 is its familiarity and extensive documentation. Many researchers have experience with MACS2 from bulk ATAC-seq or ChIP-seq projects, and the parameters are well understood. The main disadvantage is that MACS2 was not designed for single-cell data and does not account for the specific characteristics of scATAC-seq, such as the fragment size distribution that reflects nucleosome positioning or the sparsity of per-cell coverage.

When using MACS2 for scATAC-seq data, several parameters require attention. The shift and extension parameters should be adjusted based on your fragment size distribution. For ATAC-seq, the standard recommendation is to use a shift of 100 base pairs and an extension of 200 base pairs for a 50 base pair shift and 100 base pair extension on each side, but these values should be validated against your specific data. The q-value threshold determines the stringency of peak calling, and lower thresholds will produce more peaks but also more false positives.

A 2026 study constructed a generic chromatin accessibility reference by aggregating peaks from 624 high-quality bulk ATAC-seq datasets, defining about 1.4 million consensus peaks [<a href="#ref-6">6</a>]. The authors found that these consensus peaks exhibited consistent shapes across tissue types, sequencing technologies, and peak-calling methods, indicating that they represent inherent genomic features. This suggests that peak calling methods, including MACS2, are identifying real biological features instead of artifacts, but the specific boundaries and exact coordinates may vary across methods and datasets.

SnapATAC: Integrated Processing with Built-in Peak Calling

SnapATAC is a comprehensive R package for scATAC-seq analysis that includes its own peak calling approach. The method was specifically designed for single-cell data and addresses the sparsity problem by using a bin-based approach instead of calling peaks on individual cells. SnapATAC divides the genome into fixed-size bins, typically 5 kilobases, and creates a cell-by-bin matrix of accessibility counts. This matrix is then used for dimensionality reduction and clustering.

The benchmarking study mentioned earlier found that SnapATAC outperformed other methods in separating cell populations across different coverages and noise levels in both synthetic and real datasets [<a href="#ref-1">1</a>]. Notably, SnapATAC was the only method able to analyze a large dataset with more than 80,000 cells in that study. This scalability is an important consideration for modern scATAC-seq experiments that routinely generate tens of thousands of cells.

After clustering, SnapATAC can call peaks on the aggregated reads from each cluster using MACS2 or other peak callers. This cluster-based peak calling approach has the advantage of identifying peaks that are specific to particular cell types, which may be missed when calling peaks on the entire dataset. The cluster-specific peaks can then be used to create a union peak set for downstream analysis.

The bin-based approach used by SnapATAC for initial analysis has both advantages and disadvantages. Bins provide a uniform feature set that does not depend on peak calling, which can be advantageous for clustering and visualization. However, bins do not correspond to actual regulatory elements, and the resolution is limited by the bin size. A 5 kilobase bin may contain multiple regulatory elements with different accessibility patterns, and this can obscure cell-type-specific signals.

ArchR: Flexible Peak Calling with Arrow File Architecture

ArchR is another comprehensive R package for scATAC-seq analysis that has gained widespread adoption. The package uses an Arrow file format to store fragment data efficiently, allowing analysis of very large datasets. ArchR provides multiple options for peak calling, including calling peaks on the entire dataset, calling peaks per cluster, and using a pseudo-bulk approach.

A protocol for analyzing single-nucleus chromatin accessibility data during zebrafish early embryogenesis describes steps for fragment file retrieval, sample integration, quality control, Latent Semantic Indexing clustering, and peak calling via ArchR [<a href="#ref-7">7</a>]. The protocol also details procedures for species-specific motif database construction, motif enrichment, and transcription factor footprinting analysis, highlighting the flexibility of the ArchR framework for non-model organisms.

ArchR uses MACS2 as the underlying peak caller but provides additional functionality for managing and comparing peaks across conditions. The package can create a union peak set from multiple peak calling runs, merge overlapping peaks, and quantify accessibility at each peak for each cell. ArchR also provides tools for visualizing peaks, annotating peaks to genomic features, and linking peaks to genes.

One of the key advantages of ArchR is its handling of the sparsity problem through the use of pseudo-bulk replicates. ArchR can create pseudo-bulk profiles by randomly sampling cells from each cluster, calling peaks on these profiles, and then combining the results. This approach provides a more robust peak set than calling peaks on a single aggregated profile because it accounts for cell-to-cell variability within clusters.

SCARlink: Regression-Based Approach Without Peak Calling

SCARlink represents a fundamentally different approach to analyzing scATAC-seq data that avoids peak calling entirely. The method uses regularized Poisson regression on tile-level accessibility data to jointly model all regulatory effects at a gene locus, predicting single-cell gene expression from chromatin accessibility [<a href="#ref-8">8</a>]. This approach avoids the limitations of pairwise gene-peak correlations and the dependence on peak calling.

The 2024 study that introduced SCARlink found that it outperformed existing gene scoring methods for imputing gene expression from chromatin accessibility across high-coverage multi-ome datasets while giving comparable to improved performance on low-coverage datasets [<a href="#ref-8">8</a>]. Shapley value analysis on trained models identified cell-type-specific gene enhancers that were validated by promoter capture Hi-C and were enriched in fine-mapped eQTLs and genome-wide association study variants.

SCARlink is particularly relevant for multi-ome datasets where scRNA-seq and scATAC-seq are measured in the same cells. The method can link enhancers to target genes using the paired expression and accessibility data, providing insights into gene regulatory mechanisms that are not possible with peak-based approaches alone. However, SCARlink requires paired multi-ome data and is not applicable to scATAC-seq datasets without matched expression data.

CREscendo and the Shift Toward Regulatory Annotations

CREscendo, introduced in 2025, advocates for moving away from traditional peak-based quantification in scATAC-seq toward a more robust framework that relies on a standardized reference of annotated candidate regulatory elements [<a href="#ref-2">2</a>]. The method uses Tn5 cleavage frequencies and regulatory annotations to identify differential usage of candidate regulatory elements across cell types.

The motivation for this approach is the observation that conventional peak-based methods can fail to capture precise cell-type-specific regulatory signals, producing results that are difficult to interpret and lack portability across datasets [<a href="#ref-2">2</a>]. Peak boundaries vary across methods and datasets, making it difficult to compare results between studies. A standardized reference of regulatory elements would address this problem by providing a common feature set for all analyses.

This approach is supported by the finding that consensus peaks derived from bulk ATAC-seq data exhibit consistent shapes across tissue types, sequencing technologies, and peak-calling methods [<a href="#ref-6">6</a>]. If regulatory elements are inherent genomic features with consistent boundaries, then a reference-based approach may be more reproducible than de novo peak calling on each dataset.

epiConv for Batch Effect Correction and Joint Analysis

When analyzing multiple scATAC-seq datasets together, batch effects can substantially impact peak calling and downstream analyses. epiConv is an algorithm designed to perform joint analyses on scATAC-seq datasets by removing technical variations between datasets while preserving biological information [<a href="#ref-9">9</a>].

The 2022 study that introduced epiConv showed that it better corrects batch effects and is less prone to overfitting than existing methods on a collection of PBMC datasets [<a href="#ref-9">9</a>]. The method was also capable of aligning low-depth scATAC-seq from co-assay data onto high-quality ATAC-seq references, increasing the resolution of chromatin profiles. Joint analysis across multiple datasets improved the performance of clustering and differentially accessible peak calling, especially when the biological signal was weak in a single dataset.

For peak calling, batch effect correction is important because technical variations between datasets can create spurious peaks or obscure real ones. If you are combining data from multiple experiments, protocols, or sequencing runs, you should consider whether batch effect correction is needed before peak calling.

At a Glance: Peak Calling Method Comparison

Practical Workflow for Peak Calling

Step 1: Quality Control and Preprocessing

Before peak calling, you must ensure that your data meet quality standards. The preprocessing steps include fragment file retrieval, sample integration, and quality control as described in protocols for scATAC-seq analysis [<a href="#ref-7">7</a>]. Key quality metrics include the number of unique fragments per cell, the fraction of fragments in peaks, the transcription start site enrichment score, and the ratio of mitochondrial to nuclear reads.

Cells with very low fragment counts should be removed because they provide little information and can introduce noise into peak calling. Cells with very high fragment counts may represent doublets or multiplets and should also be examined. The specific thresholds depend on your protocol and data quality, but a common approach is to examine the distribution of fragments per cell and remove outliers.

The fragment size distribution provides information about nucleosome positioning. Fragments shorter than 100 base pairs represent nucleosome-free regions, while fragments of 200 to 400 base pairs represent mononucleosomes. A healthy ATAC-seq library should show a strong nucleosome-free fragment peak and a periodic pattern corresponding to nucleosome positioning.

Step 2: Choose Your Peak Calling Strategy

The choice of peak calling strategy depends on your experimental design and analysis goals. The main options are:

Whole-dataset peak calling: Aggregate all fragments from all cells and call peaks on the resulting pseudobulk profile. This approach maximizes sequencing depth and can detect peaks that are present in any cell type. However, it may miss peaks that are specific to rare cell types because the signal from those cells is diluted by the majority population.

Cluster-based peak calling: Cluster cells first using a bin-based approach or other feature definition, then call peaks on the aggregated reads from each cluster. This approach can detect cell-type-specific peaks that would be missed in whole-dataset peak calling. The tradeoff is that clustering quality depends on the initial feature definition, and errors in clustering will propagate to peak calling.

Reference-based peak calling: Use a predefined set of peaks from a reference dataset or consensus peak set. This approach ensures consistency across datasets and avoids the computational cost of de novo peak calling. The tradeoff is that reference peaks may not capture all relevant regions in your specific dataset, particularly if you are studying a rare cell type or non-model organism.

A 2026 study found that consensus peaks derived from bulk ATAC-seq data showed superior performance in scATAC-seq analyses, improving cell annotation and rare cell type identification compared to existing feature-defining methods [<a href="#ref-6">6</a>]. This suggests that reference-based approaches may be preferable for many applications, particularly when comparing across datasets.

Step 3: Parameter Selection

The specific parameters for peak calling depend on the method you choose. For MACS2, the key parameters are the q-value threshold, the shift and extension values, and the genome size. The q-value threshold controls the false discovery rate, with lower values producing more stringent peak calls. A common default is 0.05, but you may need to adjust this based on your data quality and the number of peaks you expect.

The shift and extension parameters should be set based on your fragment size distribution. For ATAC-seq, the Tn5 transposase binds as a dimer and inserts adapters with a 9 base pair gap, so the actual cut sites are offset from the read start positions. The standard correction is to shift reads by 4 base pairs on the positive strand and 5 base pairs on the negative strand, then extend reads to a fixed length. The extension length should be based on the expected fragment size, typically 100 to 200 base pairs for nucleosome-free fragments.

For ArchR, the peak calling parameters are similar to MACS2 because ArchR uses MACS2 as the underlying peak caller. However, ArchR provides additional options for creating pseudo-bulk replicates and combining peak calls across replicates. The number of pseudo-bulk replicates and the number of cells per replicate should be chosen based on your total cell count and the expected heterogeneity of your sample.

Step 4: Quality Assessment of Called Peaks

After peak calling, you should assess the quality of your peak set. Key metrics include:

Number of peaks: The expected number of peaks depends on the cell types in your sample and the sequencing depth. Human cell lines typically have 50,000 to 200,000 accessible regions, while primary tissues may have more due to cell-type heterogeneity.

Fraction of reads in peaks: This metric measures how much of your sequencing data falls within called peaks. A low fraction suggests that your peak calling was too stringent or that your data have high background. A high fraction suggests that your peak calling was too permissive or that your data are dominated by a few highly accessible regions.

Peak width distribution: The distribution of peak widths should be centered around 200 to 500 base pairs, corresponding to nucleosome-free regions. Very wide peaks may represent regions of broad accessibility or technical artifacts.

Genomic annotation distribution: Peaks should be enriched at promoters, enhancers, and other regulatory regions. If a large fraction of peaks fall in intergenic regions with no known regulatory function, you should examine whether your peak calling parameters are appropriate.

Reproducibility: If you have multiple replicates or pseudo-bulk profiles, you should assess the overlap between peak sets. High overlap indicates reproducible peak calling, while low overlap suggests that your peak calling is sensitive to sampling variation.

Step 5: Integration with Downstream Analysis

The peak set you choose will be used for all downstream analyses, including dimensionality reduction, clustering, cell-type annotation, differential accessibility, and transcription factor activity inference. The choice of peak set can substantially impact these analyses, so you should validate that your peak set produces biologically meaningful results.

For dimensionality reduction and clustering, the peak-by-cell matrix is typically normalized and transformed before applying methods like Latent Semantic Indexing or principal component analysis. The sparsity of the peak-by-cell matrix means that most entries are zero, and this sparsity must be accounted for in the normalization and transformation steps.

For differential accessibility analysis, the peak set defines the regions that will be tested for differences between conditions or cell types. If your peak set is missing relevant regions, you will miss differential accessibility at those regions. If your peak set includes spurious peaks, you will increase the multiple testing burden and reduce statistical power.

For transcription factor activity inference, the peak set defines the regions where motif enrichment will be assessed. Methods like chromVAR and SCENIC+ use the peak set to calculate motif accessibility scores and infer transcription factor activity [<a href="#ref-5">5</a>]. The choice of peak set can substantially impact these analyses, particularly for transcription factors that bind to regions that are not well represented in the peak set.

Common Failure Patterns and Troubleshooting

Failure Pattern 1: Too Few Peaks Called

If your peak calling produces very few peaks, the possible causes include:

Low sequencing depth: If your library complexity is low, you may not have enough reads to detect peaks above the background threshold. The 2024 benchmarking study found significant differences in library complexity across scATAC-seq protocols, and these differences impacted peak calling [<a href="#ref-4">4</a>]. Solutions include increasing sequencing depth, aggregating more cells, or using a reference-based peak set.

Stringent parameters: If your q-value threshold is too low or your shift and extension parameters are incorrect, you may miss real peaks. Try relaxing the q-value threshold or adjusting the shift and extension values based on your fragment size distribution.

Poor data quality: If your data have high background or high mitochondrial contamination, peak calling may be compromised. Check your quality metrics and consider whether additional filtering is needed.

Failure Pattern 2: Too Many Peaks Called

If your peak calling produces an excessive number of peaks, the possible causes include:

Permissive parameters: If your q-value threshold is too high, you will call many false positive peaks. Try a more stringent threshold and assess whether the additional peaks are biologically meaningful.

High background: If your data have high background signal, the peak caller may identify broad regions of elevated signal as peaks. This can happen with poor tagmentation specificity or high levels of ambient DNA.

Repetitive regions: Peaks in repetitive regions may represent mapping artifacts instead of true accessible chromatin. Consider filtering peaks in known repetitive regions or using a more stringent mapping quality threshold.

Failure Pattern 3: Inconsistent Peaks Across Replicates

If peak calling produces inconsistent results across replicates, the possible causes include:

Sampling variation: scATAC-seq data are sparse, and peak calling on aggregated reads is sensitive to sampling variation. Using pseudo-bulk replicates, as implemented in ArchR, can help stabilize peak calls [<a href="#ref-7">7</a>].

Batch effects: Technical variations between datasets can create spurious differences in peak calls. Methods like epiConv can correct batch effects and improve the consistency of peak calling across datasets [<a href="#ref-9">9</a>].

Biological variation: If your replicates represent different biological conditions, the peak sets may legitimately differ. In this case, you should consider whether a consensus peak set or a union peak set is more appropriate for your analysis.

Failure Pattern 4: Peaks Not Enriched at Regulatory Regions

If your peaks are not enriched at promoters and enhancers, the possible causes include:

Incorrect genome annotation: If you are using the wrong genome build or annotation, peak annotation will be incorrect. Ensure that you are using the correct genome version for your species.

Data quality issues: Poor tagmentation specificity or high background can produce peaks that are not enriched at regulatory regions. The 2024 benchmarking study found that tagmentation specificity varied across protocols and impacted downstream analyses [<a href="#ref-4">4</a>].

Cell-type composition: If your sample contains cell types with unusual chromatin landscapes, such as senescent cells or cells with large genome rearrangements, the peak distribution may differ from expectations.

Records and Measurements for Peak Calling

Documentation Requirements

For reproducible scATAC-seq analysis, you should document the following information for each peak calling run:

Input data: The fragment file or BAM file used for peak calling, including the number of cells and the number of fragments.

Software versions: The exact versions of the peak calling software and all dependencies. Version differences can produce different results, and documenting versions is essential for reproducibility.

Parameters: All parameters used for peak calling, including thresholds, shift and extension values, and genome size. These should be recorded in a configuration file or script.

Quality metrics: The quality metrics for the input data and the called peaks, including the number of peaks, the fraction of reads in peaks, and the peak width distribution.

Output files: The locations of all output files, including peak files, coverage tracks, and quality reports.

The nf-core documentation emphasizes the importance of reproducible workflows and standardized configuration for bioinformatics pipelines [<a href="#ref-10">10</a>]. Following these principles for peak calling ensures that your analysis can be reproduced and extended by others.

Quality Metrics to Track

The following quality metrics should be tracked for each peak calling run:

Number of peaks called: This should be recorded for each sample and compared to expectations based on cell types and sequencing depth.

Fraction of reads in peaks (FRiP): This metric measures the signal-to-noise ratio of your data. Low FRiP values suggest that most reads fall outside of accessible regions, which may indicate poor data quality or overly stringent peak calling.

Peak width distribution: The median and interquartile range of peak widths should be recorded. Peaks that are much wider than expected may indicate technical artifacts.

Overlap with known regulatory regions: The fraction of peaks overlapping promoters, enhancers, and other annotated regulatory regions provides a biological validation of your peak calls.

Reproducibility metrics: If you have multiple replicates or pseudo-bulk profiles, record the overlap between peak sets and the correlation of peak scores.

Data Management

scATAC-seq data are large, and peak calling produces additional files that need to be managed. The NCBI provides data resources for storing and sharing sequencing data, including raw reads and processed data [<a href="#ref-11">11</a>]. You should plan for data storage and backup before starting your analysis.

The EMBL-EBI training resources provide guidance on data management and bioinformatics analysis [<a href="#ref-12">12</a>]. Following best practices for data organization and documentation ensures that your analysis can be reproduced and shared with collaborators.

Limitations and Interpretation

Peak Calling Is Not Standardized

One of the main challenges in scATAC-seq analysis is the lack of standardization in peak calling. Different methods and parameters can produce substantially different peak sets from the same data, and this variability complicates cross-study comparisons. A 2025 study highlighted the limitations of conventional peak-based methods, noting that they can produce results that are difficult to interpret and lack portability across datasets [<a href="#ref-2">2</a>].

The development of consensus peak references, such as the cPeaks resource derived from 624 bulk ATAC-seq datasets, represents an attempt to address this problem [<a href="#ref-6">6</a>]. By providing a standardized feature set, these references enable cross-dataset consistency and improve the reproducibility of scATAC-seq analyses.

Peak Calling Is Not the Only Approach

While peak calling is the most common approach for defining features in scATAC-seq data, it is not the only approach. Methods like SCARlink use tile-level accessibility data and avoid peak calling entirely [<a href="#ref-8">8</a>]. These methods may be preferable for certain applications, particularly when paired expression data are available.

The choice between peak-based and peak-free approaches depends on your analysis goals. If you need to identify specific regulatory elements and link them to genes, peak-based approaches provide interpretable features. If you need to predict gene expression from chromatin accessibility or analyze developmental trajectories, regression-based approaches may be more appropriate.

Interpretation Limits

Peak calling identifies regions of open chromatin, but it does not tell you what is happening at those regions. A peak may represent a promoter, enhancer, insulator, or other regulatory element, and the functional significance of a peak depends on the cell type and context. Additional analyses, such as motif enrichment, transcription factor footprinting, and integration with gene expression data, are needed to interpret the biological meaning of peaks.

The sparsity of scATAC-seq data also limits the interpretation of per-cell peak accessibility. A zero count at a peak for a particular cell does not necessarily mean that the region is inaccessible in that cell, it may simply reflect the limited sampling of the assay. This missing-data problem is a fundamental limitation of scATAC-seq and should be considered when interpreting results [<a href="#ref-3">3</a>].

Safety and Regulatory Context

Data Privacy and Sharing

scATAC-seq data from human samples may contain sensitive genetic information. If you are working with human data, you must comply with applicable privacy regulations and institutional review board requirements. The NCBI provides data resources and guidelines for sharing genomic data while protecting participant privacy [<a href="#ref-11">11</a>].

Computational Resources

Peak calling and downstream analysis of scATAC-seq data require substantial computational resources. The Galaxy Training Network provides accessible workflows and training for genomic analysis that can be run on shared or cloud infrastructure [<a href="#ref-13">13</a>]. The Carpentries lessons provide foundational training in computing and data analysis that is useful for researchers who are new to bioinformatics [<a href="#ref-14">14</a>].

Reproducibility Standards

Reproducibility is a key concern in bioinformatics analysis. The nf-core documentation describes community standards for pipeline development and usage that promote reproducibility [<a href="#ref-10">10</a>]. Following these standards, including version control, containerization, and automated testing, ensures that your peak calling analysis can be reproduced by others.

Professional Escalation Criteria

You should consider seeking expert assistance or escalating to a bioinformatics specialist in the following situations:

Unusual data characteristics: If your data show unexpected patterns, such as very low library complexity, high background, or unusual fragment size distributions, you may need expert help to diagnose the problem.

Non-model organisms: If you are working with a species without a well-annotated genome or established motif databases, you may need specialized approaches for peak calling and downstream analysis. A protocol for zebrafish analysis describes the steps for constructing species-specific motif databases [<a href="#ref-7">7</a>].

Large-scale integration: If you are integrating data from multiple experiments, protocols, or studies, you may need specialized methods for batch effect correction and joint analysis. Methods like epiConv can help, but expert guidance may be needed for complex integration scenarios [<a href="#ref-9">9</a>].

Regulatory or clinical applications: If your analysis will be used for regulatory decisions or clinical applications, you should consult with experts to ensure that your methods meet the required standards.

Frequently Asked Questions

What is the difference between peak calling in bulk ATAC-seq and scATAC-seq?

Bulk ATAC-seq aggregates chromatin accessibility across millions of cells, producing a high-depth profile where peaks are readily detectable. scATAC-seq profiles individual cells, and each cell contributes only a small number of fragments. The sparsity of single-cell data means that peaks must be called on aggregated reads from many cells, and the choice of which cells to aggregate substantially impacts the results. A benchmarking study found that scATAC-seq data have only 1 to 10 percent of peaks detected per cell, compared to 10 to 45 percent of expressed genes detected per cell in scRNA-seq data [<a href="#ref-1">1</a>].

Should I call peaks on the entire dataset or per cluster?

Calling peaks on the entire dataset maximizes sequencing depth and can detect peaks present in any cell type. Calling peaks per cluster can detect cell-type-specific peaks that would be diluted in the whole-dataset approach. The choice depends on your analysis goals. If you are interested in rare cell types or cell-type-specific regulatory elements, cluster-based peak calling is recommended. If you are primarily interested in shared regulatory features, whole-dataset peak calling may be sufficient.

What is the best peak caller for scATAC-seq data?

There is no single best peak caller for all applications. MACS2 is the most widely used and well documented, but it was not designed for single-cell data. SnapATAC and ArchR provide integrated workflows that handle the specific characteristics of scATAC-seq data. A benchmarking study found that SnapATAC outperformed other methods in separating cell populations across different coverages and noise levels [<a href="#ref-1">1</a>]. The best choice depends on your data characteristics, computational resources, and analysis goals.

How many cells do I need for reliable peak calling?

The number of cells needed for reliable peak calling depends on the sequencing depth per cell and the heterogeneity of your sample. Higher sequencing depth per cell reduces the number of cells needed, while greater heterogeneity increases the number needed to capture rare cell types. The 2024 benchmarking study used human PBMCs and found that library complexity varied substantially across protocols [<a href="#ref-4">4</a>]. As a general guideline, you should have enough cells to detect peaks in the rarest cell type of interest.

What parameters should I adjust for MACS2 peak calling?

The key parameters for MACS2 are the q-value threshold, the shift and extension values, and the genome size. The q-value threshold controls the false discovery rate, with lower values producing more stringent peak calls. The shift and extension values should be based on your fragment size distribution and the Tn5 insertion pattern. For ATAC-seq, the standard correction is to shift reads by 4 base pairs on the positive strand and 5 base pairs on the negative strand, then extend reads to a fixed length.

Can I use a reference peak set instead of calling peaks de novo?

Yes, using a reference peak set is a valid approach that can improve cross-dataset consistency. A 2026 study constructed a generic chromatin accessibility reference with about 1.4 million consensus peaks from 624 bulk ATAC-seq datasets and found that it improved cell annotation and rare cell type identification in scATAC-seq analyses [<a href="#ref-6">6</a>]. Reference-based approaches are particularly useful when comparing across datasets or building cell atlases.

How do batch effects affect peak calling?

Batch effects from technical variations between datasets can create spurious peaks or obscure real ones. Methods like epiConv can correct batch effects and improve the consistency of peak calling across datasets [<a href="#ref-9">9</a>]. If you are combining data from multiple experiments, protocols, or sequencing runs, you should consider whether batch effect correction is needed before peak calling.

What is the role of peak calling in multi-ome analysis?

In multi-ome analysis where scRNA-seq and scATAC-seq are measured in the same cells, peak calling defines the chromatin features that are linked to gene expression. Methods like SCARlink use tile-level accessibility data to link enhancers to target genes, avoiding the limitations of pairwise gene-peak correlations and dependence on peak calling [<a href="#ref-8">8</a>]. The choice of peak calling approach can substantially impact the results of multi-ome integration.

Related Bioinformatics Guides

• Benchmarking Atlas-Level Data Integration in Single-Cell Genomics: Methods and Best Practices
• Longitudinal Microbiome Data Analysis: Methods and Best Practices
• Single-Cell Sequencing Workflow: From Sample Preparation to Data Analysis
• Single-Cell Sequencing Analysis Pipeline: From Raw Data to Biological Insights
• Single-Cell RNA Sequencing Depth: A Cost-Benefit Analysis for Experimental Design

Related Clinical & Scientific Guides

• A Practical Guide to Detecting Antimicrobial Resistance Genes in Shotgun Metagenomic Data
• Computational Immunology: Modeling the Immune System
• How to Set Hard Filters for Germline Variant Calling: A Practical Guide to GATK Best Practices

References and Further Reading

This article is educational and does not replace validated analysis plans, institutional policy, clinical interpretation, or specialist review.