今天继续进行基于barplot()命令的文章条形图仿制

原图展示

之前看过一篇利用肠道菌群进行疾病诊断的文章，其中涉及到一个图形，作者将筛选出的marker微生物的总丰度与肠道菌群的alpha多样性指数在同一个图像利用对称的条形图进行展示，同时给出了健康和疾病样品间的差异检验结果。

今天就来仿照绘制一下这幅图，我绘制出来的最终图像是下面这个样子。

还是照例先给出完成的绘图代码：

marker_abundance <- read.table("marker_abundance.txt",header = TRUE,row.names = 1,sep = "\t")diversity <- read.table("index.txt",header = TRUE,row.names = 1,sep = "\t")
library(ggplot2)
tiff(filename = "marker_abundance_diversity.tif",width = 5400,height = 7200,compression = "lzw",res = 600,type = "windows")par(mfrow = c(2,1))par(mar = c(5,7,1,1))par(xpd = TRUE)bar1 <- barplot(marker_abundance$Abundance,ylim = c(0,0.9),axes = FALSE,col = "grey80")axis(side = 2,at = c(0,0.2,0.4,0.6,0.8),labels = c("","","","",""),lwd = 3)text(x = 0.07,y = c(0,0.4,0.8),labels = c(0,0.4,0.8),cex = 2.5,adj = c(1,0.5))mtext(side = 2,text = "Abundance sum",line = 4,cex = 3)arrows(x0 = 0.7,y0 = 0.4478635,x1 = 0.7,y1 = 0.4478635+0.05248740,col = "black",angle = 90,length = 1,lwd = 3)arrows(x0 = 1.9,y0 = 0.4721783,x1 = 1.9,y1 = 0.4721783+0.07471587,col = "black",angle = 90,length = 1,lwd = 3)text(x = bar1,y = -0.12,labels = c("Disease","Healthy"),cex = 3.5)segments(x0 = 0.7,y0 = 0.55,x1 = 0.7,y1 = 0.75,col = "black",lwd = 3)segments(x0 = 0.7,y0 = 0.75,x1 = 1.9,y1 = 0.75,col = "black",lwd = 3)segments(x0 = 1.9,y0 = 0.75,x1 = 1.9,y1 = 0.58,col = "black",lwd = 3)text(x = 1.3,y = 0.8,labels = "NS",cex = 2.5)bar2 <- barplot(diversity$shannon,ylim = c(-12,0),axes = FALSE,col = "grey80")axis(side = 2,at = c(-12,-10,-8,-6,-4,-2,0),labels = c("","","","","","",""),lwd = 3)text(x = 0.07,y = c(-12,-8,-4,0),labels = c(-12,-8,-4,0),cex = 2.5,adj = c(1,0.5))mtext(side = 2,text = "Within-sample diversity",line = 4,cex = 3)arrows(x0 = 0.7,y0 = -8.946364,x1 = 0.7,y1 = -8.946364-0.5686171,col = "black",angle = 90,length = 1,lwd = 3)arrows(x0 = 1.9,y0 = -8.716364,x1 = 1.9,y1 = -8.716364-0.6158940,col = "black",angle = 90,length = 1,lwd = 3)segments(x0 = 0.7,y0 = -10,x1 = 0.7,y1 = -12,col = "black",lwd = 3)segments(x0 = 0.7,y0 = -12,x1 = 1.9,y1 = -12,col = "black",lwd = 3)segments(x0 = 1.9,y0 = -12,x1 = 1.9,y1 = -9.8,col = "black",lwd = 3)text(x = 1.3,y = -13,labels = "NS",cex = 2.5)dev.off()

绘制这个图像需要两个数据，一个是marker微生物总丰度，另一个是样品对应的alpha多样性指数，需要在绘图前先使用Excel分别计算实验和对照组中marker微生物总丰度和alpha多样性指数的平均值和标准偏差(sd)，具体的绘图数据如下：

特别需要注意的是：用于绘制下方条形图的数据要手动改为负值，sd不需要改为负值。

绘图过程详解

定义绘图区域

首先要制定用于绘图的区域。

