跳转至

SnapATAC2 单细胞 ATAC — 超大规模单细胞 ATAC-seq 分析框架

一句话说明

SnapATAC2 是 Python 实现的高性能单细胞 ATAC-seq 分析工具，采用谱聚类（spectral clustering）替代传统的 TF-IDF+SVD，可高效处理百万细胞级别的数据，与 Scanpy 生态完全兼容。

安装与配置

# 创建专用 conda 环境
conda create -n snapatac2 python=3.10 -y
conda activate snapatac2

# 安装 SnapATAC2（从 PyPI 安装）
pip install snapatac2

# 安装可视化依赖
pip install scanpy matplotlib seaborn

# 验证安装
python -c "import snapatac2 as snap; print(snap.__version__)"

核心用法

import snapatac2 as snap               # 导入 SnapATAC2
import anndata as ad                   # AnnData 数据格式
import numpy as np

# ── 数据导入 ──────────────────────────────────────────
# 从 Cellranger ATAC 片段文件创建数据集
# 支持一次性处理多个样本（内存高效）

# 导入片段文件（fragment 文件是 ATAC 的原始数据格式）
data = snap.pp.import_data(
    'sample_fragments.tsv.gz',         # Cellranger 输出的片段文件
    genome=snap.genome.hg38,           # 使用 hg38 基因组注释
    file='sample.h5ad',                # 输出 h5ad 文件（节省内存）
    sorted_by_barcode=False,           # 片段文件是否按 barcode 排序
    min_num_fragments=500,             # 最少片段数（过滤低质量细胞）
    tempdir='./'                       # 临时目录
)
print(f"导入完成！细胞数：{data.n_obs}，基因组 bin 数：{data.n_vars}")

# ── 质控 ──────────────────────────────────────────────
# 计算质控指标（TSS 富集、片段大小分布等）
snap.metrics.tsse(data, snap.genome.hg38)   # 计算 TSS 富集分数
snap.pl.tsse(data, min_fragment=0)           # 可视化 TSS 富集分数分布

参数详解

参数函数说明
min_num_fragmentsimport_data最少片段数过滤（建议 500-1000）
genomeimport_data基因组（snap.genome.hg38/mm10）
n_featurespp.select_features选取的高变 bin 数量
n_dimstl.spectral谱嵌入维数（类似 PCA 主成分数）
distancetl.umapUMAP 距离度量
resolutiontl.leidenLeiden 聚类分辨率

实战案例

import snapatac2 as snap
import scanpy as sc
import matplotlib.pyplot as plt

# ── 完整 SnapATAC2 分析流程 ──────────────────────────

# 1. 导入多个样本（批量处理）
sample_files = {
    'sample1': 'sample1_fragments.tsv.gz',
    'sample2': 'sample2_fragments.tsv.gz',
    'sample3': 'sample3_fragments.tsv.gz'
}

# 批量导入（SnapATAC2 支持流式处理，内存效率高）
data_list = []
for sample_name, frag_file in sample_files.items():
    d = snap.pp.import_data(
        frag_file,
        genome=snap.genome.hg38,       # 人类 hg38 基因组
        file=f'{sample_name}.h5ad',    # 每个样本保存到独立文件
        min_num_fragments=500,
        sorted_by_barcode=False
    )
    d.obs['sample'] = sample_name      # 添加样本名到元数据
    data_list.append(d)
    print(f"{sample_name}: {d.n_obs} 个细胞")

# 合并多个样本
data = snap.AnnDataSet(
    adatas=[(name, d) for name, d in zip(sample_files.keys(), data_list)],
    filename='merged_dataset.h5ads'    # 保存为 h5ads 格式（多样本数据集）
)

# 2. 质控和过滤
# 计算 TSS 富集分数
snap.metrics.tsse(data, snap.genome.hg38)

# 过滤低质量细胞
snap.pp.filter_cells(
    data,
    min_counts=1000,                   # 最少 1000 个片段
    max_counts=100000,                 # 最多 10 万片段（过滤 doublet）
    min_tsse=10,                       # TSS 富集 > 10（严格阈值）
    max_tsse=None                      # 不限制最大 TSS 富集
)
print(f"过滤后细胞数：{data.n_obs}")

