单细胞ATAC亚群分析

Posted 2023-04-03

参考技术A

上一期，跟大家简单介绍了 下单细胞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的整合分析等。
本分享更多是从知识点和分析原理来讲解和归纳总结，具体实现方法和流程脚本可以查看下面参考资料软件的官方文档，里面都写得都很详细清楚。

斯坦福Greenleaf等开发工具研究CRISPR筛选单细胞染色质开放性变化

能够在单细胞分析染色质开放性的同时分析细胞内sgRNA序列，命名为Spear-ATAC(Single-cell perturbations with an accessibility read-out using scATAC-seq)，并且用此项技术研究了肿瘤细胞系关键转录因子受到干扰后整体染色质开放性和转录因子结合模式相关变化。斯坦福William Greenleaf此项工作刚刚发布在bioRxiv。

原文链接：

https://www.biorxiv.org/content/10.1101/2020.11.02.364265v1