R包circlize绘制弦状图示例

2020/07/06 23:33
阅读数 5K
R包circlize绘制样本-物种丰度关联弦状图

鲤小白:当初学习circlize包,最大的原因是白鱼同学一直以来都很讨厌perl……

路人:等等,这俩有什么关系呢?

鲤小白:因为要画圈图啊,但是circosperl程序,不想用,所以得找个替代品,正好看到了circlize

路人:好吧你赢了……

白鱼同学和perl之间的梗,哪天有兴趣了可以当笑话给大家分享,足以写满一页纸的笑话了。不过这里得回过话题来,circlize包绘制弦状图。关于circlize的基础本篇不涉及,还请自行了解,本篇以一个也不算很复杂的弦状图为例,展示一下常规做法,circlize绘制的弦状图还是很好看滴。

circlize文档:https://jokergoo.github.io/circlize_book/book/

 

本篇所用的作图示例文件、R脚本等的百度盘链接:

https://pan.baidu.com/s/10BCdWESl6Xt7DliUyT6drw

网盘附件中的数据为某16S高通量测序所得的细菌OTU丰度、物种分类等数据。本篇将使用这些数据绘制弦状图,用以展示20种重要的OTU6个样本中的丰度关系。20OTU可划分为5个门(物种分类单元,界门纲目科属种的“门”),6个样本可划分为2个分组。


左侧为circlize做图结果,分为5圈:

第一圈,各OTU的门水平分类以及各样本的分组信息;

第二圈,OTU相对丰度的百分比信息;

第三圈,OTU及样本主区块,以不同颜色和标签区分,区块外周的刻度为OTU的绝对丰度信息;

第四圈,OTU及样本副区块,与主区块(第三圈)对应,展示了各OTU在各样本中的丰度,以及各样本所含每种OTU的丰度信息;

第五圈,与OTU及样本副区块(第四圈)相对应,连线展示OTU、样本关联信息。

右侧为图例,展示了各OTU(本示例中共计20OTU)的详细物种分类详情。

 

准备文件



网盘附件3个文件。

1)物种丰度表格文件(otu_table.txt)。每一列为样本,每一行为物种(此处为OTU),交叉区为个OTU在各样本中的丰度。


2)样本分组信息文件(group.txt)。第一列为样本ID,第二列为样本分组信息。


3)物种分类信息文件(taxonomy.txt)。第一列为OTU ID,第二列为OTU的“门水平”分类,第三列为OTU的详细物种分类。


 

circlize作图



预处理数据


首先根据OTU分类文件(taxonomy.txt)以及样本分组文件(group.txt)中的OTU、样本顺序,对OTU丰度表(otu_table.txt)中的OTU、样本进行排序,以及获取OTU和样本的名称。

#读取 taxonomy.txt 的内容,获取“OTU/分类”排序,OTU 名称
taxonomy <- read.delim('taxonomy.txt', sep = '\t', stringsAsFactors = FALSE)
tax_phylum <- unique(taxonomy$phylum)
taxonomy$phylum <- factor(taxonomy$phylum, levels = tax_phylum)
all_otu <- taxonomy$OTU_ID
taxonomy$OTU_ID <- factor(taxonomy$OTU_ID, levels = all_otu)

#读取 group.txt 的内容,获取“样本/分组”排序,样本名称
group <- read.delim('group.txt', sep = '\t', stringsAsFactors = FALSE)
all_group <- unique(group$group_ID)
group$group_ID <- factor(group$group_ID, levels = all_group)
all_sample <- group$sample_ID

#读取 otu_table.txt,排序 OTU 和样本
otu_table <- read.delim('otu_table.txt', sep = '\t')
otu_table <- merge(taxonomy, otu_table, by = 'OTU_ID')
otu_table <- otu_table[order(otu_table$phylum, otu_table$OTU_ID), ]
rownames(otu_table) <- otu_table$OTU_ID
otu_table <- otu_table[all_sample]

 

然后作些小统计,各OTU的总丰度,各样本中所有OTU的总丰度,各OTU在各样本中的丰度,以及各样本所含各OTU的丰度。顺便给各OTU及样本定义颜色。这些统计好的结果将直接用于circlize作图。

library(reshape2)

