admin管理员组

文章数量:1029193

R中单细胞RNA

引言

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

简介

RNA 速度分析最早由 La Manno 等人在 2018 年提出，他们分别在瑞典 Karolinska 研究所的 Sten Linnarsson 实验室和美国 Harvard Medical School 的 Peter Kharchenko 实验室工作。这是一种基于简单转录动态模型的分析方法。这个模型认为，一个细胞中某个基因的转录本数量是由转录、剪接和降解这几个过程共同决定的。成熟的 RNAs 反映了细胞当前的状态。如果 RNA 的降解速度与转录和剪接的速度达到平衡，细胞状态可能会保持不变；否则，它就会随时间发生变化。假设细胞状态并非稳定的，转录和 RNA 处理（比如剪接）之间的时间差就能帮助预测细胞状态未来的变化趋势，前提是知道不同转录本的降解速度。基于这一想法，La Manno 团队开发了首个算法，用外显子读取数据代表成熟的 RNA 转录本，用内含子读取数据代表还未剪接的未成熟 RNA 转录本。

在此基础上，Bergen 等人在德国 Helmholtz Center Munich 的 Alexander Wolf 实验室和 Fabian Theis 实验室改进了这个转录动态模型的估计方法，让它不再需要假设细胞状态是稳定的。他们还开发了一个 Python 工具包 scvelo，不仅功能更全面，而且计算效率更高。

接下来，将用 DS1 数据作为例子，展示如何在 R 中通过 scvelo 包进行 RNA 速度分析，方法上类似于前面提到的 PAGA。

RNA 速度分析

RNA 速度分析需要分别准备细胞的外显子和内含子计数矩阵。但遗憾的是，目前最常用的 10x Genomics scRNA-seq 平台的标准预处理工具 CellRanger 只统计外显子区域的转录本数量，所以需要额外步骤来生成所需的矩阵。通常有两种办法：

  1. 利用 CellRanger 的映射结果（BAM 格式）和基因注释信息，通过转录计数工具（比如 dropEst）生成计数矩阵。
  2. 通过伪映射技术，用工具如 kallisto 生成外显子和内含子计数矩阵。

在这个例子中，DS1 的计数矩阵是用 dropEst 生成的。

首先，需要在 R 中加载经过预处理的 DS1 Seurat 对象，以及对应的外显子和内含子计数矩阵。dropEst 生成的计数矩阵包含了所有检测到的细胞条形码，所以需要从中筛选出 Seurat 对象里记录的那些细胞。

代码语言：javascript代码运行次数：0运行复制
library(Seurat)
library(Matrix)

seurat_DS1 <- readRDS("DS1/seurat_obj_all.rds")
mats <- readRDS("DS1/mats_dropest.rds")
mats <- lapply(mats, function(mat) mat[,colnames(seurat_DS1)])

接下来，得准备一个 loom 或 h5ad 文件，里面包含运行 scvelo 所需的各种信息。这跟前面为 PAGA 准备文件的步骤很像，不过得用上已剪接和未剪接的数据矩阵。

下面是创建 loom 文件的具体方法：

代码语言：javascript代码运行次数：0运行复制
library(loomR)
cell_attrs <- list(pca = Embeddings(seurat_DS1,"pca")[,1:20],
                   umap = Embeddings(seurat_DS1,"umap"),
                   celltype = seurat_DS1@active.ident)
shared_genes <- intersect(rownames(mats$exon),rownames(mats$intron))
loom <- loomR::create("DS1/loom_obj_scvelo.loom",
                      data = mats$exon[shared_genes,],
                      layers = list(spliced = mats$exon[shared_genes,],
                                    unspliced = mats$intron[shared_genes,]),
                      cell.attrs = cell_attrs)
loom$close_all()

以下是创建H5AD文件的方法：

代码语言：javascript代码运行次数：0运行复制
library(anndata)
shared_genes <- intersect(rownames(mats$exon),rownames(mats$intron))
adata <- AnnData(X = t(mats$exon[shared_genes,]),
                 obs = data.frame(seurat_DS1@meta.data, celltype=seurat_DS1@active.ident),
                 var = NULL,
                 layers = list(spliced = t(mats$exon[shared_genes,]),
                               unspliced = t(mats$intron[shared_genes,])),
                 obsm = list(pca = Embeddings(seurat_DS1,"pca")[,1:20],
                             umap = Embeddings(seurat_DS1,"umap"))
                )
