# SnapATAC2: A Python package for single-cell epigenomics analysis

SnapATAC2 is a flexible, versatile, and scalable single-cell omics analysis framework, featuring:

- Scale to more than 10 million cells.
- Blazingly fast preprocessing tools for BAM to fragment files conversion and count matrix generation.
- Matrix-free spectral embedding algorithm that is applicable to a wide range of single-cell omics data, including single-cell ATAC-seq, single-cell RNA-seq, single-cell Hi-C, and single-cell methylation.
- Efficient and scalable co-embedding algorithm for single-cell multi-omics data integration.
- End-to-end analysis pipeline for single-cell ATAC-seq data, including preprocessing, dimension reduction, clustering, data integration, peak calling, differential analysis, motif analysis, regulatory network analysis.
- Seamless integration with other single-cell analysis packages such as Scanpy.
- Implementation of fully backed AnnData.

## API Reference

Below is a list of available functions in SnapATAC2. For detailed API usage, please use python function `help(function_name)` to retrieve the docstring.

- snapatac2.AnnData: An annotated data matrix.
- snapatac2.AnnDataSet: Similar to `AnnData`, `AnnDataSet` contains annotations of observations `obs` (`obsm`, `obsp`), variables `var` (`varm`, `varp`), and unstructured annotations `uns`. Additionally it provides lazy access to concatenated component AnnData objects, including `X`, `obs`, `obsm`, `obsp`.
- snapatac2.PyDNAMotif: Represent a DNA motif as a position weight matrix.
- snapatac2.PyDNAMotifScanner: Scan DNA sequences with a motif-specific scanner.
- snapatac2.PyDNAMotifTest: No description available.
- snapatac2.concat: Concatenates AnnData objects.
- snapatac2.datasets.Meuleman_2020: Fetch grouped transcription factor motifs from Meuleman 2020.
- snapatac2.datasets.Path: PurePath subclass that can make system calls.
- snapatac2.datasets.PyDNAMotif: Represent a DNA motif as a position weight matrix.
- snapatac2.datasets.cis_bp: Fetch CIS-BP transcription factor motifs for motif analysis.
- snapatac2.datasets.colon: Fetch five transverse colon scATAC-seq fragment datasets.
- snapatac2.datasets.cre_HEA: Fetch the curated human colon cis-regulatory element BED file.
- snapatac2.datasets.load: Fetch ``dataset`` from the registry and return whatever its loader produces.
- snapatac2.datasets.pbmc10k_multiome: Fetch the 10x Genomics 10k PBMC multiome example dataset.
- snapatac2.datasets.pbmc500: Fetch the 10x Genomics 500 PBMC scATAC-seq example dataset.
- snapatac2.datasets.pbmc5k: Fetch the 10x Genomics 5k PBMC scATAC-seq example dataset.
- snapatac2.datasets.typechecked: Instrument the target function to perform run-time type checking.
- snapatac2.export.Path: PurePath subclass that can make system calls.
- snapatac2.export.export_coverage: Export grouped genome-wide coverage tracks.
- snapatac2.export.export_fragments: Export grouped fragments to BED-like files.
- snapatac2.export.get_file_format: Infer genomic track format and compression from a file suffix.
- snapatac2.genome.Callable: No description available.
- snapatac2.genome.Decompress: Processor that decompress a file and returns the decompressed version.
- snapatac2.genome.Genome: A class that encapsulates information about a genome, including its FASTA sequence, its annotation, and chromosome sizes.
- snapatac2.genome.Path: PurePath subclass that can make system calls.
- snapatac2.genome.load: Fetch ``dataset`` from the registry and return whatever its loader produces.
- snapatac2.get_write_options: Get the current default write configuration.
- snapatac2.metrics.Genome: A class that encapsulates information about a genome, including its FASTA sequence, its annotation, and chromosome sizes.
- snapatac2.metrics.Path: PurePath subclass that can make system calls.
- snapatac2.metrics.frag_size_distr: Compute the dataset-level fragment size distribution.
- snapatac2.metrics.frip: Compute fraction of reads or insertions in selected regions.
- snapatac2.metrics.summary_by_chrom: Compute per-cell summary statistics for each chromosome.
- snapatac2.metrics.tsse: Compute transcription start site enrichment for each cell.
- snapatac2.plotting.coverage: Plot coverage tracks for grouped cells across one genomic region.
- snapatac2.plotting.frag_size_distr: Plot the fragment size distribution.
- snapatac2.plotting.motif_enrichment: Plot motif enrichment scores across groups.
- snapatac2.plotting.network_scores: Plot average network edge scores by distance-to-TSS bin.
- snapatac2.plotting.regions: Plot grouped accessibility over selected peak regions.
- snapatac2.plotting.spectral_eigenvalues: Plot spectral embedding eigenvalues and mark the elbow.
- snapatac2.plotting.tsse: Plot TSS enrichment against unique fragment counts.
- snapatac2.plotting.umap: Plot a two- or three-dimensional UMAP embedding.
- snapatac2.preprocessing.Path: PurePath subclass that can make system calls.
- snapatac2.preprocessing.add_tile_matrix: Generate a cell-by-genomic-bin count matrix.
- snapatac2.preprocessing.call_cells: Call valid cells from feature counts using the OrdMag algorithm.
- snapatac2.preprocessing.filter_cells: Filter cells by fragment-count and TSS-enrichment QC thresholds.
- snapatac2.preprocessing.filter_doublets: Remove cells classified as doublets.
- snapatac2.preprocessing.filter_kwargs: Keep only keyword arguments accepted by a target function.
- snapatac2.preprocessing.harmony: Correct batch effects in an embedding with Harmony.
- snapatac2.preprocessing.import_contacts: Import chromatin contacts into an AnnData object.
- snapatac2.preprocessing.import_fragments: Import fragment files and compute basic QC metrics.
- snapatac2.preprocessing.import_values: Import base-pair values into an AnnData object.
- snapatac2.preprocessing.knn: Build a Euclidean k-nearest-neighbor graph for observations.
- snapatac2.preprocessing.make_fragment_file: Convert a BAM file into a sorted fragment file.
- snapatac2.preprocessing.make_gene_matrix: Generate a cell-by-gene activity matrix.
- snapatac2.preprocessing.make_peak_matrix: Generate a cell-by-peak count matrix.
- snapatac2.preprocessing.mnc_correct: Correct batch effects with centroid-based mutual nearest neighbors.
- snapatac2.preprocessing.recipe_10x_metrics: Generate 10x-style ATAC QC metrics from a raw BAM file.
- snapatac2.preprocessing.scanorama_integrate: Integrate batch-specific embeddings with Scanorama.
- snapatac2.preprocessing.scrublet: Score ATAC-seq cells for doublet likelihood with Scrublet.
- snapatac2.preprocessing.select_features: Select informative genomic features for downstream analysis.
- snapatac2.read: Read `.h5ad`-formatted hdf5 file.
- snapatac2.read_10x_mtx: Read a 10x Genomics MTX directory into an AnnData object.
- snapatac2.read_dataset: Read AnnDataSet object.
- snapatac2.read_motifs: Read DNA motifs from a MEME format file.
- snapatac2.read_mtx: Read Matrix Market file.
- snapatac2.set_write_options: Set the default write configuration for all subsequent write operations.
- snapatac2.tools.LinkData: No description available.
- snapatac2.tools.NodeData: No description available.
- snapatac2.tools.add_cor_scores: Add Spearman correlation scores to network edges.
- snapatac2.tools.add_regr_scores: Add regression-based importance scores to network edges.
- snapatac2.tools.add_tf_binding: Add motif-to-region edges to a regulatory network.
- snapatac2.tools.aggregate_X: Aggregate `.X` values across cells or cell groups.
- snapatac2.tools.aggregate_cells: Assign cells to pseudo-cell groups by iterative clustering.
- snapatac2.tools.dbscan: Cluster cells with DBSCAN.
- snapatac2.tools.diff_test: Test regions for differential accessibility between two cell groups.
- snapatac2.tools.hdbscan: Cluster cells with HDBSCAN.
- snapatac2.tools.init_network_from_annotation: Build a region-to-gene network from gene annotations.
- snapatac2.tools.kmeans: Cluster cells with k-means.
- snapatac2.tools.leiden: Cluster cells with the Leiden community detection algorithm.
- snapatac2.tools.leiden_sweep: Score Leiden clustering across multiple resolutions.
- snapatac2.tools.link_tf_to_gene: Create a transcription-factor-to-gene network.
- snapatac2.tools.macs3: Call open chromatin peaks with MACS3.
- snapatac2.tools.marker_regions: Select marker regions for each group by z-score enrichment.
- snapatac2.tools.merge_peaks: Merge group-specific peak calls into a non-overlapping peak set.
- snapatac2.tools.motif_enrichment: Test transcription factor motifs for enrichment in region sets.
- snapatac2.tools.multi_spectral: Compute one spectral embedding from multiple modalities.
- snapatac2.tools.pagerank: Rank regulator nodes with personalized PageRank.
- snapatac2.tools.prune_network: Filter nodes and edges from a network.
- snapatac2.tools.smooth: Smooth cell-by-feature values over a nearest-neighbor graph.
- snapatac2.tools.spectral: Compute a spectral embedding with Laplacian Eigenmaps.
- snapatac2.tools.transfer_labels: Predict missing cell labels from labeled neighbors.
- snapatac2.tools.umap: Compute a UMAP embedding from an existing representation.


