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

数据科学工厂

发布于 2025-04-09 06:04:20

11100

代码可运行

文章被收录于专栏：数据科学（冷冻工厂）数据科学（冷冻工厂）

运行总次数：0

代码可运行

引言

本系列开启 R 中scRNA-seq数据分析教程[1]

不拆分数据重建无分支轨迹的伪时间

之前讲解了如何通过扩散图（diffusion map）为一条没有分叉的细胞轨迹重建伪时间（pseudotimes）。但如果系统很复杂，包含多种不同的终点状态，还有一些初始或中间状态被多个终点状态的轨迹共同拥有，那么要找出代表单一无分支轨迹的细胞就变得相当困难，这本身就成了一个需要解决的难题。

有意思的是，许多解决这类问题的工具会基于细胞邻居图进行随机游走（random walk），而且强烈建议把伪时间（pseudotimes）作为输入之一。这样一来，伪时间分析和轨迹分析就形成了一个相互依赖的循环。为了解决这个问题，有几种思路。比如，可以借助除基于转录组相似性的伪时间之外的其他方法来揭示细胞状态的变化，其中最常见的一种是将在下一节介绍的 RNA 速度分析（RNA velocity analysis）。另外，也可以尝试不把数据拆分成单一无分支轨迹，直接进行伪时间分析。

比如，可以直接把之前提到的扩散图（diffusion map）和扩散伪时间（DPT）应用到完整的 DS1 数据上，而不是只分析皮质部分。

library(Seurat)
library(destiny)

seurat_DS1 <- readRDS("DS1/seurat_obj_all.rds")
dm <- DiffusionMap(Embeddings(seurat_DS1, "pca")[,1:20])
dpt <- DPT(dm)
seurat_DS1$dpt <- rank(dpt$dpt)

FeaturePlot(seurat_DS1, c("dpt","SOX2","NHLH1","DCX"), ncol=4)

看起来还不错！不过有个显而易见的问题，就是重建的伪时间弄反了， SOX2+ 前体细胞的伪时间居然比 DCX+ 神经元还高。解决这个问题的简单办法当然就是手动把它调过来。

seurat_DS1$dpt <- max(seurat_DS1$dpt) - seurat_DS1$dpt
FeaturePlot(seurat_DS1, c("dpt","SOX2","NHLH1","DCX"), ncol=4)

或者，也可以直接指定随机游走的起点“末端”细胞。随便挑三个前体细胞作为末端，然后再跑一次 DPT。

set.seed(12345)
idx <- sample(which(seurat@active.ident %in% c('Dorsal telen. NPC',
                                               'G2M dorsal telen. NPC',
                                               'Dien. and midbrain NPC',
                                               'G2M Dien. and midbrain NPC')),3)
dpt2 <- DPT(dm, tips=idx)
seurat_DS1$dpt2 <- rank(dpt2$dpt)

FeaturePlot(seurat_DS1, c("dpt","dpt2"), ncol=2)

这种方法有个问题，因为 NPC 簇实际上涵盖了从前体细胞到神经元分化过程中各种不同的细胞状态，所以随机选出来的末端细胞不一定是轨迹真正的起点，很可能只是中间某个位置。为了解决这个麻烦，可以先试着在 destiny 包里用 random_root 函数，这个函数会先随便挑一个细胞，然后找出离它 DPT 值最大的那个细胞，作为末端候选者（默认的 DPT 函数也是这么干的）。之后，再从这些候选者里挑出被标注为前体细胞的那些，把它们的交集作为 DPT 的末端细胞。

tips_cand <- sapply(1:100, function(i){ random_root(dm) })
idx_NPC <- which(seurat@active.ident %in% c('Dorsal telen. NPC',
                                            'G2M dorsal telen. NPC',
                                            'Dien. and midbrain NPC',
                                            'G2M Dien. and midbrain NPC'))
tips_cand <- as.numeric(names(which.max(table(tips_cand[tips_cand %in% idx_NPC]))))
dpt3 <- DPT(dm, tips=tips_cand)
seurat_DS1$dpt3 <- rank(dpt3$dpt)

FeaturePlot(seurat_DS1, c("dpt","dpt2", "dpt3"), ncol=3)