Qiime2学习笔记之Qiime2网站示例学习笔记

官网提供了一份完整的分析教程：https://docs.qiime2.org/2018.2/tutorials/moving-pictures

这里也有牛人整理成的中文教程：https://forum.qiime2.org/t/qiime2-chinese-manual/838

就这份中文教程我们首先来了解他的数据来源和数据特征：

数据来源：

示例的的数据来自文章《Moving pictures of the human microbiome》，GenomeBiology 2011，取样来自两个人，身体四个部位，五个时间点；与Qiime1中Illumina Overview Tutorial教程数据来源相同

这篇文章下载地址：http://xueshu.baidu.com/s?wd=paperuri%3A%289788d401d2273dcf8be475b636305a51%29&filter=sc_long_sign&tn=SE_xueshusource_2kduw22v&sc_vurl=http%3A%2F%2Fwww.ncbi.nlm.nih.gov%2Fpmc%2Farticles%2FPMC3271711%2Fpdf%2Fgb-2011-12-5-r50.pdf&ie=utf-8&sc_us=13911429232176682823

数据特征：

Illumina HiSeq using the Earth MicrobiomeProject hypervariable region 4 (V4) 16S rRNA sequencing protocol：

使用Illumina HiSeq；

地球微生物组协议（了解：http://earthmicrobiome.org/）；

单端测序16sv4区域：

弄清楚数据后，我们在开始按照网站的流程运行一次：

#按照官网命令，创建相同的工作目录：

mkdirqiime2-moving-pictures-tutorial

cdqiime2-moving-pictures-tutorial

接下来下载原始数据，首先是实验分组文件：

官网给出了下载地址，但是很显然我并打不开，这里推荐使用（QIIME2中文帮助文档 (Chinese Manual)）提供的链接：wget http://bailab.genetics.ac.cn/markdown/sample-metadata.tsv

可以成功下载，这里我们可以看到和qiime1的mappig文件是一样的

下面下载序列文件和barcode文件：

#首先新建一个文件夹储存之：

mkdiremp-single-end-sequences

官网给出的链接可以使用：

#barcode文件

wget -O"emp-single-end-sequences/barcodes.fastq.gz" https://data.qiime2.org/2017.7/tutorials/moving-pictures/emp-single-end-sequences/barcodes.fastq.gz

#单端v4区域序列文件

wget -O"emp-single-end-sequences/sequences.fastq.gz" https://data.qiime2.org/2017.7/tutorials/moving-pictures/emp-single-end-sequences/sequences.fastq.gz

得到的文件存于emp-single-end-sequences文件夹中名为：

下面开始构建Qiime2要求的文件输入形式（QIIME 2 artifacts），为此，官网上解释其包含数据类型和数据来源：下面是官网给出命令：

“\”：linux下这个符号可以代表换行符，Qiime1的示例中并没有出现这个符号，Qiime2更加注重读者的阅读体验了

qiime toolsimport \

  --type EMPSingleEndSequences \

  --input-path emp-single-end-sequences \

  --output-path emp-single-end-sequences.qza

--type：定义数据类型

--input-path：输入文件路径

--output-path：输出文件路径

qza格式的文件是Qiime2独有的格式，内包含原始数据，分析过程参数和相应的图表绘制过程参数等

如图：

这里示例的数据导入形式为：地球微生物组标准混样单端数据，还有多种输入原始数据的形式；而且几乎在分析的任何步骤都可以导入文件进入Qiime2进行分析，这是导入文件官方说明页面，后面我们专门解读这篇：

https://docs.qiime2.org/2017.12/tutorials/importing/

接下来的这一步qiime1中是不需要的，但是qiime2专门拎出来做这个步骤，将所有序列按照样品拆分为单个序列文件

#下一步对数据进行拆分

qiime demuxemp-single \

  --i-seqs emp-single-end-sequences.qza \

  --m-barcodes-file sample-metadata.tsv \

  --m-barcodes-category BarcodeSequence \

  --o-per-sample-sequences demux.qza