## Tutorials

### Standard PBMC scATAC-seq Pipeline

Use this workflow to analyze the 10x Genomics 5k PBMC scATAC-seq fragment dataset with SnapATAC2. The pipeline imports fragments, computes QC metrics, filters cells, creates a tile matrix, selects features, removes doublets, computes embeddings, clusters cells, creates a gene activity matrix, and saves outputs.

#### Prerequisites

- Install SnapATAC2 and optional plotting/analysis dependencies used in the workflow.
- Use backed mode (`file="pbmc.h5ad"`) for large datasets or memory-constrained environments.
- Close backed AnnData objects before the Python process exits to avoid HDF5 file corruption.

#### Complete Script

```python
import snapatac2 as snap

# 1. Download the 10x PBMC fragment file.
fragment_file = snap.datasets.pbmc5k()

# 2. Import fragments into a backed AnnData object and compute basic QC metrics.
# Use `sorted_by_barcode=True` if fragment file is sorted by barcode names.
# Cellranger does not sort barcode names so I set to `False` here.
data = snap.pp.import_fragments(
    fragment_file,
    chrom_sizes=snap.genome.hg38,
    file="pbmc.h5ad",
    sorted_by_barcode=False,
)

# 3. Inspect fragment-size distribution and compute TSS enrichment.
snap.pl.frag_size_distr(data, interactive=False)
snap.metrics.tsse(data, snap.genome.hg38)
snap.pl.tsse(data, interactive=False, show=False, out_file="tsse.png")

# 4. Keep high-quality cells.
snap.pp.filter_cells(
    data,
    min_counts=5000,
    min_tsse=10,
    max_counts=100000,
)

# 5. Generate a 500-bp genome-wide tile matrix in data.X.
snap.pp.add_tile_matrix(data)

# 6. Select accessible features for downstream modeling.
snap.pp.select_features(data, n_features=250000)

# 7. Detect and remove doublets.
snap.pp.scrublet(data)
snap.pp.filter_doublets(data)

# 8. Compute dimensionality reduction and UMAP embedding.
snap.tl.spectral(data)
snap.tl.umap(data)

# 9. Build the KNN graph and cluster cells with Leiden.
snap.pp.knn(data)
snap.tl.leiden(data)
snap.pl.umap(data, color="leiden", interactive=False, show=False, out_file="umap.png")

# 10. Create an in-memory gene activity matrix for marker-gene annotation.
gene_matrix = snap.pp.make_gene_matrix(data, snap.genome.hg38)

# 11. Close or save outputs.
gene_matrix.write("pbmc5k_gene_mat.h5ad", compression="gzip")
data.close()
```

