单细胞ATAC高级分析

前两篇分享(单细胞ATAC概述和单细胞ATAC亚群分析)，我们分别给大家介绍了单细胞ATAC的基础知识点和亚群分析，这篇分享就是专门给大家带来一些高级分析的内容了。

主要包括以下四个方面：motifdeviation的计算、单细胞ATAC拟时间分析、单细胞RNA与ATAC的整合分析以及peak co-accessibility分析。

01motif deviation的计算

之前我们有提到过第一篇单细胞ATAC数据Greenleaf(Buenrostro, Wu et al. 2015)实验室的Nature文章，里面算是首次介绍了motif deviation score的概念。简单来说，它是通过一定的算法原理，利用公共TF motif数据库，给每个细胞的每个TF对应的motif进行一个打分值，值的高低代表该motif在该细胞的调控活跃程度。

以下是具体算法原理，感兴趣的读者可以移步去看这篇文章。

图1：motif deviation score的计算原理（Buenrostro, Wu et al. 2015）

02单细胞ATAC拟时间分析

以下内容主要以ArchR软件为例讲解，目前ArchR软件1.0.1最新版本给大家提供了3种分析策略，分别是ArchR自带的、引用monocle3和SlingShot包的。后面还给大家展示stream包的例子。

01

addMonocleTrajectory()方法

函数说明：

https://www.archrproject.com/reference/addMonocleTrajectory.html

注意这里的monocle是指monocle3（monocle3和monocle2的分析原理有较大区别）。目前monocle3官方教程里不仅提供了单细胞转录组的轨迹发育教程，也比较友好提供了单细胞ATAC的接口（大家也可以从原始peak matrix出发，直接用monocle3做拟时间分析）。

ArchR的addMonocleTrajectory函数实际上是调用monocle3包里的函数，便于大家用ArchR里降维的结果直接做拟时间分析。

接下来，让我们跟着monocle3的官方教程来了解下它是如何区别scATAC和scRNA数据来做拟时间分析的。首先，不知道大家是否还记得上篇”亚群分析“第2部分的降维分群里给大家归纳总结了二者的区别。为了加深印象，这个给大家再展示下。

Monocle3主要包括以下几个步骤来分析：preprocess_cds（数据预处理，标准化）、reduce_dimension（降维）、cluster_cells、learn_graph、order_cells、plot_cells。

结合上面表格对scATAC和scRNA的总结，我们知道二者在用monocle3做拟时间分析时主要是preprocess_cds、reduce_dimension函数的参数要注意区别：

以下是具体参数列举：

#scATAC
cds <- preprocess_cds(cds, method='LSI',norm_method='none', num_dim = 50) #ATAC, 需要先用log(TF-IDF)标准化scATAC peak matrix
cds <- reduce_dimension(cds, preprocess_method = "LSI",reduction_method = c("UMAP"))

#scRNA
cds <- preprocess_cds(cds, method='PCA',norm_method='log', num_dim = 50)
cds <- reduce_dimension(cds, preprocess_method = "PCA",reduction_method = c("UMAP"))

02

addTrajectory()方法

函数说明：

https://www.archrproject.com/reference/addTrajectory.html

为了在伪时间中对细胞进行排序，ArchR创建细胞轨迹，在ArchRProject中通过一个低n维子空间对细胞进行排序。以前，我们在二维UMAP子空间中执行过这种排序（例如monocle3），但是ArchR改进了这种方法，使其能够在n维子空间(即LSI)中对齐。

具体原理请移步官方教程，如下：

https://www.archrproject.com/bookdown/trajectory-analysis-with-archr.html

值得注意的是，addTrajectory()函数需要大家指定亚群的大体发育轨迹（至少3个点），例如trajectory =c('C2','C4','C5')，便是告诉ArchR这些群里的细胞是按照C2>C4>C5的发育轨迹的。

03

addSlingShotTrajectories()方法

函数说明：

https://www.archrproject.com/reference/addSlingShotTrajectories.html

提供用于推断低维数据中连续的分支沿路结构的函数。Slingshot的设计目的是在单细胞测序数据中模拟发育轨迹，并在降维和聚类后作为分析管道的一个组成部分。它足够灵活，可以处理任意多的分支事件，并允许通过监督图的构造将先验知识结合起来。

04

stream包

stream包(Chen, Albergante et al. 2019)同样也是一款可以做单细胞转录组和表观基因组数据发育轨迹分析的软件。它的特点是适合多个平行发育的分支。如下图，可以根据自己数据的背景选择。

图2：stream包scATAC发育轨迹结果图

03单细胞ATAC与RNA的整合分析

单细胞ATAC与RNA的整合分析主要目的是帮助scATAC亚群的命名，依据scRNA所有基因表达matrix和scATAC所有基因区域ATAC信号打分genescore matrix的相似性而整合预测。

01

利用Signac整合

Signac作为Seurat的“姐妹”包较好给大家提供了单细胞ATAC与RNA的整合分析，用的是FindTransferAnchors函数。

https://satijalab.org/signac/articles/pbmc_vignette.html#integrating-with-scrna-seq-data-1

图3：Signac包的单细胞ATAC和RNA的整合结果

02

利用ArchR包整合

ArchR包整合的工作方式是，通过比较scATAC-seq基因得分矩阵和scRNA-seq基因表达矩阵，将来自scATAC-seq的细胞与来自scRNA-seq的细胞直接对齐。在底层，这种对齐也是使用Seurat包中的FindTransferAnchors函数执行的，该函数允许跨两个数据集对齐数据。然而，为了将这个过程适当地扩展到数十万个单元，ArchR通过将整个单元划分为更小的单元组并执行单独的对齐，提供了这个过程的并行化。

图4：ArchR包利用单细胞RNA数据预测ATAC细胞类型结果

04peak co-accessibility分析

识别peak co-accessibility可以预测调控元件之间的相互作用，比如enhancer和promoter的loop。这个是利用两两peak（regulator）在多个细胞之间的相关性来预测的，利用相应的阈值（比如SCC 0.5）来筛选定义。                       

图5：ArchR包peak co-accessibility分析图

单细胞ATAC的三篇分享（基础知识点、亚群分析、高级分析）就到这里啦。有机会的话，小编后面可能还会给大家带来更多关于单细胞ATAC的分享内容。

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

假设大家对单细胞ATAC分析需求很大（欢迎读者朋友们多多点击“在看”、收藏和转给需要的朋友呀），后期可能会专门给大家带来一篇ArchR软件流程实操的经验分享。

05参考资料

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

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

3. Buenrostro JD, Wu B, Litzenburger UM, Ruff D, Gonzales ML, Snyder MP,et al. Single-cell chromatin accessibility reveals principles of regulatory variation. Nature. 2015;523:486–90. #Greenleaf(Buenrostro, Wu et al. 2015)实验室

4. 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

5. 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文章

6. 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上。

7. https://cole-trapnell-lab.github.io/monocle3/docs/introduction/ #monocle3官方教程

8. https://github.com/kstreet13/slingshot  #slingshot发育轨迹分析

9. https://nbviewer.jupyter.org/github/pinellolab/STREAM/blob/master/tutorial/2.1-STREAM_scATAC-seq_peaks.ipynb?flush_cache=true #利用stream做scATAC数据的发育轨迹分析

10. Chen, H., et al., Single-cell trajectories reconstruction, exploration and mapping of omics data with STREAM. Nat Commun, 2019. 10(1): p. 1903. #stream软件发表文章

作者：Jenny

审稿：童蒙

编辑：angelica

往 期 精 彩 看 过 来