扩增子分析QIIME2. 2分析实战Moving Pictures

欢迎点击「宏基因组」关注我们！专业干货每日推送！

声明：本文为QIIME2官方帮助文档的中文版，由中科院遗传发育所刘永鑫博士翻译并亲测有效，文档翻译己获QIIME2团队官方授权。由于QIIME2更新频繁，如使用中遇到问题请访问QIIME2官方论坛阅读最新版中文帮助。
https://forum.qiime2.org/t/qiime2-1-chinese-manual/838
如中文翻译没有急时更新，新阅读英文原版 https://docs.qiime2.org

本人只习惯使用命令行模式分析数据，图形界面和Ipython模式下使用暂不介绍。本系列的教程主要以命令行方式为大家演示。在其它方面有使用的经验的朋友，欢迎共享您的使用笔记于“宏基因组”公众号，方便大家学习和使用。

扩增子分析QIIME2. 2分析实战：人体不同部分微生物组Moving Pictures

QIIME2版本 2017.6

本文的操作的前提是完成QIIME2的安装，想安装QIIME2请阅读《扩增子分析QIIME2. 1简介和安装》阅读；

本示例的的数据来自文章《Moving pictures of the human microbiome》，Genome Biology 2011，取样来自两个人，身体四个部位，五个时间点。

启动QIIME2运行环境

对于上文提到了两种常用安装方法，我们每次在分析数据前，需要打开工作环境，根据情况选择对应的打开方式。

# 创建工作目录并进入 mkdir -p qiime2 cd qiime2 # Miniconda安装方式的请运行如下命令加载工作环境 source activate qiime2-2017.6 # 如果是docker安装方式的请运行如下命令，否则跳过此行 docker run --rm -v $(pwd):/data --name=qiime -it  qiime2/core:2017.6

准备数据

# 创建并进入工作目录 mkdir -p qiime2-moving-pictures-tutorial cd qiime2-moving-pictures-tutorial # 下载实验设计，此文件位于google服务器，服务器无法翻墙的可以用windows翻墙下载并上传到工作目录；或尝试下载我建立的备份 wget -O "sample-metadata.tsv" "https://data.qiime2.org/2017.6/tutorials/moving-pictures/sample_metadata.tsv" ## 上面一步下载失败，可尝试删除空文件并使用我建立的备份链接下载；否则跳过下面两行命令 rm sample_metadata.tsv wget http://bailab.genetics.ac.cn/markdown/sample-metadata.tsv # 下载实验测序数据 mkdir -p emp-single-end-sequences wget -O "emp-single-end-sequences/barcodes.fastq.gz" "https://data.qiime2.org/2017.6/tutorials/moving-pictures/emp-single-end-sequences/barcodes.fastq.gz" wget -O "emp-single-end-sequences/sequences.fastq.gz" "https://data.qiime2.org/2017.6/tutorials/moving-pictures/emp-single-end-sequences/sequences.fastq.gz" # 生成qiime需要的artifact文件(qiime文件格式，将原始数据格式标准化) qiime tools import \  --type EMPSingleEndSequences \  --input-path emp-single-end-sequences \  --output-path emp-single-end-sequences.qza

拆分样品

# 按barcode拆分样品 Demultiplexing sequences qiime demux emp-single \  --i-seqs emp-single-end-sequences.qza \  --m-barcodes-file sample-metadata.tsv \  --m-barcodes-category BarcodeSequence \  --o-per-sample-sequences demux.qza # 结果统计 qiime demux summarize \  --i-data demux.qza \  --o-visualization demux.qzv # 查看结果 (依赖XShell+XManager或其它ssh终端和图形界面软件) qiime tools view demux.qzv

使用qiime tools view demux.qzv或访问 https://view.qiime2.org/?src=https%3A%2F%2Fdocs.qiime2.org%2F2017.6%2Fdata%2Ftutorials%2Fmoving-pictures%2Fdemux.qzv 在线查看
下图展示拆分样品的图表结果，其中的图表示各样品测序数据柱状分布图，展示不同测序深度下样品数量分布，可下载PDF文件；表可下载CSV格式文件，即每个样品的数据量表，可用Excel打开。