adata$write_h5ad("DS1/anndata_obj_scvelo.h5ad")

接下来，跟做 PAGA 分析时差不多，可以转到安装了 scvelo 的 Python（版本 3.6 或更高）环境来进行后面的操作，或者直接用 R 的 reticulate 包。

然后，在 R 中加载 reticulate 包，安装 scvelo 包，再把它引入 R 环境里使用。

代码语言：javascript代码运行次数：0运行复制
library(reticulate)
py_install("scvelo", pip=T)
scvelo <- import("scvelo")

如果没出现任何错误提示，那就说明 scvelo 已经在你的 Python 里装好了，也顺利导入了。接下来，就可以开始进行 RNA 速度分析了。

代码语言：javascript代码运行次数：0运行复制
adata_DS1 <- scvelo$read_loom("DS1/loom_obj_scvelo.loom") # option1: load the loom file
adata_DS1 <- scvelo$read("DS1/anndata_obj_scvelo.h5ad") # option2: load the h5ad file

scvelo$pp$filter_and_normalize(adata_DS1,
                               min_shared_counts=as.integer(10),
                               n_top_genes=as.integer(3000))
scvelo$pp$moments(adata_DS1,
                  n_neighbors = as.integer(30),
                  use_rep = "pca")
scvelo$tl$velocity(adata_DS1)
scvelo$tl$velocity_graph(adata_DS1)

在 scvelo 里，你可以根据需要调整一些参数。算出速度之后，还能用 scvelo 做更多分析，比如速度伪时间估计，不过这些内容就不在这篇教程里细讲了。这一部分的最后，来把速度估计的结果可视化展示一下。

代码语言：javascript代码运行次数：0运行复制
plt <- import("matplotlib")
plt$use("Agg", force = TRUE)
scvelo$pl$velocity_embedding_stream(adata_DS1,
                                    basis="umap",
                                    color="celltype",
                                    dpi=120,
                                    figsize = c(8,8),
                                    save="DS1_scvelo_stream.png")

这些箭头展示了估算出的细胞状态变化方向，你能清楚地看到从 NPC 到神经元的转变，这反映了神经细胞分化和成熟的过程。

另外，scvelo 还能利用估算出的 RNA 速度来生成伪时间。它支持两种速度伪时间计算方式。一种叫“速度伪时间”，计算方法跟前面提到的扩散伪时间有点像。不同之处在于，扩散伪时间用的是基于扩散图嵌入算出来的对称转移矩阵来模拟随机游走，而速度伪时间用的是 RNA 速度估算出的带有方向性的转移矩阵。还有一种叫“潜在伪时间”，完全依靠速度的动态变化来生成。

代码语言：javascript代码运行次数：0运行复制
scvelo$tl$velocity_pseudotime(adata_DS1)
scvelo$tl$recover_dynamics(adata_DS1)
scvelo$tl$recover_latent_time(adata_DS1)

接下来，可以把算出来的伪时间提取出来，投影到 UMAP 上做成特征图，这样就能直观地看到结果并进行比较。

代码语言：javascript代码运行次数：0运行复制
seurat_DS1$velocity_pseudotime <- adata_DS1[['obs']]$velocity_pseudotime
seurat_DS1$latent_time <- adata_DS1[['obs']]$latent_time

FeaturePlot(seurat_DS1,
            c("velocity_pseudotime", "latent_time")) & NoAxes()

从这里可以明显看出，速度伪时间比潜在时间更能准确地反映从 NPC 到神经元的分化路径。这种情况挺常见的，但也不是绝对的。有些数据集上，潜在时间反而表现得更好。所以得具体检查和比较一下。

最后，再把处理好的 AnnData 和包含速度伪时间信息的 Seurat 对象保存起来。

代码语言：javascript代码运行次数：0运行复制
adata_DS1$write_h5ad('DS1/anndata_obj_scvelo.h5ad')
saveRDS(seurat_DS1, file='DS1/seurat_obj_all.rds')

Reference

[1] Source:

本文参与 腾讯云自媒体同步曝光计划，分享自微信公众号。原始发表：2025-04-21，如有侵权请联系 cloudcommunity@tencent 删除模型数据对象工具教程

R中单细胞RNA

引言

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

简介

