一个 pbmc5k的scATAC-seq数据基本分析：snapATAC2软件

在 scATAC-seq 专辑中，我前面已经简单的学习过了Signac软件，ArchR 软件。这两个软件都是大佬开发的，然后都是R语言生态。今天来学习一个python生态的分析软件。于 08 January 2024 发表在 Nature Methods 上面，文献标题为《A fast, scalable and versatile tool for analysis of single-cell omics data》。

我们技能树最新一期生信入门就在11月3号开课，还没上车的可以看一看瞄一瞄，详细介绍页面可以点击：生信入门&数据挖掘线上直播课11月班，完全适合0基础小白。

这个软件的网址：

https://scverse.org/SnapATAC2/index.html

https://github.com/kaizhang/single-cell-benchmark/.

https://github.com/kaizhang/SnapATAC2/.

还搜到了一个scATAC-seq的发展史，也差不多有10年了哇。简单熟悉一下：https://gitee.com/booew/SnapATAC/blob/master/notebooks/experiemnt_timeline.md

这个软件发表的文献处理了好多数据，后面要是做一些什么分析缺测试数据，可以来看看：

然后软件的原理，这里就自己下去默默啃读了。

0.安装

总体来说，不是特别好安装，然后我看到这个软件也在2024年后也没有了更新。心里不由得咯噔一下，后面分析可能会遇到很多问题，但是软件不更新！！！

开发的工具即使再怎么好，如果不维护更新，生命周期一定是非常短的。这里感叹一下 clusterProfiler 包这么多年了，还是一直更新特别快！！！你能体会到开发者的伟大了吗！

conda create -n atac python=3.9 -y
conda activate atac
pip install snapatac2
pip install jupyterlab -i https://mirrors.tuna.tsinghua.edu.cn/pypi/web/simple
conda install -y pandas=2.1.1

我这里运行的教程：https://scverse.org/SnapATAC2/tutorials/pbmc.html

1.加载模块

导入模块：

import snapatac2 as snap

snap.__version__
# 2.8.0

下载数据：'https://cf.10xgenomics.com/samples/cell-atac/2.0.0/atac_pbmc_5k_nextgem/atac_pbmc_5k_nextgem_fragments.tsv.gz'

返回的是数据路径。

# Input files
fragment_file = snap.datasets.pbmc5k()
fragment_file 
# /home/zhangj/.cache/snapatac2/atac_pbmc_5k.tsv.gz

2.预处理

过导入 fragment 文件并使用 pp.import_fragments() 函数计算基础质控指标来开始数据预处理。

该函数会将 fragment 数据压缩并存储至 AnnData 对象中以供后续使用，在此过程中，各类质控指标（如TSS富集度和每个细胞的唯一片段数量）也会被计算并同步存储至anndata对象中。

%%time
data = snap.pp.import_fragments( fragment_file, chrom_sizes=snap.genome.hg38, file="pbmc.h5ad", sorted_by_barcode=False)
data

会在当前目录下生成一个 pbmc.h5ad 文件。计算的qc指标包括：

• n_fragment：the number of unique fragments per cell
• frac_dup：fraction of duplicated reads
• frac_mito：fraction of mitochondrial read

片段长度分布：ATAC-seq数据中典型的片段长度分布呈现以下关键特征：

无核小体区域（NFR）：大部分片段较短，通常长度在80-300碱基对之间。这些短片段代表开放染色质区域，其中DNA相对可及且未与核小体结合。这些区域通常被称为核小体游离区域，对应活跃转录调控区域。

单核小体峰：在约147-200碱基对处通常会出现一个片段长度分布峰，这对应于单个核小体包裹DNA的大小。这些片段由转座酶遇到核小体时切割DNA产生，因DNA受到保护而形成相同长度的片段。

双核小体峰：除单核小体峰外，还可能在约300-400碱基对处观察到较小的峰，对应双核小体片段。这些片段产生于转座酶在两个相邻核小体之间切割DNA时。

大片段：可能会观察到一些大于500碱基对的大片段。这些片段可能来源于开放染色质的较长延伸区域或技术性伪影。

snap.pl.frag_size_distr(data, interactive=False)

TSS富集度

ATAC-seq 数据中的 TSS 富集度指标表明基因转录起始位点周围的染色质可及性或片段化程度有所增加。这通常意味着启动子区域处于开放和可及状态——这些区域通常没有核小体覆盖，并参与活跃的基因转录过程。

研究人员常将 TSS 富集度作为此类实验的质量控制指标。如果在 TSS 区域周围观察到明显且突出的 reads 或信号富集，则表明实验成功捕获了相关基因组特征，并可能产生具有生物学意义的结果。反之，若缺乏 TSS 富集则可能提示实验质量或数据处理存在问题。

单个细胞的 TSS 富集度评分可通过 metrics.tsse() 函数进行计算。

这里会用到一个 gencode.v41.basic.annotation.gff3.gz 文件，我是直接使用代码中提示的网址下载到本地然后读取的，代码自己下载会有点慢：

%%time
snap.metrics.tsse(data,  gene_anno="gencode.v41.basic.annotation.gff3.gz")
data