序列质控和生成OTU表

此步主要有DADA2和Deblur两种方法可选，推荐使用DADA2，去年发表在Nature Method上，比较同类方法结果(基于分析mock community模拟重组数据)优于其它OTU聚类结果，结果更真实；相较QIIME的UPARSE聚类方法，目前DADA2方法仅去噪去嵌合，不再按相似度聚类，相当于100%相似度的OTU。比上一代分析结果更准确。

1. DADA2
   主要作用是去除低质量序列、嵌合体；再生成OTU表，现在叫Feature表，因为不再使用聚类方法，相当于QIIME时代100%相似度的OTU表。

读者思考时间：基于上图对拆分样品的统计结果，如何设置下面生成OTU表的参数。

1. —p-trim-left 截取左端低质量序列，我们看上图中箱线图，左端质量都很高，无低质量区，设置为0；

2. —p-trunc-len 序列截取长度，也是为了去除右端低质量序列，我们看到大于120以后，质量下降极大，甚至中位数都下降至20以下，需要全部去除。

# 单端序列去噪, 去除左端0bp(--p-trim-left用于切除边缘低质量区)，序列切成120bp长；生成代表序列和OTU表；并重命名用于下游分析 qiime dada2 denoise-single \  --i-demultiplexed-seqs demux.qza \  --p-trim-left 0 \  --p-trunc-len 120 \  --o-representative-sequences rep-seqs-dada2.qza \  --o-table table-dada2.qza mv rep-seqs-dada2.qza rep-seqs.qza mv table-dada2.qza table.qza

2. Deblur
与DADA2二选一，用户可自行比较结果的差异，根据喜好选择。

# 按测序质量过滤序列 qiime quality-filter q-score \ --i-demux demux.qza \ --o-filtered-sequences demux-filtered.qza \ --o-filter-stats demux-filter-stats.qza # 去冗余生成OTU表和代表序列；结果文件名有deblur，没有用于下游分析，请读者想测试的自己尝试 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 \  --o-stats deblur-stats.qza

Feature表统计

qiime feature-table summarize \  --i-table table.qza \  --o-visualization table.qzv \  --m-sample-metadata-file sample-metadata.tsv qiime tools view table.qzv # 展打开网页版统计结果，看下图

图中展示了Feature表的统计结果

代表序列统计

qiime feature-table tabulate-seqs \  --i-data rep-seqs.qza \  --o-visualization rep-seqs.qzv qiime tools view rep-seqs.qzv # 展打开网页版统计结果，看下图

下图展示代表序列网页版详细，点"Click here"可下载fasta文件

建树：用于多样性分析

# 多序列比对 qiime alignment mafft \  --i-sequences rep-seqs.qza \  --o-alignment aligned-rep-seqs.qza # 移除高变区 qiime alignment mask \  --i-alignment aligned-rep-seqs.qza \  --o-masked-alignment masked-aligned-rep-seqs.qza # 建树 qiime phylogeny fasttree \  --i-alignment masked-aligned-rep-seqs.qza \  --o-tree unrooted-tree.qza # 无根树转换为有根树 qiime phylogeny midpoint-root \  --i-tree unrooted-tree.qza \  --o-rooted-tree rooted-tree.qza

Alpha多样性

读者思考时间：下面多样性分析，需要基于标准化的OTU表，标准化采用重抽样至序列一致，如何设计样品重抽样深度参数。—p-sampling-depth

如是数据量都很大，选最小的即可。如果有个别数据量非常小，去除最小值再选最小值。比如此分析最小值为917，我们选择1080深度重抽样，即保留了大部分样品用于分析，又去除了数据量过低的异常值。
注：本示例为454时代的测序，数据量很小。现在一般采用HiSeq PE250测序，数据量都非常大，通常可以采用3万或5万的标准筛选，仍可保留90%以上样本。过低或过高一般结果也会异常，不建议放在一起分析。

