Qiime2人体菌群分析

1. 项目介绍：

在本教程中，您将使用QIIME 2在五个时间点对四个身体部位的两个人进行人类微生物组样本的分析，其中第一个立即使用抗生素。 基于这些样品的研究最初发表于Caporaso等人。 （2011年）。 使用Earth Microbiome Project高变区4（V4）16S rRNA测序方案，在Illumina HiSeq上对本教程中使用的数据进行测序。

在开始本教程之前，创建一个新目录并切换到该目录。

---

mkdir qiime2-moving-pictures-tutorial

cd qiime2-moving-pictures-tutorial

---

2. 示例原始数据

在开始分析之前，请浏览样本元数据以熟悉本研究中使用的样本。 示例元数据以Google表格形式提供。 您可以选择“文件”>“下载为”>“制表符分隔值”，将此文件作为制表符分隔文本下载。 或者，以下命令将样本元数据下载为制表符分隔文本，并将其保存在文件sample-metadata.tsv中。 在本教程的其余部分中使用了此sample-metadata.tsv文件。

---

wget \
-O "sample-metadata.tsv" \
"

---

3. 获取数据

下载我们将在此分析中使用的序列reads。在本教程中，我们将使用一小部分完整的序列数据，以便命令能够快速运行。

---

mkdir emp-single-end-sequences

wget \
-O "emp-single-end-sequences/barcodes.fastq.gz" \
"

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

---

4. 导入数据

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

此QIIME 2工件的语义类型是EMPSingleEndSequences。 EMPSingleEndSequences QIIME 2工件包含多路复用的序列，这意味着序列尚未分配给样本（因此包含sequence.fastq.gz和barcodes.fastq.gz文件，其中barcodes.fastq.gz包含条形码 读取与sequences.fastq.gz中的每个序列相关联。）

---

qiime tools import \
--type EMPSingleEndSequences \
--input-path emp-single-end-sequences \
--output-path emp-single-end-sequences.qza

Output artifacts:

emp-single-end-sequences.qza

---

5. 解复用序列

为了解复用序列，我们需要知道哪个条形码序列与每个样本相关联。 此信息包含在示例元数据文件中。 您可以运行以下命令来解复用序列（demux emp-single命令指的是这些序列根据Earth Microbiome Project协议进行条形码编码，并且是单端读取的事实）。 demux.qza QIIME 2工件将包含解复用的序列。 第二个输出（demux-details.qza）显示Golay错误纠正详细信息，本教程不会探讨（您可以使用qiime元数据列表来显示这些数据）。

---

qiime demux emp-single \
--i-seqs emp-single-end-sequences.qza \
--m-barcodes-file sample-metadata.tsv \
--m-barcodes-column BarcodeSequence \
--o-per-sample-sequences demux.qza \
--o-error-correction-details demux-details.qza

Output artifacts:

demux-details.qza

demux.qza

---

在解复用之后，生成解复用结果的摘要是有用的。 这使您可以确定每个样本获得的序列数，还可以获得序列数据中每个位置的序列质量分布的摘要。

---

qiime demux summarize \
--i-data demux.qza \
--o-visualization demux.qzv

Output visualizations

demux.qzv

---

6. 序列质量控制和特征表构建

QIIME 2插件可用于多种质量控制方法，包括DADA2，Deblur和基于质量得分的基本过滤。 在本教程中，我们使用DADA2和Deblur演示此步骤。 这些步骤是可以互换的，因此您可以使用您喜欢的任何一种。 这两种方法的结果都是FeatureTable [Frequency] QIIME 2工件，它包含数据集中每个样本中每个唯一序列的计数（频率），以及FeatureData [Sequence] QIIME 2工件，它映射了特征标识符。 它们代表的序列的FeatureTable。

6.1 选项1:DADA2

DADA2是用于检测和校正（如果可能）Illumina扩增子序列数据的管道。 如在q2-dada2插件中实施的，该质量控制过程将另外过滤在测序数据中鉴定的任何phiX读数（通常存在于标记基因Illumina序列数据中），并将过滤嵌合序列。

dada2 denoise-single方法需要两个用于质量过滤的参数： - p-trim-left m，它修剪每个序列的前m个碱基， - p-trunc-len n，它截断每个序列 位置 这允许用户移除序列的低质量区域。 要确定要为这两个参数传递的值，您应该查看由qiime demux汇总生成的demux.qzv文件中的交互式质量图选项卡。

