This workflow describes how to process large-scale scATAC-seq fragment files into a human chromatin accessibility atlas. The pipeline follows a specific hierarchy: processing individual samples first, then combining them into a backed AnnDataSet for atlas-level analysis.
Key Workflow Steps:
- Import Fragments: Convert fragment files into backed
AnnData objects using snap.pp.import_fragments. - Sample Processing: Build tile matrices, select features, and perform doublet removal (
scrublet and filter_doublets) on a per-sample basis. - Aggregation: Combine samples into a single
snap.AnnDataSet. - Atlas Analysis: Perform feature selection, spectral embedding, batch correction (
mnc_correct), and clustering (KNN, Leiden, UMAP) on the combined dataset. - Subclustering (Optional): Convert a specific cluster to an in-memory
AnnData object using .to_adata() for higher-resolution analysis.
Best Practices for Large Datasets:
- Use Backed Files: Always use backed files for both per-sample
AnnData and the combined AnnDataSet to manage memory. - Coarse Binning: Use a larger
bin_size (e.g., 5000 bp) for atlas-scale workflows to reduce memory and runtime requirements. - Unique Barcodes: Ensure cell barcodes are unique across samples (e.g., by prefixing with the sample name) before combining them into an
AnnDataSet. - Early Filtering: Apply
min_num_fragments and min_tsse during the import stage to remove low-quality barcodes early.
import os
import numpy as np
import snapatac2 as snap
# 1. Configure input and output locations.
data_dir = "/path/to/fragment/files"
output_dir = "h5ad_output"
os.makedirs(output_dir, exist_ok=True)
fragment_files = [
os.path.join(data_dir, filename)
for filename in os.listdir(data_dir)
if filename.endswith(".bed.gz")
]
# 2. Import each fragment file into a backed AnnData object.
outputs = []
for fragment_file in fragment_files:
sample_name = os.path.basename(fragment_file).split(".bed.gz")[0]
outputs.append(os.path.join(output_dir, sample_name + ".h5ad"))
adatas = snap.pp.import_fragments(
fragment_files,
file=outputs,
genome=snap.genome.hg38,
min_num_fragments=1000,
min_tsse=7,
)
# 3. Build sample-level tile matrices, select features, and remove doublets.
snap.pp.add_tile_matrix(adatas, bin_size=5000)
snap.pp.select_features(adatas)
snap.pp.scrublet(adatas)
snap.pp.filter_doublets(adatas)
# 4. Combine samples into an AnnDataSet and make barcodes unique.
adataset = snap.AnnDataSet(
adatas=[(adata.filename.split("/")[-1].split(".h5ad")[0], adata) for adata in adatas],
filename="data.h5ads",
)
adataset.obs_names = np.array(adataset.obs["sample"]) + "+" + np.array(adataset.obs_names)
# 5. Select atlas-level features and compute the initial embedding.
snap.pp.select_features(adataset, n_features=250000)
snap.tl.spectral(adataset)
# 6. Correct sample effects while preserving broad tissue groups.
adataset.obs["tissue"] = [sample.split(":")[0] for sample in adataset.obs["sample"]
snap.pp.mnc_correct(
adataset,
batch="sample",
groupby="tissue",
key_added="X_spectral",
)
# 7. Cluster and visualize the atlas.
snap.pp.knn(adataset)
snap.tl.leiden(adataset)
snap.tl.umap(adataset, n_comps=3)
snap.pl.umap(adataset, color="leiden", interactive=False)
# 8. Optional: convert to in-memory AnnData for subclustering one major cluster.
adata = adataset.to_adata()
adata_c0 = adata[adata.obs["leiden"] == "0"].copy()
snap.pp.select_features(adata_c0, n_features=250000)
snap.tl.spectral(adata_c0)
snap.pp.mnc_correct(
adata_c0,
batch="sample",
groupby="tissue",
key_added="X_spectral",
)
snap.pp.knn(adata_c0)
snap.tl.leiden(adata_c0)
snap.tl.umap(adata_c0)
snap.pl.umap(adata_c0, color="leiden", interactive=False)
snap.pl.umap(adata_c0, color="sample", interactive=False)
