单细胞交响乐
搜索文档…
⌃K
Links

3.12 细胞轨迹推断

刘小泽写于2020.7.16

1 前言

接触单细胞数据,相信你一定听过:轨迹推断,英文名词是: Trajectory Analysis。和它相关的另一个名词是:拟时序分析(pseudotime),指的是细胞沿着这个轨迹,并且对潜在的生物活动进行量化。注意这里看字面意思就知道,并不是指真正的时间,而是指细胞与细胞之间的更替、转化的顺序或者是轨迹,可以理解为“一个连续过程的缩影”。
不同的生物过程对应的“拟时序”也是不同的:
图片来自:https://www.youtube.com/watch?v=XmHDexCtjyw
许多生物过程都伴随着细胞状态的连续性变化,比如研究发育就会经常使用到。我们可以利用单细胞数据在高维空间画一条线,贯穿于多种细胞状态。最简单是点到点的一条路径,更复杂的还有一个点出发再生成多个分支。
自2014年以后,已经开发出了超50种方法,那么选择何种方法进行分析成为了一个难题,因为我们不可能每一种都试一下,但有评测文章发现:SlingshotTSCANMonocle DDRTree这几种方法都不错
当然还有其他的几种方法值得推荐,还做成了一个流程图方便查阅 如果对评测感兴趣,可以看看dynverse

资源整理

做这个分析之前,最好先问几个问题:

  • 确定数据会体现发育轨迹吗?也就是研究的样本是不是和发育相关的?
  • 数据中的细胞会体现出中间态吗?
  • 是否认为轨迹会出现分支?

并且要注意:

  • 任何数据都可以强行画出轨迹,但不一定都有生物学意义!
  • 先要保证目前找到的HVGs和降维结果符合我们的预期,才能继续向下分析

2 学习Slingshot

它需要两个必须的输入文件:降维结果与细胞分群结果
因为它分析的基础假设就是:在低维空间上,细胞的位置是连续的并且是一个接一个的

2.1 数据准备

means <- rbind(
# non-DE genes
matrix(rep(rep(c(0.1,0.5,1,2,3), each = 300),100),
ncol = 300, byrow = TRUE),
# early deactivation
matrix(rep(exp(atan( ((300:1)-200)/50 )),50), ncol = 300, byrow = TRUE),
# late deactivation
matrix(rep(exp(atan( ((300:1)-100)/50 )),50), ncol = 300, byrow = TRUE),
# early activation
matrix(rep(exp(atan( ((1:300)-100)/50 )),50), ncol = 300, byrow = TRUE),
# late activation
matrix(rep(exp(atan( ((1:300)-200)/50 )),50), ncol = 300, byrow = TRUE),
# transient
matrix(rep(exp(atan( c((1:100)/33, rep(3,100), (100:1)/33) )),50),
ncol = 300, byrow = TRUE)
)
counts <- apply(means,2,function(cell_means){
total <- rnbinom(1, mu = 7500, size = 4)
rmultinom(1, total, cell_means)
})
rownames(counts) <- paste0('G',1:750)
colnames(counts) <- paste0('c',1:300)
> counts[1:4,1:4]
c1 c2 c3 c4
G1 0 1 2 0
G2 2 3 2 5
G3 3 8 5 4
G4 5 16 6 8
> dim(counts)
[1] 750 300
# 构建一个sce对象
sim <- SingleCellExperiment(assays = List(counts = counts))
> sim
class: SingleCellExperiment
dim: 750 300
metadata(0):
assays(1): counts
rownames(750): G1 G2 ... G749 G750
rowData names(0):
colnames(300): c1 c2 ... c299 c300
colData names(0):
reducedDimNames(0):
altExpNames(0):

2.2 基因过滤

geneFilter <- apply(assays(sim)$counts,1,function(x){
sum(x >= 3) >= 10
})
sim <- sim[geneFilter, ]
# 过滤掉11个基因
> dim(sim)
[1] 739 300

2.3 归一化

FQnorm <- function(counts){
rk <- apply(counts,2,rank,ties.method='min')
counts.sort <- apply(counts,2,sort)
refdist <- apply(counts.sort,1,median)
norm <- apply(rk,2,function(r){ refdist[r] })
rownames(norm) <- rownames(counts)
return(norm)
}
assays(sim)$norm <- FQnorm(assays(sim)$counts)

