跟着Nat Commun学作图 | 3.物种丰度堆积柱状图

原创小小木舟木舟笔记

跟着Nat Commun学作图 | 3.物种丰度堆积柱状图

今天要学习的图来自2021年10月29号发表在的Nature Communication上的一篇文章，题目是【新冠肺炎患者呼吸道菌群组成及其与宿主相互作用的临床研究】。今天先来复现其中的一幅物种丰度堆积柱状图。

DOI:10.1038/s41467-021-26500-8

读图
示例数据及作图前准备

导入数据
数据处理

绘制

结果展示

数据及代码
后记
往期内容

读图

该图为分组的物种属水平的相对丰度堆积柱状图。其中只显示了前15丰度的属，剩下的属都归于others里

示例数据及作图前准备

由于作者给开源的数据设置了权限，就用手上的Phylum水平的绝对丰度矩阵来作为示例。

导入数据

# 导入数据并查看数据集格式
rm(list = ls())
setwd("F:\\~\\mzbj\\mzbj_note\\NC\\3.stackbar")
phylum <- read.table('count_2Phylum.txt', sep="\t", header=T, row.names=1)
head(phylum)

> head(phylum)
                             KO1   KO2   KO3   KO4  KO5  KO6  OE1   OE2   OE3  OE4
Acidobacteria                 52    69    39    45   50   59  196    90   140   55
Actinobacteria              7771 13581 10242 10116 6305 8058 8130 10355 10470 8324
Armatimonadetes                2     1     1     2    0    4    2     7     2    5
Bacteroidetes                845   821  2766  1059 1889  782  975  2408  2783 2250
Candidatus_Saccharibacteria    1     1     6     0   33    4    0    10     1    0
Chlamydiae                     4     2     0     0    2    9   10     9    11   16

数据处理

# 将绝对丰度转化为百分比形式的相对丰度
phylum_per <- as.data.frame(lapply(phylum, function(x) x / sum(x)))
row.names(phylum_per) <- row.names(phylum) #加一下行名
# 计算每个门水平的平均丰度 便于后续筛选                                 
phylum.ave <- apply(phylum_per, 1, FUN=mean)
phylum.2 <- cbind(phylum_per, phylum.ave)[order(-phylum.ave),] #排个序

# 选择丰度最高的10个门 剩下的放入others里
phylum.2 <- subset(phylum.2[1:10,], select=-phylum.ave)
# 统计others丰度
phylum.2 <- rbind(phylum.2, others=apply(phylum.2, 2, function(x){1-sum(x)}))
# 加一列行名 便于后续的长宽转换
phylum.2 <- cbind(PhylumID=row.names(phylum.2), phylum.2)

# 长宽转换
library(reshape2)
# 因子排个序
phylum.2$PhylumID <- factor(phylum.2$PhylumID, levels = rev(phylum.2$PhylumID))
phylum.gg <- melt(phylum.2, id.vars="PhylumID", variable.name="SampleID", value.name="Abundance")
head(phylum.gg)

> head(phylum.gg)
        PhylumID SampleID  Abundance
1 Proteobacteria      KO1 0.66311258
2 Actinobacteria      KO1 0.25731788
3  Bacteroidetes      KO1 0.02798013
4     Firmicutes      KO1 0.01394040
5    Chloroflexi      KO1 0.01480132
6     Unassigned      KO1 0.01864238

# 添加分组信息
phylum.gg$group <- c(rep('KO', 66), rep('OE', 66), rep('WT', 66)) # 根据样本情况设置

绘制

# 为了复现文章中的图需要的颜色包
library(wesanderson)
library(colortools)
library(ggpubr)

ggbarplot(phylum.gg, x = "SampleID", y="Abundance", color="black", fill="PhylumID",
          legend="right", 
          legend.title="Top Phylum", main="Relative counts per Phylum",
          font.main = c(14,"bold", "black"), font.x = c(12, "bold"), 
          font.y=c(12,"bold")) + 
  theme_bw() +
  rotate_x_text() + 
  scale_fill_manual(values=c("gray",rev(wheel("skyblue3")[1:10]))) + # 颜色设置
  facet_grid(~ group, scales = "free_x", space='free') + 
  labs(x = "Sample", y = "Relative proportion") + 
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title = element_text(face = "bold"), 
        plot.title = element_text(face = "bold"), 
        legend.title = element_text(face = "bold")) 
ggsave(filename = "relative_counts.pdf", device="pdf", width=8, height=4)

结果展示

数据及代码

为爱发电不易~如果需要示例数据及代码(当然其实文中都已经写的很详细了)的同学，点赞&打赏 5￥以上，并在主页添加本人微信(MZBJ0001)发送截图即可。

后记

关于更详细的代码讲解、作者的原代码的一些细节以及我修改的地方会在之后的视频教程中详细讲到，有兴趣的可以关注我的B站【木舟笔记】。

往期内容

继续滑动看下一个