扩增子分析较常用的软件是QIIME2，然而其生成的图并不美观，且无法下载，是不可能够达到发表要求的，因此需要我们通过其他方法将QIIME2生成的数据结果进行可视化。本文将记录相关的分析方法和示例。

准备数据

特征表（Feature table），或者叫做OTU表，本教程中文件名为feature-table.tsv。该文件可由QIIME2生成，由biome格式转换而来。

	F01	F02	F03	F04	F05	F06	H2_1	H2_2	H2_3	H2_4	H2_5	H2_6	I1_1	I1_2	I1_3	I1_4	I1_5	I1_6
3733c766539aba5e3e8f964a5ba39048	0.0	0.0	0.0	0.0	0.0	0.0	0.0	6.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
547715c94c5b0c7ba2049658f732397b	0.0	0.0	0.0	6.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
708c69ff1252447c0709333304cd9bed	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	5.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
344ee69a9e51ffac2811294d3f04833d	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	7.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
d2e467b7867f49c25caaacdea5c7eb43	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	18.0	0.0	0.0	0.0	0.0
4adc9b2b0866f620425d07d3ff7ca7e2	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	8.0	0.0	0.0	0.0	0.0
0ffcdef70b9c1436126ba32d7a40d961	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	31.0	0.0	0.0	0.0
7170ce6c7520acc3c13189478aea54a8	53.0	0.0	0.0	0.0	0.0	0.0	37.0	25.0	0.0	0.0	0.0	0.0	23.0	21.0	0.0	0.0	0.0	0.0
b624986630f3d33fbe0a21e66c46d1a5	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	9.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
78305c3f5955043aba7de952c78859c9	0.0	0.0	22.0	57.0	73.0	55.0	0.0	0.0	0.0	82.0	73.0	0.0	0.0	0.0	43.0	132.0	82.0	95.0
60d9992b10f6dd4da4dc53826b7433f5	0.0	0.0	0.0	0.0	45.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
33468988848349138af5cb1c23b550f2	0.0	2.0	11.0	12.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0

样本元数据（Metadata），本教程中文件名为sample-metadata.tsv。该文件由分析者自己准备，用于描述各个样本的特征。

  ID	Samples	Layer	Depth	Depth2	Location	Position	Site
F01	F01	0-5	5	05 cmbsf	Turn	Upstream	F
F02	F02	5-9	10	10 cmbsf	Turn	Upstream	F
F03	F03	9-13	15	15 cmbsf	Turn	Upstream	F
F04	F04	13-17	15	15 cmbsf	Turn	Upstream	F
F05	F05	17-22	20	20 cmbsf	Turn	Upstream	F
F06	F06	22-26	25	25 cmbsf	Turn	Upstream	F
H2_1	H2_1	0-5	5	05 cmbsf	Left	Upstream	H2
H2_2	H2_2	5-10	10	10 cmbsf	Left	Upstream	H2
H2_3	H2_3	10-15	15	15 cmbsf	Left	Upstream	H2
H2_4	H2_4	15-20	20	20 cmbsf	Left	Upstream	H2
H2_5	H2_5	20-25	25	25 cmbsf	Left	Upstream	H2
H2_6	H2_6	25-29	30	30 cmbsf	Left	Upstream	H2
I1_1	I1_1	0-5	5	05 cmbsf	Left	Upstream	I1
I1_2	I1_2	5-10	10	10 cmbsf	Left	Upstream	I1
I1_3	I1_3	10-15	15	15 cmbsf	Left	Upstream	I1
I1_4	I1_4	15-20	20	20 cmbsf	Left	Upstream	I1
I1_5	I1_5	20-25	25	25 cmbsf	Left	Upstream	I1
I1_6	I1_6	25-29	30	30 cmbsf	Left	Upstream	I1

分类信息，本教程中文件名为taxonomy.tsv。该文件可由QIIME2生成。