2.4降维

方法一:PCA
pca <- prcomp(t(log1p(assays(sim)$norm)), scale. = FALSE)
rd1 <- pca$x[,1:2]
plot(rd1, col = rgb(0,0,0,.5), pch=16, asp = 1)
方法二:diffusion maps
library(destiny, quietly = TRUE)
dm <- DiffusionMap(t(log1p(assays(sim)$norm)))
rd2 <- cbind(DC1 = dm$DC1, DC2 = dm$DC2)
plot(rd2, col = rgb(0,0,0,.5), pch=16, asp = 1)
将两种结果都保存起来
reducedDims(sim) <- SimpleList(PCA = rd1, DiffMap = rd2)

2.5 聚类

方法一: Gaussian mixture modeling
library(mclust, quietly = TRUE)
#根据PCA结果
cl1 <- Mclust(rd1)$classification
colData(sim)$GMM <- cl1
library(RColorBrewer)
plot(rd1, col = brewer.pal(9,"Set1")[cl1], pch=16, asp = 1)
方法二:k-means
cl2 <- kmeans(rd1, centers = 4)$cluster
colData(sim)$kmeans <- cl2
plot(rd1, col = brewer.pal(9,"Set1")[cl2], pch=16, asp = 1)

2.6 使用slingshot

sim <- slingshot(sim, clusterLabels = 'GMM', reducedDim = 'PCA')
summary(sim$slingPseudotime_1)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 8.633 21.118 21.415 34.367 43.186
  • 如果要把slingshot的所有结果都提取出来,可以用SlingshotDataSet
  • 像是SingleCellExperiment这一类对象的结果可以用 metadata(sim)$slingshot 获取
对结果可视化
colors <- colorRampPalette(brewer.pal(11,'Spectral')[-6])(100)
plotcol <- colors[cut(sim$slingPseudotime_1, breaks=100)]
plot(reducedDims(sim)$PCA, col = plotcol, pch=16, asp = 1)
lines(SlingshotDataSet(sim), lwd=2, col='black')

3 学习TSCAN

3.1 准备数据

library(TSCAN)
data(lpsdata)
procdata <- preprocess(lpsdata)
这个preprocess函数做了三件事:
  • 对表达量进行了log2(exp+1)
  • 去掉了在超过一半细胞中表达量小于1的基因
  • 将coefficient ofcovariance小于1的基因去掉

3.2 构建拟时序

使用exprmclust函数,进行了PCA降维以及model-based的聚类
lpsmclust <- exprmclust(procdata)
# 然后看下结果
plotmclust(lpsmclust)
这个图上很乱,但不妨碍看到分了3群
接着获得排序
lpsorder <- TSCANorder(lpsmclust)

3.3 基于找到的排序检测差异基因

使用difftest函数
diffval <- difftest(procdata,lpsorder)
根据q值找差异基因
head(row.names(diffval)[diffval$qval < 0.05])
对其中某个差异基因作图
# 以STAT2基因为例
STAT2expr <- log2(lpsdata["STAT2",]+1)
singlegeneplot(STAT2expr, TSCANorder(lpsmclust,flip=TRUE,orderonly=FALSE))

4 关于monocle

monocle2的拟时序分析前期数据准备可以看:单细胞转录组学习笔记-18-scRNA包学习Monocle2
另外monocle3可以看:跟着官网学习单细胞Monocle3
以版本2为例,基本上还是分三步走:从差异分析结果选合适基因=》降维=》细胞排序

step1: 选合适基因

ordering_genes <- row.names (subset(diff_test_res, qval < 0.01))
cds <- setOrderingFilter(cds, ordering_genes)
plot_ordering_genes(cds)

step2: 降维

# 默认使用DDRTree的方法
cds <- reduceDimension(cds, max_components = 2,
method = 'DDRTree')

step3: 细胞排序

cds <- orderCells(cds)

最后可视化

plot_cell_trajectory(cds, color_by = "Biological_Condition")
这个图就可以看到细胞的发展过程
另外,plot_genes_in_pseudotime 可以对基因在不同细胞中的表达量变化进行绘图
plot_genes_in_pseudotime(cds[cg,],
color_by = "Biological_Condition")

TODO: