单细胞交响乐
  • 前言:我与《单细胞交响乐》的缘分
  • 1 准备篇:背景知识
    • 1.1 数据结构
    • 1.2 总览 | 从实验到分析
  • 2 积累篇:文献阅读
    • 2.1.1 综述 | 2019-单细胞转录组分析最佳思路
    • 2.1.2 综述 | 2018-单细胞捕获平台
    • 2.1.3 综述 | 2017-scRNA中的细胞聚类分群
    • 2.1.4 综述 | scRNA已经开发出超过1000款工具了,你用过几种?
    • 2.1.5 综述 | 2021-单细胞测序的微流控技术应用
    • 2.2.1 研究 | 2018-单细胞转录组探索癌症免疫治疗获得性抗性机理
    • 2.2.2 研究 | 2018-人类结直肠癌单细胞多组学分析
    • 2.2.3 研究 | 2020-单细胞分析揭示葡萄膜黑色素瘤新的进化复杂性
    • 2.2.4 研究 | 2020-COVID-19病人支气管免疫细胞单细胞测序分析
    • 2.2.5 研究 | 2020-原汁原味读--单细胞肿瘤免疫图谱
    • 2.2.6 研究 | 2021-多发性骨髓瘤发展过程中肿瘤和免疫细胞的共同进化
    • 2.2.7 研究 | 2021-多个组织的成纤维细胞图谱
    • 2.2.8 研究 | 2021-多组学分析肺结核队列的记忆T细胞状态
    • 2.2.9 研究 | 2021-CancerSCEM: 人类癌症单细胞表达图谱数据库
    • 2.2.10 研究| 2021-单细胞转录组分析COVID-19重症患者肺泡巨噬细胞亚型
    • 2.2.11 研究 |2021-单细胞转录组揭示肺腺癌特有的肿瘤微环境
    • 2.2.12 研究 | 2021-单细胞转录组揭示乳头状甲状腺癌起始与发展
    • 2.2.13 研究 | 2021-解析食管鳞癌化疗病人的单细胞转录组
    • 2.2.14 研究 | 2021-单细胞水平看骨髓瘤的细胞状态和基因调控
    • 2.3.1 算法|2020-BatchBench比较scRNA批次矫正方法
    • 2.3.2 算法 | 2021-scPhere——用地球仪来展示降维结果
    • 2.3.3 算法 | 2021-单细胞差异分析方法评测
    • 2.3.4 算法 | 2021-细胞分群新方法——CNA(co-varying neighborhood analysis)
    • 2.3.5 工具 | 2018-iSEE:单细胞数据可视化辅助网页工具
    • 2.3.6 工具 | 2021-MACA: 一款自动注释细胞类型的工具
    • 2.3.7 工具 | 2021-一个很有想法的工具——Ikarus,想要在单细胞水平直接鉴定肿瘤细胞
  • 3 流程篇:分析框架
    • 3.1 质控
    • 3.2 归一化
    • 3.3 挑选表达量高变化基因
    • 3.4 降维
    • 3.5 聚类
    • 3.6 Marker/标记基因检测
    • 3.7 细胞类型注释
    • 3.8 批次效应处理
    • 3.9 多样本间差异分析
    • 3.10 检测Doublet
    • 3.11 细胞周期推断
    • 3.12 细胞轨迹推断
    • 3.13 与蛋白丰度信息结合
    • 3.14 处理大型数据
    • 3.15 不同R包数据的相互转换
  • 4 实战篇:活学活用
    • 4.1 实战一 | Smart-seq2 | 小鼠骨髓
    • 4.2 实战二 | STRT-Seq | 小鼠大脑
    • 4.3 实战三 | 10X | 未过滤的PBMC
    • 4.4 实战四 | 10X | 过滤后的PBMC
    • 4.5 实战五 | CEL-seq2 | 人胰腺细胞
    • 4.6 实战六 | CEL-seq | 人胰腺细胞
    • 4.7 实战七 | SMARTer | 人胰腺细胞
    • 4.8 实战八 | Smart-seq2 | 人胰腺细胞
    • 4.9 实战九 | 不同技术数据整合 | 人胰腺细胞
    • 4.10 实战十 | CEL-seq | 小鼠造血干细胞
    • 4.11 实战十一 | Smart-seq2 | 小鼠造血干细胞
    • 4.12 实战十二 | 10X | 小鼠嵌合体胚胎
    • 4.13 实战十三 | 10X | 小鼠乳腺上皮细胞
    • 4.14 | 实战十四 | 10X | HCA计划的38万骨髓细胞
  • 5 补充篇:开拓思路
    • 5.1 10X Genomics概述
      • 5.1.1 10X Genomics 问题集锦
    • 5.2 CellRanger篇
      • 5.2.1 CellRanger实战(一)数据下载
      • 5.2.2 CellRanger实战(二) 使用前注意事项
      • 5.2.3 CellRanger实战(三) 使用初探
      • 5.2.4 CellRanger实战(四)流程概览
      • 5.2.5 CellRanger实战(五) 理解count输出的结果
    • 5.3 Seurat的使用
      • 5.3.1 Seurat V3 | 实战之2700 PBMCs分析
      • 5.3.2 Seurat V3 | 如何改造Seurat包的DoHeatmap函数?
      • 5.3.3 scRNA的3大R包对比
      • 5.3.4 Seurat两种数据比较:integrated vs RNA assay
      • 5.3.5 seurat 的几种findmaker比较
    • 5.4 Monocle的使用
      • 5.4.1 Monocle V3实战
    • 5.5 多个数据集的整合
      • 5.5.1 使用Seurat的merge功能进行整合
      • 5.5.2 如何使用sctransform去除批次效应