# 绘制每个细胞的TSS富集度评分与唯一片段数量的散点图
snap.pl.tsse(data, interactive=False)

散点图右上角的细胞代表有效或高质量细胞，而左下角的细胞则对应低质量细胞或空液滴。基于此分布图，可以设定最低 TSS 富集度为10且 最低片段数量为5,000，以此标准对细胞进行筛选。

# 过滤
# %%time
snap.pp.filter_cells(data, min_counts=5000, min_tsse=10, max_counts=100000)
data

过滤后的数据：

接下来：

1.创建一个以细胞为行、基因组区间为列的矩阵，该矩阵统计全基因组范围内500bp区间内的插入位点计数。

2.特征筛选，n_features参数决定了后续分析步骤中使用的特征（即基因组区间）数量。增加特征数量能提升分辨率并揭示更精细的细节，但也可能引入噪声。为了获得最佳结果，建议尝试不同的n_features参数值，以找到最适合您数据集的设置

# 创建矩阵
snap.pp.add_tile_matrix(data)

# 特征筛选。筛选结果将存储在data.var['selected']
snap.pp.select_features(data, n_features=250000)

3.Doublet去除

采用定制化的 Scrublet 算法来识别潜在的多重捕获细胞。调用pp.scrublet()函数将为每个细胞标注成为多重捕获细胞的概率值。随后，可以通过pp.filter_doublets函数来剔除这些多重捕获细胞。

snap.pp.scrublet(data)
# 过滤
snap.pp.filter_doublets(data) # Detected doublet rate = 2.936%
data

4.降维聚类分群

采用谱嵌入方法进行降维处理。降维结果将存储在 data.obsm['X_spectral'] 中。然后使用UMAP将细胞嵌入二维空间以实现可视化。

# 降维
snap.tl.spectral(data)
snap.tl.umap(data)

# 聚类
snap.pp.knn(data)
snap.tl.leiden(data)

# umap可视化
snap.pl.umap(data, color='leiden', interactive=False, height=500)

得到了13个clusrer，数据结构现在如下：

data

5.细胞亚群注释

创建基于基因活性的细胞矩阵

在获得细胞聚类结果后，对聚类进行注释并将其归类至已知细胞类型。为此，需要首先使用pp.make_gene_matrix函数计算每个细胞的基因活性。

与pp.import_data类似，pp.make_gene_matrix函数允许选择性地提供文件名参数用于存储AnnData对象。

# 基因活性矩阵
gene_matrix = snap.pp.make_gene_matrix(data, gene_anno="gencode.v41.basic.annotation.gff3.gz")
gene_matrix

插补处理

基因活性矩阵通常具有高度稀疏性。为便于可视化和标记基因识别，采用MAGIC算法进行数据插补与平滑处理。后续分析步骤与单细胞RNA-seq分析流程相似，可以使用scanpy工具包来完成这一过程。

首先进行基因筛选、数据标准化和数据转换，随后调用sc.external.pp.magic函数完成插补处理。执行此操作需预先安装magic-impute包：pip install magic-impute。

import scanpy as sc

sc.pp.filter_genes(gene_matrix, min_cells= 5)
sc.pp.normalize_total(gene_matrix)
sc.pp.log1p(gene_matrix)

# 插补
sc.external.pp.magic(gene_matrix, solver="approximate")

# Copy over UMAP embedding
gene_matrix.obsm["X_umap"] = data.obsm["X_umap"]

# 可视化若干标记基因的活性表达
marker_genes = ['MS4A1', 'CD3D', 'LEF1', 'NKG7', 'TREM1', 'LYZ', 'PPBP']
sc.pl.umap(gene_matrix, use_raw=False, color=["leiden"] + marker_genes)

上面的marker细胞还挺明显的！

没想到这个软件这里做亚群注释不需要使用单细胞作为参考，然后这个基因活性矩阵的方法看起来还不错！

6.保存数据

注意，在关闭Python进程前，务必关闭所有以后端模式打开的AnnData对象，以避免HDF5文件损坏！

data.close()
data

# 保存所有内存中的AnnData对象
gene_matrix.write("pbmc5k_gene_mat.h5ad", compression='gzip')

使用了这好几个软件（ArchR、Signac、snapATAC）之后，大家的分析步骤都差不多，就看后面的大数据分析速度了。

然后我发现这些软件现在基本上都不怎么维护更新了。不知道后面还有没有大佬会开发一个像scanpy或者seurat一样使用广泛且能持续更新的软件哇！

再回答一个群友的问题

再回答一个群友的问题：我们的ATAC-seq群里之前看到过一个问题，GSM8539201 这个数据没有 fragment文件怎么分析？

现在我的回答是：分析不了，我上次熬夜看了好几个软件，发现他们都需要这个 fragment 文件，snapATAC这个也是以 fragment作为的输入，然后教程中提到了：

If you do not have a fragment file for your experiment, you can make one from a BAM file, see pp.make_fragment_file(

fragement详细说明：https://support.10xgenomics.com/single-cell-atac/software/pipelines/latest/output/fragments

10X scATAC-seq的结果 fragment files 含义