桑基图在单细胞数据探索中的应用

男，

一个长大了才会遇到的帅哥，

稳健，潇洒，大方，靠谱。

一段生信缘，一棵技能树，

一枚大型测序工厂的螺丝钉，

一个随机森林中提灯觅食的津门旅客。

什么是桑基图

桑基图（Sankey diagram），即桑基能量分流图，也叫桑基能量平衡图。它是一种特定类型的流程图，图中延伸的分支的宽度对应数据流量的大小，比较适用于用户流量等数据的可视化分析。因1898年Matthew Henry Phineas Riall Sankey绘制的“蒸汽机的能源效率图”而闻名，此后便以其名字命名为“桑基图”。

用一个故事来介绍一下桑基图：

这个非常著名的图是Charles Minard在1869年所作的拿破仑东征俄国的信息图。Charles Minard是信息图表的之父，他是信息图领域的创始者。这张图描绘的是拿破仑在1812到1813年进攻俄国的情况。它的背景是一个真实的地图，西边是波兰的边境，东边是莫斯科。图上那条主线的宽度代表拿破仑军队的人数，黄色表示进攻路线，黑色表示撤退的路线：他开始于42万人，在向莫斯科进军的过程中丧失了很多人，到达莫斯科时只剩下10万人，而最后从莫斯科活着返回的只剩下1万人。

为什么说这个图好呢，因为除了主线的宽度之外，这张图还告诉了你更多的东西。画面下面的折线图告诉你当时的温度，其中最高的点是0度，最低到达过零下30度，回城的黑线周围还标注了月份，可以看出，拿破仑的军队在到达莫斯科的时候已经是将近十月分了，等到完全撤离俄国已经是12月份了，如果你仔细观察，会发现在撤退过程中他们路过了一条叫Studienska的河，军队人数在河两岸出现了剧减，原来那个时候天气寒冷，军队长促情况下淌水过河，于是在这条寒冷的河中冻死了很多人。

根据Edward Tufte所总结的信息设计原则：

• 这个图让显示出了比较关系（Show comparisons, contrasts, differences），比如军队人数的起始时候的宽度和结束时候的宽度的强烈对比，比如过那条河流的时候军队人数的剧烈的变化等等。

• 这个图解释了因果关系（Show causality, mechanism, structure, explanation），比如时间，温度和军队人数的关系。

• 这个图有多个变量（Multivariate analysis），1), 军队人数。2), 地理的位置（经度和纬度）3), 军队的行进方向。4), 温度。5), 时间。

所有的这些信息都不是独立存在的，他们结合在一起，将观众带入当时的拿破仑的旅程，同时能让人感受到无情的战争夺走人们生命的痛苦。

桑基图怎么看

• 线条的走向

• 粗细的变化

• 节点间的比较

绘制属于自己的桑基图

在单细胞数据分析中有一个关键的步骤FindClusters（分群，以启发样本中可能有的细胞类型数量），但是这个目前用的方法是非监督聚类，也就是数据驱动的，不依赖生物学背景。而且常常带来参数诅咒：如kmeans的K值不同，得到的分群数量不同；Seurat中FindClusters的不同 resolution 参数也会带来不同的分群数量。

于是，我的样本中到底有多少细胞类型？

所以只靠一个参数，往往不能满足要求，或者说启发的力度还不够。那就尝试多个分群参数吧，得到的结果可能是这样的：

> head([e-mail])
                  orig.ident nCount_RNA nFeature_RNA RNA_snn_res.0.8 letter.idents groups RNA_snn_res.1 RNA_snn_res.0.4
ATGCCAGAACGACT SeuratProject         70           47               0             A     g2             0               0
CATGGCCTGTGCAT SeuratProject         85           52               0             A     g1             0               0
GAACCTGATGAACC SeuratProject         87           50               0             B     g2             0               0
TGACTGGATTCTCA SeuratProject        127           56               0             A     g2             0               0
AGTCAGACTGCACA SeuratProject        173           53               0             A     g2             0               0
TCTGATACACGTGT SeuratProject         70           48               0             A     g1             0               0
               RNA_snn_res.1.2 RNA_snn_res.1.6 seurat_clusters RNA_snn_res.0.6 RNA_snn_res.1.4
