如题所述
上一期，跟大家简单介绍了 下单细胞ATAC的背景知识点及其10x ATAC基础数据的获取方式 。接下来就带大家从fragment.csv、singlecell.csv、peaks matrix等数据出发，做单细胞ATAC的亚群分析。
与单细胞转录组类似，单细胞ATAC的分析流程也主要包括细胞质控、peaks标准化及其降维分群、marker基因的鉴定等几个步骤。常用的单细胞ATAC分析流程软件包含 cell-ranger-atac、Signac和ArchR等。
单细胞ATAC的质控点一般包含以下几个方面：样本重复（biological replicates），bulkATAC vs scATAC的相关性、fragment length distribution、per nucleus read-depth、transcription start site (TSS) enrichment、双细胞比例等。
前面提到的样本相关性和fragments的长度分布主要是从整体水平上检查我们的单个样本数据的可靠性。
而要去掉不符合质控的细胞，我们主要从 fragments 数目 和 TSS enrichment score 这两点出发。
👉 fragments 数目：一般指单个细胞（barcode）所属的total fragments数目。这个不同的软件具体的定义不同，比如cell-ranger-atac和Signac指peaks所属区域的fragments 数目，其中singlecell.csv文件中peak_region_fragments列便是指fragments 数目，而ArchR是指全基因组所有的fragments 数（这个跟该软件的分析策略有关，后面会提到）。
备注：fragments 数目&TSS enrichment score的阈值不仅与所用软件具体的计算公式有关（不同的软件具体的参数可能不同），也与自己数据的实际情况有关。比如哺乳动物和植物的单细胞ATAC数据TSS enrichment score就不能用相同的指标cutoff来衡量，一般来说哺乳动物的TSS enrichment score值要整体偏高些。
双细胞预测几乎是所有单细胞测序技术都得考虑的一个问题，从原理来说，我们每个barcode就是一个细胞，但是因为所有的实验技术都不是100%完美的，因此往往会有一个barcode所包裹的油滴进来2个细胞。
对于10x数据来说，即使在使用标准试剂盒时，也可能有超过5%的细胞属于双细胞，这对聚类产生了重大影响。特别是在发育/轨迹分析中十分受影响，因为doublets看起来像是两种细胞类型的混合物，这可能与中间细胞类型或细胞状态混淆。
为了预测哪些“细胞”实际上是双细胞的，ArchR会从我们真实的数据中随机模拟产生混合的“双细胞”数据，这些“双细胞”数据与我们所有细胞一起做降维并UMAP可视化（"双细胞"会投影到UMAP中，并识别它们邻近的细胞），在这个过程中，ArchR会计算每个细胞的Doublet Enrichment，值越大，表示该细胞是双细胞的可能性越大。
与单细胞RNA(scRNA-seq)相比，scATAC-seq数据由于其高维度和稀疏性而更具计算分析挑战性。主要体现在标准化和降维，这两大步骤跟单细胞转录组分析所用的统计学原理完全不同，以下为归纳总结的具体内容，如下表所示：
备注：TF-IDF & LSI都是自然语言常用的统计学方法。
获得peak matrix后，跟基因类似，我们必须对其标准化。因为单细胞ATAC测的是DNA序列，对于二倍体物种来说，同一个位置最多有2套DNA序列，这便是单细胞ATAC peak matrix稀疏性的最大根源（单细胞转录组因测的是RNA，高表达的基因往往有多个转录分子）。因此，从数据实际情况出发，单细胞ATAC采取的是log(TF-IDF)( Term frequency-inverse document frequency) 标准化，简称文档频率法。
所有高维数据的分析都是采取降维的方式从多维到低纬的策略，之后还可以再次降维成2个维度并可视化（比如TSNE和UMAP）。我们对peaks是采取LSI降维的方式。
与单细胞转录组类似，降维后的单细胞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()函数用不同分辨率来分群的。
细胞分群后，我们需要知道每个cluster属于什么细胞类型，也就是细胞命名。我们知道，单细胞转录组主要是依据每个cluster的marker基因来判断细胞类型的。那么 对于单细胞ATAC，是不是也可以定义出每个cluster的特异高表达的基因集呢？
答案是肯定的，一般来说，我们是通过基因body区域加上一定范围内的上下游区域的整体ATAC信号来计算每个细胞每个基因的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中作为默认实现，有三个主要组件:
marker 基因的ATAC信号（genescore值）同样可以在umap展示，也可以用小提琴图（VlnPlot），点状图（DotPlot）展示。与单细胞转录组相比，单细胞ATAC还多了基因区域的track的可视化展示。
单细胞ATAC的亚群分析介绍就到这里，下一篇会给大家介绍单细胞ATAC的高级分析内容，比如motifdeviation、 拟时间分析、 单细胞RNA与单细胞ATAC的整合分析等。
本分享更多是从知识点和分析原理来讲解和归纳总结，具体实现方法和流程脚本可以查看下面参考资料软件的官方文档，里面都写得都很详细清楚。