par(mfrow = c(2,1))

将绘图区定义为两行一列，分别用于绘制两个互相对称的条形图。

par(mar = c(5,7,1,1))

指定第一个绘图区即上方条形图的绘图边界，注意下方和左侧要留出大一点的区域，用于样品组的名称和纵坐标，而上方和右侧的边界要尽量小。

第二个绘图区即下方条形图的绘图边界无需设置。

par(xpd = TRUE)

将绘图参数定义为可以在绘图区外显示，从而使得绘图参数的显示能够跨越上下两个绘图区。

上方条形图绘制

定义完绘图区之后，要首先进行上方条形图的绘制。

bar1 <- barplot(marker_abundance$Abundance,ylim = c(0,0.9),axes = FALSE,col = "grey80")

注意：y轴的范围要适当的大一些，为差异性检验结果的展示预留空间。

使用“axes = FALSE”隐藏默认的坐标轴信息。

axis(side = 2,at = c(0,0.2,0.4,0.6,0.8),labels = c("","","","",""),lwd = 3)

添加纵坐标，并将刻度标签设为空。

text(x = 0.07,y = c(0,0.4,0.8),labels = c(0,0.4,0.8),cex = 2.5,adj = c(1,0.5))

通过x的数值，调整标签与坐标轴刻度的间距。

添加纵坐标刻度标签，通过标签的数量和间隔实现坐标轴的主刻度和副刻度。

使用adj参数调整标签的对齐方式，使得长度不一的标签文本统一靠右居中对齐。

mtext(side = 2,text = "Abundance sum",line = 4,cex = 3)

添加纵坐标轴的标签，使用line参数调整标签与坐标轴的距离，cex调整标签字号的大小。

条形误差棒的添加

R默认的绘图命令中并不包含误差棒的绘制，这里使用arrows()命令以箭头的形式添加误差棒。

arrows(x0 = 0.7,y0 = 0.4478635,x1 = 0.7,y1 = 0.4478635+0.05248740,col = "black",angle = 90,length = 1,lwd = 3)
arrows(x0 = 1.9,y0 = 0.4721783,x1 = 1.9,y1 = 0.4721783+0.07471587,col = "black",angle = 90,length = 1,lwd = 3)

箭头的x0位对应条形在x轴的位置，因为是竖直方向的，因而x1与x0一致。

y0为对应条形的绘图数值，y1为绘图数值+sd。

angle设置为90度，使得箭头的两边处于一条水平的直线，从而形成误差棒的形式。

使用length和lwd调整误差棒的长度和线宽。

差异性检验结果展示

使用线段的形式构成一个U型框连接两个条形，之后利用文本给出组间差异性检验结果。

segments(x0 = 0.7,y0 = 0.55,x1 = 0.7,y1 = 0.75,col = "black",lwd = 3)segments(x0 = 0.7,y0 = 0.75,x1 = 1.9,y1 = 0.75,col = "black",lwd = 3)segments(x0 = 1.9,y0 = 0.75,x1 = 1.9,y1 = 0.58,col = "black",lwd = 3)text(x = 1.3,y = 0.8,labels = "NS",cex = 2.5)

具体的参数数值根据绘图数据进行调整。

此时需要添加条形对应的样本组名称，因为如果在第二个条形绘制后再添加可能需要额为设置绘图区参数。

text(x = bar1,y = -0.12,labels = c("Disease","Healthy"),cex = 3.5)

注意：在第二个条形绘制完成之后，再返回此处调整标签在y轴的位置，使得标签位于两个条形图的中间。

对称条形的绘制

对称条形的绘制方法与第一个条形图基本一致。