# 计算多样性(包括所有常用的Alpha和Beta多样性方法)，输入有根树、Feature表、样本重采样深度(一般为最小样本数据量，或覆盖绝大多数样品的数据量) qiime diversity core-metrics \  --i-phylogeny rooted-tree.qza \  --i-table table.qza \  --p-sampling-depth 1080 \  --output-dir core-metrics-results # 输出结果包括多种多样性结果，文件列表和解释如下： # beta多样性bray_curtis距离矩阵 bray_curtis_distance_matrix.qza # alpha多样性evenness(均匀度，考虑物种和丰度)指数 evenness_vector.qza # alpha多样性faith_pd(考虑物种间进化关系)指数 faith_pd_vector.qza # beta多样性jaccard距离矩阵 jaccard_distance_matrix.qza # alpha多样性observed_otus(OTU数量)指数 observed_otus_vector.qza # alpha多样性香农熵(考虑物种和丰度)指数 shannon_vector.qza # beta多样性unweighted_unifrac距离矩阵，不考虑丰度 unweighted_unifrac_distance_matrix.qza # beta多样性unweighted_unifrac距离矩阵，考虑丰度 weighted_unifrac_distance_matrix.qza # 统计faith_pd算法Alpha多样性组间差异是否显著，输入多样性值、实验设计，输出统计结果 qiime diversity alpha-group-significance \  --i-alpha-diversity core-metrics-results/faith_pd_vector.qza \  --m-metadata-file sample-metadata.tsv \  --o-visualization core-metrics-results/faith-pd-group-significance.qzv # 统计evenness组间差异是否显著 qiime diversity alpha-group-significance \  --i-alpha-diversity core-metrics-results/evenness_vector.qza \  --m-metadata-file sample-metadata.tsv \  --o-visualization core-metrics-results/evenness-group-significance.qzv # 网页展示结果，只要是qzv的文件，均可用qiime tools view查看或在线https://view.qiime2.org/查看，以后不再赘述 qiime tools view evenness-group-significance.qzv

以Evenness为例，下图展示多样性的箱线图，可以下载svg矢量图

也可点击如下网址，查看可交互结果。https://view.qiime2.org/visualization/?type=html&src=https%3A%2F%2Fdocs.qiime2.org%2F2017.6%2Fdata%2Ftutorials%2Fmoving-pictures%2Fcore-metrics-results%2Fevenness-group-significance.qzv

读者思考时间：实验设计中的那一种分组方法，与微生物群体的丰富度差异相关，这些差异显著吗？

解答：图中可按Catalogy选择分类方法，查看不同分组下箱线图间的分布与差别。图形下面的表格，详细详述组间比较的显著性和假阳性率统计。
结果我们会看到本实验设计的分组方式有Bodysite, Subject, ReportAntibioticUse，只有身体位置各组间差异明显，且下面统计结果也存在很多组间的显著性差异。

Beta多样性

# 按BodySite分组，统计unweighted_unifrace距离的组间是否有显著差异 qiime diversity beta-group-significance \  --i-distance-matrix core-metrics-results/unweighted_unifrac_distance_matrix.qza \  --m-metadata-file sample-metadata.tsv \  --m-metadata-category BodySite \  --o-visualization core-metrics-results/unweighted-unifrac-body-site-significance.qzv \  --p-pairwise # 按Subject分组，统计unweighted_unifrace距离的组间是否有显著差异 qiime diversity beta-group-significance \  --i-distance-matrix core-metrics-results/unweighted_unifrac_distance_matrix.qza \  --m-metadata-file sample-metadata.tsv \  --m-metadata-category Subject \  --o-visualization core-metrics-results/unweighted-unifrac-subject-group-significance.qzv \  --p-pairwise # 可视化三维展示unweighted-unifrac的主坐标轴分析 qiime emperor plot \  --i-pcoa core-metrics-results/unweighted_unifrac_pcoa_results.qza \  --m-metadata-file sample-metadata.tsv \  --p-custom-axis DaysSinceExperimentStart \  --o-visualization core-metrics-results/unweighted-unifrac-emperor.qzv # 可视化三维展示bray-curtis的主坐标轴分析 qiime emperor plot \  --i-pcoa core-metrics-results/bray_curtis_pcoa_results.qza \  --m-metadata-file sample-metadata.tsv \  --p-custom-axis DaysSinceExperimentStart \  --o-visualization core-metrics-results/bray-curtis-emperor.qzv # 网页展示结果，或下载在线查看 qiime tools view core-metrics-results/bray-curtis-emperor.qzv

