单细胞ATAC亚群分析

上一期，跟大家简单介绍了下单细胞ATAC的背景知识点及其10x ATAC基础数据的获取方式。接下来就带大家从fragment.csv、singlecell.csv、peaks matrix等数据出发，做单细胞ATAC的亚群分析。

与单细胞转录组类似，单细胞ATAC的分析流程也主要包括细胞质控、peaks标准化及其降维分群、marker基因的鉴定等几个步骤。常用的单细胞ATAC分析流程软件包含 cell-ranger-atac、Signac和ArchR等。

01

细胞质控

单细胞ATAC的质控点一般包含以下几个方面：样本重复（biological replicates），bulkATAC vs scATAC的相关性、fragment length distribution、per nucleus read-depth、transcription start site (TSS) enrichment、双细胞比例等。

图1：单细胞ATAC质控点

1.1 低质量细胞的过滤

前面提到的样本相关性和fragments的长度分布主要是从整体水平上检查我们的单个样本数据的可靠性。

而要去掉不符合质控的细胞，我们主要从fragments 数目&TSS enrichment score这两点出发。

👉 fragments 数目：一般指单个细胞（barcode）所属的total fragments数目。这个不同的软件具体的定义不同，比如cell-ranger-atac和Signac指peaks所属区域的fragments 数目，其中singlecell.csv文件中peak_region_fragments列便是指fragments 数目，而ArchR是指全基因组所有的fragments 数（这个跟该软件的分析策略有关，后面会提到）。