Feature ID	Taxon	Confidence
9e24eb7a415db9a59586d17d0f49cbea	Archaea;Asgardarchaeota;Lokiarchaeia	0.9999999998216045
2a1d71d9fca92f9b3fc8f3126db076d0	Bacteria;Proteobacteria;Gammaproteobacteria;Pseudomonadales;Pseudomonadaceae;Pseudomonas	0.9881304182870394
5d89868f1a1db1bc0b153eca96bbf7fd	Archaea;Crenarchaeota;Nitrososphaeria;Nitrosopumilales;Nitrosopumilaceae	0.9999999999859067
cc957f2e853493d72875a6528ca83435	Archaea;Crenarchaeota;Nitrososphaeria;Nitrosopumilales;Nitrosopumilaceae	0.9999999999854506
5179960239613d4c220e889075c6773e	Archaea;Crenarchaeota;Nitrososphaeria;Nitrosopumilales;Nitrosopumilaceae;Candidatus Nitrosopumilus	0.9558707678091154
d62ab321b324b106e11961e79384f974	Bacteria;Methylomirabilota;Methylomirabilia;Methylomirabilales;Methylomirabilaceae;wb1-A12	0.9999999975860316
920aa9728ffb783cbc5721932cadda36	Bacteria;Aerophobota;Aerophobia;Aerophobales;uncultured bacterium	0.9999939600150409
b5edcefd9280528ee377cb70dea8a6a7	Archaea;Asgardarchaeota;Lokiarchaeia;uncultured archaeon	0.99744406884555
43b02d567fe419b94ad063081c7bba58	Bacteria;Planctomycetota;Brocadiae;Brocadiales;Scalinduaceae;Candidatus Scalindua;uncultured bacterium	0.9999941829048367
0f1be1121dd04bc63b5ec4be0efc09b8	Archaea;Asgardarchaeota;Lokiarchaeia;uncultured archaeon	0.7884239553020901
430d9115ca79bad01227bcc9c7694cf1	Bacteria;Proteobacteria;Gammaproteobacteria;Xanthomonadales;Xanthomonadaceae;Thermomonas;uncultured bacterium	0.7413321650122459
0bf08973da329747d913d5c261f91944	Bacteria;Aerophobota;Aerophobia;Aerophobales	1.0000000000000024

进化树，本教程中文件名为tree.nwk。该文件可由QIIME2生成，Newick格式，是用代表序列构建的系统发育树。

Alpha diversity

Shannon、Richness、Simpson、Chao1、ACE、Pielou、PD_whole_tree指数的计算与可视化

输入文件

feature-table.tsv
sample-metadata.tsv

第一步：计算多样性

将如下代码保存到文件”compute_alpha.R“中，并将程序与输入文件放在同一目录下。在终端中运行命令”Rscript compute_alpha.R“，得到输出文件”alpha.txt“。

#!/usr/bin/env R
#定义函数
library(picante)       #picante 包加载时默认同时加载 vegan
 
alpha <- function(x, tree = NULL, base = exp(1)) {
        est <- estimateR(x)
        Richness <- est[1, ]
        Chao1 <- est[2, ]
        ACE <- est[4, ]
        Shannon <- diversity(x, index = 'shannon', base = base)
        Simpson <- diversity(x, index = 'simpson')    #Gini-Simpson 指数
        Pielou <- Shannon / log(Richness, base)
        goods_coverage <- 1 - rowSums(x == 1) / rowSums(x)
        
        result <- data.frame(Richness, Shannon, Simpson, Pielou, Chao1, ACE, goods_coverage)
        if (!is.null(tree)) {
                PD_whole_tree <- pd(x, tree, include.root = FALSE)[1]
                names(PD_whole_tree) <- 'PD_whole_tree'
                result <- cbind(result, PD_whole_tree)
        }
        result
}
 