在demux.qzv质量图中，我们看到初始碱基的质量似乎很高，因此我们不会从序列的开头修剪任何碱基。 质量似乎在位置120附近下降，因此我们将在120个碱基处截断我们的序列。 下一个命令可能需要10分钟才能运行，这是本教程中最慢的一步。

---

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 --o-denoising-stats stats-dada2.qza

Output artifacts:

stats-dada2.qza

table-dada2.qza

rep-seqs-dada2.qza

---

qiime metadata tabulate \

 --m-input-file stats-dada2.qza --o-visualization stats-dada2.qzv

Output visualizations:

stats-dada2.qzv

---

6.2 选项2：Deblur

Deblur使用序列错误概况将错误的序列读数与衍生它们的真实生物序列相关联，从而产生高质量的序列变体数据。 这分两步应用。 首先，应用基于质量分数的初始质量过滤过程。 该方法是Bokulich等人描述的质量过滤方法的实现。 （2013年）。

---

qiime quality-filter q-score \
--i-demux demux.qza \
--o-filtered-sequences demux-filtered.qza \
--o-filter-stats demux-filter-stats.qza

Output artifacts:

demux-filtered.qza

demux-filter-stats.qza

---

接下来，使用qiime deblur denoise-16S方法应用Deblur工作流程。 此方法需要一个用于质量过滤的参数， - p-trim-length n，它截断位置n处的序列。 一般而言，Deblur开发人员建议将此值设置为中位数质量得分开始降得太低的长度。 在这些数据上，质量图（质量过滤之前）表明合理的选择是在115到130序列位置范围内。 这是一种主观评估。 您可能偏离该建议的一种情况是在多个测序运行中执行元分析。 在这种类型的荟萃分析中，对于所有比较的测序运行，读取长度是相同的是至关重要的，以避免引入研究特异性偏倚。 由于我们已经为qiime dada2 denoise-single使用了120的修剪长度，并且由于120是合理的给定质量图，我们将通过-p-trim-length 120.下一个命令可能需要长达10分钟才能运行 。

---

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

Output artifacts:

deblur-stats.qza

table-deblur.qza

rep-seqs-deblur.qza

---

---

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

Output visualizations:

demux-filter-stats.qzv

deblur-stats.qzv

---

7. 功能表与功能数据摘要

质量过滤步骤完成后，您将需要浏览结果数据。 您可以使用以下两个命令执行此操作，这两个命令将创建数据的可视摘要。 feature-table summarize命令将为您提供有关每个样本和每个要素的关联序列数，这些分布的直方图以及一些相关的摘要统计信息。 feature-table tabulate-seqs命令将提供功能ID到序列的映射，并提供链接以轻松地针对NCBI nt数据库BLAST每个序列。 当您想要了解有关数据集中重要特定功能的更多信息时，后一个可视化将在本教程后期非常有用。

---

qiime feature-table summarize \
--i-table table.qza \
--o-visualization table.qzv \
--m-sample-metadata-file sample-metadata.tsv
qiime feature-table tabulate-seqs \
--i-data rep-seqs.qza \
--o-visualization rep-seqs.qzv

Output visualizations:

table.qzv

rep-seqs.qzv

---

8. 生成系统发生多样性分析树

QIIME支持多种系统发育多样性指标，包括Faith的系统发育多样性和加权和未加权的UniFrac。 除了每个样本的特征计数（即，特征表[频率] QIIME 2伪像中的数据）之外，这些度量还需要根据特征将相互关联的系统发生树相关联。 此信息将存储在Phylogeny [Rooted] QIIME 2工件中。 为了生成系统发育树，我们将使用来自q2-phylogeny插件的align-to-tree-mafft-fasttree管道。

首先，管道使用mafft程序在我们的FeatureData [Sequence]中执行序列的多序列比对，以创建FeatureData [AlignedSequence] QIIME 2工件。 接下来，管道屏蔽（或过滤）对齐以移除高度可变的位置。 通常认为这些位置会增加所产生的系统发育树的噪声。 之后，管道应用FastTree从掩蔽对齐生成系统发育树。 FastTree程序创建一个无根树，因此在本节的最后一步中，应用中点生根将树的根放置在无根树中最长的tip-to-tip距离的中点。