--i-seqs：输入文件，也就是我们上一步构建的qza文件

--m-barcodes-file：分组信息文件，也就是我们下载的第一个文件

--m-barcodes-category：指定按照分组文件中的哪一列进行拆分

--o-per-sample-sequences：输出文件

这里需要说的是Qiime2将拆分命令和质控命令分离，并没有对序列进行任何裁剪去除

拆分好的文件：（序列长度153）

对拆分结果做一个总结，简单做了每个样品序列条数的统计和全部样品的简单统计分析和整体的序列质量的一个总结；

#做一个可视化文件

qiime demuxsummarize \

  --i-data demux.qza \

  --o-visualization demux.qzv

--o-visualization：输出可视化文件，以.qzv为后缀名的文件是Qiime2中的终端文件，也是可视化文件，是Qiime2独有的文件格式同qza一样

#Qiime2提供了查看这个文件的命令

qiime tools viewdemux.qzv

总共展示了两个工作表单，第一张表格很简单，就是简单的整体的和每个样品的序列统计

第二张表格是质量的统计，这里展示一下：

这张图随机在几十万条原始序列中随机不放回的抽取一万条序列做每个碱基的质量频率分布箱线图；

下一步就是序列质控和生成feature table (OTU table),这里提到一个概念

FeatureTable[Frequency]: 频率，即Feature表(OTU表)，为每个样品中对应OTU出现频率的表格（参考）

FeatureData[Sequence]：就是代表序列文件

这里我们提到DADA2，可用于illumina数据，主要就是100%相似性聚类，对低质量的数据去燥整理。

关于这一步我们在qiime1中要去除引物序列，质控，去除嵌合体，之后才能聚类otu，然而在qiime2中并不是这样；

下面是命令

qiime dada2denoise-single \

  --i-demultiplexed-seqs demux.qza \

  --p-trim-left 0 \

  --p-trunc-len 120 \

  --o-representative-sequencesrep-seqs-dada2.qza \

  --o-table table-dada2.qza

--i-demultiplexed-seqs：输入上一步文件

--p-trim-left:切除左端边缘低质量区

--p-trunc-len：序列裁切为120bp

--o-representative-sequences：生成代表序列文件

--o-table：这就是我们的FeatureTable文件

将两个后续分析的文件改名，用于后续分析

mvrep-seqs-dada2.qza rep-seqs.qza

mvtable-dada2.qza table.qza

下面我们进行简单的结果统计

#对feature表进行统计，每个样品有多少天序列，每条序列出现多少次，就相当于Qiime1的otu.table的一个简单统计

qiime feature-tablesummarize \

  --i-table table.qza \

  --o-visualization table.qzv \

  --m-sample-metadata-file sample-metadata.tsv

#就是相当于咱们一般认识的代表序列文件，只不过序列名没那么人性化方便识别

qiimefeature-table tabulate-seqs \

  --i-data rep-seqs.qza \

  --o-visualization rep-seqs.qzv

#可视化总结文件

qiime tools viewtable.qzv

qiime tools viewrep-seqs.qzv

跟换另一种算法作分析：

deblur：前些日子地球微生物组就采用这种方法做的分析，据说还被他们推荐使用，比DADA2好；文章地址在这里：http://msystems.asm.org/content/2/2/e00191-16

#去过滤质量的序列

qiimequality-filter q-score \

 --i-demux demux.qza \

 --o-filtered-sequences demux-filtered.qza \

 --o-filter-stats demux-filter-stats.qza

--o-filtered-sequences：输出质控后的文件

--o-filter-stats：输出一个总结性质文件，详细记录每个样品多少条序列，质控完成还剩多少条之类的，之后我们会可视化这个文件

#生成otu表格和代表序列文件

