3.12 细胞轨迹推断

刘小泽写于2020.7.16

1 前言
接触单细胞数据，相信你一定听过：轨迹推断，英文名词是： Trajectory Analysis。和它相关的另一个名词是：拟时序分析（pseudotime），指的是细胞沿着这个轨迹，并且对潜在的生物活动进行量化。注意这里看字面意思就知道，并不是指真正的时间，而是指细胞与细胞之间的更替、转化的顺序或者是轨迹，可以理解为“一个连续过程的缩影”。
不同的生物过程对应的“拟时序”也是不同的：
图片来自：https://www.youtube.com/watch?v=XmHDexCtjyw
许多生物过程都伴随着细胞状态的连续性变化，比如研究发育就会经常使用到。我们可以利用单细胞数据在高维空间画一条线，贯穿于多种细胞状态。最简单是点到点的一条路径，更复杂的还有一个点出发再生成多个分支。
看一下现在做相关分析的工具：
自2014年以后，已经开发出了超50种方法，那么选择何种方法进行分析成为了一个难题，因为我们不可能每一种都试一下，但有评测文章发现：Slingshot、TSCAN、Monocle DDRTree这几种方法都不错
当然还有其他的几种方法值得推荐，还做成了一个流程图方便查阅 如果对评测感兴趣，可以看看dynverse​
资源整理
ELIXIR-SE 公开课视频推荐：Trajectory inference analysis of scRNA-seq data from ELIXIR-SE 
他们也有线上课程资料​
Anthony Gitter整理了一份轨迹推断工具的清单​
做这个分析之前，最好先问几个问题：
确定数据会体现发育轨迹吗？也就是研究的样本是不是和发育相关的？
数据中的细胞会体现出中间态吗？
是否认为轨迹会出现分支？
并且要注意：
任何数据都可以强行画出轨迹，但不一定都有生物学意义！
先要保证目前找到的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
说明文档在：https://bioconductor.org/packages/3.11/bioc/vignettes/TSCAN/inst/doc/TSCAN.pdf​
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：
4 ”关于monocle“章节中的单细胞转录组学习笔记-18-scRNA包学习Monocle2和
​跟着官网学习单细胞Monocle3 放到补充篇

最近更新 2yr ago