---

qiime phylogeny align-to-tree-mafft-fasttree \

--i-sequences rep-seqs.qza \
--o-alignment aligned-rep-seqs.qza \
--o-masked-alignment masked-aligned-rep-seqs.qza \
--o-tree unrooted-tree.qza \
--o-rooted-tree rooted-tree.qza

Output artifacts:

aligned-rep-seqs.qza

masked-aligned-rep-seqs.qza

rooted-tree.qza

unrooted-tree.qza

---

9. α和β多样性分析

QIIME 2的多样性分析可通过q2-diversity插件获得，该插件支持计算alpha和beta多样性指标，应用相关统计测试以及生成交互式可视化。 我们将首先应用核心 - 度量 - 系统发育方法，该方法将FeatureTable [频率]稀疏到用户指定的深度，计算几个alpha和beta多样性度量，并使用Emperor生成每个的主要坐标分析（PCoA）图。 beta多样性指标。 默认计算的指标是：

Alpha多样性

香农的多样性指数（社区丰富度的量化指标）

观察到的OTU（社区丰富度的定性测量）

信仰的系统发育多样性（社区丰富性的质量测量，包含特征之间的系统发育关系）

均匀度（或Pielou的均匀度;社区均匀度的衡量标准）

Beta多样性

Jaccard距离（社区差异的定性测量）

Bray-Curtis距离（社区差异的量化指标）

未加权的UniFrac距离（社区差异的定性测量，包含特征之间的系统发育关系）

加权UniFrac距离（社区差异的定量测量，包含特征之间的系统发育关系）

---

qiime diversity core-metrics-phylogenetic \
--i-phylogeny rooted-tree.qza \
--i-table table.qza \
--p-sampling-depth 1103 \
--m-metadata-file sample-metadata.tsv \
--output-dir core-metrics-results

Output artifacts:

core-metrics-results/faith_pd_vector.qza

core-metrics-results/unweighted_unifrac_distance_matrix.qza

core-metrics-results/bray_curtis_pcoa_results.qza

core-metrics-results/shannon_vector.qza

core-metrics-results/rarefied_table.qza

core-metrics-results/weighted_unifrac_distance_matrix.qza

core-metrics-results/jaccard_pcoa_results.qza

core-metrics-results/observed_otus_vector.qza

core-metrics-results/weighted_unifrac_pcoa_results.qza

core-metrics-results/jaccard_distance_matrix.qza

core-metrics-results/evenness_vector.qza

core-metrics-results/bray_curtis_distance_matrix.qza

core-metrics-results/unweighted_unifrac_pcoa_results.qza

Output visualizations:

core-metrics-results/unweighted_unifrac_emperor.qzv

core-metrics-results/jaccard_emperor.qzv

core-metrics-results/bray_curtis_emperor.qzv

core-metrics-results/weighted_unifrac_emperor.qzv

---

在计算多样性度量之后，我们可以开始在样本元数据的上下文中探索样本的微生物组成。 此信息存在于先前下载的示例元数据文件中。我们将首先测试分类元数据列和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

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

Output visualizations:

core-metrics-results/faith-pd-group-significance.qzv

core-metrics-results/evenness-group-significance.qzv

---

接下来，我们将使用PERMANOVA（首次在Anderson（2001）中描述）使用beta-group-significance命令在分类元数据的上下文中分析样本组成。 以下命令将测试组内样本之间的距离（例如来自同一身体部位（例如，肠道）的样本）是否彼此更相似，然后它们是来自其他组的样本（例如，舌头，左手掌， 和右掌）。 如果你使用--p-pairwise参数调用这个命令，就像我们在这里做的那样，它也会执行成对测试，这将允许你确定哪些特定的组对（例如，舌头和肠道）彼此不同， 如果有的话。 此命令可能运行缓慢，尤其是在传递--p-pairwise时，因为它基于置换测试。 因此，与以前的命令不同，我们将在我们有兴趣探索的特定元数据列上运行beta-group-significant，而不是适用于它的所有元数据列。 在这里，我们将使用两个样本元数据列将其应用于未加权的UniFrac距离，如下所示。

---

qiime diversity beta-group-significance \
--i-distance-matrix core-metrics-results/unweighted_unifrac_distance_matrix.qza \
--m-metadata-file sample-metadata.tsv \
--m-metadata-column BodySite \
--o-visualization core-metrics-results/unweighted-unifrac-body-site-significance.qzv \
--p-pairwise

