scRNA-小鼠发育学习笔记-5-发育谱系推断及可视化

Oct 24, 2019 4 min read cnposts

刘小泽写于19.10.24 笔记目的：根据生信技能树的单细胞转录组课程探索Smartseq2技术及发育相关的分析 课程链接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=55 这次会介绍如何利用diffusion-map和Slingshot进行发育谱系推断，并结合作者包装的代码进行可视化。对应视频第三单元12-13讲

前言

细胞的变化是连续性的，它们从一个时间到另一个时间的变化轨迹是非常需要了解的，这也就是为何谱系推断这么重要的原因。

有了表达矩阵、高变化基因、分群信息和发育时期信息，就能进行谱系推断，有很多方法可以构建发育谱系，比如DiffusionMap、Slingshot

image-20191023110515767

1 进行DiffusionMap

1.1 准备表达矩阵、HVGs、分群信息

# 表达矩阵
load('../female_rpkm.Rdata')
dim(females)

# HVGs
load('../step1-female-RPKM-tSNE/females_hvg_matrix.Rdata')
dim(females_data)

# 6个发育时期获取
head(colnames(females))
female_stages <- sapply(strsplit(colnames(females), "_"), `[`, 1)
names(female_stages) <- colnames(females)
table(female_stages)

# 4个cluster获取
cluster <- read.csv('../step1-female-RPKM-tSNE/female_clustering.csv')
female_clustering=cluster[,2];names(female_clustering)=cluster[,1]
table(female_clustering)

1.2 进行DiffusionMap

包装的代码很简单

female_dm <- run_diffMap(
  females_data, 
  female_clustering,
  sigma=15
)
# 这个包装的函数其实做了下面几行代码的事情
data=females_data
condition=female_clustering
sigma=15
destinyObj <- as.ExpressionSet(as.data.frame(t(data)))
destinyObj$condition <- factor(condition)
dm <- DiffusionMap(destinyObj, sigma, rotate = TRUE)

save(female_dm,females_data,female_clustering,female_stages,
     file = 'diffusionMap_output.Rdata')

1.3 作图探索

画出特征值，这个很像PCA的碎石图screeplots或者elbowplot，也是看拐点
plot_eigenVal(
  dm=female_dm
)

image-20191024212306320

# 探索4个分群
female_clusterPalette <- c(
  C1="#560047",
  C2="#a53bad", 
  C3="#eb6bac", 
  C4="#ffa8a0"
)

plot_dm_3D(
  dm=female_dm, 
  dc=c(1:3),
  condition=female_clustering, 
  colour=female_clusterPalette
)
# 探索6个发育时间
female_stagePalette=c(
  E10.5="#2754b5", 
  E11.5="#8a00b0", 
  E12.5="#d20e0f", 
  E13.5="#f77f05", 
  E16.5="#f9db21",
  P6="#43f14b"
)
plot_dm_3D(
  dm=female_dm, 
  dc=c(1:3),
  condition=female_stages, 
  colour=female_stagePalette
)

可以看到，这个函数主要就是选取了前3个DC成分，做了三维空间的映射，然后把点的颜色分别按照cluster和stage两种不同的属性上色

2 进行Slingshot

将会利用前面DiffusionMap的结果

看上面的plot_eigenVal结果，看到DC4处是一个拐点，于是可以认为前4个DC是重要的

dm=female_dm
dim=c(1:4)
condition=factor(female_clustering)

data <- data.frame(
  dm@eigenvectors[,dim]
)

female_lineage <- slingshot(
  data, 
  condition, 
  start.clus = "C1", 
  end.clus=c("C2", "C4"),
  maxit=100000,
  shrink.method="cosine"
  # shrink.method="tricube"
)
# 看下结果
> female_lineage
class: SlingshotDataSet 

 Samples Dimensions
     563          4

lineages: 2 
Lineage1: C1  C3  C4  
Lineage2: C1  C2  

curves: 2 
Curve1: Length: 1.3739	Samples: 453.62
Curve2: Length: 0.74646	Samples: 312.73

它推断的细胞发育谱系结果在：

female_pseudotime <- get_pseudotime(female_lineage, wthres=0.9)
rownames(female_pseudotime) <- colnames(females)

image-20191024215602738

从这个结果可以看出：行名是细胞，curve1是第一条推断的发育轨迹，curve2是第二条；每个细胞在不同轨迹中所处的位置不同，并且有的细胞只在第一条轨迹中存在，在第二条中就是NA

再回来看这个发育轨迹图，

  • 最初是由E10.5作为起点发育的，然后分化成两个方向

  • 看到P6这个时期是在两条轨迹中都存在的，说明这个时期的细胞存在异质性。也就是说，虽然都是P6时期取的细胞，但是它们实际的发育方向是不用的

  • 这个算法就是帮助我们认识到不同时期内部的各个细胞，它们依然还存在着差异

    image-20191024215816009

小结

构建发育谱系，先走一下DiffusionMap的流程，得到几个重要的DC；接着走一下Slingshot函数，就会得到谱系结果

3 谱系发育相关基因可视化