#现在直接使用定义好的命令 alpha()，一步得到多种 Alpha 多样性指数
#加载 OTU 丰度表和进化树文件
otu <- read.delim('feature-table.tsv', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
otu <- t(otu)
tree <- read.tree('tree.nwk')
 
#不包含谱系多样性，无需指定进化树；Shannon 公式的 log 底数我们使用 2
##alpha_all <- alpha(otu, base = 2)
#包含谱系多样性时，指定进化树文件；Shannon 公式的 log 底数我们使用 2
alpha_all <- alpha(otu, tree, base = 2)
 
#输出保存在本地
write.table(alpha_all, 'alpha.txt', sep="\t", quote = FALSE)

第二步：制图

这一步需要在Rstudio中运行如下命令，在终端里直接运行会报错（莫名其妙），最终输出7个PDF图形文件。

#!/usr/bin/env R
# 基于github安装包，需要devtools，检测是否存在，不存在则安装
if (!requireNamespace("devtools", quietly = TRUE))
    install.packages("devtools")

# 加载github包安装工具
library(devtools)

# 检测amplicon包是否安装，没有从源码安装
if (!requireNamespace("amplicon", quietly = TRUE))
    install_github("liaochenlanruo/amplicon")
# 提示升级，选择3 None不升级；升级会容易出现报错

# library加载包，suppress不显示消息和警告信息
suppressWarnings(suppressMessages(library(amplicon)))

# 设置工作目录，根据实际情况设定
setwd("用户自己的工作目录")

# 读取元数据，参数指定包括标题行(TRUE)，列名为1列，制表符分隔，无注释行，不转换为因子类型
metadata = read.table("sample-metadata.tsv", header=T, row.names=1, sep="\t", comment.char="", stringsAsFactors = F)

# 读取第一步中vegan计算的6种alpha多样性指数
#alpha_div = read.table("alpha.txt", header=T, row.names=1, sep="\t", comment.char="")
alpha_div <- read.delim("alpha.txt", quote="", stringsAsFactors=FALSE)

# 绘制各组香农指数分布，外层()可对保存的图形同时预览
(p1 = alpha_boxplot(alpha_div, index = "Shannon", metadata, groupID = "Depth2"))

##物种丰富度 Richness 指数
(p2 = alpha_boxplot(alpha_div, index = "Richness", metadata, groupID = "Depth2"))

#Gini-Simpson 指数（我们平时常用的 Simpson 指数即为 Gini-Simpson 指数）
(p3 = alpha_boxplot(alpha_div, index = "Simpson", metadata, groupID = "Depth2"))

#Chao1 指数
(p4 = alpha_boxplot(alpha_div, index = "Chao1", metadata, groupID = "Depth2"))

#ACE 指数
(p5 = alpha_boxplot(alpha_div, index = "ACE", metadata, groupID = "Depth2"))

#Shannon 均匀度（Pielou 均匀度）
(p6 = alpha_boxplot(alpha_div, index = "Pielou", metadata, groupID = "Depth2"))

##谱系多样性
(p7 = alpha_boxplot(alpha_div, index = "PD_whole_tree", metadata, groupID = "Depth2"))


# 保存图片，指定图片为pdf格式方便后期修改，图片宽89毫米，高75毫米
ggsave(paste0("alpha_boxplot_shannon.pdf"), p1, width=89, height=75, units="mm")
ggsave(paste0("alpha_boxplot_Richness.pdf"), p2, width=89, height=75, units="mm")
ggsave(paste0("alpha_boxplot_Simpson.pdf"), p3, width=89, height=75, units="mm")
ggsave(paste0("alpha_boxplot_Chao1.pdf"), p4, width=89, height=75, units="mm")
ggsave(paste0("alpha_boxplot_ACE.pdf"), p5, width=89, height=75, units="mm")
ggsave(paste0("alpha_boxplot_Pielou.pdf"), p6, width=89, height=75, units="mm")
ggsave(paste0("alpha_boxplot_PD_whole_tree.pdf"), p7, width=89, height=75, units="mm")

演示视频

物种多样性条形图的绘制

输入文件

处理后的Feature Table

从QIIME2可视化网站上下载的各level的物种丰度表，此处以level-2（Phylum水平）为例，命名为“level-2.csv”。该文件可以用excel或WPS打开，然后复制所有内容至一个新建的名字为“level-2.txt”的文本文档中，删除无关的行或列（如路径），只保留丰度数据，并保存。格式如下：

index	species1	species2	species3	species4
Sample1	0	8152	11	65
Sample2	12	97	1097	1315
...
Samplex	15	46	2928	1011

sample-metadata.tsv
taxonomy.tsv

绘图

以下代码在Rstudio中逐步输入运行

setwd("用户自己的路径")

if (!requireNamespace("devtools", quietly=TRUE))
    install.packages("devtools")
library(devtools)

if (!requireNamespace("amplicon", quietly=TRUE))
    install_github("microbiota/amplicon")
suppressWarnings(suppressMessages(library(amplicon)))

# 从文件读取元数据，特征表和物种注释
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)

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

#堆叠柱状图。以分组均值绘制，展示丰度最高的19个门，其余归类为其他(Other)
# 门水平物种组成表和元数据作为输入，分组列名为Group，默认显示前19个分类，按丰度排序
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

演示视频

Beta多样性

NMDS分析及可视化

输入文件

