R中单细胞RNA-seq分析教程 (14)

引言

本系列开启 R 中scRNA-seq数据分析教程，持续更新，欢迎关注，转发！想要获取更多教程内容或者生信分析服务可以添加文末的学习交流群或客服QQ:941844452。

高级分析选择

前面提到的分析大多是关于 scRNA-seq 数据的预处理，比如 normalization（标准化）、dimension reduction（降维）和 data integration（数据整合），以及最基础的分析，例如 clustering（聚类）和 marker identification（标记识别）。根据 scRNA-seq 数据反映的具体生物系统，还可以进一步运用一些分析方法，来挖掘更深层次的生物学信息。这些方法包括但不限于 pseudotime analysis（伪时间分析，前面已提到）、不同条件下基因表达差异的 differential expression analysis、RNA velocity analysis（RNA速度分析）、branching point analysis（分支点分析）、基于 ligand-receptor pairing（配体-受体配对）的 cell-cell communication analysis（细胞间通信分析），还有 gene regulatory network inferences（基因调控网络推断）。在下面的内容中，会简单介绍其中一些针对 scRNA-seq 数据的高级分析方法。

用 PAGA 分析群间的联系

之前讲到的 clustering analysis（聚类分析）是整理数据中细胞类型或状态多样性的常用手段。不过，基础的聚类分析无法告诉各个细胞集群之间是如何关联起来的。而这对于理解像发育或再生这样的动态过程至关重要。为此，有几种方法可以弥补这个不足。一种方法是 branching analysis（分支分析），它不直接划分集群，而是把细胞类型或状态的变化描绘成一棵树：以最初的细胞状态为树根，每当两种不同的细胞类型或状态从同一个祖先状态分化出来时，就标记一个 branching point（分支点）。目前有不少工具可以用来做这种分支分析，最受欢迎的有 monocle/monocle2/monocle3、URD 和 Slingshot。实际上，分支分析通常会和 pseudotime reconstruction（伪时间重构）一起进行，作为 trajectory analysis（轨迹分析）的一部分。所以，这些分支分析工具大多内置了计算 pseudotimes（伪时间）的功能；而许多专门用来生成伪时间的工具，比如前面提到的 diffusion map（扩散图），也能识别分支点。

除了分支分析，还有一种办法是利用已有的聚类结果，通过统计方法评估每两个集群之间的联系强度。通常，这种联系是看一个集群里的细胞有多大可能成为另一个集群里细胞的“最近邻”。做这类分析最常用的工具是 PAGA。

library(Seurat)
library(Matrix)
seurat_DS1 <- readRDS("DS1/seurat_obj_all.rds")

因为 scanpy 是 Python 的一个工具包，无法直接处理 Seurat 对象，所以得把这些信息转换成 scanpy 能识别的格式，比如 h5ad 或 loom。下面是用 loomR 包生成一个 loom 文件的例子，里面包含了必要的数据，像 PCA、UMAP 以及细胞类型的标注，方便后续分析。

library(loomR)
cell_attrs <- list(pca = Embeddings(seurat_DS1,"pca")[,1:20],
                   umap = Embeddings(seurat_DS1,"umap"),
                   celltype = [e-mail])
loom <- loomR::create("DS1/loom_obj.loom",
                      data = seurat_DS1[['RNA']]@data,
                      layers = list(counts = seurat[['RNA']]@counts),
                      cell.attrs = cell_attrs)
loom$close_all()

这些步骤会生成一个新文件（DS1/loom_obj.loom），之后就可以直接用它来运行 scanpy。
另外，也可以用 anndata 包生成一个 h5ad 文件，里面包含类似的数据。现在这种方式其实更受欢迎，因为 anndata 已经被很多 Python 单细胞 RNA 测序分析工具（比如 scanpy、scvelo 和 cellrank）认可为默认格式。