以Bray_curtis可视化结果为例，下图展示多样性的三维散点图

读者思考时间：按subject分组有显著区别吗？按body-site分组有显著区别吗？那些body-site组间存在区别？

大家仔细查看按subject进行主坐标轴分析的结果，如下链接，我是没看到区别；
https://view.qiime2.org/?src=https%3A%2F%2Fdocs.qiime2.org%2F2017.6%2Fdata%2Ftutorials%2Fmoving-pictures%2Fcore-metrics-results%2Funweighted-unifrac-subject-group-significance.qzv
大家仔细查看按body-site进行主坐标轴分析的结果，如下链接，只有左右手间无明显差别，其它各组间均有显著差别。
https://view.qiime2.org/?src=https%3A%2F%2Fdocs.qiime2.org%2F2017.6%2Fdata%2Ftutorials%2Fmoving-pictures%2Fcore-metrics-results%2Funweighted-unifrac-body-site-significance.qzv

按其它距离计算的结果，读者可以仔细看看不同距离矩阵计算结果的区别。个人感觉，一般比较好解释科学问题的方法就是适合的方法。

物种分类

# 下载物种注释 wget -O "gg-13-8-99-515-806-nb-classifier.qza" "https://data.qiime2.org/2017.6/common/gg-13-8-99-515-806-nb-classifier.qza" # 物种分类 qiime feature-classifier classify-sklearn \  --i-classifier gg-13-8-99-515-806-nb-classifier.qza \  --i-reads rep-seqs.qza \  --o-classification taxonomy.qza # 物种结果转换表格，可用于查看 qiime taxa tabulate \  --i-data taxonomy.qza \  --o-visualization taxonomy.qzv # 展示taxonomy.qzv结果如下： #Feature ID    Taxonomy #d12759fe8dda1d65fe9077cc1ca9cf28    k__Bacteria; p__Bacteroidetes; c__Flavobacteriia; o__Flavobacteriales; f__[Weeksellaceae]; g__Chryseobacterium; s__ #5ada68b9a081358e1a7d5f1d351e656a    k__Bacteria; p__Fusobacteria; c__Fusobacteriia; o__Fusobacteriales; f__Leptotrichiaceae; g__Leptotrichia; s__ #d9095748835ade1b8914c5f57b6acbcf    k__Bacteria; p__Proteobacteria; c__Gammaproteobacteria; o__Aeromonadales; f__Aeromonadaceae; g__Oceanisphaera; s__ # 物种分类柱状图 qiime taxa barplot \  --i-table table.qza \  --i-taxonomy taxonomy.qza \  --m-metadata-file sample-metadata.tsv \  --o-visualization taxa-bar-plots.qzv qiime tools view taxa-bar-plots.qzv

以堆叠柱状图，展示物种各组、各分类级别下的相对丰度，图形可交互查看每个柱的详细信息，也可以下载svg矢量图。

读者思考时间1：代表序列文件rep-seqs.qzv可视化结果中，可以下载fasta文件采用NCBI进行blast注释物种信息，与我们目前的结果比较，看看有什么不同，各分类级别的注释定义的相似程度是什么？
读者思考时间2：查看门水平(level2)分类结果柱状图，看每一类body-site中主要丰度的门类是什么？

差异丰度分析

差异丰度分析采用ANCOM (analysis of composition of microbiomes)，是2015年发布在Microb Ecol Health Dis上的方法，文章称在微生物组方面更专业，但不接受零值(零在二代测序结果表中很常见)。我个人一直用edgeR，感觉靠谱，因为高通量测序本质上是相同的。