##生成作图数据
#circlize 外圈属性数据
all_ID <- c(all_otu, all_sample)
accum_otu <- rowSums(otu_table)
accum_sample <- colSums(otu_table)
all_ID_xlim <- cbind(rep(0, length(all_ID)),data.frame(c(accum_otu, accum_sample)))

#circlize 内圈连线数据
otu_table$otu_ID <- all_otu
plot_data <- melt(otu_table, id = 'otu_ID') #此处使用了reshape2包中的melt()命令
colnames(plot_data)[2] <- 'sample_ID'
plot_data$otu_ID <- factor(plot_data$otu_ID, levels = all_otu)
plot_data$sample_ID <- factor(plot_data$sample_ID, levels = all_sample)
plot_data <- plot_data[order(plot_data$otu_ID, plot_data$sample_ID), ]
plot_data <- plot_data[c(2, 1, 3, 3)]

#颜色设置
color_otu <- c('#8DD3C7', '#FFFFB3', '#BEBADA', '#FB8072', '#80B1D3', '#FDB462', '#B3DE69', '#FCCDE5', '#BC80BD', '#CCEBC5', '#FFED6F', '#E41A1C', '#377EB8', '#4DAF4A', '#984EA3', '#FF7F00', '#FFFF33', '#A65628', '#F781BF', '#66C2A5')
color_sample <- c('#6181BD', '#F34800', '#64A10E', '#FF00FF', '#c7475b', '#049a0b')
color_phylum <- c('#BEAED4', '#FDC086', '#FFFF99', '#386CB0', '#F0027F')
color_group <- c('#4253ff', '#ff4308')

names(color_otu) <- all_otu
names(color_sample) <- all_sample

 

定义整体布局


然后接下来就是使用circlize包绘制弦状图了。

注:绘图细节部分的调整很多也很繁琐,此处默认使用的各细节参数设置(如字体大小等)均根据示例文件而来。若是大家后续更换为自己的数据时,还需多加调试参数了。


首先定义整体的布局,这一步生成空白画板。

library(circlize)

#整体画板布局
gap_size <- c(rep(3, length(all_otu) - 1), 6, rep(3, length(all_sample) - 1), 6)
circos.par(cell.padding = c(0, 0, 0, 0), start.degree = 270, gap.degree = gap_size)
circos.initialize(factors = factor(all_ID, levels = all_ID), xlim = all_ID_xlim)

gap_size将用于设定circlize图中,各OTU或样本区块之间的间距(如下图所示);circos.par就是画板设置了;circos.initialize定义绘图因子,此处使用了前述得到的circlize外圈属性数据框all_ID_xlim

 

第一圈,OTU分类、样本分组区块


然后开始绘图。Rcirclize也是由外圈往内圈逐个添加。

#绘制 OTU 分类、样本分组区块(第一圈)
circos.trackPlotRegion(
ylim = c(0, 1), track.height = 0.03, bg.border = NA,
panel.fun = function(x, y) {
sector.index = get.cell.meta.data('sector.index')
xlim = get.cell.meta.data('xlim')
ylim = get.cell.meta.data('ylim')
} )

for (i in 1:length(tax_phylum)) {
tax_OTU <- {subset(taxonomy, phylum == tax_phylum[i])}$OTU_ID
highlight.sector(tax_OTU, track.index = 1, col = color_phylum[i], text = tax_phylum[i], cex = 0.5, text.col = 'black', niceFacing = FALSE)
}

for (i in 1:length(all_group)) {
group_sample <- {subset(group, group_ID == all_group[i])}$sample_ID
highlight.sector(group_sample, track.index = 1, col = color_group[i], text = all_group[i], cex = 0.7, text.col = 'black', niceFacing = FALSE)
}

第一条命令用于绘制第一圈。ylim控制y轴范围,track.height控制高度,bg.border控制边框颜色;定义的函数panel.fun()以循环的方式,读取all_ID_xlim中的内容,逐一绘制外圈区块。运行完上述的第一条命令之后,发现绘图区域啥也没显示。但事实上,并非没有绘制任何元素,只是将绘制元素的填充颜色及边框颜色设置为NA了,因为此处不需要任何颜色作为添加。

