多种方法绘制酷炫的桑基图

希望所有学员都可以站在生信技能树的舞台上发光发热！

下面是粉丝随机投稿

就因为在群里问了句怎么做桑基图，今天被健明老师“强行”安排了写教程的任务，并获赠高阶加密的桑基图手札一份，既然无法拒绝，那就硬着头皮解密下。（对于我这样的菜鸟，这份手札的亮点就不是炫酷的图，而是“dplyr”函数，简直从零开始被教育）

1、前置数据处理

开局拿到一个包含单细胞分群和marker的sce.markers（data.frame）：

如果你没有自己的单细胞数据分析结果，可以参考下面的代码，一个简单的测试数据 ：

library(Seurat)
data("pbmc_small")
table(pbmc_small$RNA_snn_res.1)
pbmc.markers <- FindAllMarkers(pbmc_small,
                               only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
sce.markers=pbmc.markers
table(sce.markers$cluster)
length(unique(sce.markers$gene))


我使用的是我自己的实战项目的单细胞分群和marker的sce.markers（data.frame），如下所示：

数据的结构很简单，关键信息也只有两列，即“cluster”和“gene”。其中包含1673个基因，分属于0~16这些cluster（同一个基因可以属于不同的cluster）。

library(plyr)
library(dplyr)
library(ggalluvial)
library(viridis)

table(sce.markers$cluster)
length(unique(sce.markers$gene))


需要简单加工下数据，用管道符把dplyr函数串起来，大意就是：赋值data_gga %>% 挑出“cluster”和“gene”俩列 %>% 按“gene”分组（变化不能直接看出来，能便利后续操作） %>% 给每个gene计数并添加列 “n” %>% 挑选 n>=6 的基因 %>% 按 n 给数据框排序。

## filter gene count > 6
data_gga <- sce.markers %>%
dplyr::select(cluster,gene) %>%
group_by(gene) %>%
add_count(gene) %>%
filter(n >= 6) %>%
arrange(n)


得到的data_gga（grouped_df）结构也比较简单，后续绘图用的就是他（不得不感叹自己学艺不精，从这段代码学了很多dplyr函数，我处理数据的速度又要加快了）。

2、ggalluvial自定义桑基图

当前有很多优秀的包可以绘制桑基图，操作亦友好，但基于“授人以渔”的理念，健明老师提供了一段充满魔幻色彩的ggplot代码，给我的感受就是：

很炫酷。。。炫酷。。。酷。。。。。。。。。~炫目~ 还头晕~，先感受下：

首先自定义一个数据处理函数sankeyData，即利用“ggalluvial”包的“to_lodes_form”函数将输入数据转变为用于绘图的长数据。

sankeyData <- function(data) {
data <- data %>%
dplyr::select(cluster, gene, n)
res <- to_lodes_form(
data,
key = "Class",
axes = 1:2)
return(res)
}

sankey_data <- sankeyData(data_gga)


用函数处理data_gga，得到的sankey_data（tbl_df）结构也比较清晰：“n”仍然是计数，“alluvium”代表gene和cluster的两两配对关系，“class”类别，“stratum”是gene和cluster的值。

下面是绘图代码，用诺兰的话说：Don’t try to understand it, feel it. 大意就是：这样，这样，然后再这样，图就画好了（开个玩笑，不理解怎么能画得好呢）。绘图的主体是“ggalluvial”包的geom_stratum和geom_flow函数，填色和文字标签比较讲究，参数无须过多解释，自定义空间太足了。

p <- ggplot(data = sankey_data,
aes(x = Class, y = n,
stratum = stratum, alluvium = alluvium)) +
geom_flow(aes(fill = stratum),
curve_type = "xspline", show.legend = F,aes.flow = "forward",alpha = 0.7) +
geom_stratum(aes(fill = as.character(stratum)),
stat = "stratum",width = 1/4,color = NA) +
# geom_stratum(aes(fill = ifelse(as.numeric(Class) == 2, as.character(stratum), NA)),
# stat = "stratum",width = 1/4,color = NA) +
# geom_text(aes(label = ifelse(as.numeric(Class) == 1, as.character(stratum), NA)),stat = "stratum", size = 3) +
scale_fill_manual(values = c(viridis::viridis(n),rainbow(ng)), name = "Cluster", na.translate = F) +
geom_text(aes(label = as.character(stratum)),stat = "stratum", size = 3) +
scale_x_discrete(limits = c("cluster","gene"), expand = c(.05, .05)) +
labs(x = NULL, y = "Freq")+
theme(legend.position="none",
axis.ticks = element_line(linetype = "blank"),
panel.grid.minor = element_line(linetype = "blank"),
axis.title = element_text(family = "serif",
size = 13, face = "bold", colour = "chocolate4"),
axis.text = element_text(family = "serif"),
axis.text.x = element_text(size = 12,
colour = "black"),
panel.background = element_rect(fill = NA),
plot.background = element_rect(fill = "aliceblue"))
p


用“Cairo”做渲染并保存图片（我第一次见Cairo，以前只会ggsave，看了下介绍好像除色泽外差别不大）

library(Cairo)
Cairo(8000, 8000, file="plotggalluvial.pdf",
type="pdf", bg="white",dpi = 300,units = "px")
p
dev.off()


3、sankeywheel

用“sankeywheel”包会比较简便，只需要将data_gga的列名“cluster”、“gene”和“n”分别改成“from”、“to”和“weight”。绘图时“type”有环形“dependencywheel”和经典桑基“sankey”两种选择，可以调整theme为喜欢的风格（没有我喜欢的），其他细节需要大家自己探索。

library(sankeywheel)
library(dplyr)

col_names <- c("from", "to", "weight")
sankeydf <- data_gga %>% dplyr::select(cluster, gene, n) # 挑选需要的3列，前面已挑过的可忽略这句
colnames(sankeydf) <- col_names

wheel <- sankeywheel(from = sankeydf$from,
to = sankeydf$to,
weight = sankeydf$weight,
type = "dependencywheel",
width = "100%",
theme = "sandsignika")
wheel

sankey <- sankeywheel(from = sankeydf$from,
to = sankeydf$to,
weight = sankeydf$weight,
type = "sankey",
width = "100%",
theme = "gridlight")
sankey


此外比较方便的是可以生成交互式的HTML文件，便于分享和查看细节。

htmlwidgets::saveWidget(widget = wheel, file = "plotsankeywheel_wheel.html")

htmlwidgets::saveWidget(widget = sankey, file = "plotsankeywheel_sankey.html")


4、highcharter

类似的还有“highcharter”包，highchart函数可将图形渲染为可交互的HTML文件（这个包可绘制的不止桑基图，其核心是交互，比较适合需动态展示的场景）。

library(highcharter)

plot <- highchart() %>%
hc_title(text = "桑基图") %>%
hc_add_series(data = sankeydf,type = "sankey",hcaes(from = from,to = to,weight = weight)) %>%
hc_add_theme(hc_theme_google())

htmlwidgets::saveWidget(widget = plot, file = "plotHighcharter.html")
webshot::webshot(url = "plotHighcharter.html",
file = "plotHighcharter.pdf")


5、circlize

“circlize”包大家比较熟悉，毕竟是绘制弦图的优等生。需要注意的是需用reshape2::dcast把长矩阵sankeydf变回宽矩阵（mat），并改为matrix类型。

library(circlize) # 加载包
library(textshape)

mat <- reshape2::dcast(sankeydf,from~to)

mat[is.na(mat)] <- 0
mat <- column_to_rownames(mat,1)
mat <- as.matrix(mat)


需要熟悉绘图函数“chordDiagram”中那些命名和文档解释都比较晦涩的参数，此外，颜色和标签的自定义相对繁琐，需要通过“circos.track”函数附加内容（有点看不懂）。于我而言，circlize不如ggplot2易于理解和操作，但为了一张漂亮的图，这点困难也不算什么，啥问题是熬夜不能解决的呢？

grid.col <- setNames(c(c('#8dd3c7','#ffffb3','#bebada','#fb8072','#80b1d3','#fdb462', '#b3de69','#fccde5','#d9d9d9','#bc80bd','#ccebc5','#ffed6f'),
viridis::viridis(n-12),
rainbow(ng)), union(rownames(mat), colnames(mat)))
par(mar = c(0, 0, 0, 0))

# original image

chordDiagram(mat, grid.col = grid.col, annotationTrack = "grid",
link.sort = TRUE, link.decreasing = TRUE, directional = -1,
transparency = 0, link.zindex = rank(mat),
preAllocateTracks = list(track.height = max(strwidth(unlist(dimnames(mat))))))

# we go back to the first track and customize sector labels
circos.track(track.index = 1, panel.fun = function(x, y) {
circos.text(CELL_META$xcenter, CELL_META$ylim[1], CELL_META$sector.index,
facing = "clockwise", niceFacing = TRUE, adj = c(0, 0.5))
}, bg.border = NA) # here set bg.border to NA is important

circos.clear()


用“Cairo”渲染和保存文件

# save in pdf -------------------------------------------------------------
library(Cairo)
Cairo(8000, 8000, file="plotCirclize.pdf", type="pdf", bg="white",dpi = 300,units = "px")
chordDiagram(mat, grid.col = grid.col, annotationTrack = "grid",
link.sort = TRUE, link.decreasing = TRUE, directional = -1,
transparency = 0, link.zindex = rank(mat),
preAllocateTracks = list(track.height = max(strwidth(unlist(dimnames(mat))))))
# we go back to the first track and customize sector labels
circos.track(track.index = 1, panel.fun = function(x, y) {
circos.text(CELL_META$xcenter, CELL_META$ylim[1], CELL_META$sector.index,
facing = "clockwise", niceFacing = TRUE, adj = c(0, 0.5))
}, bg.border = NA) # here set bg.border to NA is important

circos.clear()
dev.off()


还感叹下，在解密这份手札的过程中，除熟悉绘制桑基图，我最大的收获竟然是学习dplyr函数。再次印证了健明老师那句话：基础不牢，地动山摇。想在生信路上走得远，打好基础很重要，须多动手、多思考，忌眼高手低。感谢健明老师和生信技能树提供的交流平台，及丰富而成体系的教程资料。

文末友情推荐

做教学我们是认真的，如果你对我们的马拉松授课（直播一个月互动教学）有疑问，可以看完我们从2000多个提问互动交流里面精选的200个问答！ 2021第二期_生信入门班_微信群答疑整理，以及 2021第二期_数据挖掘班_微信群答疑笔记

与十万人一起学生信，你值得拥有下面的学习班：

• 生信入门课-2021第6期

• 数据挖掘（GEO,TCGA,单细胞）2021第7期