Peak Calling Pseudo Bulk Aggregation

peak-calling-pseudo-bulk-aggregation

Summary

Aggregate single-cell ATAC-seq fragments into pseudo-bulk samples grouped by cluster assignment, then call peaks using MACS3 to identify chromatin-accessible regions that are consistent within cell populations. This approach recovers peaks from sparse single-cell data by leveraging population-level signal.

When to use

After clustering single-cell ATAC-seq data (e.g., via Leiden clustering on spectral embeddings), use this skill to identify peaks within each cluster. Triggering conditions: (1) you have sparse, per-cell insertion counts organized in a tile matrix; (2) cells have been assigned to discrete clusters or cell-type annotations; (3) you want to call peaks that are reproducible within cell populations rather than in individual cells; (4) downstream analyses require cluster-specific peak lists for differential accessibility or motif enrichment.

When NOT to use

• Input is already a pre-computed peak set or peak matrix; peak calling has been performed elsewhere.
• Cluster assignments are unreliable or highly fragmented (many clusters with <100 cells each), reducing pseudo-bulk signal below MACS3 detection threshold.
• You need single-cell resolution peaks (e.g., to study peak heterogeneity within a cluster); pseudo-bulk aggregation will obscure cell-to-cell variation.

Inputs

• Single-cell ATAC-seq fragment file (imported and stored as paired-end insertions in AnnData)
• Cluster assignments for each cell (from tl.leiden or equivalent clustering method)
• Tile matrix or other sparse count matrix representing chromatin accessibility

Outputs

• Per-cluster peak list (BED format or equivalent, called by MACS3)
• Merged peak set (union or consensus of all cluster-specific peaks)
• Peak matrix (annotation of cells by peak presence/absence, via pp.make_peak_matrix)

How to apply

Group fragments by cluster assignment and aggregate them into a single pseudo-bulk BAM or fragment file per cluster. Pass each pseudo-bulk file to tl.macs3 with appropriate parameters (e.g., genome size, p-value threshold, minimum read count). MACS3 will perform peak calling on the aggregated signal from all cells within a cluster, improving signal-to-noise compared to single-cell calling. After calling peaks independently for each cluster, merge overlapping peaks across clusters using tl.merge_peaks to create a unified peak set. This strategy exploits the population-level chromatin accessibility pattern while avoiding the sparsity problem of individual-cell peak calling. Validate peak calls by comparing cluster assignments and UMAP coordinates before and after peak calling to ensure clustering quality is preserved.

Related tools

• SnapATAC2 (Provides tl.macs3 for pseudo-bulk peak calling and tl.merge_peaks for peak merging; orchestrates fragment import, tile matrix generation, clustering, and peak calling in a unified pipeline.) — https://github.com/scverse/SnapATAC2
• MACS3 (Peak calling algorithm applied to pseudo-bulk aggregated fragments; detects chromatin-accessible regions from enriched read pileups.)
• Leiden (Clustering algorithm used to assign cells to populations before pseudo-bulk aggregation.)
• UMAP (Visualization of cell embeddings; used to validate that cluster quality and structure are preserved after peak calling.)

Examples

import snapatac2 as snap; adata = snap.pp.import_fragments(fragments_file); snap.pp.add_tile_matrix(adata); snap.tl.spectral(adata); snap.tl.leiden(adata); snap.tl.macs3(adata, groupby='leiden'); snap.tl.merge_peaks(adata)

Evaluation signals

• Cluster assignments and UMAP coordinates are stable and reproducible before and after peak calling (consistent with reference tutorial outcomes).
• Number of peaks per cluster is consistent with expected chromatin complexity (typically 10k–100k peaks for PBMC populations).
• Overlapping peaks across clusters are successfully merged; merged peak set size is reasonable relative to individual cluster peak counts.
• Peak calls are enriched at known regulatory regions (e.g., transcription start sites, enhancers); can be validated by motif enrichment or comparison to published ATAC-seq datasets.
• Per-cluster peak matrices can be generated successfully via pp.make_peak_matrix and have expected dimensions (cells × merged peaks).

Limitations

• Pseudo-bulk aggregation requires sufficient cell counts per cluster; clusters with very few cells (<50–100) may have insufficient signal for robust peak calling.
• Peak merging strategy (overlap threshold, consensus vs. union) can affect downstream interpretation; choice depends on whether you want permissive or strict peak sets.
• MACS3 parameters (p-value, minimum read count, genome size) are critical; suboptimal tuning may yield false positives or false negatives. Default parameters may not suit all tissues or species.
• Pseudo-bulk averaging obscures within-cluster heterogeneity; cells with divergent chromatin states will be homogenized in the aggregated signal.

Evidence

• [intro] End-to-end analysis pipeline for single-cell ATAC-seq data including preprocessing, dimension reduction, clustering, data integration, and peak calling: "End-to-end analysis pipeline for single-cell ATAC-seq data, including preprocessing, dimension reduction, clustering, data integration, peak calling"
• [other] Perform clustering using tl.leiden with default resolution parameter. Call peaks using tl.macs3 in pseudo-bulk mode grouped by cluster assignment. Merge overlapping peaks using tl.merge_peaks.: "Perform clustering using tl.leiden with default resolution parameter. Call peaks using tl.macs3 in pseudo-bulk mode grouped by cluster assignment. Merge overlapping peaks using tl.merge_peaks"
• [methods] Peak calling tools in SnapATAC2 support MACS3 and peak merging for single-cell ATAC-seq analysis: "tl.macs3, tl.merge_peaks for peak calling"
• [readme] Blazingly fast preprocessing tools for BAM to fragment files conversion and count matrix generation: "Blazingly fast preprocessing tools for BAM to fragment files conversion and count matrix generation"