library(anndata)
adata <- AnnData(X = t(seurat_DS1[['RNA']]@data),
                 obs = data.frame(celltype = [e-mail], row.names = colnames(seurat_DS1)),
                 var = seurat_DS1[['RNA']]@meta.features,
                 layers = list(counts = t(seurat_DS1[['RNA']]@counts)),
                 obsm = list(pca = Embeddings(seurat_DS1,"pca")[,1:20],
                             umap = Embeddings(seurat_DS1,"umap"))
                )
adata$write_h5ad("DS1/anndata_obj.h5ad")

这些步骤会生成一个新文件 DS1/anndata_obj.h5ad，同样可以用它来运行 scanpy。

接下来，你可以选择转到 Python（版本至少 3.6）环境并装好 scanpy 来继续操作。如果不想离开 R，就得用一个叫 reticulate 的 R 包。这个包是 RStudio 开发的（详情见 https://rstudio.github.io/reticulate/），能让 R 和 Python 打通，方便在 R 里直接跑 Python 代码。一般来说，装 Seurat 的时候它就应该已经装好了，但你可以通过运行 install.packages("reticulate") 再确认一下。如果想要用最新开发版，可以用 remotes::install_github("rstudio/reticulate") 来安装。

不过得注意，reticulate 只是个桥梁，不是 Python 本身，所以你还是得装一个 Python（版本 3.6 或以上）让它调用。之后，在 R 里加载 reticulate 包，装好 scanpy，然后把它引入 R 环境用。

library(reticulate)
py_install("scanpy", pip=T)
sc <- import("scanpy")

如果没报错，那就说明 scanpy 已经在你的 Python 环境里装好了，也顺利导入了。如前面提到的，scanpy 是个很全面的 scRNA-seq 数据分析工具包，跟 Seurat 差不多，但它是基于 Python 开发的，所以功能特别丰富，不仅仅限于 PAGA。这次只是用它来做 PAGA 分析，所需的输入数据还是靠 Seurat 的分析结果提供的。

adata_DS1 <- sc$read_loom("DS1/loom_obj.loom") # option1: load the loom file
adata_DS1 <- sc$read("DS1/anndata_obj.h5ad") # option2: load the h5ad file

sc$pp$neighbors(adata_DS1, n_neighbors=20L, use_rep='pca')
sc$tl$paga(adata_DS1, groups='celltype')
adata_DS1$write_h5ad("DS1/anndata_obj.h5ad")

PAGA 的核心思路是先用给定的降维矩阵（比如例子里的 'pca'）构建一个细胞间的 k 近邻网络，再数一数连接两个细胞簇之间的边有多少条。然后，把这个数字跟随机网络里预期的边数比一比，就能估算出两个簇之间的 PAGA 连通性。这个连通性是个 0 到 1 的值，0 表示完全没联系，1 表示联系很紧密。所以不难看出，除了聚类结果本身，影响 PAGA 结果的关键因素还有：1) 用哪种降维方式来找细胞的邻居（也就是 neighbors 函数里的 use_rep 参数）；2) 每个细胞找几个邻居（n_neighbors 参数）。另外值得一提的是，neighbors 函数里还有些其他参数，比如 method 参数决定怎么衡量细胞间的连通性，metric 参数决定用什么方式算细胞间的距离，这些都可能会改变邻居网络，进而影响到 PAGA 的结果。

接下来，可以把 PAGA 估算出来的簇间连通性画出来看看。

plt <- import("matplotlib")
plt$use("Agg", force = TRUE)
sc$pl$paga(adata_DS1,
           color='celltype',
           fontsize=7,
           frameon=FALSE,
           save="DS1_paga.png")

虽然不算完美，但从簇级别的总结图上还是能看出，它在一定程度上反映了数据集整体的多样性。那些连通性最强、画得最粗的边，往往确实对应着细胞的分化过程，比如背侧端脑中间神经元（dorsal telen. IP）和背侧端脑神经元（dorsal telen. neurons）之间那种紧密的联系。

本文由mdnice多平台发布