然后接下来使用两个循环,将预先设定的OTU门分类颜色以及样本分组颜色添加在刚才未定义颜色的第一圈。highlight.sector()实现这个功能,并可以将处于同一分组的样本合并在一起。track.index = 1指定了将颜色添加在circlize图的第一圈;coltext参数分别指定颜色和标签内容;cex指定标签字体大小,text.col指定标签字体颜色,niceFacing用于展示标签字体的方向,默认为TRUE但个人推荐使用FALSE。更多参数信息可?highlight.sector


 

第二圈,OTU丰度的百分比信息


继续往内圈绘制。

#添加 OTU 百分比注释(第二圈)
circos.trackPlotRegion(
ylim = c(0, 1), track.height = 0.05, bg.border = NA,
panel.fun = function(x, y) {
sector.index = get.cell.meta.data('sector.index')
xlim = get.cell.meta.data('xlim')
ylim = get.cell.meta.data('ylim')
} )

circos.track(
track.index = 2, bg.border = NA,
panel.fun = function(x, y) {
xlim = get.cell.meta.data('xlim')
ylim = get.cell.meta.data('ylim')
sector.name = get.cell.meta.data('sector.index')
xplot = get.cell.meta.data('xplot')

by = ifelse(abs(xplot[2] - xplot[1]) > 30, 0.25, 1)
for (p in c(0, seq(by, 1, by = by))) circos.text(p*(xlim[2] - xlim[1]) + xlim[1], mean(ylim) + 0.4, paste0(p*100, '%'), cex = 0.4, adj = c(0.5, 0), niceFacing = FALSE)

circos.lines(xlim, c(mean(ylim), mean(ylim)), lty = 3)
} )

首先我们使用circos.trackPlotRegion()绘制一圈无边框、无颜色的空白区域,目的为添加的刻度标签等内容留出展示的空间。

然后使用circos.track()添加百分比刻度标签。其中,track.index = 2意为将刻度标签展示在circlize图的第二圈;定义的函数panel.fun()使用循环的方式,读取all_ID_xlim中的内容(丰度信息)并在读取后进行计算得到百分比数值,同时额外使用if判断OTU丰度的绝对数值,若数值足够大则考虑将展示更多的刻度(以25%为一刻度),若数值过小则仅展示0100%这两个刻度;circos.lines()添加了外周的点状虚线。


  

第三、四圈, OTU及样本名称区块


其中第三圈为主区块,第四圈为副区块。

#绘制 OTU、样本主区块(第三圈)
circos.trackPlotRegion(
ylim = c(0, 1), track.height = 0.03, bg.col = c(color_otu, color_sample), bg.border = NA, track.margin = c(0, 0.01),
panel.fun = function(x, y) {
xlim = get.cell.meta.data('xlim')
sector.name = get.cell.meta.data('sector.index')
circos.axis(h = 'top', labels.cex = 0.4, major.tick.percentage = 0.4, labels.niceFacing = FALSE)
circos.text(mean(xlim), 0.2, sector.name, cex = 0.4, niceFacing = FALSE, adj = c(0.5, 0))
} )

#绘制 OTU、样本副区块(第四圈)
circos.trackPlotRegion(ylim = c(0, 1), track.height = 0.03, track.margin = c(0, 0.01))

继续使用circos.trackPlotRegion()循环读取all_ID_xlim中的内容并绘制,此时可以将预先定义的OTU颜色及样本颜色添加上。同时在内部使用circos.axis()circos.text()函数,在循环绘图的同时,将刻度线和标签添加上。此处,sector.name为在循环过程中每次读取的OTU或样本名称;circos.text()中的mean(xlim)0.2分别意为将给定的标签文字添加在每个区块的X轴中间位置以及Y=0.2的位置。其余参数项类似上述不再多说,更细致的参数可参阅帮助。

第三圈主圈绘制完成,第四圈副圈暂未添加颜色(将在后续操作中添加)。


   

最内圈,连线区


最内圈为样本-OTU丰度关联信息,以连线的方式展示。