qiime diversity beta-group-significance \
--i-distance-matrix core-metrics-results/unweighted_unifrac_distance_matrix.qza \
--m-metadata-file sample-metadata.tsv \
--m-metadata-column Subject \
--o-visualization core-metrics-results/unweighted-unifrac-subject-group-significance.qzv \
--p-pairwise

Output visualizations:

core-metrics-results/unweighted-unifrac-body-site-significance.qzv

core-metrics-results/unweighted-unifrac-subject-group-significance.qzv

---

同样，我们对此数据集的连续样本元数据都没有与样本组成相关联，因此我们不会在此测试这些关联。 如果您对执行这些测试感兴趣，可以将qiime元数据距离矩阵与qiime Diversity mantel和qiime diversity bioenv命令结合使用。

最后，排序是在样本元数据的背景下探索微生物群落组成的流行方法。 我们可以使用Emperor工具在样本元数据的上下文中探索主坐标（PCoA）图。 虽然我们的core-metrics-phylogenetic命令已经生成了一些Emperor图，但我们想要传递一个可选参数--p-custom-axes，这对于探索时间序列数据非常有用。 核心度量 - 系统发育中使用的PCoA结果也是可用的，这使得使用Emperor生成新的可视化变得容易。 我们将为未加权的UniFrac和Bray-Curtis生成Emperor图，以便生成的图将包含主坐标1，主坐标2和自实验开始以来的天数的轴。 我们将使用最后一个轴来探索这些样本如何随时间变化。

---

qiime emperor plot \
--i-pcoa core-metrics-results/unweighted_unifrac_pcoa_results.qza \
--m-metadata-file sample-metadata.tsv \
--p-custom-axes DaysSinceExperimentStart \
--o-visualization core-metrics-results/unweighted-unifrac-emperor-DaysSinceExperimentStart.qzv

qiime emperor plot \
--i-pcoa core-metrics-results/bray_curtis_pcoa_results.qza \
--m-metadata-file sample-metadata.tsv \
--p-custom-axes DaysSinceExperimentStart \
--o-visualization core-metrics-results/bray-curtis-emperor-DaysSinceExperimentStart.qzv

Output visualizations:

core-metrics-results/unweighted-unifrac-emperor-DaysSinceExperimentStart.qzv

core-metrics-results/bray-curtis-emperor-DaysSinceExperimentStart.qzv

---

10. 阿尔法绘制稀疏

在本节中，我们将使用qiime多样性alpha-rarefaction可视化器探索alpha多样性作为采样深度的函数。 该可视化器在多个采样深度处计算一个或多个α分集度量，步长在1之间（可选地用-p-min-depth控制）和作为-p-max-depth提供的值。 在每个采样深度步骤，将生成10个稀疏表，并且将为表中的所有样本计算分集度量。 可以使用-p-iterations控制迭代次数（在每个采样深度计算的稀疏表）。 将在每个均匀采样深度处为每个样本绘制平均多样性值，并且如果样本元数据与--m-metadata-file参数一起提供，则可以基于结果可视化中的元数据对样本进行分组。

---

qiime diversity alpha-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

Output visualizations:

alpha-rarefaction.qzv

---

可视化将有两个图。 顶部图是α稀疏图，主要用于确定样品的丰富程度是否已被完全观察或测序。 如果图中的线条在沿x轴的某个采样深度处看起来“平整”（即接近零斜率），则表明收集超出该采样深度的其他序列将不太可能导致观察 其他功能。 如果图中的线没有达到平衡，这可能是因为尚未完全观察到样品的丰富程度（因为收集的序列太少），或者它可能表明仍存在大量测序错误 在数据中（被误认为是新颖的多样性）。

在按元数据对样本进行分组时，此可视化中的底部图非常重要。 它说明了当特征表稀疏到每个采样深度时每个组中保留的样本数。 如果给定的采样深度d大于样本s的总频率（即，针对样本s获得的序列的数量），则不可能在采样深度d处计算样本s的分集度量。 如果一组中的许多样本的总频率低于d，那么在顶部图中d处给出的该组的平均多样性将是不可靠的，因为它将在相对较少的样本上计算。 因此，在按元数据对样本进行分组时，必须查看底部图，以确保顶部图中显示的数据是可靠的。