qiime deblurdenoise-16S \

  --i-demultiplexed-seqs demux-filtered.qza \

  --p-trim-length 120 \

  --o-representative-sequencesrep-seqs-deblur.qza \

  --o-table table-deblur.qza \

  --p-sample-stats \

  --o-stats deblur-stats.qza

--o-stats：输出一个类似于总结性质的文件，之后我们会可视化这个文件

#制作输出文件，总结性质的文件

qiime metadatatabulate \

  --m-input-file demux-filter-stats.qza \

  --o-visualization demux-filter-stats.qzv

#

qiime deblurvisualize-stats \

  --i-deblur-stats deblur-stats.qza \

  --o-visualization deblur-stats.qzv

#可视化，可以提取出来的

qiimetools view demux-filter-stats.qzv

qiime tools viewdeblur-stats.qzv

#修改代表序列文件

mvrep-seqs-deblur.qza rep-seqs.qza

#修改otu表格文件

mvtable-deblur.qza table.qza

两种算法，结果用于后续多样性分析；

正如官网上所言，Qiime2的确支持多种系统发育多样性分析，分析矩阵主要是Faith’s Phylogenetic Diversity and weighted and unweighted UniFrac，这些矩阵的构建需要有根进化树（rooted phylogenetic tree）下面我们来构建有根发育树

#首先我们对代表序列文件进行多序列比对

qiime alignmentmafft \

  --i-sequences rep-seqs.qza \

  --o-alignment aligned-rep-seqs.qza

接下来我们移除高变区（positions that are highly variable），Qiime2对高变区的去除是认为这些高边区域通常会给系统发生树增加噪音。

qiime alignmentmask \

  --i-alignment aligned-rep-seqs.qza \

  --o-masked-alignmentmasked-aligned-rep-seqs.qza

我们开始构建无根树

qiime phylogenyfasttree \

  --i-alignment masked-aligned-rep-seqs.qza \

  --o-tree unrooted-tree.qza

无根树转化为有根树

qiime phylogenymidpoint-root \

  --i-tree unrooted-tree.qza \

  --o-rooted-tree rooted-tree.qza

这里我们来对比qiime1，这里要说的就是qiime1中做的是无根树，也没有转化为有根树

同样是构建发育树在qiime1中我们可以这样：

#对齐，这是qiime1跑进化树的方法

align_seqs.py -irep_seqs.fa -o pynast_aligned_seqs

#过滤

filter_alignment.py-o pynast_aligned_seqs/ -i pynast_aligned_seqs/rep_seqs_aligned.fasta

#做树

make_phylogeny.py-i pynast_aligned_seqs/rep_seqs_aligned_pfiltered.fasta -o rep_set.tree

#或者改进一下安装多序列比对软件clustalo，一样的结果

sudo apt-getinstall clustalo

# clusto多序列比对

clustalo -i tree_tomato/rep_seqs_k5.fa-o tree_tomato/rep_seqs_k5_clus.fa --seqtype=DNA --full --force --threads=5

# 使用qiime的脚本调用fastree建树

make_phylogeny.py-i tree_tomato/rep_seqs_k5_clus.fa -o tree_tomato/rep_seqs_k5.tree

下面来做多样性分析，这里官网给出的命令同时计算alpha和beta多样性：

qiime diversitycore-metrics-phylogenetic \

  --i-phylogeny rooted-tree.qza \

  --i-table table.qza \

  --p-sampling-depth 1109 \

  --m-metadata-file sample-metadata.tsv \

  --output-dir core-metrics-results

--i-phylogeny：输入有根进化树

--p-sampling-depth：抽样深度，这里没有使用最小的样品深度，选择了1109，就等于去掉了两个样品；使用的重抽样计算alpha多样性和beta多样性，不是唯一的选择，这里我们今后会谈到其他的的保准化办法来分析多样性，结果也会略有不同；

--m-metadata-file：分组文件，现在是咱们第一次使用，都快把它忘记了；

这里会计算很多多样性指标，包括数据文件和图形展示

图形主要展示pcoa结果图形：