由 GitBook 提供支持
在本页
  • 1 前言
  • 资源整理
  • 做这个分析之前,最好先问几个问题:
  • 并且要注意:
  • 2 学习Slingshot
  • 2.1 数据准备
  • 2.2 基因过滤
  • 2.3 归一化
  • 2.4降维
  • 2.5 聚类
  • 2.6 使用slingshot
  • 3 学习TSCAN
  • 3.1 准备数据
  • 3.2 构建拟时序
  • 3.3 基于找到的排序检测差异基因
  • 4 关于monocle
  • step1: 选合适基因
  • step2: 降维
  • step3: 细胞排序
  • 最后可视化
  • TODO:

这有帮助吗?

  1. 3 流程篇:分析框架

3.12 细胞轨迹推断

刘小泽写于2020.7.16

上一页3.11 细胞周期推断下一页3.13 与蛋白丰度信息结合

最后更新于4年前

这有帮助吗?

1 前言

接触单细胞数据,相信你一定听过:轨迹推断,英文名词是: Trajectory Analysis。和它相关的另一个名词是:拟时序分析(pseudotime),指的是细胞沿着这个轨迹,并且对潜在的生物活动进行量化。注意这里看字面意思就知道,并不是指真正的时间,而是指细胞与细胞之间的更替、转化的顺序或者是轨迹,可以理解为“一个连续过程的缩影”。

不同的生物过程对应的“拟时序”也是不同的:

图片来自:https://www.youtube.com/watch?v=XmHDexCtjyw

许多生物过程都伴随着细胞状态的连续性变化,比如研究发育就会经常使用到。我们可以利用单细胞数据在高维空间画一条线,贯穿于多种细胞状态。最简单是点到点的一条路径,更复杂的还有一个点出发再生成多个分支。

资源整理

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

  • 确定数据会体现发育轨迹吗?也就是研究的样本是不是和发育相关的?

  • 数据中的细胞会体现出中间态吗?

  • 是否认为轨迹会出现分支?

并且要注意:

  • 任何数据都可以强行画出轨迹,但不一定都有生物学意义!

  • 先要保证目前找到的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

以版本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:

看一下:

自2014年以后,已经开发出了超50种方法,那么选择何种方法进行分析成为了一个难题,因为我们不可能每一种都试一下,但有评测文章发现:、、这几种方法都不错

当然还有其他的几种方法值得推荐,还做成了一个流程图方便查阅 如果对评测感兴趣,可以看看

ELIXIR-SE 公开课视频推荐: 他们也有

Anthony Gitter整理了一份

说明文档在:

monocle2的拟时序分析前期数据准备可以看:

另外monocle3可以看:

4 ”关于monocle“章节中的和

放到补充篇

现在做相关分析的工具
Slingshot
TSCAN
Monocle DDRTree
dynverse
Trajectory inference analysis of scRNA-seq data from ELIXIR-SE
线上课程资料
轨迹推断工具的清单
https://bioconductor.org/packages/3.11/bioc/vignettes/TSCAN/inst/doc/TSCAN.pdf
单细胞转录组学习笔记-18-scRNA包学习Monocle2
跟着官网学习单细胞Monocle3
单细胞转录组学习笔记-18-scRNA包学习Monocle2
跟着官网学习单细胞Monocle3