bar2 <- barplot(diversity$shannon,ylim = c(-12,0),axes = FALSE,col = "grey80")axis(side = 2,at = c(-12,-10,-8,-6,-4,-2,0),labels = c("","","","","","",""),lwd = 3)text(x = 0.07,y = c(-12,-8,-4,0),labels = c(12,8,4,0),cex = 2.5,adj = c(1,0.5))mtext(side = 2,text = "Within-sample diversity",line = 4,cex = 3)arrows(x0 = 0.7,y0 = -8.946364,x1 = 0.7,y1 = -8.946364-0.5686171,col = "black",angle = 90,length = 1,lwd = 3)arrows(x0 = 1.9,y0 = -8.716364,x1 = 1.9,y1 = -8.716364-0.6158940,col = "black",angle = 90,length = 1,lwd = 3)segments(x0 = 0.7,y0 = -10,x1 = 0.7,y1 = -12,col = "black",lwd = 3)segments(x0 = 0.7,y0 = -12,x1 = 1.9,y1 = -12,col = "black",lwd = 3)segments(x0 = 1.9,y0 = -12,x1 = 1.9,y1 = -9.8,col = "black",lwd = 3)text(x = 1.3,y = -13,labels = "NS",cex = 2.5)

需要注意的是：

1. 对称条形图因为是在下方展示，所以所有的数值均是负数；

2. 在绘制下方条形图纵坐标标签时，要将标签设置为对应的正数；

3. 由于此时误差棒的方向是向下，因而y1的数值为mean-sd。

原图展示

还是上一篇推文中介绍的对称条形图的文献，其中包含一个对疾病和健康marker微生物的丰度进行相关性分析，进而使用条形图展示其相关系数分布的图形。

今天就来仿照绘制一下这幅图，我绘制出来的最终图像是下面这个样子。

还是照例先给出完成的绘图代码：

disease <- read.table("Disease.txt",header = TRUE,row.names = 1,sep = "\t")healthy <- read.table("Healthy.txt",header = TRUE,row.names = 1,sep = "\t")disease <- t(disease)healthy <- t(healthy)
cor.d <- cor(disease,method = "spearman")cor.h <- cor(healthy,method = "spearman")write.table(cor.d,"correlation_disease.xls",sep="\t",quote=F,col.names=NA)write.table(cor.h,"correlation_healthy.xls",sep="\t",quote=F,col.names=NA)
library(reshape2)cor.d <- melt(cor.d)cor.h <- melt(cor.h)
d.value <- cor.d$valued.value <- as.vector(d.value)d.value <- as.numeric(d.value)d.value <- d.value[d.value > 0]h.value <- cor.h$valueh.value <- as.vector(h.value)h.value <- as.numeric(h.value)h.value <- h.value[h.value > 0]
hist1 <- hist(d.value,breaks = c(0,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1),plot = FALSE)hist2 <- hist(h.value,breaks = c(0,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1),plot = FALSE)
tiff(filename = "correlation_density.tif",width = 5400,height = 7200,compression = "lzw",res = 600,type = "windows")par(mfrow = c(1,2))par(mar = c(3,8,6,2))par(xpd = TRUE)bar1 <- barplot(hist1$density,space = 0,horiz = TRUE,axes = FALSE,col = "grey80")axis(side = 1,at = c(0,0.5,1.0,1.5,2.0),labels = c("","","","",""),lwd = 3)text(x = c(0,1.0,2.0),y = -1.5,labels = c(0,1.0,2.0),cex = 2.5)axis(side = 2,line = 1,at = c(0,4,8,12,16,20),labels = c("","","","","",""),lwd = 3)text(x = -0.28,y = c(0,4,8,12,16,20),labels = c("0.0","0.2","0.4","0.6","0.8","1.0"),cex = 2.5,adj = c(1,0.5))mtext(side = 2,text = "Occurrence rate (n = 22)",line = 5.8,cex = 3)text(x = 0.8,y = 21,labels = "Disease-enriched\nmarkers",cex = 2.5)par(mar = c(3,0.5,6,9.5))bar2 <- barplot(hist2$density,space = 0,horiz = TRUE,axes = FALSE,col = "grey80")axis(side = 1,at = c(0,0.5,1.0,1.5,2.0),labels = c("","","","",""),lwd = 3)text(x = c(0,1.0,2.0),y = -1.5,labels = c(0,1.0,2.0),cex = 2.5)text(x = 1.2,y = 21,labels = "Healthy-enriched\nmarkers",cex = 2.5)text(x = 1.5,y = 0.5,labels = "Density",cex = 3)dev.off()

