可视化：两种方法绘制桑基图
今天继续分享生信分析中不那么常见的图形 -- Sankey diagram。桑基图即桑基能量分流图，是一种可视化特定流程图的技术，图中延伸的分支宽度对应数据流量的大小(即宽度与流的重要性成正比)。所以，想要用桑基图展示数据，必须满足一个条件，即数据之间有方向性的流动。因此，该图常用于展示数据的动态变化，如各地移民情况，数据源的流向情况。下面来分别使用networkD3和ggalluvial包绘制桑基图。
networkD3
1、准备数据
首先，准备好数据，包含三列的数据框，第一列是源数据名，第二列是目标数据名，第三列是流量值：
library(reshape2)
data <- read.table("https://raw.githubusercontent.com/holtzy/data_to_viz/master/Example_dataset/13_AdjacencyDirectedWeighted.csv", header=T, check.name=F)
data$source <- rownames(data)
data_long <- melt(data,id.vars='source', variable.name='target')
data_long <- subset(data_long, value > 0)
data_long$target <- paste(data_long$target, " ", sep="")
head(data_long)
source target value
1 Africa Africa 3.142471
2 East Asia East Asia 1.630997
3 South Asia East Asia 0.525881
4 South East Asia East Asia 0.145264
5 Africa Europe 2.107883
6 East Asia Europe 0.601265
......准备数据的时候需要特别注意，源数据名和目标数据名有相同的，但二者从数据处理的角度需要加以区别。例如，这里源数据的`Africa`会流向很多目标，当然也包括目标数据的 Africa，这时为了区分源数据名和目标数据名，可以采用一个很简单巧妙的方法来解决，在目标数据名的后面加一个空格，这样后续绘图时，数据处理上既能区分二者，而绘图时名称显示上又可以保持一致。真是妙啊！
......
2、绘图
使用networkD3绘图，只用上面的数据还不行，那只是半成品，还需要准备节点名数据框。同时，还要根据节点名数据框给上一步的数据框添加索引列。因为该包在画图时，使用的并不是节点的名称，而是节点的索引。准备好这些，画图就水到渠成了：
library(networkD3)
nodes <- data.frame(name=unique(c(as.character(data_long$source), as.character(data_long$target))))
data_long$IDsource <- match(data_long$source, nodes$name)-1
data_long$IDtarget <- match(data_long$target, nodes$name)-1
p <- sankeyNetwork(Links = data_long, Nodes = nodes,
Source = "IDsource", Target = "IDtarget",
Value = "value", NodeID = "name",
sinksRight=FALSE,, nodeWidth=40, fontSize=13, nodePadding=20)
p
saveNetwork(p, "sankey.html")
结果如下：
注意，节点名的索引是 0-base 的。其实，从 networkD3 包名可以顾名思义，该包可以用来画网络图，桑基图只是一个功能而已。其他功能大家自行摸索吧！
3、格式转化
虽然networkD3画桑基图挺简单的，但该包有一个硬伤就是其基于shiny的JS框架，可以保存为html格式，但没法直接保存为图片格式。若想要保存为pdf、png，需要借助webshot包，而该包又需要依赖phantomjs，格式转化前先安装好这两个东西：
library(webshot)
if(!is_phantomjs_installed()){
install_phantomjs()
}
webshot("sankey.html" , "sankey.pdf")
结果如下：
ggalluvial
该包也可以用来绘制桑基图，同时该包基于ggplot2，可扩展性就非常强了。绘图代码如下：
data <- read.table("https://raw.githubusercontent.com/holtzy/data_to_viz/master/Example_dataset/13_AdjacencyDirectedWeighted.csv", header=T, check.name=F)
data$source <- rownames(data)
data_long <- melt(data,id.vars='source', variable.name='target')
data_long <- subset(data_long, value > 0)
data_long <- melt(data_long,id.vars='value', variable.name='subject',value.name='type')
head(data_long)
value subject type
1 3.142471 source Africa
2 1.630997 source East Asia
3 0.525881 source South Asia
4 0.145264 source South East Asia
5 2.107883 source Africa
6 0.601265 source East Asia
p <- ggplot(data_long , aes(x = subject, stratum = type, alluvium = value, label =type)) +
geom_stratum(aes(fill = type),alpha= 0.5) +
geom_text(stat = "stratum", size = 2.5) +
geom_flow() +
theme_void() +
theme(legend.position= 'none')
p
结果如下：
......与 networkD3 包相比，数据格式有所不同，相当于在其基础上再将数据变长一些，由于 ggalluvial 有单独一列信息来区分 source 和 target，所以这里的 type 变量里面就不用在最后添加空格以示区分了。不得不说，还是 ggalluvial 包更具吸引力！
......
--------- The End ---------
基因组 | 转录组 | 表观组 | 单细胞
寻找灵魂的工具人
长按扫码加关注