#绘制 OTU-样本关联连线(最内圈)
for (i in seq_len(nrow(plot_data))) {
circos.link(
plot_data[i,2], c(accum_otu[plot_data[i,2]], accum_otu[plot_data[i,2]] - plot_data[i,4]),
plot_data[i,1], c(accum_sample[plot_data[i,1]], accum_sample[plot_data[i,1]] - plot_data[i,3]),
col = paste0(color_otu[plot_data[i,2]], '70'), border = NA )

circos.rect(accum_otu[plot_data[i,2]], 0, accum_otu[plot_data[i,2]] - plot_data[i,4], 1, sector.index = plot_data[i,2], col = color_sample[plot_data[i,1]], border = NA)
circos.rect(accum_sample[plot_data[i,1]], 0, accum_sample[plot_data[i,1]] - plot_data[i,3], 1, sector.index = plot_data[i,1], col = color_otu[plot_data[i,2]], border = NA)

accum_otu[plot_data[i,2]] = accum_otu[plot_data[i,2]] - plot_data[i,4]
accum_sample[plot_data[i,1]] = accum_sample[plot_data[i,1]] - plot_data[i,3]
}

此处以循环的方式添加最内圈连线,连线颜色以OTU的预设颜色为准(绘制的同时将颜色的不透明度设置为70%)。先前绘制外圈第四圈的时候未添加颜色,也在此处将颜色添加上,即命令中使用的两个circos.rect()函数分别添加。

运行完毕后,我们的circlize初图就完成了。


 

可选,添加图例


其实到了这一步就已经完成了,不过示例图中还包含了图例信息。circlize包中是没有绘制图例的函数的,因此还需额外调用其它的包绘制图例。

circlize的官方说明文档中,提到可通过ComplexHeatmapLegend()函数,为circlize图添加图例。


备注:为了更好的体验,请使用3.6版本及以上的R。经过多次踩坑,发现像3.33.4等版本的R中,ComplexHeatmap包的很多功能无法适用,函数残缺且不兼容等问题…….

 

library(ComplexHeatmap) #可用此包添加图例
library(grid) #可用此包调整画板

#添加图例
otu_legend <- Legend(
at = all_otu, labels = taxonomy$detail, labels_gp = gpar(fontsize = 8),
grid_height = unit(0.5, 'cm'), grid_width = unit(0.5, 'cm'), type = 'points', pch = NA, background = color_otu)

pushViewport(viewport(x = 0.85, y = 0.5))
grid.draw(otu_legend)
upViewport()

#推荐作图完成后,清除目前的 circlize 样式,以免影响继续作图
circos.clear()

图例添加完毕后,最终结果如下所示。


 

其他


如果图片比例失调,除了更改作图参数外,还可通过在本地输出为pdf/png时设置保存图形的长宽比列来解决。

#比方说先打开本地 pdf,设置合适的长宽比
pdf('circlize_plot.pdf', width = 20, height = 8)
circle_size = unit(1, 'snpc')

#然后作图
#整体布局
gap_size <- c(rep(3, length(all_otu) - 1), 6, rep(3, length(all_sample) - 1), 6)
circos.par(cell.padding = c(0, 0, 0, 0), start.degree = 270, gap.degree = gap_size)
circos.initialize(factors = factor(all_ID, levels = all_ID), xlim = all_ID_xlim)

#中间那些作图步骤(略)
……

##添加图例
otu_legend <- Legend(
at = all_otu, labels = taxonomy$detail, labels_gp = gpar(fontsize = 8),
grid_height = unit(0.5, 'cm'), grid_width = unit(0.5, 'cm'), type = 'points', pch = NA, background = color_otu)

pushViewport(viewport(x = 0.85, y = 0.5))
grid.draw(otu_legend)
upViewport()

##清除 circlize 样式并关闭画板/pdf 文件
circos.clear()
dev.off()

 

 

链接

R语言绘制曼哈顿图示例

R语言绘制差异火山图示例

R语言绘制三元图(Ternary Plot)示例

突破韦恩图数量限制,R包UpSetR集合可视化

R语言绘制箱线图

R语言绘制双向柱状图

R语言绘制堆叠面积图

R语言绘制堆叠柱状图

R语言绘制星形图示例

R语言绘制饼图(扇形图)

R语言绘制花瓣图示例

R语言绘制韦恩图的几种方法示例



本文分享自微信公众号 - 小明的数据分析笔记本()。
如有侵权,请联系 support@oschina.cn 删除。
本文参与“OSC源创计划”,欢迎正在阅读的你也加入,一起分享。

展开阅读全文
加载中
点击引领话题📣 发布并加入讨论🔥
打赏
0 评论
1 收藏
0
分享
返回顶部
顶部