最初我们知道细胞有6个时期（就是取样的6个时间点）；然后进行聚类发现这些细胞能分成4个cluster（意思就是虽然是一个时间点取的细胞，依然可能属于不同类型）；后来进行谱系推断，又增加了一个细胞属性（就是2条不同的发育轨迹）

3.1 载入之前结果

rm(list = ls()) 
options(warn=-1) 
options(stringsAsFactors = F)
source("../analysis_functions.R")

load('../female_rpkm.Rdata')
load(file = 'step4.1-diffusionMap_output.Rdata')
load(file = 'step4.2-female_pseudotime.Rdata')

3.2 对谱系推断结果进行归一化

目的就是让两条轨迹可以比较，采用的方法就是每条轨迹的每个值分别除以各自的最大值

## 第一条
pseudotime_lin <- female_pseudotime[,"curve1"]
max_pseudotime <- max(pseudotime_lin, na.rm = TRUE)
pseudotime_lin1_percent <- (pseudotime_lin*100)/max_pseudotime
## 第二条
pseudotime_lin <- female_pseudotime[,"curve2"]
max_pseudotime <- max(pseudotime_lin, na.rm = TRUE)
pseudotime_lin2_percent <- (pseudotime_lin*100)/max_pseudotime
# 现在female_pseudotime中的两条轨迹结果都在0-100之间了
female_pseudotime[,"curve1"] <- pseudotime_lin1_percent
female_pseudotime[,"curve2"] <- pseudotime_lin2_percent

不同于stage和cluster两种细胞属性，这个发育谱系属性是连续型的 。既然细胞是按照一定顺序排列的，那么就会有一些基因表达量会跟着这个连续变量进行变化

这也正是单细胞数据的优势所在，过去只有离散型的分类变量，因此只能先通过差异分析得到结果，然后对结果去注释。

3.3 可视化

将感兴趣的基因在感兴趣的谱系中进行展示

作者包装的代码非常复杂，一个包装好的函数就有140多行代码

## 先给一个颜色
# 6个stage颜色
female_stagePalette <- c(
  E10.5="#2754b5", 
  E11.5="#8a00b0", 
  E12.5="#d20e0f", 
  E13.5="#f77f05", 
  E16.5="#f9db21",
  P6="#43f14b"
)
# 4个cluster颜色
female_clusterPalette <- c(
  C1="#560047",
  C2="#a53bad", 
  C3="#eb6bac", 
  C4="#ffa8a0"
)
# 2个发育谱系颜色
female_clusterPalette2 <- c(
  "#ff6663", 
  "#3b3561"
)

## 做第一个谱系的一个基因（以Amhr2基因为例）
plot_smoothed_gene_per_lineage(
  rpkm_matrix=females, # RPKM表达矩阵
  pseudotime=female_pseudotime,  #谱系推断结果
  lin=c(1), # 对第一个谱系操作
  gene="Amhr2",  #画Amhr2基因变化
  stages=female_stages, # 发育时间点分类
  clusters=female_clustering, # cluster分类
  stage_colors=female_stagePalette,
  cluster_colors=female_clusterPalette,
  lineage_colors=female_clusterPalette2
)

从得到的图可以看出，这个加入了散点geom_point、平滑线geom_smooth、地毯线等等

image-20191024223743209

# 做某个基因在两个谱系中的变化（以Amhr2基因为例）
plot_smoothed_gene_per_lineage(
  rpkm_matrix=females, 
  pseudotime=female_pseudotime, 
  lin=c(1,2),
  gene="Amhr2", 
  stages=female_stages, 
  clusters=female_clustering, 
  stage_colors=female_stagePalette,
  cluster_colors=female_clusterPalette,
  lineage_colors=female_clusterPalette2
)

image-20191024223901393

看到这个Amhr2基因在第一个谱系中变化很大，尤其是到后期；而在第二个谱系中基本保持平衡，这就说明这个基因就是第一个谱系中重要的基因

最后就是对多个基因批量作图
gene_list <- c(
  "Sall4",
  "Sox11",
  "Gata4",
  "Lgr5",
  "Runx1",
  "Foxl2",
  "Hey2",
  "Wnt5a",
  "Pdgfra",
  "Nr2f2",
)

plot_smoothed_genes <- function(genes, lin){
  female_clusterPalette2 <- c("#ff6663", "#3b3561")
  for (gene in genes){
    plot_smoothed_gene_per_lineage(
      rpkm_matrix=females, 
      pseudotime=female_pseudotime, 
      lin=lin,
      gene=gene, 
      stages=female_stages, 
      clusters=female_clustering, 
      stage_colors=female_stagePalette,
      cluster_colors=female_clusterPalette,
      lineage_colors=female_clusterPalette2
    )
  }
}

pdf("interesting_genes_in_lineage.pdf", width=4, height=4)
plot_smoothed_genes(gene_list, 1) # plot only lineage 1
plot_smoothed_genes(gene_list, 2) # plot only lineage 2
plot_smoothed_genes(gene_list, c(1,2)) # plot the two moleages in the same graph to see the divergence
dev.off()
scRNA
Yunze Liu
Yunze Liu
Bioinformatics Sharer

Co-founder of Bioinfoplanet(生信星球)

Next
Previous

Related