什么是闪酉瘢基圖

桑基圖（Sankey diagram），即桑基能量分流圖锦募，也叫杀记常基能量平衡圖馆纳。它是一種特定類型的流程圖，圖中延伸的分支的寬度對(duì)應(yīng)數(shù)據(jù)流量的大小汹桦，比較適用于用戶流量等數(shù)據(jù)的可視化分析鲁驶。因1898年Matthew Henry Phineas Riall Sankey繪制的“蒸汽機(jī)的能源效率圖”而聞名，此后便以其名字命名為“晌杪妫基圖”钥弯。

用一個(gè)故事來介紹一下桑基圖：

這個(gè)非常著名的圖是Charles Minard在1869年所作的拿破侖東征俄國的信息圖葛作。Charles Minard是信息圖表的之父寿羞，他是信息圖領(lǐng)域的創(chuàng)始者。這張圖描繪的是拿破侖在1812到1813年進(jìn)攻俄國的情況赂蠢。它的背景是一個(gè)真實(shí)的地圖绪穆，西邊是波蘭的邊境，東邊是莫斯科虱岂。圖上那條主線的寬度代表拿破侖軍隊(duì)的人數(shù)玖院，黃色表示進(jìn)攻路線，黑色表示撤退的路線： 他開始于42萬人第岖，在向莫斯科進(jìn)軍的過程中喪失了很多人，到達(dá)莫斯科時(shí)只剩下10萬人郊酒，而最后從莫斯科活著返回的只剩下1萬人键袱。

為什么說這個(gè)圖好呢燎窘，因?yàn)槌酥骶€的寬度之外褐健，這張圖還告訴了你更多的東西。畫面下面的折線圖告訴你當(dāng)時(shí)的溫度澜汤，其中最高的點(diǎn)是0度，最低到達(dá)過零下30度谁不，回城的黑線周圍嗨標(biāo)注了月份烛缔，可以看出轩拨，拿破侖的軍隊(duì)在達(dá)打到莫斯科的時(shí)候已經(jīng)是將近十月分了，等到完全撤離俄國已經(jīng)是12月份了砍濒，如果你仔細(xì)觀察爸邢，會(huì)發(fā)現(xiàn)在撤退過程中他們路過了一條叫Studienska的河杠河，軍隊(duì)人數(shù)在河兩岸出現(xiàn)了劇減，原來那個(gè)時(shí)候天氣寒冷待诅，軍隊(duì)長促情況下淌水過河卑雁，于是在這條寒冷的河中凍死了很多人测蹲。

根據(jù)Edward Tufte所總結(jié)的信息設(shè)計(jì)原則：

• 這個(gè)圖讓顯示出了比較關(guān)系（Show comparisons, contrasts, differences）弛房，比如軍隊(duì)人數(shù)的起始時(shí)候的寬度和結(jié)束時(shí)候的寬度的強(qiáng)烈對(duì)比，比如過那條河流的時(shí)候軍隊(duì)人數(shù)的劇烈的變化等等粹排。
• 這個(gè)圖解釋了因果關(guān)系（Show causality, mechanism, structure, explanation）顽耳，比如時(shí)間膝迎，溫度和軍隊(duì)人數(shù)的關(guān)系限次。
• 這個(gè)圖有多個(gè)變量（Multivariate analysis）卖漫，1), 軍隊(duì)人數(shù)。 2), 地理的位置（經(jīng)度和緯度）3), 軍隊(duì)的行進(jìn)方向突委。 4), 溫度。 5), 時(shí)間钧唐。

所有的這些信息都不是獨(dú)立存在的钝侠，他們結(jié)合在一起，將觀眾帶入當(dāng)時(shí)的拿破侖的旅程忽舟，同時(shí)能讓人感受到無情的戰(zhàn)爭(zhēng)奪走人們生命的痛苦叮阅。

桑基圖怎么看

• 線條的走向
• 粗細(xì)的變化
• 節(jié)點(diǎn)間的比較

繪制屬于自己的衫盏基圖

在單細(xì)胞數(shù)據(jù)分析中有一個(gè)關(guān)鍵的步驟FindClusters（分群拌汇，以啟發(fā)樣本中可能有的細(xì)胞類型數(shù)量）担猛，但是這個(gè)目前用的方法是非監(jiān)督聚類，也就是數(shù)據(jù)驅(qū)動(dòng)的，不依賴生物學(xué)背景貌嫡。而且常常帶來參數(shù)詛咒：如kmeans的K值不同岛抄，得到的分群數(shù)量不同；Seurat中FindClusters的不同 resolution 參數(shù)也會(huì)帶來不同的分群數(shù)量蹭秋。

于是，我的樣本中到底有多少細(xì)胞類型洞豁？

所以只靠一個(gè)參數(shù)丈挟，往往不能滿足要求，或者說啟發(fā)的力度還不夠桐绒。那就嘗試多個(gè)分群參數(shù)吧咧叭，得到的結(jié)果可能是這樣的：

> 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.下的分群結(jié)果可能會(huì)比較困難派撕。不過终吼，借助R方便地看出某一分群下商佛，每個(gè)群的細(xì)胞數(shù)量：

> 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 

但是有了桑基圖情況就不一樣了玛追，變得一目了然起來：

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

繪制細(xì)胞分群的桑基圖：

#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每個(gè)群的來源，也可以啟發(fā)在RNA_snn_res.1.2時(shí)cluster0可能有三個(gè)亞群妇汗，cluster4有兩個(gè)亞群杨箭。這不僅解析了resolution 參數(shù)（其他的分群算法也一樣）捣郊，同時(shí)啟發(fā)了樣本的異質(zhì)性。

有了這個(gè)赡锢基圖的框架琐旁，其實(shí)很多我們想在這個(gè)圖上展示的metadata信息就變得容易了，比如我們可以看一下某一細(xì)胞類型流向或者樣本的流向，只需要在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")

clustertree

在聚類分析中麻削，由于它的啟發(fā)性本質(zhì)呛哟，經(jīng)常需要比較不同分群的結(jié)果。下面提供另一種(簡(jiǎn)單直白的)“杀罟拢基圖“淌铐，供大家參考：

clustree(pbmc_small@meta.data, prefix = "RNA_snn_res.")