feature-table.tsv
sample-metadata.tsv

分析过程

在Rstudio中运行如下代码：

setwd("用户的工作目录")
library("vegan")

# 读取数据，一份otu.table文件和一份分组信息文件
otu <- read.table("feature-table.tsv",row.names=1, header=T,sep="\t",check.names=F)

design <- read.table("sample-metadata.tsv",header=T,sep="\t",row.names=1,check.names=F)

# 数据调整为列名是OTU，行名是样本名
otu=t(otu)

# 标准化，当然method有很多，可以通过?decostand查看其它method
vare.hel<-decostand(otu,method="hellinger")

# 计算Bray-curtis距离矩阵
vare.dis <- vegdist(vare.hel,method="bray")

# 使用NMDS的方法
vare.mds <- metaMDS(vare.hel,distance = "bray")

# 首先提取前两轴坐标
point = scores(vare.mds)

# 将分组文件和得分文件合并
index = merge(design, point,by="row.names",all=F)

# 查看Stress值
# Stress值是反映模型合适程度的指标，NMDS会多次打乱数据计算Stress值，直到找到最合适的模型，也就是最低的Stress值；理想状况下，Stress值为0，一般Stress值低于0.1较为合理。
vare.mds

# Stress:     0.07498046 #记住这个数字，后面会打印在图中

# 显着性检验；anosim本质是基于排名的算法，更加适合NMDS，这里基于元数据中的深度（Depth）计算
anosim.result<-anosim(vare.dis, design$Depth,permutations =999)

summary(anosim.result)

# ANOSIM statistic R: 0.8036 #记住这个数字，后面会打印在图中
#       Significance: 0.001 #记住这个数字，后面会打印在图中

# tiff输出图形，适合大部分出版刊物，入门级别分辩率300,18*14的长宽；
#tiff(file="beta_bray_NMDS.tiff", res = 300, compression ="none", width=180,height=140,units= "mm")

# 在这里我选择输出为PDF，后续可以进一步的编辑
pdf(file="beta_bray_NMDS.pdf", width=180,height=140)

#开始出图，代码如下：
library("ggplot2")

p = ggplot(index, aes(x=NMDS1, y=NMDS2, color=as.factor(Depth))) +
      geom_point(size=0) +
     #scale_colour_manual(values = c("red","blue", "green")) +#这行代码是用于自定义颜色的，有几个组就用几个颜色，少于3个组的话不支持自定义
     labs(x=paste("NMDS1"), y=paste("NMDS2"), title="")

#置信区间当然要加上，有三种方式，线条类型也可以更改。将上面得到的三个指标在图中更换Stress，R，p。

p+stat_ellipse(type = "t", linetype = 2, level = 0.95, show.legend = TRUE)+
    annotate("text",x=-2.5,y=-1.25,parse=TRUE,size=4,label="'R: '*0.8036",family="serif",fontface="italic",colour="darkred",hjust = 0)+
    annotate("text",x=-2.5,y=-1.35,parse=TRUE,size=4,label="'p: '*0.001",family="serif",fontface="italic",colour="darkred",hjust = 0)+
    annotate("text",x=-2.5,y=-1.45,parse=TRUE,size=4,label="'Stress: '*0.075",family="serif",fontface="italic",colour="darkred",hjust = 0) + geom_point(aes(x=NMDS1, y=NMDS2, shape=Site), size=2)+
	#scale_shape_manual(values = c(15, 19, 17))+#这行代码是用于自定义形状的，有几个组就用几个形状，少于3个组的话不支持自定义
    theme_classic()

dev.off()

演示视频

参考

https://github.com/YongxinLiu/MicrobiomeStatPlot

敬告：使用文中脚本请引用本文网址，请尊重本人的劳动成果，谢谢！Notice: When you use the scripts in this article, please cite the link of this webpage. Thank you!

了尘兰若的小坑

分享微生物生物信息学分析方法，欢迎加入QQ群交流945751012，不接受群内广告！

准备数据

Alpha diversity

Shannon、Richness、Simpson、Chao1、ACE、Pielou、PD_whole_tree指数的计算与可视化

输入文件

第一步：计算多样性

第二步：制图

演示视频

物种多样性条形图的绘制

输入文件

绘图

演示视频

Beta多样性

NMDS分析及可视化

输入文件

分析过程

演示视频

参考

你的赏识是我前进的动力