- 1.所有的宏基因组/扩增子测序的分析的起点是原始序列数据。原始数据多数为fastq格式，其包含有DNA序列数据和每个碱基的质量值。
- 2.我们必须进行混池的样本拆分（demultiplex），以便确定每条序列来自于哪个样本。
- 3.然后进行序列去噪（denoised）以获得扩增子序列变异（amplicon sequence variants, ASVs），或聚类为可操作分类单元（operational taxonomic units, OTUs），这样做目的有二个：⑴降低测序错误；⑵序列去冗余。
- 4.特征表和代表性序列是去噪获得的关键结果，是下游分析的核心数据。一个特征表简单说就是一张Excel表，行名为ASV或OTU名称，列名为样本名。
- 5.我们可以基于特征表做很多事，常用分析包括：
- 序列的物种分类。比如，这里面有什么物种？
- Alpha(α)和Beta(β)多样性分析，即分别描述样本内或样本间的多样性。比如我们可以了解样本间有多少物种是一样的，即相似性如何？
- 许多的多样性分析依赖于个体特征的进化相似性。如果你测序的是系统发育的标记基因，如16S rRNA基因，你可以采用多序列比对方式评估特征间的系统发育关系。
- 不同实验组间特征的差异丰度分析，确定哪些OTUs或ASVs显著的多或少。
2.双端数据分析流程¶
2.样本拆分¶
接下来，我们按Barcode序列信息进行样品拆分。这需要样本元数据文件，必须指明该文件中的哪个列包含每个样本的条形码。 在该样本元数据文件中，该列名称是条形码序列BarcodeSequence。
在数据集中，条形码读长是样本元数据文件中包含的条形码读长的反向互补序列，因此我们还包括--p-rev-comp-mapping-barcodes参数。 在样品拆分之后，我们可以生成并查看每个样本获得多少序列的摘要。
QIIME 2样本拆分过程见下面的流程图（暂时忽略此流程图的右侧）：
样本拆分，EMP双端序列类型，指定元数据，barcode所在列，¶
并序列取反向互补(barcode加有右端时使用)，i为输出文件，o为输出文件和统计¶
time qiime demux emp-paired \ --m-barcodes-file Mapping_file_second_run.txt \ --m-barcodes-column BarcodeSequence \ --p-rev-comp-mapping-barcodes \ --i-seqs emp-paired-end-sequences.qza \ --o-per-sample-sequences demux.qza \ --o-error-correction-details demux-details.qza
输出对象：
- demux.qza: 样品拆分结果文件。
- demux-details.qza：样品拆分统计文件。
摘要可视化为表格¶
qiime demux summarize \ --i-data demux.qza \ --o-visualization demux.qzv
输出可视化：
- demux.qzv: 样本拆分结果可视化。
网页中交互式图形可以查看每个碱基位置的详细信息。质量分析后，我们根据上图结果和相关表格来确定下步denoise分析参数。
3. DADA2降噪¶
由样本拆分结果可视化demux.qzv的结果，可以确定降噪的参数--p-trim-left-f和--p-trim-left-r以及--p-trunc-len-f和--p-trunc-len-r为0，0，250，250。
time qiime dada2 denoise-paired \ --i-demultiplexed-seqs demux.qza \ --p-trim-left-f 0 \ --p-trim-left-r 0 \ --p-trunc-len-f 250 \ --p-trunc-len-r 250 \ --o-table table.qza \ --o-representative-sequences rep-seqs.qza \ --o-denoising-stats denoising-stats.qza \ --p-n-threads 0
输出结果为：
- table.qza: 特征表。
- rep-seqs.qza: 代表序列。
- denoising-stats.qza: 去噪过程统计。
Feature表可视化¶
time qiime feature-table summarize \ --i-table table.qza \ --o-visualization table.qzv \ --m-sample-metadata-file Mapping_file_second_run.txt
代表序列统计可视化¶
time qiime feature-table tabulate-seqs \ --i-data rep-seqs.qza \ --o-visualization rep-seqs.qzv
去噪过程统计数据可视化¶
time qiime metadata tabulate \ --m-input-file denoising-stats.qza \ --o-visualization denoising-stats.qzv
可视化结果：
- table.qzv：特征表可视化。
- rep-seqs.qzv：代表序列可视化。
- denoising-stats.qzv：去噪过程统计数据可视化。
table.qzv：特征表可视化
rep-seqs.qzv：代表序列可视化
denoising-stats.qzv：去噪过程统计数据可视化。
导出table.qza的步骤
1.将table.qza文件导出为biom格式
qiime tools export \ --input-path table.qza \ --output-path exported-feature-table
导出的biom文件位于exported-feature-table文件夹中， 名为feature-table.biom，可用biom程序对文件进行格式转换和分析
2.将taxonomy.qza文件导出为tsv格式
qiime tools export \ --input-path taxonomy.qza \ --output-path exported-feature-table
3.append taxonomy classification to feature table and convert to biom format
biom add-metadata -i feature-table.biom -o table.w_omd.biom --observation-metadata-fp taxonomy.tsv --observation-header OTUID,taxonomy --sc-separated taxonomy
4.将table.w_omd.biom转换为table.from_biom_w_taxonomy.txt格式
biom convert -i table.w_omd.biom -o table.from_biom_w_taxonomy.txt --to-tsv --header-key taxonomy
4.构建进化树¶
time 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
输出结果：
- aligned-rep-seqs.qza
- masked-aligned-rep-seqs.qza
- unrooted-tree.qza
- rooted-tree.qza
5. 多样性分析¶
划重点：理解下面4种alpha和beta多样性指数的所代表的生物学意义至关重要。
α多样性
- [x] 香农(Shannon’s)多样性指数（群落丰富度的定量度量，即包括丰富度richness和均匀度evenness两个层面）
- [x] 可观测的OTU(Observed OTUs，群落丰富度的定性度量，只包括丰富度）
- [x] Faith’s系统发育多样性（包含特征之间的系统发育关系的群落丰富度的定性度量）
- [x] 均匀度Evenness（或 Pielou’s均匀度；群落均匀度的度量）
β多样性
- [x] Jaccard距离（群落差异的定性度量，即只考虑种类，不考虑丰度）
- [x] Bray-Curtis距离（群落差异的定量度量，较常用）
- [x] 非加权UniFrac距离（包含特征之间的系统发育关系的群落差异定性度量）
- [x] 加权UniFrac距离（包含特征之间的系统发育关系的群落差异定量度量）
计算多样性(包括所有常用的Alpha和Beta多样性方法)，输入有根树、Feature表、样本重采样深度¶
取样深度看table.qzv文件确定（一般为样本最小的sequence count，或覆盖绝大多数样品的sequence count）¶
time qiime diversity core-metrics-phylogenetic \ --i-phylogeny rooted-tree.qza \ --i-table table.qza \ --p-sampling-depth 15000 \ --m-metadata-file Mapping_file_second_run.txt \ --output-dir core-metrics-results
输出结果包括多种多样性结果，文件列表和解释如下：¶
输出对象(13个数据文件):
core-metrics-results/faith_pd_vector.qza: 考虑进化的faith指数(Alpha多样性)。
core-metrics-results/observed_otus_vector.qza: observed otus指数(Alpha多样性)。
core-metrics-results/evenness_vector.qza: 均匀度指数(Alpha多样性)。
core-metrics-results/shannon_vector.qza: Alpha多样性香农指数(Alpha多样性)。
core-metrics-results/bray_curtis_distance_matrix.qza: Bray-Curtis距离矩阵(Beta多样性)。
core-metrics-results/jaccard_distance_matrix.qza: jaccard距离矩阵(Beta多样性)。
core-metric-results/unweighted_unifrac_distance_matrix.qza: 无权重unifrac距离矩阵(Beta多样性)。
core-metrics-results/weighted_unifrac_distance_matrix.qza: 有权重的unifrac距离矩阵(Beta多样性)。
- core-metrics-results/rarefied_table.qza: 等量重采样后的特征表。
core-metrics-results/bray_curtis_pcoa_results.qza: 基于Bray-Curtis距离PCoA的结果。
core-metrics-results/jaccard_pcoa_results.qza: jaccard距离PCoA结果。
core-metrics-results/unweighted_unifrac_pcoa_results.qza: 无权重的unifrac距离的PCoA结果。
core-metrics-results/weighted_unifrac_pcoa_results.qza: 基于有权重的unifrac距离的PCoA结果。
输出对象(4种可视化结果):
core-metrics-results/unweighted_unifrac_emperor.qzv： 无权重的unifrac距离PCoA结果采用emperor可视化。
core-metrics-results/jaccard_emperor.qzv： jaccard距离PCoA结果采用emperor可视化。
core-metrics-results/bray_curtis_emperor.qzv： Bray-Curtis距离PCoA结果采用emperor可视化。
core-metrics-results/weighted_unifrac_emperor.qzv： 有权重的unifrac距离PCoA结果采用emperor可视化。
我们对元数据文件Mapping_file_second_run.txt 进行了处理，增加了三列source(样本来源：flower、leaf、seed)，species(ZS4,ZS9,没有该指标的样本，该列设为NA),position(Upper,Lower，没有该指标的样本，该列设为NA)，此时样本元数据文件命名为sample_metadata.txt 。
5.1 Alpha多样性分析¶
使用alpha-group-significance测试alpha多样性数据和分类元数据列之间的关联¶
查看帮助¶
qiime diversity alpha-group-significance --help
输入多样性值、sample_medata.txt，输出统计结果¶
统计evenness组间差异是否显著¶
time qiime diversity alpha-group-significance \ --i-alpha-diversity core-metrics-results/evenness_vector.qza \ --m-metadata-file sample_metadata.txt \ --o-visualization core-metrics-results/evenness-group-significance.qzv
统计faith_pd算法Alpha多样性组间差异是否显著¶
time qiime diversity alpha-group-significance \ --i-alpha-diversity core-metrics-results/faith_pd_vector.qza \ --m-metadata-file sample_metadata.txt \ --o-visualization core-metrics-results/faith-pd-group-significance.qzv
统计observed_otus指数组间差异是否显著¶
time qiime diversity alpha-group-significance \ --i-alpha-diversity core-metrics-results/observed_otus_vector.qza \ --m-metadata-file sample_metadata.txt \ --o-visualization core-metrics-results/observed_otus_group-significance.qzv
统计shannon熵指数组间差异是否显著¶
time qiime diversity alpha-group-significance \ --i-alpha-diversity core-metrics-results/shannon_vector.qza \ --m-metadata-file sample_metadata.txt \ --o-visualization core-metrics-results/shannon_group-significance.qzv
- evenness-group-significance.qzv
Kruskal-Wallis(KW)检验是一个关于三组或更多数据的非参数性检验。它是用来检验总体函数分布的一致性原假设和其替代假设，关于至少两个样本之间存在差异的假设。
- faith-pd-group-significance.qzv
- observed_otus_group-significance.qzv
- shannon_group-significance.qzv
5.2 Beta多样性分析¶
分别按source、species、position分组，统计unweighted_unifrace距离的组间是否有显著差异¶
time qiime diversity beta-group-significance \ --i-distance-matrix core-metrics-results/unweighted_unifrac_distance_matrix.qza \ --m-metadata-file sample_metadata.txt \ --m-metadata-column source \ --o-visualization core-metrics-results/unweighted-unifrac-source-significance.qzv \ --p-pairwise
time qiime diversity beta-group-significance \ --i-distance-matrix core-metrics-results/unweighted_unifrac_distance_matrix.qza \ --m-metadata-file sample_metadata.txt \ --m-metadata-column species \ --o-visualization core-metrics-results/unweighted-unifrac-species-significance.qzv \ --p-pairwise
time qiime diversity beta-group-significance \ --i-distance-matrix core-metrics-results/unweighted_unifrac_distance_matrix.qza \ --m-metadata-file sample_metadata.txt \ --m-metadata-column position \ --o-visualization core-metrics-results/unweighted-unifrac-position-significance.qzv \ --p-pairwise
输出对象：
- unweighted-unifrac-source-significance.qzv
- unweighted-unifrac-species-significance.qzv
- unweighted-unifrac-position-significance.qzv
置换多元方差分析（Permutational multivariate analysis of variance，PERMANOVA），又称非参数多因素方差分析（nonparametric multivariate analysis of variance）、或者ADONIS分析，其本质是基于F统计的方差分析，依据距离矩阵对总方差进行分解的非参数多元方差分析方法。使用PERMANOVA可分析不同分组因素对样品差异的解释度，并使用置换检验进行显著性统计。
PERMANOVA也常应用在物种多度数据的β多样性分析中。例如在PCA、PCoA、NMDS等排序分析中，通过排序图显示了不同样品之间群落结构区分明显。但由于这些非约束排序仅为探索性分析，不涉及统计检验，无法提供一个明确的指示值，告诉我们所观察到的差异是否是有效的，或者差异的程度。这时，即可通过PERMANOVA确定群落结构差异的显著性。同样地，对于无法在排序图中区分的分组，也可通过PERMANOVA评估潜在的差异。
unweighted-unifrac-species-significance.qzv
unweighted-unifrac-position-significance.qzv
- alpha-rarefaction.qzv
生成物种可视化，即每个Feature对应的物种注释和可信度¶
time qiime metadata tabulate \ --m-input-file taxonomy.qza \ --o-visualization taxonomy.qzv
9. ANCOM差异度分析¶
按属比较，需要先合并¶
time qiime taxa collapse \ --i-table table.qza \ --i-taxonomy taxonomy.qza \ --p-level 2 \ --o-collapsed-table table-l2.qza
格式转换¶
time qiime composition add-pseudocount \ --i-table table-l2.qza \ --o-composition-table comp-table-l2.qza
time qiime composition ancom \ --i-table comp-table-l2.qza \ --m-metadata-file sample_metadata.txt \ --m-metadata-column source \ --o-visualization l6-ancom-source.qzv
分类学差异直接有名称，不用feature再对应物种注释¶
time qiime composition ancom \ --i-table comp-table-l2.qza \ --m-metadata-file sample_metadata.txt \ --m-metadata-column species \ --o-visualization l6-ancom-species.qzv
time qiime composition ancom \ --i-table comp-table-l2.qza \ --m-metadata-file sample_metadata.txt \ --m-metadata-column position \ --o-visualization l6-ancom-position.qzv
FDR（False Discovery Rate）
- 在差异表达分析过程中，采用统计学方法对原有假设检验得到的p值进行校正，并最终采用校正后的p值，即为FDR。
FC（Fold Change）
- 是两个/组样品间基因表达水平的比值，是表示差异倍数的变量。
一般差异表达分析中会同时控制这两个参数来筛选显著差异表达基因，如取FDR<0.05和FC≧2。
火山图可以直观地展现所有基因的FDR与Fold Change之间的关系，以便快速查看基因在两组样品间表达水平的差异程度和其统计学显著性。火山图的本质其实就是散点图。
横坐标表示某一个基因在两样品中表达水平差异倍数的对数值，即log2(FC)；横坐标绝对值越大，说明基因表达量在两样品间的倍数差异越大。纵坐标表示FDR的负对数值，即-log10(FDR)。纵坐标值越大，表明差异表达越显著，筛选得到的差异表达基因越可靠。
l6-ancom-species.qzv
l6-ancom-position.qzv
元数据¶
1.使用QIIME2对象作为元数据
微生物组分析的Alpha、Beta多样性值，在本质上，也是样本的一种属性，可加入元数据表中，方便后续分析。
将此对象faith_pd_vector.qza视为元数据，只需将其传递给希望看到元数据的任何方法或可视化工具（例如，metadata tabulate或emperor plot）：
qiime metadata tabulate \ --m-input-file faith_pd_vector.qza \ --o-visualization tabulated-faith-pd-metadata.qzv
可视化对象
- tabulated-faith-pd-metadata.qzv：多样性元数据可视化。
2.合并原数据
由于元数据可以来自许多不同的源，因此在运行命令时，QIIME 2支持元数据合并。基于上面的示例，只需多次传递--m-input-file就可以将元数据列组合到指定的文件中：
qiime metadata tabulate \ --m-input-file sample_metadata.txt \ --m-input-file core-metrics-results/faith_pd_vector.qza \ --m-input-file core-metrics-results/evenness_vector.qza \ --o-visualization core-metrics-results/tabulated-combined-metadata.qzv
- tabulated-combined-metadata.qzv：多样性元数据可视化。
合并后生成的元数据将包含所有指定文件中标识符的交集。换句话说，合并的元数据将只包含在所有提供的元数据文件中共享的标识符。这是一个使用数据库术语的内部联接。
在QIIME 2中接受元数据的任何地方都支持元数据合并。例如，根据研究元数据或样本alpha多样性为Emperor plot着色可能很有意思。这可以通过提供样本元数据文件和SampleData[AlphaDiversity]对象来实现：
组合Alpha和Beta多样性¶
qiime emperor plot \ --i-pcoa core-metrics-results/unweighted_unifrac_pcoa_results.qza \ --m-metadata-file sample_metadata.txt \ --m-metadata-file core-metrics-results/faith_pd_vector.qza \ --o-visualization core-metrics-results/unweighted-unifrac-emperor-with-alpha.qzv
可视化对象
- unweighted-unifrac-emperor-with-alpha.qzv：Beta多样性可视化添加Alpha多样性。
3.探索特征元数据
qiime metadata tabulate \ --m-input-file rep-seqs.qza \ --m-input-file taxonomy.qza \ --o-visualization tabulated-feature-metadata.qzv
可视化对象
- tabulated-feature-metadata.qzv：特征元数据，包括序列、物种注释和置信度。
QIIME1 用户注意： QIIME 2对象FeatureTable[Frequency]等价于QIIME 1 OTU或BIOM表，QIIME 2对象FeatureData[Sequence]等价于QIIME 1代表序列文件。由于DADA2和Deblur产生的“OTU”是通过对唯一序列进行分组而创建的，因此这些OTU相当于来自QIIME 1的100%相似度的OTU，通常称为序列变体。在QIIME 2中，这些OTU比QIIME 1默认的97%相似度聚类的OTU具有更高的分辨率，并且它们具有更高的质量，因为这些质量控制步骤比QIIME 1中实现更好。因此，与QIIME 1相比，可以对样本的多样性和分类组成进行更准确的估计。