👉 TSS enrichment score：相当于计算每个细胞的信噪比（signal-to-background ratio），ENCODE项目已经定义了一个ATAC-seq目标评分，该评分基于TSS中心的片段与TSS侧翼区域的片段的比例(见https://www.encodeproject.org/data-standards/terms/)。较差的ATAC-seq实验通常会有较低的TSS浓缩分数。Signac软件可以用TSSEnrichment()函数为每个细胞计算TSS enrichment score，而ArchR包也是利用类似的原理createArrowFiles()函数在读取基础数据时就为每个细胞计算了该指标。

备注：fragments 数目&TSS enrichment score的阈值不仅与所用软件具体的计算公式有关（不同的软件具体的参数可能不同），也与自己数据的实际情况有关。比如哺乳动物和植物的单细胞ATAC数据TSS enrichment score就不能用相同的指标cutoff来衡量，一般来说哺乳动物的TSS enrichment score值要整体偏高些。

1.2 双细胞的去除

双细胞预测几乎是所有单细胞测序技术都得考虑的一个问题，从原理来说，我们每个barcode就是一个细胞，但是因为所有的实验技术都不是100%完美的，因此往往会有一个barcode所包裹的油滴进来2个细胞。

对于10x数据来说，即使在使用标准试剂盒时，也可能有超过5%的细胞属于双细胞，这对聚类产生了重大影响。特别是在发育/轨迹分析中十分受影响，因为doublets看起来像是两种细胞类型的混合物，这可能与中间细胞类型或细胞状态混淆。

为了预测哪些“细胞”实际上是双细胞的，ArchR会从我们真实的数据中随机模拟产生混合的“双细胞”数据，这些“双细胞”数据与我们所有细胞一起做降维并UMAP可视化（"双细胞"会投影到UMAP中，并识别它们邻近的细胞），在这个过程中，ArchR会计算每个细胞的Doublet Enrichment，值越大，表示该细胞是双细胞的可能性越大。                                                                     

图2：ArchR包双细胞预测原理

02

降维分群

与单细胞RNA(scRNA-seq)相比，scATAC-seq数据由于其高维度和稀疏性而更具计算分析挑战性。主要体现在标准化和降维，这两大步骤跟单细胞转录组分析所用的统计学原理完全不同，以下为归纳总结的具体内容，如下表所示：

备注：TF-IDF & LSI都是自然语言常用的统计学方法。

2.1 peaks标准化

获得peak matrix后，跟基因类似，我们必须对其标准化。因为单细胞ATAC测的是DNA序列，对于二倍体物种来说，同一个位置最多有2套DNA序列，这便是单细胞ATAC peak matrix稀疏性的最大根源（单细胞转录组因测的是RNA，高表达的基因往往有多个转录分子）。因此，从数据实际情况出发，单细胞ATAC采取的是log(TF-IDF)( Term frequency-inverse document frequency) 标准化，简称文档频率法。

TF-IDF：是一种用于信息检索与数据挖掘的常用加权技术。TF是词频(Term Frequency)，IDF是逆文本频率指数(Inverse Document Frequency)。单个词汇在一篇文章中出现的次数越多，越重要。但是在语料库多次出现，重要性越来越低。IDF : 计算A term 出现稀少度。越稀少，越重要。

2.2 peaks降维

所有高维数据的分析都是采取降维的方式从多维到低纬的策略，之后还可以再次降维成2个维度并可视化（比如TSNE和UMAP）。我们对peaks是采取LSI降维的方式。

LSI：潜在语义索引(Latent Semantic Indexing,以下简称LSI)，有的文章也叫Latent Semantic  Analysis（LSA）。其实是一个东西，后面我们统称LSI，它是一种简单实用的主题模型。LSI是基于奇异值分解（SVD）的方法来得到文本的主题的。

2.3 细胞分群

与单细胞转录组类似，降维后的单细胞ATAC数据也同样可以采取graph-based clustering的分群方法。Graph-based图聚类算法包括两步：首先用降维（PCA或者LSI）的数据构建一个细胞间的k近邻稀疏矩阵，即将一个细胞与其欧式距离上最近的k个细胞聚为一类，然后在此基础上用Louvain算法进行模块优化(Blondel, Guillaume, Lambiotte, & Lefebvre, 2008)，旨在找到图中高度连接的模块。最后通过层次聚类将位于同一区域内没有差异表达基因(B-H adjusted p-value 低于0.05)的cluster进一步融合，重复该过程直到没有clusters可以合并。

备注：Signac和ArchR都是直接调用Seurat包的FindClusters()函数用不同分辨率来分群的。

03

marker基因的鉴定

细胞分群后，我们需要知道每个cluster属于什么细胞类型，也就是细胞命名。我们知道，单细胞转录组主要是依据每个cluster的marker基因来判断细胞类型的。那么对于单细胞ATAC，是不是也可以定义出每个cluster的特异高表达的基因集呢？

答案是肯定的，一般来说，我们是通过基因body区域加上一定范围内的上下游区域的整体ATAC信号来计算每个细胞每个基因的genescore值。

3.1 Genescore的计算

1）Signac是通过GeneActivity()函数https://satijalab.org/signac/reference/geneactivit来实现的，默认参数是基因上游2kb到TES区域。

2）而ArchR是通过addGeneScoreMatrix()函数https://www.archrproject.com/reference/addGeneScoreMatrix.html来实现的（createArrowFiles函数也会用默认参数得到genescore matrix矩阵），注意其计算原理稍微复杂，ArchR考虑到远端调控元件对基因活性的影响，因此默认的upstream和downstream范围更广。

在ArchR作者的发表文章中，他们测试了50多个不同的基因评分模型，并确定了一类在各种测试条件下表现始终优于其他模型的模型。这个模型类，在ArchR中作为默认实现，有三个主要组件:

a.整个基因体内的可及性有助于基因得分。

b.一种指数加权函数，以一种距离依赖的方式来解释假定的远端调节元件的活动。

c.施加基因边界，使不相关的调控元素对基因得分的贡献最小化。

图3：ArchR包Gene Scores计算原理

3.2 marker基因的可视化

marker 基因的ATAC信号（genescore值）同样可以在umap展示，也可以用小提琴图（VlnPlot），点状图（DotPlot）展示。与单细胞转录组相比，单细胞ATAC还多了基因区域的track的可视化展示。

1）以下为Signac包里单细胞ATAC marker 基因的ATAC信号（genescore值）结果展示图，Signac包与Seurat包一样，都是satijalab实验室团队开发的，因此该包继承了很多我们熟悉的Seurat包的方法函数。以下为例图举例。

 图4：Signac包单细胞ATAC可视化展示例图

2）ArchR包同样也可以做很多可视化的图。左边是CD14基因genescore值umap可视化展示，右边是track可视化图。

图5：ArchR包单细胞ATAC可视化展示例图

单细胞ATAC的亚群分析介绍就到这里，下一篇会给大家介绍单细胞ATAC的高级分析内容，比如motifdeviation、 拟时间分析、 单细胞RNA与单细胞ATAC的整合分析等。

本分享更多是从知识点和分析原理来讲解和归纳总结，具体实现方法和流程脚本可以查看下面参考资料软件的官方文档，里面都写得都很详细清楚。

如果大家对单细胞ATAC分析需求很大（欢迎读者朋友们多多点赞、在看、收藏和分享呀），后期给大家专门安排一期ArchR软件流程实操的经验分享。

04

参考资料

1.https://support.10xgenomics.com/single-cell-atac/software/pipelines/latest/what-is-cell-ranger-atac  # cell-ranger-atac

2.https://satijalab.org/signac/articles/pbmc_vignette.html #Signac官方教程

3.https://github.com/GreenleafLab/ArchR/ #ArchR github网站

4.https://www.archrproject.com/bookdown/creating-arrow-files.html #ArchR官方教程

5.Granja, J. M., et al. (2021). "ArchR is a scalable software package for integrative single-cell chromatin accessibility analysis." Nat Genet 53(3): 403-411.https://www.nature.com/articles/s41588-021-00790-6

6.Single-cell multiomic analysis identifies regulatory programs in mixed-phenotype acute leukemia. Nature Biotechnology (Granja JM, Klemm SK, McGinnis LM*, et al. 2019)https://www.nature.com/articles/s41587-019-0332-7 #人scATAC文章

7.A cis-regulatory atlas in maize at single-cell resolution. https://www.cell.com/cell/fulltext/S0092-8674(21)00493-1#articleInformation #玉米scATAC文章。5月7日发表在Cell上。

作者：Jenny

审稿：童蒙

编辑：angelica

往期精彩