# 3. 特征选择（选高变 bin）
snap.pp.select_features(
    data,
    n_features=500000,                 # 选择 50 万个高变 genomic bin
    filter_lower_quantile=0.005,       # 过滤底部 0.5% 低频 bin
    filter_upper_quantile=0.005        # 过滤顶部 0.5% 高频 bin（重复区域）
)
print(f"选择的特征数：{sum(data.var.selected)}")

# 4. 谱嵌入（SnapATAC2 的核心降维算法）
# 比 TF-IDF + SVD 更能保留单细胞 ATAC 的拓扑结构
snap.tl.spectral(
    data,
    n_comps=50,                        # 计算 50 个谱嵌入维度
    features='selected'                # 使用上一步选择的特征
)

# 5. 批次校正（使用 Harmony）
snap.pp.harmony(
    data,
    batch='sample',                    # 批次列名
    max_iter_harmony=20                # Harmony 最大迭代次数
)

# 6. 构建邻居图
snap.pp.knn(
    data,
    n_neighbors=15,                    # KNN 近邻数
    use_dims=list(range(1, 30)),       # 使用第 2-30 个谱嵌入维度（跳过第1个）
    use_rep='X_spectral_harmony'       # 使用 Harmony 校正后的嵌入
)

# 7. 聚类
snap.tl.leiden(
    data,
    resolution=0.5,                    # 聚类分辨率
    random_state=42
)

# 8. UMAP 可视化
snap.tl.umap(data, random_state=42)    # 计算 UMAP 坐标

# 可视化聚类和样本分布
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
snap.pl.umap(
    data,
    color='leiden',                    # 按聚类着色
    ax=axes[0], show=False,
    title='Leiden 聚类'
)
snap.pl.umap(
    data,
    color='sample',                    # 按样本着色（检查批次效应）
    ax=axes[1], show=False,
    title='样本来源'
)
plt.tight_layout()
plt.savefig('snapatac2_umap.pdf', dpi=150)

# 9. Peak Calling（导出合并后的覆盖度并调用 MACS3）
# 将覆盖度导出为 fragment 文件（按 cluster 分组）
snap.ex.export_coverage(
    data,
    groupby='leiden',                  # 按 cluster 导出伪 bulk 覆盖度
    suffix='.bw'                       # 输出 bigwig 格式
)

# 10. 将谱嵌入转换为 Scanpy AnnData 对象（方便后续分析）
adata = data.to_adata()                # 转换为 Scanpy AnnData
sc.pp.neighbors(adata, use_rep='X_spectral')
sc.tl.umap(adata)
sc.pl.umap(adata, color='leiden')      # Scanpy 可视化

# 11. 保存结果
adata.write_h5ad('snapatac2_result.h5ad')
data.close()                           # 关闭 h5ads 文件
print("SnapATAC2 分析完成！")

常见报错与解决

报错原因解决方法
genome not found基因组注释未正确加载检查 snap.genome.hg38 是否可用
fragment file not sorted片段文件未排序sort -k1,1 -k2,2n 重新排序
谱嵌入很慢特征维度太多减少 n_features=100000
内存不足数据集太大使用 h5ads 格式流式处理
knn 步骤报错维度数不够确认 n_comps >= 30

速查表

# 安装
# pip install snapatac2

# 标准流程
data = snap.pp.import_data('fragments.tsv.gz', genome=snap.genome.hg38, file='out.h5ad')
snap.metrics.tsse(data, snap.genome.hg38)
snap.pp.filter_cells(data, min_counts=1000, min_tsse=10)
snap.pp.select_features(data, n_features=500000)
snap.tl.spectral(data, n_comps=50)
snap.pp.knn(data, n_neighbors=15, use_dims=list(range(1,30)))
snap.tl.leiden(data, resolution=0.5)
snap.tl.umap(data)

# SnapATAC2 vs ArchR vs Signac 对比
# SnapATAC2  → Python生态，超大数据集（>100万细胞），谱聚类
# ArchR      → R生态，功能最全，中大型数据集（<100万细胞）
# Signac     → R生态，Seurat整合，多组学分析首选