qiime tools viewcore-metrics-results/unweighted_unifrac_emperor.qzv

我们通过调节坐标轴，可以展示pcoa结果，可以展示任何轴，可以三维展示（默认）

下一步我们来做显著性分析：

官网例子选择的是faith_pd指标和evenness指标，使用非参数检验Kruskal-Wallis首先对多组数据进行检验，之后又使用非参数检验两组之间的显著性

#alpha多样性

qiime diversityalpha-group-significance \

  --i-alpha-diversitycore-metrics-results/faith_pd_vector.qza \

  --m-metadata-file sample-metadata.tsv \

  --o-visualizationcore-metrics-results/faith-pd-group-significance.qzv

qiime diversityalpha-group-significance \

  --i-alpha-diversitycore-metrics-results/evenness_vector.qza \

  --m-metadata-file sample-metadata.tsv \

  --o-visualizationcore-metrics-results/evenness-group-significance.qzv

--i-alpha-diversity：上一步计算的到的alpha指数文件

--m-metadata-file：分组文件

--o-visualization：输出文件

#这里我们做一个展示

qiime tools viewcore-metrics-results/ faith-pd-group-significance.qzv

qiime tools viewcore-metrics-results/evenness-group-significance.qzv

接下来对beta多样性指标进行检置换验

qiime diversitybeta-group-significance \

  --i-distance-matrixcore-metrics-results/unweighted_unifrac_distance_matrix.qza \

  --m-metadata-file sample-metadata.tsv \

  --m-metadata-column BodySite \

  --o-visualizationcore-metrics-results/unweighted-unifrac-body-site-significance.qzv \

  --p-pairwise

qiime diversitybeta-group-significance \

  --i-distance-matrix core-metrics-results/unweighted_unifrac_distance_matrix.qza\

  --m-metadata-file sample-metadata.tsv \

  --m-metadata-column Subject \

  --o-visualizationcore-metrics-results/unweighted-unifrac-subject-group-significance.qzv \

  --p-pairwise

--i-distance-matrix：输入距离矩阵文件

--m-metadata-column：指定分组的列

--p-pairwise：两组之间检验

#这里我们做一个展示

qiime tools viewcore-metrics-results/unweighted-unifrac-body-site-significance.qzv

qiime tools viewcore-metrics-results/unweighted-unifrac-subject-group-significance.qzv

做一个可视化

qiime emperorplot \

  --i-pcoacore-metrics-results/unweighted_unifrac_pcoa_results.qza \

  --m-metadata-file sample-metadata.tsv \

  --p-custom-axis DaysSinceExperimentStart \

  --o-visualizationcore-metrics-results/unweighted-unifrac-emperor-DaysSinceExperimentStart.qzv

qiime emperorplot \

  --i-pcoacore-metrics-results/bray_curtis_pcoa_results.qza \

  --m-metadata-file sample-metadata.tsv \

  --p-custom-axis DaysSinceExperimentStart \

  --o-visualizationcore-metrics-results/bray-curtis-emperor-DaysSinceExperimentStart.qzv

#自行查看一下

qiime tools viewcore-metrics-results/bray-curtis-emperor-DaysSinceExperimentStart.qzv

注意一下官网上示例在实际运行过程中会有不一样的地方

#按照—help的实际参数调整即可

qiime emperorplot --help

下一步绘制稀释曲线：

qiime diversityalpha-rarefaction \

  --i-table table.qza \

  --i-phylogeny rooted-tree.qza \

  --p-max-depth 4000 \

  --m-metadata-file sample-metadata.tsv \

  --o-visualization alpha-rarefaction.qzv

--i-table：使用前面otu.table数据

--p-max-depth：设置抽样深度为4000

--i-phylogeny：使用进化树数据

#自行查看一下

qiime tools viewalpha-rarefaction.qzv

下面我们就开始进行物种注释了：

点击官网链接下载注释数据库，还挺大的八百多兆

