内部培训资料（2022.10.29）

QIIME 2分析实例--人类肠癌不同部位微生物多样性分析

“Moving Pictures” tutorial

熊金波实验室出品 ⚠️⚠️⚠️内部交流培训使用

宁波大学海洋学院 Larry 陆

本节1.6万字，14张图。阅读时间大约50分钟。

在本次分享中，我们将使用QIIME2软件进行微生物数据分析。在五个时间点对来自两个人四个身体部位的微生物组样本进行分析，第一个时间点取样之后使用了抗生素处理。基于这批数据的研究可以参考2011年发表在Genome Biology的《Moving pictures of the human microbiome》（10分+的文章哦）。本次分享将会完整完成文章中的所有扩增子数据分析。本教程中使用的数据基于Illumina HiSeq产出，使用地球微生物组计划扩增16S rRNA基因高变区4（V4）测序的方法。

安装

Qiime一直以来被诟病最多的就是非常难以安装，目前据我了解有三种主流安装方法，根据环境任选其一即可，实际的安装涉及到服务器运维工作，与我们实验室日常工作相关性不大，所以此处不做过多介绍。我们205实验室的的T430服务器和Biostar服务器都会部署好QIIME1代和2代，可以直接使用。如果对于个人PC有安装需求的同学，可以以后向我了解相关材料。

启动QIIME2运行环境

我们在服务器上进行生物信息数据分析的时候非常重要的就是有良好的运行环境习惯，有了良好的环境管理习惯可以避免大量报错和数据分析问题。

我们在每次分析开始前，必须先进入工作目录，除非你是一个把什么东西都放在桌面上还很工作更有效率的人😂。

# 定义工作目录变量，方便以后多次使用
wd=~/github/QIIME2ChineseManual/2019.7
mkdir -p $wd
# 进入工作目录，是不是很简介，这样无论你在什么位置就可以快速回到项目文件夹
cd $wd

# 方法1\. 进入QIIME 2 conda工作环境
conda activate qiime2-2019.7
# 这时我们的命令行前面出现 (qiime2-2019.7) 表示成功进入工作环境

# 方法2\. conda版本较老用户，使用source进入QIIME 2
source activate qiime2-2019.7

# 方法3\. 如果是docker安装的请运行如下命令，默认加载当前目录至/data目录
docker run --rm -v $(pwd):/data --name=qiime -it  qiime2/core:2019.7

# 创建本节学习目录
mkdir qiime2-moving-pictures-tutorial
cd qiime2-moving-pictures-tutorial

激活分析环境

qiime2

创建个人文件夹

mkdir -p larry

首先获取分析的数据

cp -r raw-data larry
cp TruSeq3-PE.fa larry

首先获取分析样本元数据信息表格

cp  mapping_file.txt larry

正式开始分析数据要先对数据进行指控流程

激活软件所在的环境

cd larry
conda activate /biostack/bioconda

创建质控文件夹

mkdir -p data_analysis/trimming/reads

参数解析 :