绘制图像所需的数据就是用于进行相关性的marker微生物在各样品中的丰度数据。

绘图过程详解

相关性分析

首先将绘图数据导入。

disease <- read.table("Disease.txt",header = TRUE,row.names = 1,sep = "\t")
disease <- t(disease)

相关性分析要求行为样本，进行相关性统计的变量为列，因而有些数据可能需要转置处理。

cor.d <- cor(disease,method = "spearman")
write.table(cor.d,"correlation_disease.xls",sep="\t",quote=F,col.names=NA)

使用cor()命令进行相关性分析，其可选方法包括spearman、pearson和kendall，输出结果为相关系数统计矩阵。

使用write.table()将相关性分析结果输出为xls表格，方便日后阅读和使用。

⚠️ cor()只能输出相关系数，而不同输出显著性p值，如果需要进行显著性检验，需要使用其它命令进行操作，相关内容会在之后的推文中介绍。

相关系数频率分布统计

首先要将“宽格式”的相关系数数据转化为“长格式”的数据。

library(reshape2)cor.d <- melt(cor.d)

之后将其中相关系数数值列提出来变为向量。

d.value <- cor.d$valued.value <- as.vector(d.value)d.value <- as.numeric(d.value)d.value <- d.value[d.value > 0]

⚠️ 因为本次绘图只是统计正相关的数据，因而使用[d.value > 0]值保留大于0的相关系数。

之后统计相关系数在不同范围内的分布频数。

hist1 <- hist(d.value,breaks = c(0,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1),plot = FALSE)

使用hist()进行频数的统计，应用breaks定义统计的范围，治理以0.05为间隔统计相关系数的分布频数。

图像绘制

最后对统计的相关系数分布频数进行图像绘制。

par(mfrow = c(1,2))par(mar = c(3,8,6,2))par(xpd = TRUE)bar1 <- barplot(hist1$density,space = 0,horiz = TRUE,axes = FALSE,col = "grey80")axis(side = 1,at = c(0,0.5,1.0,1.5,2.0),labels = c("","","","",""),lwd = 3)text(x = c(0,1.0,2.0),y = -1.5,labels = c(0,1.0,2.0),cex = 2.5)axis(side = 2,line = 1,at = c(0,4,8,12,16,20),labels = c("","","","","",""),lwd = 3)text(x = -0.28,y = c(0,4,8,12,16,20),labels = c("0.0","0.2","0.4","0.6","0.8","1.0"),cex = 2.5,adj = c(1,0.5))mtext(side = 2,text = "Occurrence rate (n = 22)",line = 5.8,cex = 3)text(x = 0.8,y = 21,labels = "Disease-enriched\nmarkers",cex = 2.5)par(mar = c(3,0.5,6,9.5))bar2 <- barplot(hist2$density,space = 0,horiz = TRUE,axes = FALSE,col = "grey80")axis(side = 1,at = c(0,0.5,1.0,1.5,2.0),labels = c("","","","",""),lwd = 3)text(x = c(0,1.0,2.0),y = -1.5,labels = c(0,1.0,2.0),cex = 2.5)text(x = 1.2,y = 21,labels = "Healthy-enriched\nmarkers",cex = 2.5)text(x = 1.5,y = 0.5,labels = "Density",cex = 3)

图像绘制的命令与正常的条形图绘制基本一致，这里就不再赘述了，详细的参数用法可以参见之前的几篇推文。

⚠️ 因为图像的条形是横向分布，因而在绘制时要使用horiz = TRUE参数。

• 声明： 本号旨在传播、传递、交流，对相关文章内容观点保持中立态度。涉及内容如有侵权或其他问题，请与本号联系，第一时间做出撤回。

  推荐课程（点击进入联系详细了解）：

  基于Ｒ语言的非参数与半参数回归模型

  R语言回归及混合效应（多水平/层次/嵌套）模型及贝叶斯实现实践技术

  基于R语言与Octave的Kriging、贝叶斯优化及高斯过程模型优化

  基于Citespace和vosviewer的文献信息可视化分析