ATGCCAGAACGACT               0               3               3               0               4
CATGGCCTGTGCAT               5              10              10               0               8
GAACCTGATGAACC               5               9               9               0               6
TGACTGGATTCTCA               0               7               7               0               1
AGTCAGACTGCACA               0               3               3               0               4
TCTGATACACGTGT               0               3               3               0               4


如果用人类的肉眼来比较不同RNA_snn_res.下的分群结果可能会比较困难。不过，借助R方便地看出某一分群下，每个群的细胞数量：

> table([e-mail]$RNA_snn_res.1.6)

0  1 10 11  2  3  4  5  6  7  8  9 
17 14  2  2  9  6  7  5  5  3  6  4 


但是有了桑基图情况就不一样了，变得一目了然起来：

#先执行不同resolution 下的分群
library(Seurat)
pbmc_small <- FindClusters(
  object = pbmc_small,
  resolution = c(seq(.4,1.6,.2))
)


绘制细胞分群的桑基图：

#install.packages("ggalluvial")
library(ggalluvial)
library(tidyverse)

head([e-mail])
ggplot(data = [e-mail],
       aes(axis1 = RNA_snn_res.0.4, axis2 = RNA_snn_res.0.6,axis3 = RNA_snn_res.0.8,axis4 = RNA_snn_res.1,
           axis5 = RNA_snn_res.1.2,axis6 = RNA_snn_res.1.4,axis7 = RNA_snn_res.1.6)) +
  scale_x_discrete(limits = c(paste0("RNA_snn_res.",seq(.4,1.6,.2))), expand = c(.01, .05)) +
  geom_alluvium(aes(fill = RNA_snn_res.1.6)) +
  geom_stratum() + geom_text(stat = "stratum", infer.label = TRUE) +
  #coord_polar()+
  theme(axis.text.x = element_text(angle = 90, hjust = 1))+
  ggtitle("cell number in each cluster")


可以清晰地看出，RNA_snn_res.1.6每个群的来源，也可以启发在RNA_snn_res.1.2时cluster0可能有三个亚群，cluster4有两个亚群。这不仅解析了resolution 参数（其他的分群算法也一样），同时启发了样本的异质性。

有了这个桑基图的框架，其实很多我们想在这个图上展示的metadata信息就变得容易了，比如我们可以看一下某一细胞类型流向或者样本的流向，只需要在metadata中加上一列即可。

ggplot(data = [e-mail],
       aes(axis1 = RNA_snn_res.0.4, axis2 = RNA_snn_res.0.6,axis3 = RNA_snn_res.0.8,axis4 = RNA_snn_res.1,
           axis5 = RNA_snn_res.1.2,axis6 = RNA_snn_res.1.4,axis7 = RNA_snn_res.1.6)) +
  scale_x_discrete(limits = c(paste0("RNA_snn_res.",seq(.4,1.6,.2))), expand = c(.01, .05)) +
  geom_alluvium(aes(fill = groups)) +
  geom_stratum() + geom_text(stat = "stratum", infer.label = TRUE) +
  #coord_polar()+
  theme(axis.text.x = element_text(angle = 90, hjust = 1))+
  ggtitle("cell number in each cluster")


或者看nFeature_RNA 的变化，即是不是nFeature_RNA 高的分到一群呢？

ggplot(data = [e-mail],
       aes(axis1 = RNA_snn_res.0.4, axis2 = RNA_snn_res.0.6,axis3 = RNA_snn_res.0.8,axis4 = RNA_snn_res.1,
           axis5 = RNA_snn_res.1.2,axis6 = RNA_snn_res.1.4,axis7 = RNA_snn_res.1.6)) +
  scale_x_discrete(limits = c(paste0("RNA_snn_res.",seq(.4,1.6,.2))), expand = c(.01, .05)) +
  geom_alluvium(aes(fill = nFeature_RNA )) +
  geom_stratum() + geom_text(stat = "stratum", infer.label = TRUE) +
  #coord_polar()+
  theme(axis.text.x = element_text(angle = 90, hjust = 1))+
  ggtitle("cell number in each cluster")


image

clustertree

在聚类分析中，由于它的启发性本质，经常需要比较不同分群的结果。下面提供另一种(简单直白的)“桑基图“，供大家参考：

clustree([e-mail], prefix = "RNA_snn_res.")