差异Features/OTUs分析

# OTU表添加假count，因为ANCOM不允许有零 qiime composition add-pseudocount \  --i-table table.qza \  --o-composition-table comp-table.qza # 采用ancon，按BodySite分组进行差异统计 qiime composition ancom \  --i-table comp-table.qza \  --m-metadata-file sample-metadata.tsv \  --m-metadata-category BodySite \  --o-visualization ancom-BodySite.qzv # 查看结果 qiime tools view ancom-BodySite.qzv

读者思考时间：不同身体部分有那些Features存在丰度差异？那一组是最高或最低丰度？这此差异的Features属那些分类单元？

使用qiime tools view ancom-BodySite.qzv查看结果，或下载结果用qiime2网页在线打开查看，回答上面的问题。

差异分类学级别分析：以按门水平合并再统计差异

# 按门水平进行合并，统计各门的总reads qiime taxa collapse \  --i-table table.qza \  --i-taxonomy taxonomy.qza \  --p-level 2 \  --o-collapsed-table table-l2.qza # 同理去除零 qiime composition add-pseudocount \  --i-table table-l2.qza \  --o-composition-table comp-table-l2.qza # 在门水平按取样部分分析 qiime composition ancom \  --i-table comp-table-l2.qza \  --m-metadata-file sample-metadata.tsv \  --m-metadata-category BodySite \  --o-visualization l2-ancom-BodySite.qzv

可视化结果包括统计、各组丰度分位数和交互式火山图

读者思考时间：不同身体部分有那些Features存在丰度差异？那一组是最高或最低丰度？这此差异的Features属那些分类单元？
使用qiime tools view ancom-BodySite.qzv查看结果，或点如下链接查看详细，回答上面的问题。
https://view.qiime2.org/?src=https%3A%2F%2Fdocs.qiime2.org%2F2017.6%2Fdata%2Ftutorials%2Fmoving-pictures%2Fl2-ancom-BodySite.qzv

结果描述：结果的可视化(Visual)页面，一共分为三部分。

1. 第一个表为ANCOM statistical results，只列出组间存在显著差异的门，其统计值W的计算及解释尚不清楚，查原始文章也没有找到。有待更新版中解释。

2. 第二个表为各组的丰度分位数，就是箱线图的原始数据，为什么作者没有直接出图，我将与作者沟通讨论；目前可以比较各组的分布，来具体看组间的差异，但不够直观；

3. 统计各门类的火山图，坐标轴还没有详细解释，但其意思是越靠上越显著差异。此图采用Python的bokeh库生成的交互式图形，可以点击图中的点来查看具体的详细，如具体的分类学信息。相当于表1的可视化。

结果的网页还有其它页面，如peek页面可以查看此文件的基本信息，Provenance页面显示当前结果的生成过程图，点击过程中的点可以查看具体的程序和参数；链接按扭可以生成共享链接；下载按扭可以下载原始文件。

Reference

1. https://docs.qiime2.org/2017.6/tutorials/moving-pictures/

2. Microb Ecol Health Dis. 2015 May 29;26:27663. doi: 10.3402/mehd.v26.27663. eCollection 2015. https://www.ncbi.nlm.nih.gov/pubmed/26028277

3. Nat Methods. 2016 Jul;13(7):581-3. doi: 10.1038/nmeth.3869. Epub 2016 May 23. DADA2: High-resolution sample inference from Illumina amplicon data. https://www.ncbi.nlm.nih.gov/pubmed/27214047

4. Deblur Rapidly Resolves Single-Nucleotide Community Sequence Patterns http://msystems.asm.org/content/2/2/e00191-16

5. Python交互可视化图形库 http://bokeh.pydata.org/en/latest/

科研经验系列文章

扩增子图表解读第一季

扩增子分析

想了解更多宏基因组、16S分析相关文章，

快关注“宏基因组”公众号，干货第一时间推送。

系统学习生物信息，快关注“生信宝典”，

那里有几千志同道合的小伙伴一起学习。