wget -O"gg-13-8-99-515-806-nb-classifier.qza""https://data.qiime2.org/2018.2/common/gg-13-8-99-515-806-nb-classifier.qza"

qiimefeature-classifier classify-sklearn \

  --i-classifiergg-13-8-99-515-806-nb-classifier.qza \

  --i-reads rep-seqs.qza \

  --o-classification taxonomy.qza

#凡是以qzv结尾的都是输出文件，供人们查看

qiime metadatatabulate \

  --m-input-file taxonomy.qza \

  --o-visualization taxonomy.qzv

--i-classifier：参考数据库文件

--i-reads：代表序列文件

#自行查看一下

qiime tools viewtaxonomy.qzv

#能鉴定到中的数量不多

做一个柱状图展示一下门类情况

qiime taxabarplot \

  --i-table table.qza \

  --i-taxonomy taxonomy.qza \

  --m-metadata-file sample-metadata.tsv \

  --o-visualization taxa-bar-plots.qzv

#自行查看一下

qiime tools viewtaxa-bar-plots.qzv

#与qiime1不同，最明显的便就是出现了level7，种水平的数据

qiime2给我们的惊喜还是挺多的，现在我们可以做差异分析，不多qiime2使用的方法是ANCOM，我之前没有接触过，但其官方声称更加适合微生物数据

#这是qiime2中过滤数据的命令，这里只选择了gut相关的样本

qiimefeature-table filter-samples \

  --i-table table.qza \

  --m-metadata-file sample-metadata.tsv \

  --p-where "BodySite='gut'" \

  --o-filtered-table gut-table.qza

ANCOM不允许丰度值为0，就需要将0改为最小的数字

qiimecomposition add-pseudocount \

  --i-table gut-table.qza \

  --o-composition-table comp-gut-table.qza

开始做差异分析，这里注意，官网教程中的参数选项出现错误，大家按照我们命令修改就好

qiimecomposition ancom \

  --i-table comp-gut-table.qza \

  --m-metadata-file sample-metadata.tsv \

  --m-metadata-category Subject \

  --o-visualization ancom-Subject.qzv

#查看结果

qiime tools viewancom-Subject.qzv

在这里可以对任何一个分类等级进行显著性分析

#这里是对

qiime taxacollapse \

  --i-table gut-table.qza \

  --i-taxonomy taxonomy.qza \

  --p-level 6 \

  --o-collapsed-table gut-table-l6.qza

qiimecomposition add-pseudocount \

  --i-table gut-table-l6.qza \

  --o-composition-table comp-gut-table-l6.qza

qiimecomposition ancom \

  --i-table comp-gut-table-l6.qza \

  --m-metadata-file sample-metadata.tsv \

  --m-metadata-category Subject \

  --o-visualization l6-ancom-Subject.qzv

#查看结果

qiime tools viewl6-ancom-Subject.qzv

这里的结果部分比较新，我以前也没有见过，一共分为三个部分，第一个就是分析结果统计，这里提到了W这个统计量，应该是大于某个阈值的被认定为显著差异，但是有些值也挺大，却不是被判定为差异；在第一张表格中，这里官网示例采用分组变量下，只有一个显著的，其他不显著的我们可以通过查看蓝色的csv文件；

第二张表格展示的丰度信息，在四个百分比下对每个组的分类信息做一个统计表格，可能是因为考虑到差异物种可能很多，所以使用表格更加方便一下，箱线图这种很好的表示方法就不用了；

第三章图是一个火山图，展示所用物种信息，相当一第一张表格csv文件的一个可视化；

到此，qiime2的示例文件就算是学习完成了，里面有很多新的东西，有很多我们不一样的分析习惯，还有许多我们不了解的概念和观点，至于今后我们应该如何分析，选择怎样的一个流程，我们还得根据自己的实际情况和个人习惯选择合适的流程

所以接下来我主要就是根据自己的分析需求和实际情况，结合qiime2，做一套适合自己的分析流程

学习永无止境，分享永不停歇！