专栏首页单细胞天地桑基图在单细胞数据探索中的应用

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

分享是一种态度

作者 | 周运来

男,

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

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

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

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

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

什么是桑基图

桑基图(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(pbmc_small@meta.data)
                  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(pbmc_small@meta.data$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(pbmc_small@meta.data)
ggplot(data = pbmc_small@meta.data,
       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 = pbmc_small@meta.data,
       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 = pbmc_small@meta.data,
       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(pbmc_small@meta.data, prefix = "RNA_snn_res.")

本文分享自微信公众号 - 单细胞天地(sc-ngs),作者:周运来

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-03-18

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 你到底想要什么样的umap/tsne图?

    随着生物学背景知识的增加,单细胞图谱的可视化直接用10X的Loup或者seurat的Dimplot函数直接绘制的umap/tsne图往往很难达到要求了,这就要求...

    生信技能树jimmy
  • 单细胞转录组高级分析四:scRNA数据推断CNV

    inferCNV是大名鼎鼎的broad研究所开发的,可以使用单细胞转录组数据分析肿瘤细胞CNV。相关文章2014年就发表在了Science上,之后算法不断优化,...

    生信技能树jimmy
  • 统计细胞检测的基因数量

    (g) number of exon mapping reads. Cutoff: 10000 (8 cells failed).

    生信技能树jimmy
  • 让你热血沸腾的电竞直播,有一群深藏不露的技术极客

    前有雅加达亚运会电竞赛事上,中国代表队交出两金一银这一近乎完美的成绩单;后有iG战队八年磨一剑,斩获LPL历史上首个英雄联盟全球总决赛冠军。

    量子位
  • 高通量测序后的实验验证手段——转录组篇(上)

    对于高通量技术大家一定都不陌生,但生物信息分析之后该做什么呢?接下来的日子,小编会和大家探讨并分享高通量测序后的实验验证,即该用什么技术做什么验证!

    生信交流平台
  • 怎样判断你提取的RNA质量好坏?

    与DNA不同,RNA是极为脆弱的,由于其单链结构,RNA的碱基和氢键全都暴露在环境中,极易被环境中的各种化学物质和RNA酶降解。一旦降解,后期的实验纯属浪费时间...

    Mark Chen
  • AI系统可以学习设计用于研究的RNA分子

    核糖核酸(RNA)存在于所有活细胞中。它充当信使,携带来自DNA(脱氧核糖核酸)的指令。当它不能正常工作时,它会严重影响神经系统、心血管和肌肉的调节过程,导致肿...

    AiTechYun
  • ios 百度地图 获取拖动或缩放手势

    在项目中遇到一个问题,在拖动或者缩放百度地图的时候要请求数据。但是百度地图SDK中没有明确如何获取拖动和缩放手势 官方推荐使用如下两个方法,通过判断状态来获取...

    赵哥窟
  • QB 期刊 | 转录组层面上研究蛋白质-RNA相互作用的技术及方法综述

    两周前,小编给大家推荐了清华大学张强峰教授在QB期刊上发表的关于转录组层面上研究RNA-RNA相互作用的干湿实验方法的综述文章后,引起了许多小伙...

    生信宝典
  • 为什么 wrk 和 ab, locust 压测的结果差异这么大?

    本地使用 wrk 压测 nginx 的一个 hello,world!接口。wrk 的 qps 是 20 多万,全部 200 状态吗,无异常。ab 结果是 1 万...

    随心助手

扫码关注云+社区

领取腾讯云代金券