11. 分类学分析

在接下来的部分中，我们将开始探索样本的分类组成，并再次将其与样本元数据联系起来。 此过程的第一步是将分类法分配给FeatureData [Sequence] QIIME 2工件中的序列。 我们将使用预先训练的Naive Bayes分类器和q2-feature-classifier插件来实现。 该分类器在Greengenes 13_8 99％OTU上进行训练，其中序列已被修剪为仅包括来自在该分析中测序的16S区域的250个碱基（V4区域，由515F / 806R引物对结合）。 我们将这个分类器应用于我们的序列，并且我们可以生成从序列到分类的结果映射的可视化。

---

wget \

-O "gg-13-8-99-515-806-nb-classifier.qza" \
"https://data.qiime2.org/2019.4/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 metadata tabulate \
--m-input-file taxonomy.qza \
--o-visualization taxonomy.qzv

Output artifacts:

taxonomy.qza

gg-13-8-99-515-806-nb-classifier.qza

Output visualizations:

taxonomy.qzv

---

接下来，我们可以使用交互式条形图查看样本的分类组成。 使用以下命令生成这些图，然后打开可视化。

---

qiime taxa barplot \
--i-table table.qza \
--i-taxonomy taxonomy.qza \
--m-metadata-file sample-metadata.tsv \
--o-visualization taxa-bar-plots.qzv

Output visualizations:

taxa-bar-plots.qzv

---

12. ANCOM差分丰度

ANCOM差分丰度测试ANCOM可用于识别样本组间差异丰富的特征(即存在于不同的丰度中)。与任何生物信息学方法一样，在使用ANCOM之前，您应该了解它的假设和限制。我们建议在使用此方法之前先阅读ANCOM论文。

ANCOM在q2-composition插件中实现。 ANCOM假设很少（少于约25％）的功能在组之间发生变化。 如果您希望组之间的更多功能发生变化，则不应使用ANCOM，因为它更容易出错（可能会增加类型I和类型II错误）。 因为我们期望很多特征在身体部位中大量变化，所以在本教程中我们将过滤我们的完整特征表以仅包含肠道样本。 然后我们将应用ANCOM来确定我们两个受试者的肠道样本中哪些（如果有的话）序列变体和属是差异丰富的。

---

qiime feature-table filter-samples \

--i-table table.qza \
--m-metadata-file sample-metadata.tsv \
--p-where "BodySite='gut'" \
--o-filtered-table gut-table.qza

Output artifacts:

gut-table.qza

---

ANCOM在FeatureTable [Composition] QIIME 2工件上运行，该工件基于每个样本的特征频率，但不能容忍零频率。 要构建合成工件，必须提供FeatureTable [Frequency]工件来添加伪计数（插补方法），这将生成FeatureTable [合成]工件

---

qiime composition add-pseudocount \
--i-table gut-table.qza \
--o-composition-table comp-gut-table.qz

Output artifacts:

comp-gut-table.qza

---

然后我们可以在“主题”列上运行ANCOM，以确定两个受试者的肠道样本中丰富度的特征。

---

qiime composition ancom \
--i-table comp-gut-table.qza \
--m-metadata-file sample-metadata.tsv \
--m-metadata-column Subject \
--o-visualization ancom-Subject.qzv

Output visualizations:

ancom-Subject.qzv

---

我们还经常对在特定的分类水平上执行差异丰度测试感兴趣。为此，我们可以在感兴趣的分类级别上折叠我们的FeatureTable[频率]中的特性，然后重新运行上面的步骤。在本教程中，我们将折叠属级别的特性表(即Greengenes分类法的第6级)。

---

qiime taxa collapse \
--i-table gut-table.qza \
--i-taxonomy taxonomy.qza \
--p-level 6 \
--o-collapsed-table gut-table-l6.qza

qiime composition add-pseudocount \
--i-table gut-table-l6.qza \
--o-composition-table comp-gut-table-l6.qza

qiime composition ancom \
--i-table comp-gut-table-l6.qza \
--m-metadata-file sample-metadata.tsv \
--m-metadata-column Subject \
--o-visualization l6-ancom-Subject.qz

Output artifacts:

gut-table-l6.qza

comp-gut-table-l6.qza

Output visualizations:

l6-ancom-Subject.qzv

---