Alpha多样性绘图
EasyAmplicon 流程提供生成各种Alpha多样性图的功能,可深入了解每个样本中微生物群落的丰富度和均匀度。这些图包括Alpha柱状图(通常显示观察到的OTU数量)、Alpha箱线图(可视化OTU计数分布)和Alpha稀释曲线(估计同一测序深度的多样性)。都是通过数据文件(例如alpha.txt、vegan.txt、alpha_rare.txt、otu_group_exist.txt)与ggplot2和vegan等R包来生成的(图 4A)。
箱线图
运行以下命令绘制Alpha箱线图:
Rscript ${db}/script/alpha_boxplot.R -h
Rscript ${db}/script/alpha_boxplot.R --alpha_index richness \
--input result/alpha/vegan.txt --design result/metadata.txt \
--group Group --output result/alpha/ \
--width 89 --height 59
使用循环绘制6种常用指数:
for i in `head -n1 result/alpha/vegan.txt|cut -f 2-`;do
Rscript ${db}/script/alpha_boxplot.R --alpha_index ${i} \
--input result/alpha/vegan.txt --design result/metadata.txt \
--group Group --output result/alpha/ \
--width 89 --height 59
done
mv alpha_boxplot_TukeyHSD.txt result/alpha/
此外,要绘制Alpha多样性柱状图并计算标准差,请运行以下命令:
Rscript ${db}/script/alpha_barplot.R --alpha_index richness \
--input result/alpha/vegan.txt --design result/metadata.txt \
--group Group --output result/alpha/ \
--width 89 --height 59
稀释曲线
稀释曲线说明了组测序深度平均值与特征标准误差(OTU或ASV)计数。输入文件alpha_rare.txt,指定实验设计和分组,设置输出目录,并定义图像尺寸(宽度和高度)。为了绘制曲线,请运行以下代码:
Rscript ${db}/script/alpha_rare_curve.R \
--input result/alpha/alpha_rare.txt --design result/metadata.txt \
--group Group --output result/alpha/ \
--width 120 --height 59
维恩图
方法1:要绘制Alpha多样性维恩图,需输入文件otu_group_exist.txt,该文件包含两列:ID和group。命令行如下:其中“-f”用于输入文件,“-a/b/c/d/g”用于组名,“-w/u”用于宽度和高度(单位为英寸),“-p”用于输出文件名(默认值与输入相同)(图 4A)。运行以下代码比较三个组:
bash ${db}/script/sp_vennDiagram.sh \
-f result/alpha/otu_group_exist.txt \
-a WT -b KO -c OE \
-w 3 -u 3 \
-p WT_KO_OE
运行以下代码比较四个组:
bash ${db}/script/sp_vennDiagram.sh \
-f result/alpha/otu_group_exist.txt \
-a WT -b KO -c OE -d All \
-w 3 -u 3 \
-p WT_KO_OE_All
方法2(使用R Markdown(Rmd)绘制Alpha箱线图):本节将引导您使用位于EasyAmplicon主文件夹的“result”目录或EasyMicrobiome的脚本目录。在RStudio中打开Diversity.Rmd文件以运行代码并生成可视化。在Rmd文件的实验设计部分,仔细检查与您的实验相关的元数据。这涉及验证分组变量(group)、图像尺寸(width和height),以及确认用于Alpha多样性计算的文件和索引类型,默认为“Richness”。有两种方法可以执行代码:首先,您可以逐行运行代码,这样您就可以在代码执行过程中监控每个变量的创建。或者,您可以使用RStudio中的右箭头按钮按顺序执行所有代码块,或者“Knit”整个文档,这将生成包含所有结果和计算的综合网页报告。生成的输出可以在“alpha”目录中找到。Alpha多样性(丰富度)的箱线图位于文件“alpha/alpha_boxplot_Richness.pdf”中。用户还可以在“alpha”目录中找到稀释曲线(如果包含),具体在“alpha/alpha_rarefaction_curve.pdf ”中。“#稀释曲线”部分中的说明将提供特定可视化图输入和输出文件位置的详细信息。
方法3:此外,用户可以访问ImageGP(https://www.bic.ac.cn/BIC/#/),在线绘制箱线图(图 S1-A)。首先在Excel中打开数据文件(metadata.txt和vegan.txt),然后将“group”和“richness”列复制到新工作表中。将这些组合数据粘贴到网站文本框中,并选择分组(group)和丰富度(Y 轴)。单击“PLOT ”以生成箱线图。同样,用户可以使用此链接绘制维恩图和稀释曲线(图 S1-A)。
Beta多样性
距离矩阵热图
距离矩阵样本的微生物组成来表示样本之间的差异。它在微生物组分析中起着至关重要的作用,例如对相似样本进行聚类。要使用Bray-Curtis作为示例绘制距离矩阵热图,您可以运行以下代码:
bash ${db}/script/sp_pheatmap.sh \
-f result/beta/bray_curtis.txt \
-H 'TRUE' -u 6 -v 5
要在第2列和第4列添加分组注释(例如基因型和位置),请运行以下代码:
cut -f 1-2 result/metadata.txt > temp/group.txt
接下来绘制热图:
bash ${db}/script/sp_pheatmap.sh \
-f result/beta/bray_curtis.txt \
-H 'TRUE' -u 6.9 -v 5.6 \
-P temp/group.txt -Q temp/group.txt
其中“-f”是输入文件,“-H”是聚类(TRUE / FALSE),“-u / v”是宽度和高度(英寸),“-P”为文件添加行注释,“-Q”添加列注释(图 4A)。
主坐标分析PCoA
本节介绍如何执行主坐标分析(PCoA)来可视化数据中的Beta多样性。主要有两种方法:
方法1:此方法使用R脚本”beta_pcoa.R”,利用bray_curtis.txt和metadata.txt作为输入文件并输出bray_curtis.cpcoa.pdf。
Rscript ${db}/script/beta_pcoa.R \
--input result/beta/bray_curtis.txt --design result/metadata.txt \
--group Group --label FALSE --width 89 --height 59 \
--output result/beta/bray_curtis.pcoa.pdf
要使用参数“--label”添加示例标签,请运行以下脚本:
Rscript ${db}/script/beta_pcoa.R \
--input result/beta/bray_curtis.txt --design result/metadata.txt \
--group Group --label TRUE --width 89 --height 59 \
--output result/beta/bray_curtis.pcoa.label.pdf
mv beta_pcoa_stat.txt result/beta/
若要查看组间两两比较P值,请运行以下代码(概率事件每次略有不同):
cat result/beta/beta_pcoa_stat.txt
方法2:使用“Diversity.Rmd”绘制PCoA:使用R studio打开“Diversity.Rmd”文件,检查并运行“#β diversity多样性——##PCoA主坐标轴分析”。Knit输出文档生成结果和计算过程,包含beta/pcoa_bray_curtis.pdf和beta_pcoa_stat.txt(包含组间两两比较P值)两个文件。
方法3:同样,用户也可以使用在线工具http://www.bic.ac.cn/ImageGP/或https://www.bic.ac.cn/BIC/#/在线绘制PCoA图和热图。(图 S1-A)。
限制性主坐标分析CPCoA
本节介绍限制性主坐标分析(CPCoA),可探究环境因素或其他变量如何影响Beta多样性。执行CPCoA有两种主要方法:
方法1:此方法利用提供的R脚本“beta_cpcoa.R”。脚本参数与“beta_pcoa.R”脚本类似:
Rscript ${db}/script/beta_cpcoa.R \
--input result/beta/bray_curtis.txt --design result/metadata.txt \
--group Group --output result/beta/bray_curtis.cpcoa.pdf \
--width 89 --height 59
使用--label TRUE参数添加样本标签:
Rscript ${db}/script/beta_cpcoa.R \
--input result/beta/bray_curtis.txt --design result/metadata.txt \
--group Group --label TRUE --width 89 --height 59 \
--output result/beta/bray_curtis.cpcoa.label.pdf
方法2:使用“Diversity.Rmd”脚本可视化CPCoA图。打开“Diversity.Rmd”文件,找到“#β diversity多样性——## Constrained有监督PCoA/CCA”部分。修改脚本中的文件路径为实际Beta多样性矩阵位置(可能位于“result/beta”目录中)。文件路径调整后,“knit”R Markdown文档生成一份全面的网页报告,包含PCoA图和其他相关结果,以及一个名为“beta/cpcoa_bray_curtis.pdf”的附加文件,展示含有Bray-Curtis距离的CPCoA图。
方法3:虽然提供的脚本支持Bray-Curtis距离,但如果用户想要使用其他距离算法,请考虑使用在线工具http://www.bic.ac.cn/ImageGP/(图 S1-A)。此工具提供了十多个非进化距离选项,但无法计算与Unifrac相关的距离。
物种组成分类
用户可以通过多种可视化方式(如堆叠条形图stackplots、圈图chord diagrams和树状图)(图 4A)探索不同层级微生物群落的分类组成。每种图都有两种绘制方法。
堆叠柱状图
方法1:以门(p)水平为例,结果包括.sample.pdf和.group.pdf两种文件。输入文件是先前生成的sum_p.txt。
Rscript ${db}/script/tax_stackplot.R \
--input result/tax/sum_p.txt --design result/metadata.txt \
--group Group --color ggplot --legend 7 --width 89 --height 59 \
--output result/tax/sum_p.stackplot
可以使用ggplot manual(22)、Paired(12)或Set3(12)修改颜色。
Rscript ${db}/script/tax_stackplot.R \
--input result/tax/sum_p.txt --design result/metadata.txt \
--group Group --color Paired --legend 12 --width 181 --height 119 \
--output result/tax/sum_p.stackplotPaired
批量绘制包含门(p),纲(c),目(o),科(f),属(g)五个分类水平的输入数据:
for i in p c o f g; do
Rscript ${db}/script/tax_stackplot.R \
--input result/tax/sum_${i}.txt --design result/metadata.txt \
--group Group --output result/tax/sum_${i}.stackplot \
--legend 8 --width 89 --height 59; done
方法2(可选):通过.Rmd文件绘制堆叠直方图:用R studio打开“Diversity.Rmd”,检查“#物种组成——堆叠柱状图”段落参数,调整分类水平、输入输出文件位置和图片大小,点击命令即可运行。
弦/圈图
以纲(c)为例,运行以下代码绘制出前5组:
i=c
Rscript ${db}/script/tax_circlize.R \
--input result/tax/sum_${i}.txt --design result/metadata.txt \
--group Group --legend 5
结果位于当前目录“circlize_pdf”(随机颜色)和“circlize_legend.pdf”(指定颜色和图例)。要移动文件并改名使之与分类级一致:
mv circlize.pdf result/tax/sum_${i}.circlize.pdf
mv circlize_legend.pdf result/tax/sum_${i}.circlize_legend.pdf
方法2:运行.Rmd文件绘制弦/圈图。打开R studio的“Diversity.Rmd”并检查“## Circlize plot弦图”段落参数,调整分类水平、图例数量和分组列名。运行段落并在文件“circlize.pdf”和“circlize_legend.pdf”中查看结果。
树形图可视化
要可视化描述物种关系的分层树状图,请运行以下代码。通过输入特征表和物种注释来生成树状图,一定要指定特征的数量和图像的尺寸。
Rscript ${db}/script/tax_maptree.R \
--input result/otutab.txt --taxonomy result/taxonomy.txt \
--output result/tax/tax_maptree.pdf \
--topN 100 --width 183 --height 118
R差异分析
视频教程链接:https://youtu.be/FsdYNQ5DtHg
差异比较
这一步的重点是识别组间的显著不同的特征。为了分析差异,用户必须指定输入数据,定义比较组,设置参数,进行统计检验并生成火山图、热图、曼哈顿图、特征图和三元图(图 4B)。为此,需要在result文件夹中创建一个名为“compare”的目录,运行以下代码:mkdir -p result/compare/。使用输入特征表(otutab.txt)和元数据(metadata.txt)指定分组列名、比较组和丰度。选择edgeR方法,p值阈值为0.05,FDR阈值为0.2并制定输出目录,运行以下代码:
compare="KO-WT"
Rscript ${db}/script/compare.R \
--input result/otutab.txt --design result/metadata.txt \
--group Group --compare ${compare} --threshold 0.1 \
--method edgeR --pvalue 0.05 --fdr 0.2 \
--output result/compare/
故障提示:Error in file(file, ifelse(append, "a", "w"))。解决方法:提示输出目录不存在,新建目录即可。
差异丰度结果的可视化与分析
EasyAmplicon可通过绘制各种可视化图像来展示不同的丰度结果。
火山图
为了绘制火山图(图 4B)并可视化差异倍数和统计显著性之间的关系,运行以下代码:
Rscript ${db}/script/compare_volcano.R \
--input result/compare/${compare}.txt \
--output result/compare/${compare}.volcano.pdf \
--width 89 --height 59
热图
要绘制显示样本和组间特征丰度模式的热图,请运行以下代码,该代码筛选列数、指定元数据和分组、物种注释、图形大小和字体大小(图 4B)。
bash ${db}/script/compare_heatmap.sh -i result/compare/${compare}.txt -l 7 \
-d result/metadata.txt -A Group \
-t result/taxonomy.txt \
-w 8 -h 5 -s 7 \
-o result/compare/${compare}
曼哈顿图
本节介绍如何生成曼哈顿图,以可视化不同组间物种丰度的显著差异(图 4B)。这些图有助于识别区分样本的关键特征:
bash ${db}/script/compare_manhattan.sh -i result/compare/${compare}.txt \
-t result/taxonomy.txt \
-p result/tax/sum_p.txt \
-w 183 -v 59 -s 7 -l 10 \
-o result/compare/${compare}.manhattan.p.pdf
其中,“-i”表示差异比较结果,“-t”表示物种注释,“-p”表示图例,“-w”表示宽度,“-v”表示高度,“-s”表示大小,“-l”表示图例最大值。若要显示详细信息,请使用以下代码切换到c类和L类:
bash ${db}/script/compare_manhattan.sh -i result/compare/${compare}.txt \
-t result/taxonomy.txt \
-p result/tax/sum_c.txt \
-w 183 -v 59 -s 7 -l 10 -L Class \
-o result/compare/${compare}.manhattan.c.pdf
要显示完整的图例,请运行以下代码:
bash ${db}/script/compare_manhattan.sh -i result/compare/${compare}.txt \
-t result/taxonomy.txt \
-p result/tax/sum_c.txt \
-w 183 -v 149 -s 7 -l 10 -L Class \
-o result/compare/${compare}.manhattan.c.legend.pdf
注意:或者,用户可以使用ImageGP(http://www.bic.ac.cn/ImageGP/)绘制火山图、热图和曼哈顿图(图 S1-A)。
单个特征的绘制
这一步的重点是可视化和分析在前面步骤中识别的差异丰富序列变体(ASV)和操作分类单元(OTU)。运行下列代码,筛选差异ASV,按KO组丰度降序排列,并展示前10名ID:
awk '$4<0.05' result/compare/KO-WT.txt | sort -k7,7nr | cut -f1 | head
故障提示:如果使用awk命令后没有返回任何结果(空输出),则可能表明数据中缺乏差异丰富的ASV(p值> 0.05)。您可以通过调整显著性阈值(0.05)或比较组的潜在问题以调查数据。
差异OTU分析:为了显示差异OTU详细信息,运行以下命令:
Rscript ${db}/script/alpha_boxplot.R --alpha_index ASV_2 --input result/otutab.txt --design result/metadata.txt --transpose TRUE --scale TRUE --width 89 --height 59 --group Group --output result/compare/feature_
代码Rscript ${db}/script/alpha_boxplot.R将启动${db}/script目录下一个名为alpha_boxplot.R的脚本。该脚本将生成一个箱线图,比较“metadata.txt”文件中定义的不同组之间特定Alpha多样性指数(ASV_2)的丰度,脚本中的其他参数设置输出格式和尺寸。
故障提示:Error in data.frame(..., check.names = FALSE) :此错误表明箱线图数据中的行数有问题。如果“otutab.txt”或“metadata.txt”文件针对每个样品的行数不同,则可能出现这种情况。检查数据处理或过滤步骤中的不一致处,确保两个文件的样本信息互相匹配。
属级分析:要按平均属丰度降序对列进行排序,请运行以下命令:
csvtk -t sort -k All:nr result/tax/sum_g.txt | head
此代码利用“csvtk”工具对位于“result/tax”目录中的“sum_g.txt”文件进行排序。“-t”标志指定制表符分隔的文件,“-k All:nr”按降序对所有列进行数字排序,“head”显示最上面的行(显示最丰富的属)。要进一步探索属性级别差异的细节,请运行以下代码:
Rscript ${db}/script/alpha_boxplot.R --alpha_index Lysobacter \
--input result/tax/sum_g.txt --design result/metadata.txt \
--transpose TRUE \
--width 89 --height 59 \
--group Group --output result/compare/feature_
故障提示:如果根据属丰度数据生成的箱线图不能提供有效信息(细节不足),可能是由于属的数量太多。可以考虑过滤数据,将重点放在前几个最丰富的属上,或者使用柱状图等其他可视化方法。
使用STAMP进行差异丰度分析
视频教程链接:https://youtu.be/PptGLAI93eE
STAMP是一个用于分析微生物组数据的独立软件,本节概述了如何使用STAMP进行差异丰度分析(图4C)。
生成输入文件
该步骤涉及STAMP数据的格式化,以及生成差异丰度结果的可视化。
Rscript ${db}/script/format2stamp.R -h
mkdir -p result/stamp
Rscript ${db}/script/format2stamp.R --input result/otutab.txt --taxonomy result/taxonomy.txt --threshold 0.01 --output result/stamp/tax
输出结果位于results/stamp目录中。用户可以在STAMP软件中打开这些文件。结果包括名为tax_1至tax_8的文件。tax_1至tax_7是界、门、纲、目、科、属和种的分类摘要,tax_8是过滤后的OTU表。其中,例如0.01代表平均过滤丰度。此流程还提供了通过format2stamp.Rmd生成输入文件的替代方法,位于EasyAmplicon主文件夹的结果目录中或EasyMicrobiome数据库文件夹的脚本目录中。要使用format2stamp.Rmd,用户必须在结果目录中拥有“otutab.txt”和“taxonomy.txt”文件。接下来,使用R studio打开“format2stamp.Rmd”。在运行代码之前,请检查并调整代码参数。单击RStudio中的“Knit”按钮执行代码。默认情况下,结果将保存在result/stamp目录中。
绘制扩展柱状图和表格
要进行两组比较(例如KO-WT),建议使用Welch 's t检验+扩展误差条/热图。将ASV (result/otutab.txt)替换为属(result/tax/sum_g.txt),并运行如下命令:compare="KO-WT"
Rscript ${db}/script/compare_stamp.R \
--input result/stamp/tax_5Family.txt --metadata result/metadata.txt \
--group Group --compare ${compare} --threshold 0.1 \
--method "t.test" --pvalue 0.05 --fdr "none" \
--width 189 --height 159 \
--output result/stamp/${compare}
还可以通过“compareStamp.Rmd”作为可选方法,该方法可以在EasyAmplicon主文件夹的结果目录中或EasyMicrobiome文件夹的脚本目录中找到。此文件包含用户可以根据需要查看和调整的代码,包括输入和输出文件的名称和位置或特定阈值。对设置满意后,单击RStudio中的“Knit”按钮执行代码,生成所需的输入文件。默认情况下,结果将放置在“result/compare”目录中。
注意:多组比较建议采用PCA、ANOVA + FDR +单特征箱线图。
生成LEfSe输入文件
视频教程链接:https://youtu.be/VXJvdidyFlU
LEfSe分析用于识别组间存在显著差异的微生物特征。本节概述了用于生成LEfSe分析输入文件的两种方法(图 4D)。
方法1:在此方法中,运行以下命令生成输入文件。其中,threshold参数调节筛选丰度,控制图中分支数量。
mkdir -p result/lefse
Rscript ${db}/script/format2lefse.R --input result/otutab.txt --taxonomy result/taxonomy.txt --design result/metadata.txt --group Group --threshold 0.4 --output result/lefse/LEfSe
方法2(可选):或者, 用户可以通过使用format2lefse.Rmd文件生成输入文件。结果包含三个文件:otutab.txt、metadata.txt、taxonomy.txt。要继续分析的话请到db/script目录,找到format2lefse.Rmd文件,然后使用RStudio打开它。运行代码将生成两个输入文件:LEfSe.txt和LEfSe2.txt。为了进行全面的LEfSe分析和可视化,我们推荐两个在线工具。首先,使用ImageGP在线服务器(https://www.bic.ac.cn/BIC/#/analysis?page=b%27MzY%3D%27)。只需提交您的LEfSe输入文件(LEfSe.txt),大约3分钟内,您就会收到高质量的可视化图(图 S1-B)。或者,LEfSe官网(https://huttenhower.sph.harvard.edu/lefse/)也提供在线LEfSe分析平台。
方法3:Linux下的LEfSe本地分析(可选),该部分将在下文“Linux下的分析”中进行说明。
特征预测分析
本节介绍如何使用PICRUSt和Bugbase分别预测微生物群落的功能和形态特征。此处的特征预测是指从扩增子测序数据生成的序列中预测或识别功能的过程。
PICRUSt特征预测
视频教程链接:https://youtu.be/bq9IkEJVAoQ
PICRUSt是一种专门设计用于根据标记基因数据(通常是16S rDNA基因序列)预测微生物群落功能能力(例如代谢途径)的分析工具(图 4E)。
PICRUSt 1.0
方法1(Windows本地运行):第一步是文件准备,即通过比对Greengenes数据库生成OTU表。这已经在上面完成(请参阅pipeline.sh/pipeline_en.sh,步骤9以供参考),生成“otutab.txt”文件。第二步是使用在线网络服务器ImageGP(https://www.bic.ac.cn/BIC/#/)上传gg/otutab.txt文件。上传前请确保将OUT ID 修改为“OTUID”。生成的结果有四个级别(图 S1-C),下载所有输出文件并保存在工作目录路径EasyAmplicon/result/picrust中(提前在结果目录中生成picrust文件夹)。获得的文件可以直接在STAMP中进行可视化,因为STAMP是整合和显示各级分析结果的强大工具。
现在,运行以下代码以从ImageGP下载的文件中删除多余的细节,以便使文件格式与后续分析保持一致。
cd ${wd}/result/picrust/
find . -type f -name "*all_level*" -exec bash -c 'mv "$0" "${0/*all_level/all_level}"' {} \;
cd ${wd}
接下来运行下面给出的代码:
l=L2
sed '/# Const/d;s/OTU //' result/picrust/all_level.ko.${l}.txt > result/picrust/${l}.txt
num=`head -n1 result/picrust/${l}.txt|wc -w`
paste
> result/picrust/${l}.spf
cut -f 2- result/picrust/${l}.spf > result/picrust/${l}.mat.txt
awk 'BEGIN{FS=OFS="\t"} {print $2,$1}' result/picrust/${l}.spf | sed 's/;/\t/' | sed '1 s/ID/Pathway\tCategory/' \
> result/picrust/${l}.anno.txt
为了可视化KEGG各层级柱状图,执行以下代码生成输入文件。compare="KO-WT"
Rscript ${db}/script/compare.R --input result/picrust/${l}.mat.txt --design result/metadata.txt --group Group --compare ${compare} --threshold 0 --method wilcox --pvalue 0.05 --fdr 0.2 --output result/picrust/
然后,基于KEGG预测函数绘制微生物群落丰度分布图,将结果${compare}过滤成“.txt 文件”。执行以下代码生成柱状图:
Rscript ${db}/script/compare_hierarchy_facet.R --input result/picrust/${compare}.txt --data MeanA --annotation result/picrust/${l}.anno.txt --output result/picrust/${compare}.MeanA.bar.pdf
接下来,运行以下代码来生成描绘两组之间显著差异的柱状图,并根据更高的分类级别进行面划分:
Rscript ${db}/script/compare_hierarchy_facet2.R \
--input result/picrust/${compare}.txt \
--pvalue 0.05 --fdr 0.1 \
--annotation result/picrust/${l}.anno.txt \
--output result/picrust/${compare}.bar.pdf
方法2(本地操作-仅限于Linux系统,需安装picrust环境):首先启动Linux子系统并安装Conda。使用conda install picrust安装特征预测需要的软件环境。为了与下游分析兼容,请将您的OTU表转换为通用格式。此方法将在“Linux系统下的分析”部分详细说明。
FAPROTAX
视频教程链接:https://youtu.be/kAiMq08mLNs
FAPROTAX数据库(http://www.loucalab.com/archive/FAPROTAX)可基于16S rDNA基因扩增子数据预测潜在功能。它提供了一个免费的Python脚本(collapse_table.py),可将分类学概况(如OTU表)转换为潜在的功能概况。FAPROTAX十分全面,包含4,600多个分类群的7,600多个功能注释,涵盖80多种功能,以及代谢途径和生态相关过程,如硝化、反硝化和发酵。该资源将特定的原核生物分类群(例如属或种)映射到功能上,帮助研究人员揭示微生物群落的代谢潜力。
方法1:ImageGP作为一键式的在线网络平台,可使用FAPROTAX对微生物群落进行功能注释。将OTU表上传到ImageGP在线网络服务器(https://www.bic.ac.cn/BIC/#/),它将可以进行FAPROTAX分析(图 S1-D),提供全面的结果和可视化。
方法2(Linux系统下分析,可选):详细信息可在“Linux系统下的分析”部分中找到。这种方法为Linux经验丰富的用户提供了灵活选择。
Bugbase细菌表型预测
视频教程链接:https://youtu.be/OvkZxjOjwmE
BugBase是一个分析微生物组样本表型的工具。该网站可以根据OTU表和映射文件(Mapping files)预测和比较大量的信息,包括以下七个方面:革兰氏阳性、革兰氏阴性、生物膜形成、致病潜力、移动元件含量、氧利用和氧化胁迫耐受。研究中通常使用BugBase基于16S预测样品中革兰氏阳性/阴性、胁迫和致病菌的相对丰度。
方法1(Window本地运行):基于R 4.x更新代码和包。在pipeline.sh中指定软件目录作为Bugbase变量,输入gg OTU表、元数据,指定分组列名和输出目录。
cd ${wd}/result
bugbase=${db}/script/BugBase
rm -rf bugbase/
Rscript ${bugbase}/bin/run.bugbase.r -L ${bugbase} -i gg/otutab.txt -m metadata.txt -c Group -o bugbase/
方法2:或者,用户可以使用在线网络服务器通过ImageGP进行Bugbase分析(http://www.bic.ac.cn/ImageGP/index.php/Home/Index/BugBase.html) (图 S1-E)。此外,用户还可以访问Bugbase网站https://bugbase.cs.umn.edu/进行分析。但是,由于其容易出错,我们并不建议这样做。
方法3(Linux系统下分析,可选):Linux系统下Bugbase的安装和使用详细说明请参见“Linux系统下分析”部分。
高级系统发育树的构建
视频教程链接:https://youtu.be/c3cM-zwbChs
本节介绍如何使用各种方案构建和可视化系统发育树。通过执行以下代码,在“result”文件夹中创建一个名为“tree”的新文件夹来存储生成的文件(图 4F):
cd ${wd}
mkdir -p result/tree
cd ${wd}/result/tree
筛选高丰度/特定特征
方法1:基于丰度筛选特征。通常选择0.001或0.005,且OTU计数在30-150的范围内。使用代码tail -n+2 ../otutab_rare.txt | wc -l计算特征表中ASV的个数。然后运行usearch -otutab_trim ../otutab_rare.txt -min_otu_freq 0.002 -output otutab.txt命令,根据相对丰度阈值0.2%筛选高丰度的OTU,最后使用 tail -n+2 otutab.txt | wc -l计算筛选后的OTU表中的特征数量。
方法2:按数量/计数筛选。首先按丰度排序,默认为降序排序(从大到小):
usearch -otutab_sortotus ../otutab_rare.txt \
-output otutab_sort.txt
运行以下代码,提取指定数量的高丰度OTU,例如前100个:
sed '1 s/#OTU ID/OTUID/' otutab_sort.txt \
| head -n101 > otutab.txt
筛选完成后执行sed -i '1 s/#OTU ID/OTUID/' otutab.txt修改特征ID的列名。运行代码cut -f 1 otutab.txt > otutab_high.id提取用于序列检索的ID。然后,通过运行以下代码筛选高丰度细菌/指定差异细菌对应的OTU序列:
usearch -fastx_getseqs ../otus.fa -labels otutab_high.id -fastaout otus.fa
head -n 2 otus.fa
基于物种注释筛选OTU:
awk 'NR==FNR{a[$1]=$0} NR>FNR{print a[$1]}' ../taxonomy.txt \
otutab_high.id > otutab_high.tax
获得与OTU相对应的组平均值,用于生成样本热图(依赖于之前otu_mean.R计算过的按Group分组的均值):
awk 'NR==FNR{a[$1]=$0} NR>FNR{print a[$1]}' ../otutab_mean.txt otutab_high.id \
| sed 's/#OTU ID/OTUID/' > otutab_high.mean
将物种注释和丰度合并到注释文件中:
cut -f 2- otutab_high.mean > temp
paste otutab_high.tax temp > annotation.txt
head -n 3 annotation.txt
构建进化树
起始文件是result/tree目录下的otus.fa(序列),annotation.txt(物种和相对丰度)。执行命令muscle -in otus.fa -out otus_aligned.fas调用软件Muscle进行序列比对。
方法1:使用IQ-TREE快速构建ML进化树(耗时2m)
rm -rf iqtree
mkdir -p iqtree
iqtree -s otus_aligned.fas -bb 1000 -redo -alrt 1000 -nt AUTO -pre iqtree/otus
方法2:在Ubuntu上使用apt install fasttree安装FastTree,非常适合大型系统发育树(数百个OTU)的构建。通过执行代码fasttree -gtr -nt otus_aligned.fas > otus.nwk构建树。
进化树美化
本节概述了如何自定义系统发育树的外观,重点关注物种注释和丰度信息的外圈着色以及塑造。利用之前我们生成的annotation.txt文件包含与OTU对应的物种注释和丰度值,可以使用以下几个方案来进行美化:
方案1:外圈的颜色、形状和丰度注释文件
cd ${wd}/result/tree
Rscript ${db}/script/table2itol.R -a -c double -D plan1 -i OTUID -l Genus -t %s -w 0.5 annotation.txt
方案2:丰度柱状图注释文件
Rscript ${db}/script/table2itol.R -a -d -c none -D plan2 -b Phylum -i OTUID -l Genus -t %s -w 0.5 annotation.txt
方案3:热图注释文件
Rscript ${db}/script/table2itol.R -c keep -D plan3 -i OTUID -t %s otutab.txt
方案4:将整数转换为因子以生成注释文件
Rscript ${db}/script/table2itol.R -a -c factor -D plan4 -i OTUID -l Genus -t %s -w 0 annotation.txt
完成方案文件之后访问网站http://itol.embl.de/,在分析结果的“tree”文件夹中找到需要美化的进化树文件(通常命名为“otus.contree”)并将该文件拖放到iTOL上传区域。同时上传方案文件(plan1/2/3/4)进一步定制树形(例如,突出显示特定的分支)。iTOL可通过各种选项来调整树的外观,用户可根据喜好更改选项从而美化进化树。美化完成后,将进化树保存在iTOL中(图 4F)。或者用户可前往另一个在线可视化平台https://www.bic.ac.cn/BIC/#/,同样用户友好。
Linux下的附加分析(可选)
EasyAmplicon的这一部分提供了在Linux服务器上运行LEfSe、FAPROTAX和PICRUSt函数等分析的说明。熟悉命令行操作的用户会发现这些步骤很简单。但是,需要安装特定的软件,包括Linux系统、Conda包和QIIME2。虽然PICRUSt2是预测宏基因组功能内容的可选工具,但需要注意的是它资源较大。在开始之前,请确保您的终端从“Git Bash”切换到“在Linux下的Windows子系统(Windows Subsystem for Linux)”。请记住,在Linux上运行分析可能需要具有足够内存的服务器环境(>16GB)。为了可视化LEfSe分析,建议使用前面提到的在线网络服务器。