读图
Snipaste_2021-11-15_21-09-21
该图为分组的物种属水平的相对丰度堆积柱状图。其中只显示了前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)
结果展示