RNA 速度分析最早由 La Manno 等人在 2018 年提出，他们分别在瑞典 Karolinska 研究所的 Sten Linnarsson 实验室和美国 Harvard Medical School 的 Peter Kharchenko 实验室工作。这是一种基于简单转录动态模型的分析方法。这个模型认为，一个细胞中某个基因的转录本数量是由转录、剪接和降解这几个过程共同决定的。成熟的 RNAs 反映了细胞当前的状态。如果 RNA 的降解速度与转录和剪接的速度达到平衡，细胞状态可能会保持不变；否则，它就会随时间发生变化。假设细胞状态并非稳定的，转录和 RNA 处理（比如剪接）之间的时间差就能帮助预测细胞状态未来的变化趋势，前提是知道不同转录本的降解速度。基于这一想法，La Manno 团队开发了首个算法，用外显子读取数据代表成熟的 RNA 转录本，用内含子读取数据代表还未剪接的未成熟 RNA 转录本。

在此基础上，Bergen 等人在德国 Helmholtz Center Munich 的 Alexander Wolf 实验室和 Fabian Theis 实验室改进了这个转录动态模型的估计方法，让它不再需要假设细胞状态是稳定的。他们还开发了一个 Python 工具包 scvelo，不仅功能更全面，而且计算效率更高。

接下来，将用 DS1 数据作为例子，展示如何在 R 中通过 scvelo 包进行 RNA 速度分析，方法上类似于前面提到的 PAGA。

RNA 速度分析

RNA 速度分析需要分别准备细胞的外显子和内含子计数矩阵。但遗憾的是，目前最常用的 10x Genomics scRNA-seq 平台的标准预处理工具 CellRanger 只统计外显子区域的转录本数量，所以需要额外步骤来生成所需的矩阵。通常有两种办法：

  1. 利用 CellRanger 的映射结果（BAM 格式）和基因注释信息，通过转录计数工具（比如 dropEst）生成计数矩阵。
  2. 通过伪映射技术，用工具如 kallisto 生成外显子和内含子计数矩阵。

在这个例子中，DS1 的计数矩阵是用 dropEst 生成的。

首先，需要在 R 中加载经过预处理的 DS1 Seurat 对象，以及对应的外显子和内含子计数矩阵。dropEst 生成的计数矩阵包含了所有检测到的细胞条形码，所以需要从中筛选出 Seurat 对象里记录的那些细胞。

代码语言：javascript代码运行次数：0运行复制
library(Seurat)
library(Matrix)

seurat_DS1 <- readRDS("DS1/seurat_obj_all.rds")
mats <- readRDS("DS1/mats_dropest.rds")
mats <- lapply(mats, function(mat) mat[,colnames(seurat_DS1)])

接下来，得准备一个 loom 或 h5ad 文件，里面包含运行 scvelo 所需的各种信息。这跟前面为 PAGA 准备文件的步骤很像，不过得用上已剪接和未剪接的数据矩阵。

下面是创建 loom 文件的具体方法：

代码语言：javascript代码运行次数：0运行复制
library(loomR)
cell_attrs <- list(pca = Embeddings(seurat_DS1,"pca")[,1:20],
                   umap = Embeddings(seurat_DS1,"umap"),
                   celltype = seurat_DS1@active.ident)
shared_genes <- intersect(rownames(mats$exon),rownames(mats$intron))
loom <- loomR::create("DS1/loom_obj_scvelo.loom",
                      data = mats$exon[shared_genes,],
                      layers = list(spliced = mats$exon[shared_genes,],
                                    unspliced = mats$intron[shared_genes,]),
                      cell.attrs = cell_attrs)
loom$close_all()

以下是创建H5AD文件的方法：

代码语言：javascript代码运行次数：0运行复制
library(anndata)
shared_genes <- intersect(rownames(mats$exon),rownames(mats$intron))
adata <- AnnData(X = t(mats$exon[shared_genes,]),
                 obs = data.frame(seurat_DS1@meta.data, celltype=seurat_DS1@active.ident),
                 var = NULL,
                 layers = list(spliced = t(mats$exon[shared_genes,]),
                               unspliced = t(mats$intron[shared_genes,])),
                 obsm = list(pca = Embeddings(seurat_DS1,"pca")[,1:20],
                             umap = Embeddings(seurat_DS1,"umap"))
                )
adata$write_h5ad("DS1/anndata_obj_scvelo.h5ad")