使用 [trimmomatic](http://www.usadellab.org/cms/?page=trimmomatic) 去除拆分序列中可能可能存在的接头和切掉低质量的碱基， 针对一些小片段扩增产物，测序数据可能包含接头序列。

命令模式:  SE指定单端数据 ，PE指定双端数据，对于PE模式: 2个输入文件（正向和反向reads）和4个输出文件（正向配对、正向未配对、反向配对和反向未配对）
threads： 设定线程数
phred<qual>： 设定碱基质量编码模式，两种模式 `-phred33` 或 `-phred64`， 如果未指定质量编码，则将自动识别, `fastq` 质量值编码可以参考：[https://wiki.bits.vib.be/index.php/FASTQ](https://wiki.bits.vib.be/index.php/FASTQ)

ILUMINACLIP:<fastaWithAdaptersEtc>:<seed mismatches>:<palindrome clip threshold>:<simple clip threshold>:<minAdapterLength>:<keepBothReads>： 该步骤用于寻找并去除Illumina接头    

    <fastaWithAdaptersEtc1>:     指定包含所有接头、PCR序列等的fasta文件的路径;
    <seed mismatches>：          指定允许执行完全匹配的最大不匹配数；
    <palindrome clip threshold>：指定两个成对接头reads之间的匹配对于双端回文read比对的精度；
    <Simple clip threshold>：    指定任何接头等序列与read之间的匹配精度阈值；
    <minAdapterLength>:          设定接头的最小长度。如果未指定，默认为8个bp;
    <keepBothReads>：            在回文模式检测到read测穿并删除接头序列后，反向读取包含与正向读取相同的序列信息。因此，默认行为是完全删除反向读取。通过为该参数指定true，反向读取也将被保留

SLIDINGWINDOW:<windowSize>:<requiredQuality>: 滑窗修剪, 它从5'端开始扫描，当窗口内的平均质量低于阈值时，剔除该窗口内的所有碱基;

    <windowSize>:       设定窗口大小, 覆盖碱基数量；
    <requiredQuality>:  设定平均质量。

HEADCROP:<length> ：切除read起始端低于阈值的碱基

    <length> 从`read` 起始端开始要切除的长度

TRAILING:<quality>：从末端移除低质量的碱基。只要碱基的质量值低于阈值，则切除该碱基，并调查下一个碱基（因为Trimmomatic从3'primeend开始，将是位于刚切除碱基之前的碱基）。 此方法可用于去除Illumina低质量段区域（质量分数标记为2）, 
 
    <quality>: 指定保留碱基所需的最低质量
 
MINLEN:<length>: 设置保留reads的最小长度。
    <length>:可被保留的最短 read 长度。

trimmomatic PE -threads 4 -phred33 \
raw-data/A1_1.fastq \
raw-data/A1_2.fastq \
data_analysis/trimming/reads/A1_1.fastq \
data_analysis/trimming/reads/A1_1.singleton.fastq \
data_analysis/trimming/reads/A1_2.fastq \
data_analysis/trimming/reads/A1_2.singleton.fastq \
ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 SLIDINGWINDOW:4:15 LEADING:3 TRAILING:3 MINLEN:36

trimmomatic PE -threads 4 -phred33 \
raw-data/A2_1.fastq \
raw-data/A2_2.fastq \
data_analysis/trimming/reads/A2_1.fastq \
data_analysis/trimming/reads/A2_1.singleton.fastq \
data_analysis/trimming/reads/A2_2.fastq \
data_analysis/trimming/reads/A2_2.singleton.fastq \
ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 SLIDINGWINDOW:4:15 LEADING:3 TRAILING:3 MINLEN:36

trimmomatic PE -threads 4 -phred33 \
raw-data/D1_1.fastq \
raw-data/D1_2.fastq \
data_analysis/trimming/reads/D1_1.fastq \
data_analysis/trimming/reads/D1_1.singleton.fastq \
data_analysis/trimming/reads/D1_2.fastq \
data_analysis/trimming/reads/D1_2.singleton.fastq \
ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 SLIDINGWINDOW:4:15 LEADING:3 TRAILING:3 MINLEN:36

trimmomatic PE -threads 4 -phred33 \
raw-data/D2_1.fastq \
raw-data/D2_2.fastq \
data_analysis/trimming/reads/D2_1.fastq \
data_analysis/trimming/reads/D2_1.singleton.fastq \
data_analysis/trimming/reads/D2_2.fastq \
data_analysis/trimming/reads/D2_2.singleton.fastq \
ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 SLIDINGWINDOW:4:15 LEADING:3 TRAILING:3 

trimmomatic PE -threads 4 -phred33 \
raw-data/H1_1.fastq \
raw-data/H1_2.fastq \
data_analysis/trimming/reads/H1_1.fastq \
data_analysis/trimming/reads/H1_1.singleton.fastq \
data_analysis/trimming/reads/H1_2.fastq \
data_analysis/trimming/reads/H1_2.singleton.fastq \
ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 SLIDINGWINDOW:4:15 LEADING:3 TRAILING:3 MINLEN:36

trimmomatic PE -threads 4 -phred33 \
raw-data/H2_1.fastq \
raw-data/H2_2.fastq \
data_analysis/trimming/reads/H2_1.fastq \
data_analysis/trimming/reads/H2_1.singleton.fastq \
data_analysis/trimming/reads/H2_2.fastq \
data_analysis/trimming/reads/H2_2.singleton.fastq \
ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 SLIDINGWINDOW:4:15 LEADING:3 TRAILING:3 MINLEN:36

将双端数据进行合并

创建目录 :

mkdir -p data_analysis/mergepairs/reads

合并双端序列 :

使用usearch [-fastq_mergepairs] 命令合并A1的双端序列，并输出结果到data_analysis/mergepairs/reads/目录下:

usearch -fastq_mergepairs参数解析

-fastq_mergepairs       正向FASTQ文件名（输入）
-reverse                反向FASTQ文件名（输入）
-fastqout               合并后的FASTQ文件名（输出）
-fastqout_notmerged_fwd 未合并的正向序列FASTQ文件名（输出）
-fastqout_notmerged_rev 未合并的反向序列FASTQ文件名（输出）
-log                    log文件名（输出）
-report                 总结报告文件名（输出）
-threads                线程数
-fastq_minmergelen      合并序列的最小长度
-fastq_maxmergelen      合并序列的最大长度  
-fastq_maxdiffs         最大错配数，默认为5
-fastq_pctid            最小序列质量，默认为90%
-fastq_minovlen         合并后的序列如果低于指定值则过滤，默认为16
-fastq_trunctail        在第一个小于指定Q值处截断序列，默认为
-fastq_minlen           如果-fastq_trunctail截断后的序列小于指定值，则过滤掉序列

usearch -fastq_mergepairs data_analysis/trimming/reads/A1_1.fastq \
-reverse data_analysis/trimming/reads/A1_2.fastq \
-fastqout data_analysis/mergepairs/reads/A1.fastq \
-fastqout_notmerged_fwd  data_analysis/mergepairs/reads/A1.notmerged_fwd.fastq \
-fastqout_notmerged_rev  data_analysis/mergepairs/reads/A1.notmerged_rev.fastq \
-log  data_analysis/mergepairs/reads/A1.log \
-report data_analysis/mergepairs/reads/A1.report \
-threads 1 \
-fastq_minmergelen 250 \
-fastq_maxmergelen 500 \
-fastq_maxdiffs 10  \
-fastq_pctid 90  \
-fastq_minovlen 16   \
-fastq_trunctail 2 \
-fastq_minlen 64

usearch -fastq_mergepairs data_analysis/trimming/reads/A2_1.fastq \
-reverse data_analysis/trimming/reads/A2_2.fastq \
-fastqout data_analysis/mergepairs/reads/A2.fastq \
-fastqout_notmerged_fwd  data_analysis/mergepairs/reads/A2.notmerged_fwd.fastq \
-fastqout_notmerged_rev  data_analysis/mergepairs/reads/A2.notmerged_rev.fastq \
-log  data_analysis/mergepairs/reads/A2.log \
-report data_analysis/mergepairs/reads/A2.report \
-threads 1 \
-fastq_minmergelen 250 \
-fastq_maxmergelen 500 \
-fastq_maxdiffs 10  \
-fastq_pctid 90  \
-fastq_minovlen 16   \
-fastq_trunctail 2 \
-fastq_minlen 64

usearch -fastq_mergepairs data_analysis/trimming/reads/D1_1.fastq \
-reverse data_analysis/trimming/reads/D1_2.fastq \
-fastqout data_analysis/mergepairs/reads/D1.fastq \
-fastqout_notmerged_fwd  data_analysis/mergepairs/reads/D1.notmerged_fwd.fastq \
-fastqout_notmerged_rev  data_analysis/mergepairs/reads/D1.notmerged_rev.fastq \
-log  data_analysis/mergepairs/reads/D1.log \
-report data_analysis/mergepairs/reads/D1.report \
-threads 1 \
-fastq_minmergelen 250 \
-fastq_maxmergelen 500 \
-fastq_maxdiffs 10  \
-fastq_pctid 90  \
-fastq_minovlen 16   \
-fastq_trunctail 2 \
-fastq_minlen 64

usearch -fastq_mergepairs data_analysis/trimming/reads/D2_1.fastq \
-reverse data_analysis/trimming/reads/D2_2.fastq \
-fastqout data_analysis/mergepairs/reads/D2.fastq \
-fastqout_notmerged_fwd  data_analysis/mergepairs/reads/D2.notmerged_fwd.fastq \
-fastqout_notmerged_rev  data_analysis/mergepairs/reads/D2.notmerged_rev.fastq \
-log  data_analysis/mergepairs/reads/D2.log \
-report data_analysis/mergepairs/reads/D2.report \
-threads 1 \
-fastq_minmergelen 250 \
-fastq_maxmergelen 500 \
-fastq_maxdiffs 10  \
-fastq_pctid 90  \
-fastq_minovlen 16   \
-fastq_trunctail 2 \
-fastq_minlen 64

usearch -fastq_mergepairs data_analysis/trimming/reads/H1_1.fastq \
-reverse data_analysis/trimming/reads/H1_2.fastq \
-fastqout data_analysis/mergepairs/reads/H1.fastq \
-fastqout_notmerged_fwd  data_analysis/mergepairs/reads/H1.notmerged_fwd.fastq \
-fastqout_notmerged_rev  data_analysis/mergepairs/reads/H1.notmerged_rev.fastq \
-log  data_analysis/mergepairs/reads/H1.log \
-report data_analysis/mergepairs/reads/H1.report \
-threads 1 \
-fastq_minmergelen 250 \
-fastq_maxmergelen 500 \
-fastq_maxdiffs 10  \
-fastq_pctid 90  \
-fastq_minovlen 16   \
-fastq_trunctail 2 \
-fastq_minlen 64

usearch -fastq_mergepairs data_analysis/trimming/reads/H2_1.fastq \
-reverse data_analysis/trimming/reads/H2_2.fastq \
-fastqout data_analysis/mergepairs/reads/H2.fastq \
-fastqout_notmerged_fwd  data_analysis/mergepairs/reads/H2.notmerged_fwd.fastq \
-fastqout_notmerged_rev  data_analysis/mergepairs/reads/H2.notmerged_rev.fastq \
-log  data_analysis/mergepairs/reads/H2.log \
-report data_analysis/mergepairs/reads/H2.report \
-threads 1 \
-fastq_minmergelen 250 \
-fastq_maxmergelen 500 \
-fastq_maxdiffs 10  \
-fastq_pctid 90  \
-fastq_minovlen 16   \
-fastq_trunctail 2 \
-fastq_minlen 64

** 统计序列信息 :**

使用atlas-utils [fqchk] 统计序列质量

cat data_analysis/mergepairs/reads/A1.fastq \
|atlas-utils \
fqchk -q 33 -l A1 \
- >data_analysis/mergepairs/reads/A1.txt

cat data_analysis/mergepairs/reads/A2.fastq \
|atlas-utils \
fqchk -q 33 -l A2 \
- >data_analysis/mergepairs/reads/A2.txt

cat data_analysis/mergepairs/reads/D1.fastq \
|atlas-utils \
fqchk -q 33 -l D1 \
- >data_analysis/mergepairs/reads/D1.txt

cat data_analysis/mergepairs/reads/D2.fastq \
|atlas-utils \
fqchk -q 33 -l D2 \
- >data_analysis/mergepairs/reads/D2.txt

cat data_analysis/mergepairs/reads/H1.fastq \
|atlas-utils \
fqchk -q 33 -l H1 \
- >data_analysis/mergepairs/reads/H1.txt

cat data_analysis/mergepairs/reads/H2.fastq \
|atlas-utils \
fqchk -q 33 -l H2 \
- >data_analysis/mergepairs/reads/H2.txt

统计信息 :

合并所有mergepairs/reads/目录下后缀为.txt的文件并传递给tsv-utils view 去除中间的注释行, 使用usearch-utils mergepairs统计log文件

mkdir -p data_analysis/mergepairs/report

cat data_analysis/mergepairs/reads/*.txt | tsv-utils view - >data_analysis/mergepairs/report/sequencing.stats.txt

比对引物信息并去除

创建引物信息目录

mkdir -p data_analysis/primer_strip/reads

usearch -search_pcr2 data_analysis/mergepairs/reads/A1.fastq \
-fwdprimer ACTCCTACGGGAGGCAGCA \
-revprimer GGACTACHVGGGTWTCTAAT \
-minamp 250 \
-maxamp 480 \
-strand both -maxdiffs 2 \
-tabbedout data_analysis/primer_strip/reads/A1.txt \
-fastqout  data_analysis/primer_strip/reads/A1.fastq \
-log data_analysis/primer_strip/reads/A1.log

usearch -search_pcr2 data_analysis/mergepairs/reads/A2.fastq \
-fwdprimer ACTCCTACGGGAGGCAGCA \
-revprimer GGACTACHVGGGTWTCTAAT \
-minamp 250 \
-maxamp 480 \
-strand both -maxdiffs 2 \
-tabbedout data_analysis/primer_strip/reads/A2.txt \
-fastqout  data_analysis/primer_strip/reads/A2.fastq \
-log data_analysis/primer_strip/reads/A2.log

usearch -search_pcr2 data_analysis/mergepairs/reads/D1.fastq \
-fwdprimer ACTCCTACGGGAGGCAGCA \
-revprimer GGACTACHVGGGTWTCTAAT \
-minamp 250 \
-maxamp 480 \
-strand both -maxdiffs 2 \
-tabbedout data_analysis/primer_strip/reads/D1.txt \
-fastqout  data_analysis/primer_strip/reads/D1.fastq \
-log data_analysis/primer_strip/reads/D1.log

usearch -search_pcr2 data_analysis/mergepairs/reads/H1.fastq \
-fwdprimer ACTCCTACGGGAGGCAGCA \
-revprimer GGACTACHVGGGTWTCTAAT \
-minamp 250 \
-maxamp 480 \
-strand both -maxdiffs 2 \
-tabbedout data_analysis/primer_strip/reads/H1.txt \
-fastqout  data_analysis/primer_strip/reads/H1.fastq \
-log data_analysis/primer_strip/reads/H1.log

usearch -search_pcr2 data_analysis/mergepairs/reads/D2.fastq \
-fwdprimer ACTCCTACGGGAGGCAGCA \
-revprimer GGACTACHVGGGTWTCTAAT \
-minamp 250 \
-maxamp 480 \
-strand both -maxdiffs 2 \
-tabbedout data_analysis/primer_strip/reads/D2.txt \
-fastqout  data_analysis/primer_strip/reads/D2.fastq \
-log data_analysis/primer_strip/reads/D2.log

usearch -search_pcr2 data_analysis/mergepairs/reads/H2.fastq \
-fwdprimer ACTCCTACGGGAGGCAGCA \
-revprimer GGACTACHVGGGTWTCTAAT \
-minamp 250 \
-maxamp 480 \
-strand both -maxdiffs 2 \
-tabbedout data_analysis/primer_strip/reads/H2.txt \
-fastqout  data_analysis/primer_strip/reads/H2.fastq \
-log data_analysis/primer_strip/reads/H2.log

导入数据

首先需要构建数据信息表格

mkdir -p  data_analysis/dada2/report
mkdir -p data_analysis/dada2/qza/

用于输入到QIIME 2的所有数据都以QIIME 2对象的形式出现，其中包含有关数据类型和数据源的信息。因此，我们需要做的第一件事是将这些序列数据文件导入到QIIME 2对象中。
开始导入数据

qiim-utils manifest mapping_file.txt \
/project/qiime-tk/larry/data_analysis/primer_strip/reads \
> data_analysis/dada2/report/manifest.txt

这个QIIME 2对象的语义类型是SequencesWithQuality。 QIIME 对象SequencesWithQuality是带有质量序列的样本，。要导入其他格式的序列数据，请参阅导入数据教程。

time qiime tools import --type 'SampleData[SequencesWithQuality]' \
--input-path  data_analysis/dada2/report/manifest.txt \
--output-path data_analysis/dada2/qza/demux.qza \
--input-format SingleEndFastqManifestPhred33V2

结果统计

time qiime demux summarize \
  --i-data data_analysis/dada2/qza/demux.qza \
  --o-visualization data_analysis/dada2/qza/demux.qzv

image

2019.7.4.01.jpg918×1877](http://bailab.genetics.ac.cn/markdown/qiime2/fig/2019.7.4.01.jpg "2019.7.4.01.jpg")

图1. 样本拆分结果统计结果——样本数据量可视化图表。

主要分为三部分：上部为摘要；中部为样本不同数据量分布频率柱状图，可下载PDF，下部为每个样本的测序量。上方面板还可切换至交互式质量图Interactive Qaulity Plot页面。如下图2。（此处再次点赞👍QIIME 2的强大交互性，几乎所有的可视化文件都可以最快速度最低成本形成发表级可视化文件，大大节约时间）

image

2019.7.4.02.gif838×1412](http://bailab.genetics.ac.cn/markdown/qiime2/fig/2019.7.4.02.gif "2019.7.4.02.gif")

图2. 交互式质量图Interactive Qaulity Plot查看页面。

同样为三部分：上部为每个位置碱基的质量分布交互式箱线图，鼠标悬停在上面，即可在下面(中部)文字和表格中显示鼠标所在位置碱基质量的详细信息；下部为拆分样本的长度摘要(一般等长测序无差别)。

qiime tools view demux.qzv

直接在服务器操作的话可以直接利用这条命令可视化，这条命令的显示需要图形界面的支持，如在有图型界面的Linux上，但仅使用SSH登陆方式无法显示图形。

版本区别提示：相比QIIME 1.9.1
在QIIME 1.9.1中，我们一般建议通过QIIME执行样本拆分（例如，使用split_libraries.py或split_libraries_fastq.py这两个脚本），因为这个步骤还执行序列的质量控制👍。现在我们将样本拆分和质量控制步骤分开，因此你可以使用混合多样本序列（如我们在此所做的）或拆分后的序列开始QIIME 2分析。

在样本拆分之后，生成拆分结果的统计信息非常重要。统计结果允许我们确定每个样本获得多少序列，并且还可以获得序列数据中每个位置处序列质量分布的摘要。

推荐使用 https://view.qiime2.org 网址显示结果

可选使用XShell+XManager支持SSH方式的图型界面、虚拟机图形界面下或服务器远程桌面方式支持上面命令的图形结果。我们上学期的Linux相关知识介绍里面有讲到过这个相关知识。以后有机会的话，会在宏基因组流程化数据分析里面再讲一下。

目前命令行方式想要查看结果可能很多使用服务器人员无法实现 (即依赖服务器安装了桌面，本地依赖XShell+XManager或其它ssh终端和图形界面软件)

本地查看可解压.qzv，目录中的data目录包括详细的图表文件，主要关注 pdf 和 html 文件，目录结构如下。

── demux
   └── 8743ab13-72ca-4adf-9b6c-d97e2dbe8ee3
       ├── checksums.md5
       ├── data
       │   ├── data.jsonp
       │   ├── demultiplex-summary.pdf
       │   ├── demultiplex-summary.png
       │   ├── dist
       │   │   ├── bundle.js
       │   │   ├── d3-license.txt
       │   │   └── vendor.bundle.js
       │   ├── forward-seven-number-summaries.csv
       │   ├── index.html
       │   ├── overview.html
       │   ├── per-sample-fastq-counts.csv
       │   ├── q2templateassets
       │   │   ├── css
       │   │   │   ├── bootstrap.min.css
       │   │   │   ├── normalize.css
       │   │   │   └── tab-parent.css
       │   │   ├── fonts
       │   │   │   ├── glyphicons-halflings-regular.eot
       │   │   │   ├── glyphicons-halflings-regular.svg
       │   │   │   ├── glyphicons-halflings-regular.ttf
       │   │   │   ├── glyphicons-halflings-regular.woff
       │   │   │   └── glyphicons-halflings-regular.woff2
       │   │   ├── img
       │   │   │   └── qiime2-rect-200.png
       │   │   └── js
       │   │       ├── bootstrap.min.js
       │   │       ├── child.js
       │   │       ├── jquery-3.2.0.min.js
       │   │       └── parent.js
       │   └── quality-plot.html
       ├── metadata.yaml
       ├── provenance
       │   ├── action
       │   │   └── action.yaml
       │   ├── artifacts
       │   │   ├── 9594ef07-c414-4658-9345-c726de100d8d
       │   │   │   ├── action
       │   │   │   │   └── action.yaml
       │   │   │   ├── citations.bib
       │   │   │   ├── metadata.yaml
       │   │   │   └── VERSION
       │   │   └── a7a882f3-5e4f-4b5e-8a35-6a1098d21608
       │   │       ├── action
       │   │       │   ├── action.yaml
       │   │       │   └── barcodes.tsv
       │   │       ├── citations.bib
       │   │       ├── metadata.yaml
       │   │       └── VERSION
       │   ├── citations.bib
       │   ├── metadata.yaml
       │   └── VERSION
       └── VERSION

qzv文件解压后文件详细，可直接访问data/index.html打开结果报告式网页。里面的重要结果，全部可以通过此网页进行索引。总而言之最重要的文件其实就是pdf格式文件和html格式文件，这集中文件是我们日常工作中真正会使用到的文件。

序列质控和生成特征表

Sequence quality control and feature table construction

QIIME 2插件多种质量控制方法可选，包括DADA2、Deblur和基于基本质量分数的过滤。在本教程中，我们使用DADA2和Deblur两种方法分别介绍这个步骤。这些步骤是可互相替换的，因此你可以使用自己喜欢的方法。这两种方法的结果将是一个QIIME 2特征表FeatureTable[Frequency]和一个代表性序列FeatureData[Sequence]对象，Frequency对象包含数据集中每个样本中每个唯一序列的计数（频率），Sequence对象将FeatureTable中的特征ID与序列对应。

提醒⏰：此步主要有DADA2和Deblur两种方法可选，推荐使用DADA2，2016年发表在Nature Method上，在阴道菌群研究中比OTU聚类结果看到更多细节；相较USEARCH的UPARSE算法，目前DADA2方法仅去噪去嵌合，不再按相似度聚类，结果与真实物种的序列更接近。

注意⚠️：本节中此次存在两种可选方法时，将创建具有特定方法名称的对象（例如，使用dada2去噪生成的特性表将被称为table.qza）。在创建这些对象之后，你将把两个选项之一的对象重命名为更通用的文件名（例如，table.qza）。为对象创建特定名称，然后对其进行重命名的过程仅允许你选择在本步骤中使用的两个选项中之一完成案例，而不必再次关注该选项。需要注意的是，在这个步骤或QIIME 2中的任何步骤中，你给对象或可视化的文件命名并不重要。

方法1. DADA2

Option 1: DADA2

DADA2是用于检测和校正（如果有可能的话）Illumina扩增序列数据的工作流程。正如在q2-dada2插件中实现的，这个质量控制过程将过滤掉在测序数据中鉴定的任何phiX序列（通常存在于标记基因Illumina测序数据中，用于提高扩增子测序质量），并同时过滤嵌合序列。

DADA2是Susan P. Holmes团队于2016年发表于Nature Methods的文章，截止18年12月22号Google学术统计引用483次；DADA2自身也是一套在R语言中完整的扩增子分析流程

dada2 denoise-single方法需要两个用于质量过滤的参数：--p-trim-left m，它去除每个序列的前m个碱基(如引物、标签序列barcode)；--p-trunc-len n，它在位置n截断每个序列。这允许用户去除序列的低质量区域、引物或标签序列等。为了确定要为这两个参数传递什么值，你应该查看上面由qiime demux summarize生成的demux.qzv文件中的交互质量图选项卡。

单端序列去噪, 输入样本拆分后结果；去除左端 0 bp (–p-trim-left，有时用于切除低质量序列、barocde或引物)，序列切成 120 bp 长(–p-trunc-len)；生成代表序列、特征表和去噪过程统计。

下面的步骤计算量较大。

time qiime dada2 denoise-single --i-demultiplexed-seqs data_analysis/dada2/qza/demux.qza \
--p-trim-left 0 \
--p-trunc-len 0 \
--p-n-threads 0 \
--o-representative-sequences data_analysis/dada2/qza/rep-seqs.qza \
--o-table data_analysis/dada2/qza/table.qza \
--o-denoising-stats data_analysis/dada2/qza/stats.qza
 # 实际计算时间，即受服务器配置影响，还受同台服务器上任务量影响

生成三个输出文件：

• /stats.qza: dada2计算统计结果。
• table.qza: 特征表。
• rep-seqs.qza: 代表序列。

对特征表统计进行进行可视化

qiime metadata tabulate \
  --m-input-file data_analysis/dada2/qza/stats.qza \
  --o-visualization data_analysis/dada2/qza/stats.qzv

qiime tools export --input-path data_analysis/dada2/qza/table.qza \
--output-path data_analysis/dada2/denoise ;
qiime tools export --input-path data_analysis/dada2/qza/rep-seqs.qza \
--output-path data_analysis/dada2/denoise ;
qiime tools export --input-path  data_analysis/dada2/qza/stats.qza \
--output-path data_analysis/dada2/denoise ;
biom convert \
--input-fp  data_analysis/dada2/denoise/feature-table.biom \
--output-fp data_analysis/dada2/denoise/biom-table.tsv --to-tsv;
tail -n +2 data_analysis/dada2/denoise/biom-table.tsv \
|tsv-utils view -r \
- > data_analysis/dada2/denoise/feature-table.tsv

方法2. Deblur（个人不是很推荐，内存占用较大，时间慢，所以不演示了，看一下流程吧）

Deblur使用序列错误配置文件将错误的序列与从其来源的真实生物序列相关联，从而得到高质量的序列变异数据，主要为两个步骤。首先，应用基于质量分数的初始质量过滤过程，是Bokulich等人2013年发表的质量过滤方法。

Deblur是本软件作者作为通讯作者2013发表于Nature Methods的重要扩增子代表序列鉴定方法，截止19年8月25号Google学术统计引用1259次，

按测序碱基质量过滤序列

time qiime quality-filter q-score \
 --i-demux demux.qza \
 --o-filtered-sequences demux-filtered.qza \
 --o-filter-stats demux-filter-stats.qza
# 用时40s

输出对象:

• demux-filtered.qza: 序列质量过滤后结果。 查看 | 下载
• demux-filter-stats.qza: 序列质量过滤后结果统计。 查看 | 下载

接下来，使qiime deblur denoise-16S方法应用于Deblur工作流程。此方法需要一个用于质量过滤的参数，即截断位置n长度的序列的--p-trim-length n。通常，Deblur开发人员建议将该值设置为质量分数中位数开始下降至低质量区时的长度。在本次数据上，质量图（在质量过滤之前）表明合理的选择是在115至130序列位置范围内。这是一个主观的评估。你可能不采用该建议的一种原因是存在多个批次测序的元分析。在这种情况的元分析中，比较所有批次的序列长度是否相同，以避免人为引入特定的偏差，全局考虑这些是非常重要的。由于我们已经使用修剪长度为120 bp用于qiime dada2 denoise-single分析，并且由于120 bp是基于质量图的结果，这里我们将使用--p-trim-length 120参数。下一个命令可能需要10分钟才能运行完成。

deblur去噪16S过程，输入文件为质控后的序列，设置截取长度参数，生成结果文件有代表序列、特征表、样本统计。

time qiime deblur denoise-16S \
  --i-demultiplexed-seqs demux-filtered.qza \
  --p-trim-length 120 \
  --o-representative-sequences rep-seqs-deblur.qza \
  --o-table table-deblur.qza \
  --p-sample-stats \
  --o-stats deblur-stats.qza

• deblur-stats.qza: 过程统计。 查看 | 下载
• table-deblur.qza: 特征表。 查看 | 下载
• rep-seqs-deblur.qza: 代表序列。 查看 1 | 下载

time qiime metadata tabulate \
  --m-input-file demux-filter-stats.qza \
  --o-visualization demux-filter-stats.qzv
time qiime deblur visualize-stats \
  --i-deblur-stats deblur-stats.qza \
  --o-visualization deblur-stats.qzv

输出结果:

• demux-filter-stats.qzv: 质量过程统计表，同上面提到的stats-dada2.qzv统计表类似。 查看| 下载 1

示例如下：包括6列，第一列为样本名称，2-6列分别为总输入读长、总保留高读长、截断的读长、截断后太短的读长和超过最大模糊碱基的读长的数量统计。我们通常只关注2，3列数量即可，其它列常用于异常的输助判断。

• deblur-stats.qzv: deblur分析统计表，有分析中每个步骤的统计表 查看 | 下载

image

2019.7.4.03.jpg1870×760](http://bailab.genetics.ac.cn/markdown/qiime2/fig/2019.7.4.03.jpg "2019.7.4.03.jpg")

图3. deblur去噪和鉴定ASV处理过程统计结果

如果你想用此处结果下游分析，可以改名为下游分析的起始名称：

这处演示不运行下面两行代码，前面添加"#"号代表注释，需要运行请自行删除行首的“#”

#mv rep-seqs-deblur.qza rep-seqs.qza
#mv table-deblur.qza table.qza

💕统一完整的数据命名也是非常好的服务器运行习惯。

物种组成分析

Taxonomic analysis

在这一节中，我们将开始探索样本的物种组成，并将其与样本元数据再次组合。这个过程的第一步是为FeatureData[Sequence]的序列进行物种注释。我们将使用经过Naive Bayes分类器预训练的，并由q2-feature-classifier插件来完成这项工作。这个分类器是在silva-138-99上训练的，其中序列被修剪到仅包括来自16S区域的250个碱基，该16S区域在该分析中采用V4区域的515F/806R引物扩增并测序。我们将把这个分类器应用到序列中，并且可以生成从序列到物种注释结果关联的可视化。

注意：物种分类器根据你特定的样品制备和测序参数进行训练时表现最好，包括用于扩增的引物和测序序列的长度。因此，一般来说，你应该按照使用q2-feature-classifier的训练特征分类器的说明来训练自己的物种分类器。我们在数据资源页面上提供了一些通用的分类器，包括基于Silva的16S分类器，不过将来我们可能会停止提供这些分类器，而让用户训练他们自己的分类器，这将与他们的序列数据最相关。

mkdir -p data_analysis/classify/sklearn

提醒⏰：这一步会占用大量线程所以请大家不要操作，我这边演示后会提供在上层目录中

qiime feature-classifier classify-sklearn --i-classifier ../db/silva-138-99-nb-classifier.qza \
--i-reads data_analysis/dada2/qza/rep-seqs.qza \
--p-n-jobs 4 \
--o-classification data_analysis/classify/sklearn/taxonomy.qza 

qiime metadata tabulate \
--m-input-file data_analysis/classify/sklearn/taxonomy.qza \
--o-visualization data_analysis/classify/sklearn/taxonomy.qzv

我们可以在网页端口导入可视化文件来看一下具体的结果

导出物种注释信息

qiime tools export --input-path  data_analysis/classify/sklearn/taxonomy.qzv \
--output-path data_analysis/classify/taxonomy

基于此我们就得到了样品中的微生物信息可以了解样品中究竟有哪些物种在其中

qiime taxa barplot \
--i-table data_analysis/dada2/qza/table.qza \
--i-taxonomy data_analysis/classify/sklearn/taxonomy.qza \
--m-metadata-file mapping_file.txt \
--o-visualization data_analysis/classify/sklearn/taxa-bar-plots.qzv

qiime tools export \
--input-path data_analysis/classify/sklearn/taxa-bar-plots.qzv \
--output-path data_analysis/classify/barplot

构建进化树用于多样性分析

Generate a tree for phylogenetic diversity analyses

QIIME 2支持几种系统发育多样性度量方法，包括Faith’s Phylogenetic Diversity、weighted和unweighted UniFrac。除了每个样本的特征计数（即QIIME2对象FeatureTable[Frequency]）之外，这些度量还需要将特征彼此关联结合有根进化树。此信息将存储在一个QIIME 2对象的有根系统发育对象Phylogeny[Rooted]中。为了生成系统发育树🌲，我们将使用q2-phylogeny插件中的align-to-tree-mafft-fasttree工作流程。

首先，工作流程使用mafft程序执行对FeatureData[Sequence]中的序列进行多序列比对，以创建QIIME 2对象FeatureData[AlignedSequence]。接下来，流程屏蔽（mask或过滤）对齐的的高度可变区(高变区)，这些位置通常被认为会增加系统发育树的噪声。随后，流程应用FastTree基于过滤后的比对结果生成系统发育树。FastTree程序创建的是一个无根树，因此在本节的最后一步中，应用根中点法将树的根放置在无根树中最长端到端距离的中点，从而形成有根树。

mkdir -p data_analysis/phylogeny/qza

qiime phylogeny align-to-tree-mafft-fasttree --i-sequences data_analysis/dada2/qza/rep-seqs.qza \
--o-alignment data_analysis/phylogeny/qza/aligned-rep-seqs.qza \
--o-masked-alignment data_analysis/phylogeny/qza/masked-aligned-rep-seqs.qza \
--o-tree data_analysis/phylogeny/qza/unrooted-tree.qza \
--o-rooted-tree data_analysis/phylogeny/qza/rooted-tree.qza

qiime tools export --input-path data_analysis/phylogeny/qza/rooted-tree.qza \
--output-path data_analysis/phylogeny/tree

多样性分析

根据最小测序深度计算物种多样性信息

rm -rf data_analysis/metrics/diversity 
qiime diversity core-metrics-phylogenetic \
--i-phylogeny data_analysis/phylogeny/qza/rooted-tree.qza \
--i-table data_analysis/dada2/qza/table.qza \
--p-sampling-depth 6000 \
--m-metadata-file mapping_file.txt \
--output-dir data_analysis/metrics/diversity

获得特征序列表

qiime tools export --input-path data_analysis/metrics/diversity/rarefied_table.qza \
--output-path data_analysis/metrics/report ;
biom convert --input-fp  data_analysis/metrics/report/feature-table.biom \
--output-fp data_analysis/metrics/report/biom-table.tsv --to-tsv;
tail -n +2 data_analysis/metrics/report/biom-table.tsv \
|tsv-utils view -r \
- >data_analysis/metrics/report/feature_table.tsv 

Alpha和beta多样性分析

Alpha and beta diversity analysis
此处多样性分析的基础知识和相关方法众多，详细参考我在公众号平台上的总结，此处不论述。QIIME 2的多样性分析使用q2-diversity插件，该插件支持计算α和β多样性指数、并应用相关的统计检验以及生成交互式可视化图表。我们将首先应用core-metrics-phylogenetic方法，该方法将FeatureTable[Frequency](特征表[频率])抽平到用户指定的测序深度，然后计算几种常用的α和β多样性指数，并使用Emperor为每个β多样性指数生成主坐标分析（PCoA）图。默认情况下计算的方法有：

划重点：理解下面4种alpha和beta多样性指数的所代表的生物学意义至关重要⚠️。

• α多样性
  • 香农(Shannon’s)多样性指数（群落丰富度的定量度量，即包括丰富度richness和均匀度evenness两个层面）
  • 可观测的OTU(Observed OTUs，群落丰富度的定性度量，只包括丰富度）
  • Faith’s系统发育多样性（包含特征之间的系统发育关系的群落丰富度的定性度量）
  • 均匀度Evenness（或 Pielou’s均匀度；群落均匀度的度量）
• β多样性
  • Jaccard距离（群落差异的定性度量，即只考虑种类，不考虑丰度）
  • Bray-Curtis距离（群落差异的定量度量，较常用）
  • 非加权UniFrac距离（包含特征之间的系统发育关系的群落差异定性度量）
  • 加权UniFrac距离（包含特征之间的系统发育关系的群落差异定量度量）

需要提供给这个脚本的一个重要参数是--p-sampling-depth，它是指定重采样（即稀疏/稀疏rarefaction）深度。因为大多数多样指数对不同样本的不同测序深度敏感，所以这个脚本将随机地将每个样本的测序量重新采样至该参数值。例如，提供--p-sampling-depth 500，则此步骤将对每个样本中的计数进行无放回抽样，从而使得结果表中的每个样本的总计数为500。如果任何样本的总计数小于该值，那么这些样本将从多样性分析中删除。选择这个值很棘手。建议通过查看上面创建的表table.qzv文件中呈现的信息并选择一个尽可能高的值（因此每个样本保留更多的序列）同时尽可能少地排除样本来进行选择。

下面多样性分析，需要基于重采样/抽平(rarefaction)标准化的特征表，标准化采用无放回重抽样至序列一致，如何设计样品重采样深度参数--p-sampling-depth呢？
如是数据量都很大，选最小的即可。如果有个别数据量非常小，去除最小值再选最小值。比如此分析最小值为917，我们选择1109深度重采样，即保留了大部分样品用于分析，又去除了数据量过低的异常值。本示例为近10年前测序技术的通量水平，454测序时代抽平至1000条即可，现在看来数据量很小。目录一般采用HiSeq2500或NovaSeq6000的 PE250模式测序，数据量都非常大，通常可以采用3万或5万的标准抽平，仍可保留90%以上样本。过低或过高一般结果也会波动较大，不建议放在一起分析。

qiime tools export \
--input-path data_analysis/metrics/diversity/shannon_vector.qza \
--output-path data_analysis/alpha/diversity ;
mv data_analysis/alpha/diversity/alpha-diversity.tsv data_analysis/alpha/diversity/shannon.tsv ;
qiime tools export \
--input-path  data_analysis/metrics/diversity/faith_pd_vector.qza \
--output-path data_analysis/alpha/diversity ;

mv data_analysis/alpha/diversity/alpha-diversity.tsv data_analysis/alpha/diversity/faith_pd.tsv ;
qiime tools export \
--input-path  data_analysis/metrics/diversity/observed_features_vector.qza \
--output-path data_analysis/alpha/diversity ;
mv data_analysis/alpha/diversity/alpha-diversity.tsv data_analysis/alpha/diversity/observed_features.tsv ;

qiime tools export \
--input-path data_analysis/metrics/diversity/evenness_vector.qza \
--output-path data_analysis/alpha/diversity ;
mv data_analysis/alpha/diversity/alpha-diversity.tsv data_analysis/alpha/diversity/evenness.tsv

计算稀释度曲线

mkdir -p data_analysis/alpha/qza

qiime diversity alpha-rarefaction \
--p-min-depth 1 \
--p-steps 10  \
--p-iterations 10 \
--p-metrics shannon simpson chao1 observed_features goods_coverage faith_pd \
--i-table data_analysis/dada2/qza/table.qza \
--i-phylogeny data_analysis/phylogeny/qza/rooted-tree.qza \
--p-max-depth 7500 \
--m-metadata-file mapping_file.txt \
--o-visualization data_analysis/alpha/qza/alpha-rarefaction.qzv ;
qiime tools export \
--input-path  data_analysis/alpha/qza/alpha-rarefaction.qzv \
--output-path data_analysis/alpha/rarefaction ;

beta多样性

mkdir -p data_analysis/beta/matrix
qiime tools export \
--input-path data_analysis/metrics/diversity/jaccard_distance_matrix.qza \
--output-path data_analysis/beta/matrix ;
mv data_analysis/beta/matrix/distance-matrix.tsv data_analysis/beta/matrix/jaccard.tsv ;

qiime tools export --input-path  data_analysis/metrics/diversity/bray_curtis_distance_matrix.qza \
--output-path data_analysis/beta/matrix ;
mv data_analysis/beta/matrix/distance-matrix.tsv data_analysis/beta/matrix/bray_curtis.tsv ;

qiime tools export \
--input-path data_analysis/metrics/diversity/weighted_unifrac_distance_matrix.qza \
--output-path data_analysis/beta/matrix ;
mv data_analysis/beta/matrix/distance-matrix.tsv data_analysis/beta/matrix/weighted_unifrac.tsv ;

qiime tools export \
--input-path  data_analysis/metrics/diversity/unweighted_unifrac_distance_matrix.qza \
--output-path data_analysis/beta/matrix ;
mv data_analysis/beta/matrix/distance-matrix.tsv data_analysis/beta/matrix/unweighted_unifrac.tsv

mkdir -p data_analysis/beta/emperor/jaccard ;
qiime tools export \
--input-path data_analysis/metrics/diversity/jaccard_emperor.qzv \
--output-path data_analysis/beta/emperor/jaccard ;
mkdir -p data_analysis/beta/emperor/bray_curtis ;
qiime tools export \
--input-path  data_analysis/metrics/diversity/bray_curtis_emperor.qzv \
--output-path data_analysis/beta/emperor/bray_curtis ;
mkdir -p data_analysis/beta/emperor/weighted_unifrac ;
qiime tools export \
--input-path  data_analysis/metrics/diversity/weighted_unifrac_emperor.qzv \
--output-path data_analysis/beta/emperor/weighted_unifrac ;
mkdir -p data_analysis/beta/emperor/unweighted_unifrac ;
qiime tools export --input-path data_analysis/metrics/diversity/unweighted_unifrac_emperor.qzv \
--output-path data_analysis/beta/emperor/unweighted_unifrac

PCoA多样性分析

mkdir -p data_analysis/beta/pcoa/jaccard ;
qiime tools export \
--input-path  data_analysis/metrics/diversity/jaccard_pcoa_results.qza \
--output-path data_analysis/beta/pcoa/jaccard ;

mkdir -p data_analysis/beta/pcoa/bray_curtis ;
qiime tools export \
--input-path  data_analysis/metrics/diversity/bray_curtis_pcoa_results.qza \
--output-path data_analysis/beta/pcoa/bray_curtis ;

mkdir -p data_analysis/beta/pcoa/weighted_unifrac ;
qiime tools export \
--input-path  data_analysis/metrics/diversity/weighted_unifrac_pcoa_results.qza \
--output-path data_analysis/beta/pcoa/weighted_unifrac ;

mkdir -p data_analysis/beta/pcoa/unweighted_unifrac ;
qiime tools export \
--input-path data_analysis/metrics/diversity/unweighted_unifrac_pcoa_results.qza \
--output-path data_analysis/beta/pcoa/unweighted_unifrac

rank_abundance

mkdir -p data_analysis/rank_abundance/report/
rank_abundance.R data_analysis/metrics/report/feature_table.tsv \
data_analysis/rank_abundance/report/rank_abundance.pdf ;
pdf2png data_analysis/rank_abundance/report/rank_abundance.pdf ;

specaccum_curve

mkdir -p data_analysis/specaccum_curve/report/
specaccum_curve.R data_analysis/metrics/report/feature_table.tsv \
data_analysis/specaccum_curve/report/specaccum_curve.pdf
pdf2png data_analysis/specaccum_curve/report/specaccum_curve.pdf

picrust2 功能预测

这一步也需要占据大量的计算资源所以请大家不要自行分析，请看我进行演示

deactive
conda activate /biostack/bioconda/envs/picrust2
rm -rf data_analysis/picrust2/pipeline;
picrust2_pipeline.py \
-s  data_analysis/dada2/denoise/dna-sequences.fasta \
-i  data_analysis/metrics/report/feature_table.tsv \
-o  data_analysis/picrust2/pipeline \
-p  64 --no_pathways \
--stratified \
--wide_table \
--per_sequence_contrib 

mkdir -p data_analysis/picrust2/prediction
mkdir -p  data_analysis/picrust2/report
tsv-utils strip mapping_file.txt  | cut -f1,2 | grep -v '#' >data_analysis/picrust2/report/metadata.txt;
tsv-utils definition \
-d " " /biostack/tools/protocols/qiim-tk-0.0.4/bins/../db/picrust2/ko_info.tsv data_analysis/picrust2/pipeline/KO_metagenome_out/pred_metagenome_unstrat.tsv.gz > data_analysis/picrust2/prediction/ko.txt ;
tsv-utils definition \
-d " " /biostack/tools/protocols/qiim-tk-0.0.4/bins/../db/picrust2/ec_level4_info.tsv data_analysis/picrust2/pipeline/EC_metagenome_out/pred_metagenome_unstrat.tsv.gz \
>data_analysis/picrust2/prediction/enzyme.txt

功能预测基础可视化

qiime2
mkdir -p data_analysis/picrust2/prediction/ko ;
mkdir -p data_analysis/picrust2/prediction/enzyme
PCA.R -g T data_analysis/picrust2/prediction/ko.txt data_analysis/picrust2/report/metadata.txt data_analysis/picrust2/prediction/ko/ko.pca.pdf;
pdf2png data_analysis/picrust2/prediction/ko/ko.pca.pdf ;
PCoA.R -g T -m bray data_analysis/picrust2/prediction/ko.txt data_analysis/picrust2/report/metadata.txt /project/qiime-tk/data_analysis/picrust2/prediction/ko/ko.pcoa.pdf;
pdf2png data_analysis/picrust2/prediction/ko/ko.pcoa.pdf ;
NMDS.R -g T -m bray data_analysis/picrust2/prediction/ko.txt data_analysis/picrust2/report/metadata.txt /project/qiime-tk/data_analysis/picrust2/prediction/ko/ko.nmds.pdf;
pdf2png data_analysis/picrust2/prediction/ko/ko.nmds.pdf ;
mkdir -p data_analysis/picrust2/prediction/enzyme ;
PCA.R -g T data_analysis/picrust2/prediction/enzyme.txt data_analysis/picrust2/report/metadata.txt data_analysis/picrust2/prediction/enzyme/enzyme.pca.pdf;
pdf2png data_analysis/picrust2/prediction/enzyme/enzyme.pca.pdf ;
PCoA.R -g T -m bray data_analysis/picrust2/prediction/enzyme.txt data_analysis/picrust2/report/metadata.txt data_analysis/picrust2/prediction/enzyme/enzyme.pcoa.pdf;
pdf2png data_analysis/picrust2/prediction/enzyme/enzyme.pcoa.pdf ;
NMDS.R -g T -m bray data_analysis/picrust2/prediction/enzyme.txt data_analysis/picrust2/report/metadata.txt data_analysis/picrust2/prediction/enzyme/enzyme.nmds.pdf;
pdf2png data_analysis/picrust2/prediction/enzyme/enzyme.nmds.pdf

######################################################