#### Workflow Notes

- `snap.pp.import_fragments(..., file="pbmc.h5ad")` creates a backed AnnData object. This streams data to disk and is preferred for large fragment files.
- `snap.metrics.tsse` adds per-cell TSS enrichment to `data.obs["tsse"]` and library-level TSS metrics to `data.uns`.
- `snap.pp.filter_cells` subsets the AnnData object in place by default. In this tutorial, cells are retained when they have at least 5000 fragments, TSS enrichment of at least 10, and at most 100000 fragments. The best threshold varies depends on the protocol and sample.
- `snap.pp.add_tile_matrix` writes a cell-by-bin count matrix to `data.X`.
- `snap.pp.select_features` stores the feature mask in `data.var["selected"]`; it does not subset the matrix directly.
- `snap.pp.scrublet` stores doublet scores and probabilities in `data.obs`; `snap.pp.filter_doublets` removes predicted doublets.
- `snap.tl.spectral` writes the spectral embedding to `data.obsm["X_spectral"]`; `snap.tl.umap` writes UMAP coordinates to `data.obsm["X_umap"]`.
- `snap.tl.leiden` stores cluster labels in `data.obs["leiden"]`.
- `snap.pp.make_gene_matrix` returns a new AnnData object with cells as observations and genes as variables.

#### Key Outputs

- `pbmc.h5ad`: backed AnnData file containing imported fragments, QC metrics, tile matrix, selected features, doublet annotations, embeddings, graph, and clusters.
- `pbmc5k_gene_mat.h5ad`: Gene activity matrix.
- `data.obs`: expected to contain fields such as `n_fragment`, `frac_dup`, `frac_mito`, `tsse`, `doublet_probability`, `doublet_score`, and `leiden`.
- `data.var["selected"]`: boolean mask for selected genomic bins.
- `data.obsm["X_spectral"]` and `data.obsm["X_umap"]`: dimensionality-reduction results.

### Other Tutorials

- "examples/atlas.md": Atlas-Scale Human Chromatin Accessibility Analysis
- "examples/modality.md": Joint Embedding of PBMC Multiome RNA and ATAC Data
