# 首先需要安装依赖包 “amplicon”，如果没有安装 “devtools” 需要先安装，命令如下：

if (!requireNamespace("devtools", quietly=TRUE))

    install.packages("devtools")

library(devtools)

# 安装依赖包 “amplicon”，我对原作者的程序做了些许的优化，所以，可以从我的 GitHub 进行安装，命令如下：

if (!requireNamespace("amplicon", quietly=TRUE))

    install_github("liaochenlanruo/amplicon")

suppressWarnings(suppressMessages(library(amplicon)))

# 设置工作目录，视自己的具体情况而定，工作目录中当包含输入文件：

setwd("E:/Researches/Xiaqian/NGS/CleanData/ALL/细菌V6V8-1/Analysis_20210624/barplot")

# 从文件读取元数据和物种注释

metadata=read.table("sample-metadata.tsv", header=T, row.names=1, sep="\t", comment.char="", stringsAsFactors=F)

taxonomy=read.table("taxonomy.tsv", header=T, row.names=1, sep="\t", comment.char="", stringsAsFactors=F)

# 绘制堆叠柱状图，颜色自定义。界水平物种组成表和元数据作为输入，分组列名为 Group，显示前 4 个分类，其余归类为其他 (Other)，按丰度排序。

otutab1=read.table("level-1.csv", header=T, row.names=1, sep=",", comment.char="", stringsAsFactors=F)

otutab1=t(otutab1)

#按照站位分组，通过 “groupID” 参数设置分组信息，这里我选择根据站位（Site）分组，用户当根据自己的实际情况来设置。

(ds=tax_stackplot(otutab1, metadata, groupID="Site", topN=4, style="sample", sorted="abundance"))

ds2 = ds + scale_fill_manual(values = rev(c("#9467bd", "#ff7f0e", "#1f77b4", "#17becf"))) + theme(axis.text.x = element_text(angle = 90, hjust = 0.5, vjust = 0.5))

# 出图

ds2

#按照深度分组

(dd=tax_stackplot(otutab1, metadata, groupID="Depth", topN=4, style="sample", sorted="abundance"))

dd2 = dd + scale_fill_manual(values = rev(c("#9467bd", "#ff7f0e", "#1f77b4", "#17becf"))) + theme(axis.text.x = element_text(angle = 90, hjust = 0.5, vjust = 0.5))

dd2

# 门水平物种组成表和元数据作为输入，分组列名为 Group，显示前 19 个分类，按丰度排序

otutab2=read.table("level-2.csv", header=T, row.names=1, sep=",", comment.char="", stringsAsFactors=F)

otutab2=t(otutab2)

#按照站位分组

(ps=tax_stackplot(otutab2, metadata, groupID="Site", topN=20, style="sample", sorted="abundance"))

ps2 = ps + scale_fill_manual(values = rev(c("#1f77b4", "#aec7e8", "#ff7f0e", "#ffbb78", "#2ca02c", "#98df8a", "#d62728", "#ff9896", "#9467bd", "#c5b0d5", "#8c564b", "#c49c94", "#e377c2", "#f7b6d2", "#7f7f7f", "#c7c7c7", "#bcbd22", "#dbdb8d", "#17becf", "#9edae5"))) + theme(axis.text.x = element_text(angle = 90, hjust = 0.5, vjust = 0.5))

ps2

#按照深度分组

(pd=tax_stackplot(otutab2, metadata, groupID="Depth", topN=20, style="sample", sorted="abundance"))

pd2 = pd + scale_fill_manual(values = rev(c("#1f77b4", "#aec7e8", "#ff7f0e", "#ffbb78", "#2ca02c", "#98df8a", "#d62728", "#ff9896", "#9467bd", "#c5b0d5", "#8c564b", "#c49c94", "#e377c2", "#f7b6d2", "#7f7f7f", "#c7c7c7", "#bcbd22", "#dbdb8d", "#17becf", "#9edae5"))) + theme(axis.text.x = element_text(angle = 90, hjust = 0.5, vjust = 0.5))

pd2

# 纲水平物种组成表和元数据作为输入，分组列名为 Group，显示前 19 个分类，按丰度排序

otutab3=read.table("level-3.csv", header=T, row.names=1, sep=",", comment.char="", stringsAsFactors=F, quote="\"")

otutab3=t(otutab3)

#按照站位分组

(cs=tax_stackplot(otutab3, metadata, groupID="Site", topN=20, style="sample", sorted="abundance"))

cs2 = cs + scale_fill_manual(values = rev(c("#1f77b4", "#aec7e8", "#ff7f0e", "#ffbb78", "#2ca02c", "#98df8a", "#d62728", "#ff9896", "#9467bd", "#c5b0d5", "#8c564b", "#c49c94", "#e377c2", "#f7b6d2", "#7f7f7f", "#c7c7c7", "#bcbd22", "#dbdb8d", "#17becf", "#9edae5"))) + theme(axis.text.x = element_text(angle = 90, hjust = 0.5, vjust = 0.5))

cs2

#按照深度分组

(cd=tax_stackplot(otutab3, metadata, groupID="Depth", topN=20, style="sample", sorted="abundance"))

cd2 = cd + scale_fill_manual(values = rev(c("#1f77b4", "#aec7e8", "#ff7f0e", "#ffbb78", "#2ca02c", "#98df8a", "#d62728", "#ff9896", "#9467bd", "#c5b0d5", "#8c564b", "#c49c94", "#e377c2", "#f7b6d2", "#7f7f7f", "#c7c7c7", "#bcbd22", "#dbdb8d", "#17becf", "#9edae5"))) + theme(axis.text.x = element_text(angle = 90, hjust = 0.5, vjust = 0.5))

cd2

# 属水平物种组成表和元数据作为输入，分组列名为 Site，显示前 19 个分类，按丰度排序

otutab6=read.table("level-6.csv", header=T, row.names=1, sep=",", comment.char="", stringsAsFactors=F)

otutab6=t(otutab6)

#按照站位分组

(gs=tax_stackplot(otutab6, metadata, groupID="Site", topN=20, style="sample", sorted="abundance"))

gs2 = gs + scale_fill_manual(values = rev(c("#1f77b4", "#aec7e8", "#ff7f0e", "#ffbb78", "#2ca02c", "#98df8a", "#d62728", "#ff9896", "#9467bd", "#c5b0d5", "#8c564b", "#c49c94", "#e377c2", "#f7b6d2", "#7f7f7f", "#c7c7c7", "#bcbd22", "#dbdb8d", "#17becf", "#9edae5"))) + theme(axis.text.x = element_text(angle = 90, hjust = 0.5, vjust = 0.5))

gs2

#按照深度分组

(gd=tax_stackplot(otutab6, metadata, groupID="Depth", topN=20, style="sample", sorted="abundance"))

gd2 = gd + scale_fill_manual(values = rev(c("#1f77b4", "#aec7e8", "#ff7f0e", "#ffbb78", "#2ca02c", "#98df8a", "#d62728", "#ff9896", "#9467bd", "#c5b0d5", "#8c564b", "#c49c94", "#e377c2", "#f7b6d2", "#7f7f7f", "#c7c7c7", "#bcbd22", "#dbdb8d", "#17becf", "#9edae5"))) + theme(axis.text.x = element_text(angle = 90, hjust = 0.5, vjust = 0.5))

gd2