专栏首页小白鱼的生统笔记R包circlize绘制弦状图示例

R包circlize绘制弦状图示例

R包circlize绘制样本-物种丰度关联弦状图

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

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

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

路人:好吧你赢了……

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

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

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

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

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

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

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

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

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

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

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

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

准备文件

网盘附件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分类、样本分组区块

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

#绘制 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图的第一圈;col和text参数分别指定颜色和标签内容;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%为一刻度),若数值过小则仅展示0和100%这两个刻度;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的官方说明文档中,提到可通过ComplexHeatmap包Legend()函数,为circlize图添加图例。

备注:为了更好的体验,请使用3.6版本及以上的R。经过多次踩坑,发现像3.3、3.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()

本文分享自微信公众号 - 小白鱼的生统笔记(gh_5f751e893315),作者:生信小白鱼 鲤小白

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

原始发表时间:2019-11-07

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

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 关于这种“网络模块”和“模块饼图”的可视化方法

    知道好多同学都很好奇这种网络图,特别是右边这个是怎么得到的,本篇就作个简单操作演示。

    用户7585161
  • ​R语言计算群落多样性分析中常见Alpha多样性指数

    本篇以某16S高通量测序所得细菌群落数据为例,简介使用R计算几种在微生物群落分析中常用的Alpha多样性指数。

    用户7585161
  • R包randomForest的随机森林回归模型以及对重要变量的选择

    关于随机森林(random forest),前文“随机森林分类以及对重要变量的选择”中已经对其基本原理作了简单概括。在前文中,响应变量是一组类别变量(代表了样本...

    用户7585161
  • 关于too many files open的解决办法

    当我们用一些大的测试程序时,有时可能会报错,too many files open之类的错误,系统默认的同时打开文件数是1024,可以用这个命令查看: #uli...

    三杯水Plus
  • 微信网页开发页面上滑效果

    我记得在之前我写过两篇关于微信网页开发上滑效果的文章,在那两篇文章中滑动是全部页面都滑动,但是会使页面的机动性变差,如果说我这个页面想滑动,但是那个页面又不想滑...

    无邪Z
  • 100行代码爬取全国所有必胜客餐厅信息

    极客猴,热衷于 Python,目前擅长利用 Python 制作网络爬虫以及 Django 框架。

    Python中文社区
  • Spring Boot---(2)SpringBoot多环境配置和使用

    我们在开发Spring Boot应用时,通常同一套程序会被应用和安装到几个不同的环境,比如:开发、测试、生产等。其中每个环境的数据库地址、服务器端口等等配置都会...

    IT云清
  • python笔记19-yaml文件写入(ruamel.yaml)

    yaml作为配置文件是非常友好的一种格式,前面一篇讲了yaml的一些基础语法和读取方法,本篇继续讲yaml文件写入方法 用yaml模块写入字典嵌套字典这种复杂的...

    上海-悠悠
  • SDS/超融合哪家强?来看看WhatMatrix怎么说!

    近期,市场分析和咨询公司WhatMatrix出了一份软件定义存储/超融合(SDS/HCI)的研究报告。

    大数据在线
  • SpringBoot 笔记(十一):Servlet容器

    lwen

扫码关注云+社区

领取腾讯云代金券