scATAC 基本分析思路介绍

基因组中的大部分 DNA 以染色体的形式紧密盘绕在细胞核内，基因进行复制和转录时，高度折叠的染色质结构需要暴露出 DNA 序列，这段暴露的区域就是染色质开放区域（也叫开放染色质 open chromatin），这个区域可以供转录因子和其他调控元件结合。暴露的染色质开放区域会包含启动子、增强子、绝缘子、沉默子等顺式调控元件，其可以被反式作用因子接近结合的特性，为染色质的可及性（也叫染色质开放性 chromatin accessibility ）。

DNA 开放区域与基因表达类似，在不同的组织细胞中、在发育或细胞的不同时期动态变化，呈现时空变化特异性，通过 ATAC-seq（Assay for Transposase-Accessible Chromatin using sequencing）技术可实现 bulk 样本染色质可及性的分析从而进行相关的表观研究。

随着近几年单细胞水平技术平台的发展，采用 10x Genomics Chromium Single Cell ATAC 等产品平台可以得到单个细胞中染色质开放区域图谱数据，也可采用 10x Genomics Single Cell Multiome ATAC + Gene Expression 类似的多组学平台同时获得单细胞水平的基因表达数据和表观调控数据，从而进行特定细胞类型下更加精细化基因表达的调控机制研究。

今天为大家分享 scATAC 单个平台数据的分析内容和分析思路，希望对各位老师有所帮助！

scATAC 分析思路

图1 scATAC 常见分析思路总结

1、质控

TileMatrix：fragments 可以理解为 scATAC 测序得到 DNA 开放区域片段，根据 fragments 的基因组定位信息，如以 5kb 为一个滑动窗口 bin，依次统计单个细胞在每个 bin 中 fragments 的数目，最终形成 TileMatrix，详见 ArchR分 析软件[1]。

PeakMatrix：对全部细胞的 fragments 进行 peak calling，过滤掉低可信度的 DNA 开放区域数据，保留在多个细胞中保守的染色质开放区域，即为 peak，最终生成 Peak-Barcode 可及性矩阵。大部分文献和软件（如 cellranger-atac[2] 和 Signac[3]）都是默认基于 PeakMatirx 进行分析，后续观察 marker gene 上下游区域尤其是 TSS 区域在不同细胞群间的分布差异，指导细胞类型定义。单细胞转录组数据的基因 feature 数目一般在2~3万，而 scATAC 数据 peak calling 结果一般会包含~10万左右的 peaks。

图2 博奥晶典 cellranger-atac 有效 barcode 识别

图中淡黄色部分有有效的细胞 barcodes，蓝色的点表示非细胞 barcodes，横坐标表示每个 barcode fragments 的数目，纵坐标为单个 barcode 内属于 peak 部分的 fragments 占全部 fragments 的比例。从上图可以看出 cellranger-atac 算法认为只有包含 fragmengts 要超过一定的阈值，且保守的可信 peaks 区的 fragments 的比例大于一定的阈值的 barcodes 才是有效的细胞。

图3 博奥晶典 TSS enrichment 结合 fragments 识别细胞 barcode 示例

TSSenrichment score： 也是文献中常见的一个细胞识别过滤参数，细胞保持活性需要一定数量基因的表达翻译，而基因转录时就需要开放 TSS 区域，好结合转录因子、DNA 聚合酶等转录起始复合物，一般情况下 TSS 区域的 fragments 能占到细胞全部 fragments 的25%~35%，如上图所示，一般要求有效细胞 barcode 的 TSS 富集得分>4且 barcode 包含 fragments>1000，也可设置更为严格的过滤参数，如 TSS 富集得分 >6 或者 >8。

scATAC 数据过滤的其它参数有：

nucleosome_signal：反映 DNA 开放区域的长度信息，DNA 开放区域的长度通常会在 1~2 个核小体长度（147bp）范围内，所以有些文献会过滤掉 nucleosome_signal>4的 barcodes。

blacklist regions ratio：ENCODE 数据库中收集了 blacklist 区域（基因组总反常的或者无论在二代测序的哪个实验中都是高信号的区域，排除掉这些区域可以更好的分析功能基因组数据），因此要求 barcode 中比对到 blacklist 的 fragments 的比例要低于一定的阈值。

图4 博奥晶典 scATAC 双细胞推断示例

与单细胞转录组数据双细胞过滤分析类似，如果一个 barcode 包含的 fragments 太多，则此 barcode 可能包含了两个或多个细胞，除了设定 fragments 数目阈值进行简单过滤外，同样可采用 demuxlet[7] 进行双细胞过滤，注意（图3）中 Doublet Enrichment 得分并不等同于细胞数目，仅按照此得分从大到小排序，然后按比例进行过滤。

2、转录因子可及性分析

TF-motif Matrix：一个转录因子（Transcription Factor）通常包含一个可结合 DNA 基序（motif），少数情况下一个转录因子会有两个或多个 motif，通过对近 10 万 peaks 的反向搜索，统计每个 TF-motif 在每个细胞中可以结合的 fragments 数目，从而得到单细胞层面转录因子水平的可及性矩阵 TF-motifMatrix，常见的分析软件有 chromVAR[4]，常用的JASPAR[5] 核心库包含的 TF~motif 在 600 左右，因此从 PeakMatrix 转化为 TF-motif Matrix 也可以看作为一个数据降维的过程。

图5 博奥晶典 chromVAR 转录因子变化排秩示例

