发育谱系推断及可视化
作者 | 单细胞天地小编 刘小泽
课程链接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=55
这次会介绍如何利用diffusion-map和Slingshot进行发育谱系推断，并结合作者包装的代码进行可视化。对应视频第三单元12-13讲
前言
细胞的变化是连续性的，它们从一个时间到另一个时间的变化轨迹是非常需要了解的，这也就是为何谱系推断这么重要的原因。
有了表达矩阵、高变化基因、分群信息和发育时期信息，就能进行谱系推断，有很多方法可以构建发育谱系，比如DiffusionMap、Slingshot
1 进行DiffusionMap
1.1 准备表达矩阵、HVGs、分群信息
1# 表达矩阵
2load('../female_rpkm.Rdata')
3dim(females)
4
5# HVGs
6load('../step1-female-RPKM-tSNE/females_hvg_matrix.Rdata')
7dim(females_data)
8
9# 6个发育时期获取
10head(colnames(females))
11female_stages <- sapply(strsplit(colnames(females), "_"), `[`, 1)
12names(female_stages) <- colnames(females)
13table(female_stages)
14
15# 4个cluster获取
16cluster <- read.csv('../step1-female-RPKM-tSNE/female_clustering.csv')
17female_clustering=cluster[,2];names(female_clustering)=cluster[,1]
18table(female_clustering)
1.2 进行DiffusionMap
包装的代码很简单
1female_dm <- run_diffMap(
2 females_data,
3 female_clustering,
4 sigma=15
5)
6# 这个包装的函数其实做了下面几行代码的事情
7data=females_data
8condition=female_clustering
9sigma=15
10destinyObj <- as.ExpressionSet(as.data.frame(t(data)))
11destinyObj$condition <- factor(condition)
12dm <- DiffusionMap(destinyObj, sigma, rotate = TRUE)
13
14save(female_dm,females_data,female_clustering,female_stages,
15 file = 'diffusionMap_output.Rdata')
1.3 作图探索
画出特征值，这个很像PCA的碎石图screeplots或者elbowplot，也是看拐点
1plot_eigenVal(
2 dm=female_dm
3)
1# 探索4个分群
2female_clusterPalette <- c(
3 C1="#560047",
4 C2="#a53bad",
5 C3="#eb6bac",
6 C4="#ffa8a0"
7)
8
9plot_dm_3D(
10 dm=female_dm,
11 dc=c(1:3),
12 condition=female_clustering,
13 colour=female_clusterPalette
14)
15# 探索6个发育时间
16female_stagePalette=c(
17 E10.5="#2754b5",
18 E11.5="#8a00b0",
19 E12.5="#d20e0f",
20 E13.5="#f77f05",
21 E16.5="#f9db21",
22 P6="#43f14b"
23)
24plot_dm_3D(
25 dm=female_dm,
26 dc=c(1:3),
27 condition=female_stages,
28 colour=female_stagePalette
29)
可以看到，这个函数主要就是选取了前3个DC成分，做了三维空间的映射，然后把点的颜色分别按照cluster和stage两种不同的属性上色
2 进行Slingshot
将会利用前面DiffusionMap的结果
看上面的plot_eigenVal结果，看到DC4处是一个拐点，于是可以认为前4个DC是重要的
1dm=female_dm
2dim=c(1:4)
3condition=factor(female_clustering)
4
5data <- data.frame(
6 dm@eigenvectors[,dim]
7)
8
9female_lineage <- slingshot(
10 data,
11 condition,
12 start.clus = "C1",
13 end.clus=c("C2", "C4"),
14 maxit=100000,
15 shrink.method="cosine"
16 # shrink.method="tricube"
17)
18# 看下结果
19> female_lineage
20class: SlingshotDataSet
21
22 Samples Dimensions
23 563 4
24
25lineages: 2
26Lineage1: C1 C3 C4
27Lineage2: C1 C2
28
29curves: 2
30Curve1: Length: 1.3739 Samples: 453.62
31Curve2: Length: 0.74646 Samples: 312.73
它推断的细胞发育谱系结果在：
1female_pseudotime <- get_pseudotime(female_lineage, wthres=0.9)
2rownames(female_pseudotime) <- colnames(females)
从这个结果可以看出：行名是细胞，curve1是第一条推断的发育轨迹，curve2是第二条；每个细胞在不同轨迹中所处的位置不同，并且有的细胞只在第一条轨迹中存在，在第二条中就是NA
再回来看这个发育轨迹图，
最初是由E10.5作为起点发育的，然后分化成两个方向
看到P6这个时期是在两条轨迹中都存在的，说明这个时期的细胞存在异质性。也就是说，虽然都是P6时期取的细胞，但是它们实际的发育方向是不用的
这个算法就是帮助我们认识到不同时期内部的各个细胞，它们依然还存在着差异
小结
构建发育谱系，先走一下DiffusionMap的流程，得到几个重要的DC；接着走一下Slingshot函数，就会得到谱系结果
3 谱系发育相关基因可视化
最初我们知道细胞有6个时期（就是取样的6个时间点）；然后进行聚类发现这些细胞能分成4个cluster（意思就是虽然是一个时间点取的细胞，依然可能属于不同类型）；后来进行谱系推断，又增加了一个细胞属性（就是2条不同的发育轨迹）
3.1 载入之前结果
1rm(list = ls())
2options(warn=-1)
3options(stringsAsFactors = F)
4source("../analysis_functions.R")
5
6load('../female_rpkm.Rdata')
7load(file = 'step4.1-diffusionMap_output.Rdata')
8load(file = 'step4.2-female_pseudotime.Rdata')
3.2 对谱系推断结果进行归一化
目的就是让两条轨迹可以比较，采用的方法就是每条轨迹的每个值分别除以各自的最大值
1## 第一条
2pseudotime_lin <- female_pseudotime[,"curve1"]
3max_pseudotime <- max(pseudotime_lin, na.rm = TRUE)
4pseudotime_lin1_percent <- (pseudotime_lin*100)/max_pseudotime
5## 第二条
6pseudotime_lin <- female_pseudotime[,"curve2"]
7max_pseudotime <- max(pseudotime_lin, na.rm = TRUE)
8pseudotime_lin2_percent <- (pseudotime_lin*100)/max_pseudotime
9# 现在female_pseudotime中的两条轨迹结果都在0-100之间了
10female_pseudotime[,"curve1"] <- pseudotime_lin1_percent
11female_pseudotime[,"curve2"] <- pseudotime_lin2_percent
不同于stage和cluster两种细胞属性，这个发育谱系属性是连续型的 。既然细胞是按照一定顺序排列的，那么就会有一些基因表达量会跟着这个连续变量进行变化
这也正是单细胞数据的优势所在，过去只有离散型的分类变量，因此只能先通过差异分析得到结果，然后对结果去注释。
3.3 可视化
将感兴趣的基因在感兴趣的谱系中进行展示
作者包装的代码非常复杂，一个包装好的函数就有140多行代码
1## 先给一个颜色
2# 6个stage颜色
3female_stagePalette <- c(
4 E10.5="#2754b5",
5 E11.5="#8a00b0",
6 E12.5="#d20e0f",
7 E13.5="#f77f05",
8 E16.5="#f9db21",
9 P6="#43f14b"
10)
11# 4个cluster颜色
12female_clusterPalette <- c(
13 C1="#560047",
14 C2="#a53bad",
15 C3="#eb6bac",
16 C4="#ffa8a0"
17)
18# 2个发育谱系颜色
19female_clusterPalette2 <- c(
20 "#ff6663",
21 "#3b3561"
22)
23
24## 做第一个谱系的一个基因（以Amhr2基因为例）
25plot_smoothed_gene_per_lineage(
26 rpkm_matrix=females, # RPKM表达矩阵
27 pseudotime=female_pseudotime, #谱系推断结果
28 lin=c(1), # 对第一个谱系操作
29 gene="Amhr2", #画Amhr2基因变化
30 stages=female_stages, # 发育时间点分类
31 clusters=female_clustering, # cluster分类
32 stage_colors=female_stagePalette,
33 cluster_colors=female_clusterPalette,
34 lineage_colors=female_clusterPalette2
35)
从得到的图可以看出，这个加入了散点geom_point、平滑线geom_smooth、地毯线等等
1# 做某个基因在两个谱系中的变化（以Amhr2基因为例）
2plot_smoothed_gene_per_lineage(
3 rpkm_matrix=females,
4 pseudotime=female_pseudotime,
5 lin=c(1,2),
6 gene="Amhr2",
7 stages=female_stages,
8 clusters=female_clustering,
9 stage_colors=female_stagePalette,
10 cluster_colors=female_clusterPalette,
11 lineage_colors=female_clusterPalette2
12)
看到这个Amhr2基因在第一个谱系中变化很大，尤其是到后期；而在第二个谱系中基本保持平衡，这就说明这个基因就是第一个谱系中重要的基因
最后就是对多个基因批量作图
1gene_list <- c(
2 "Sall4",
3 "Sox11",
4 "Gata4",
5 "Lgr5",
6 "Runx1",
7 "Foxl2",
8 "Hey2",
9 "Wnt5a",
10 "Pdgfra",
11 "Nr2f2",
12)
13
14plot_smoothed_genes <- function(genes, lin){
15 female_clusterPalette2 <- c("#ff6663", "#3b3561")
16 for (gene in genes){
17 plot_smoothed_gene_per_lineage(
18 rpkm_matrix=females,
19 pseudotime=female_pseudotime,
20 lin=lin,
21 gene=gene,
22 stages=female_stages,
23 clusters=female_clustering,
24 stage_colors=female_stagePalette,
25 cluster_colors=female_clusterPalette,
26 lineage_colors=female_clusterPalette2
27 )
28 }
29}
30
31pdf("interesting_genes_in_lineage.pdf", width=4, height=4)
32plot_smoothed_genes(gene_list, 1) # plot only lineage 1
33plot_smoothed_genes(gene_list, 2) # plot only lineage 2
34plot_smoothed_genes(gene_list, c(1,2)) # plot the two moleages in the same graph to see the divergence
35dev.off()