接下来，跟做 PAGA 分析时差不多，可以转到安装了 scvelo 的 Python（版本 3.6 或更高）环境来进行后面的操作，或者直接用 R 的 reticulate 包。

然后，在 R 中加载 reticulate 包，安装 scvelo 包，再把它引入 R 环境里使用。

代码语言：javascript代码运行次数：0运行复制
library(reticulate)
py_install("scvelo", pip=T)
scvelo <- import("scvelo")

如果没出现任何错误提示，那就说明 scvelo 已经在你的 Python 里装好了，也顺利导入了。接下来，就可以开始进行 RNA 速度分析了。

代码语言：javascript代码运行次数：0运行复制
adata_DS1 <- scvelo$read_loom("DS1/loom_obj_scvelo.loom") # option1: load the loom file
adata_DS1 <- scvelo$read("DS1/anndata_obj_scvelo.h5ad") # option2: load the h5ad file

scvelo$pp$filter_and_normalize(adata_DS1,
                               min_shared_counts=as.integer(10),
                               n_top_genes=as.integer(3000))
scvelo$pp$moments(adata_DS1,
                  n_neighbors = as.integer(30),
                  use_rep = "pca")
scvelo$tl$velocity(adata_DS1)
scvelo$tl$velocity_graph(adata_DS1)

在 scvelo 里，你可以根据需要调整一些参数。算出速度之后，还能用 scvelo 做更多分析，比如速度伪时间估计，不过这些内容就不在这篇教程里细讲了。这一部分的最后，来把速度估计的结果可视化展示一下。

代码语言：javascript代码运行次数：0运行复制
plt <- import("matplotlib")
plt$use("Agg", force = TRUE)
scvelo$pl$velocity_embedding_stream(adata_DS1,
                                    basis="umap",
                                    color="celltype",
                                    dpi=120,
                                    figsize = c(8,8),
                                    save="DS1_scvelo_stream.png")

这些箭头展示了估算出的细胞状态变化方向，你能清楚地看到从 NPC 到神经元的转变，这反映了神经细胞分化和成熟的过程。

另外，scvelo 还能利用估算出的 RNA 速度来生成伪时间。它支持两种速度伪时间计算方式。一种叫“速度伪时间”，计算方法跟前面提到的扩散伪时间有点像。不同之处在于，扩散伪时间用的是基于扩散图嵌入算出来的对称转移矩阵来模拟随机游走，而速度伪时间用的是 RNA 速度估算出的带有方向性的转移矩阵。还有一种叫“潜在伪时间”，完全依靠速度的动态变化来生成。

代码语言：javascript代码运行次数：0运行复制
scvelo$tl$velocity_pseudotime(adata_DS1)
scvelo$tl$recover_dynamics(adata_DS1)
scvelo$tl$recover_latent_time(adata_DS1)

接下来，可以把算出来的伪时间提取出来，投影到 UMAP 上做成特征图，这样就能直观地看到结果并进行比较。

代码语言：javascript代码运行次数：0运行复制
seurat_DS1$velocity_pseudotime <- adata_DS1[['obs']]$velocity_pseudotime
seurat_DS1$latent_time <- adata_DS1[['obs']]$latent_time

FeaturePlot(seurat_DS1,
            c("velocity_pseudotime", "latent_time")) & NoAxes()

从这里可以明显看出，速度伪时间比潜在时间更能准确地反映从 NPC 到神经元的分化路径。这种情况挺常见的，但也不是绝对的。有些数据集上，潜在时间反而表现得更好。所以得具体检查和比较一下。

最后，再把处理好的 AnnData 和包含速度伪时间信息的 Seurat 对象保存起来。

代码语言：javascript代码运行次数：0运行复制
adata_DS1$write_h5ad('DS1/anndata_obj_scvelo.h5ad')
saveRDS(seurat_DS1, file='DS1/seurat_obj_all.rds')

Reference

[1] Source:

本文参与 腾讯云自媒体同步曝光计划，分享自微信公众号。原始发表：2025-04-21，如有侵权请联系 cloudcommunity@tencent 删除模型数据对象工具教程

本文标签： R中单细胞RNA

版权声明：本文标题：R中单细胞RNA 内容由热心网友自发贡献，该文观点仅代表作者本人， 转载请联系作者并注明出处：http://it.en369.cn/jiaocheng/1747576316a2180527.html， 本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌抄袭侵权/违法违规的内容，一经查实，本站将立刻删除。