注：纵坐标表示每个 TF motif 在全部细胞中的可及性变化程度，横坐标为排秩后的转录因子。

3、细胞聚类

除了常规的基于 TileMatrix、PeakMatrix 进行降维，细胞聚类分群外， Jason D. Buenrostro 2018 年在 Cell 发表的关于人类造血细胞分化中表观调控的文献，采用了 TF motif Matirx 数据进行细胞聚类分析（如下图所示），当然基于 GeneScore Matirx 参考单细胞转录分析流程进行聚类分析也是一种思路。

图6 基于 TF-motif Matrix 的细胞聚类

图7 聚类 -peak calling- 再聚类模式

初始聚类->peak calling->再聚类模式：上图为 Ansuman T. Satpathy 在 NATURE BIOTECHNOLOGY 发表的关于人类免疫细胞发育的表观图谱文献的分析思路。

简介如下：

1. 保留 Unique fragments >1000且 TSS 区域富集得分>8的 barcode 为有效细胞；

2. 基因组位置上每 2.5kb 为一个 bin，生成 TileMatrix，选择变化最大最高的 top2 万 bins 进行初始聚类；

3. 依据初始分类结果，每 200 个细胞合并成一个 pseudo bulk ATAC 数据；

4. 对每个 Cluster 的 pseudo bulk ATAC 数据分别 call peaks，最后合并全部的 peaks 数据；

5. 生成最终的 PeakMatrix、降维、聚类、细胞类型注释及后续分析。

采用类似思路三轮聚类作为最终结果的还有 Alexandro E. Trevino 2020 年在 Cell 上发表的关于人类大脑皮层发育的表观基因表达调控研究的文献。此分析思路总结如下：一方面，合并单细胞数据生成 pseudo bulk ATAC 数据，可以有效减少由于单个细胞 ATAC 数据捕获不足造成的数据缺失问题；另一方面，每个初始 cluster 重新 call peak，可以得到细胞群间特异的 peaks，便于发现细胞特异的表观调控机制。

图8 博奥晶典初始聚类 -peak calling 结果示例

4、细胞类型注释

Marker 基因方法：在 scATAC 细胞分群结果的基础上，进行类间差异 peak 可及性分析，根据 marker 基因 TSS 区域或整个基因区域的可及性类间差异，对细胞群进行人工注释。

基因活性方法：GeneScoreMatrix 基因活性得分矩阵，考虑 peak 与 TSS/genebody 的距离、测序深度等因素采用广义线性模型将 PeakMatrix 矩阵转换为 GeneScore Matrix，常用的软件有 Cicero[6]、Signac[3] 和 ArchR[1]。基因活性得分可以理解为基于表观数据预测出的基因表达值，可以通过观察 markergene 在细胞群中活性得分值的特异分布，进行细胞类型的人工注释。也可采用 SingleR 和 Cellassign 等软件进行细胞类型注释。

图9 博奥晶典 Gene Score SingleR 细胞类型注释示例

联合 scRNA 转移类标签方法：根据实验设计寻找最接近且已经注释好细胞类型的 scRNA 数据，基于 Gene Score 数据，采用 Seurat TransferData 的方法，获得 scATAC 的细胞类型注释。

5、共可及性分析

peak-peak 共可及性分析，常用来研究基因组上的顺势作用元件间（如启动子和增强子）的相互作用，可采用 Cicero[6]包实现。

图10 博奥晶典 peak-peak 共可及性分析示例

如上图所示，横坐标是预测到的共可及性网络区间内的 peaks 和注释到的基因，纵坐标是共可及性值，图中 0.2 处的虚线是存在共可及性关系的阈值线，紫色弧线是 peak 之间的相互作用连接线。

6、拟时间分析

基于基因活性得分矩阵，采用 Monocle[8] 软件可实现 scATAC 数据的拟时间分析（又称趋势分析）。下图为 scATAC 的细胞在 Monocle 降维结果上 pseudotime 的分布图。

图11 博奥晶典 scATAC 拟时间分析示例

参考文献

1. Granja JM, Corces MR et al., ArchR is a scalable software package for integrative single-cell chromatin accessibility analysis. Nature Genetics (2021).
2. https://support.10xgenomics.com/single-cell-atac/software/pipelines/latest/algorithms/overview.
3. Tim Stuart, Avi Srivastava, et al. Single-cell chromatin state analysis with Signac. Nat Methods . 2021 Nov;18(11):1333-1341. HelenaTodorov. et al.
4. Schep AN, Wu B, Buenrostro JD and Greenleaf WJ. chromVAR: inferring transcription-factor-associated accessibility from single-cell epigenomic data. Nature Methods. doi: 10.1038/nmeth.4401.
5. Khan, A. JASPAR 2018: update of the open-access database of transcription factor binding profiles and its web framework. Nucleic Acids Res. 2018; 46: D260–D266.
6. Hannah A. Pliner, Jonathan S. Packer, et al. Cicero Predicts cis-Regulatory DNA Interactions from Single-Cell Chromatin Accessibility Data. Molecular Cell. TECHNOLOGY| VOLUME 71, ISSUE 5, P858-871. E8, SEPTEMBER 06, 2018
7. Kang et al.,Multiplexed droplet single-cell RNA-sequencing using natural genetic variation. Nature Biotechnology 2017.
8. Cole Trapnell, Davide Cacchiarelli. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nature Biotechnology 2014.

来源：博奥晶典科研服务