扩增子分析流程

生信菜鸟团

发布于 2020-06-17 11:14:34

4.3K0

发布于 2020-06-17 11:14:34

文章被收录于专栏：生信菜鸟团

本文翻译自 https://github.com/LangilleLab/microbiome_helper/wiki/Amplicon-SOP-v2-(qiime2-2020.2) 并适当加入自己的理解。

该 SOP 基于 QIIME2 2020.2，学习之前建议先过一遍 QIIME2 “Moving Pictures” tutorial[1]。

准备工作

安装 QIIME 2，FASTQC 和 MultiQC。

# qiime2
wget https://data.qiime2.org/distro/core/qiime2-2020.2-py36-osx-conda.yml
conda env create -n qiime2-2020.2 --file qiime2-2020.2-py36-osx-conda.yml
conda activate qiime2-2020.2

# fastqc
conda install -c bioconda fastqc -y

# multiqc
pip install multiqc

或者也可以直接 Pull 我创建的 docker 镜像，包含了这三个软件：

docker pull zwbao/ampliconsop:2020.2

1. 第一步

1.1 格式化 metadata 文件

根据 QIIME2 所需的格式定义 metadata（样本信息文件），后续 QIIME 2 将根据这些分组信息进行分析。示例如下：

sample-id    barcode-sequence    body-site    year    month    day    subject    reported-antibiotic-usage    days-since-experiment-start
#q2:types    categorical    categorical    numeric    numeric    numeric    categorical    categorical    numeric
L1S8    AGCTGACTAGTC    gut    2008    10    28    subject-1    Yes    0
L1S57    ACACACTATGGC    gut    2009    1    20    subject-1    No    84
L1S76    ACTACGTGTGGT    gut    2009    2    17    subject-1    No    112
L1S105    AGTGCGATGCGT    gut    2009    3    17    subject-1    No    140
L2S155    ACGATGCGACCA    left palm    2009    1    20    subject-1    No    84
L2S175    AGCTATCCACGA    left palm    2009    2    17    subject-1    No    112
L2S204    ATGCAGCTCAGT    left palm    2009    3    17    subject-1    No    140
L2S222    CACGTGACATGT    left palm    2009    4    14    subject-1    No    168

以 # 开头的为注释行，会被 QIIME 2 忽略。关于 metadata 详细说明请查阅 QIIME 2 文档：https://docs.qiime2.org/2020.2/tutorials/metadata/

设置环境变量

METADATA="/home/user/metadata.txt"

1.2 设置运行核数

NCORES=1

1.3 检查 reads 质量

可视化原始 reads 的序列质量，测序深度等其他指标对后续的质控非常重要。虽然 QIIME 2 自带可视化测序质量的插件，但一般来说直接用 FASTQC + MultiQC 的组合更为方便且详细。

mkdir fastqc_out
fastqc -t $NCORES raw_data/*.fastq.gz -o fastqc_out

FASTQC 为每个样本单独生成一个报告，用 MultiQC 对这些网页报告进行合并：

cd fastqc_out
multiqc .
cd ..

完整的报告为 multiqc_report.html ，直接用浏览器即可查看。生成可视化报告最重要的原因是要确保样品是否高质量（大多数 reads 中每碱基的质量是否大于 30），是否存在异常的样本。

1.4 激活 QIIME 2 conda 环境

在接下来的步骤中，我们就需要用到 QIIME2 了。若之前用 conda 安装了 QIIME2，可用以下命令激活 conda 环境：

conda activate qiime2-2020.2
# 退出环境：conda deactivate

若使用 docker ，可用如下命令创建容器：

docker run -t -i -v $(pwd):/data zwbao/ampliconsop:2020.2 /bin/bash

1.5 导入 fastq 文件

为了标准化 QIIME 2 分析并跟踪分析过程，所有 QIIME 2 输入和输出文件都用到了一种特殊的格式，即 artifact（扩展名为.qza）。所以第一步就是将 .fastq 导入为 .qza。

mkdir reads_qza

qiime tools import \
   --type SampleData[PairedEndSequencesWithQuality] \
   --input-path raw_data/ \
   --output-path reads_qza/reads.qza \
   --input-format CasavaOneEightSingleLanePerSampleDirFmt

注：以上命令仅适用于 paired-end MiSeq 数据（FASTQ 格式，位于 raw_data 文件夹中）导入，同时每个测序文件都需命名为特定格式，例如：105CHE6WT_S325_L001_R1_001.fastq.gz，包含五个字段，用 _ 字符分隔，才能被命令所识别。

1.The sample name ( 105CHE6WT )2.The sample number on the run ( S325 )3.The lane number ( L001 )4.The read number ( R1 )5.The set number ( 001 )

若命名稍有不同，也可用其他 QIIME2 命令来导入数据。另一种常用的导入数据方式是写一个包含样本名和测序数据绝对路径的文本（代码中命名为 pe-33-manifest），然后用该文本导入。文本内容如下：

sample-id     forward-absolute-filepath       reverse-absolute-filepath
sample-1      $PWD/some/filepath/sample0_R1.fastq.gz  $PWD/some/filepath/sample1_R2.fastq.gz
sample-2      $PWD/some/filepath/sample2_R1.fastq.gz  $PWD/some/filepath/sample2_R2.fastq.gz
sample-3      $PWD/some/filepath/sample3_R1.fastq.gz  $PWD/some/filepath/sample3_R2.fastq.gz
sample-4      $PWD/some/filepath/sample4_R1.fastq.gz  $PWD/some/filepath/sample4_R2.fastq.gz

用代码批量导入：

qiime tools import \
  --type 'SampleData[PairedEndSequencesWithQuality]' \
  --input-path pe-33-manifest \
  --output-path reads_qza/reads.qza \
  --input-format PairedEndFastqManifestPhred33V2

所有 .fastq 测序数据都将保存在 reads_qza/reads.qza 文件中。QIIME 2 支持多种数据导入方式，详见：https://docs.qiime2.org/2020.2/tutorials/importing/ 。

1.6 用 cutadapt 插件去除引物序列

下面的引物对应于16S V6-V8 区域，用的时候应修改为适合自己实验的引物：

qiime cutadapt trim-paired \
   --i-demultiplexed-sequences reads_qza/reads.qza \
   --p-cores $NCORES \
   --p-front-f ACGCGHNRAACCTTACC \
   --p-front-r ACGGGCRGTGWGTRCAA \
   --p-discard-untrimmed \
   --p-no-indels \
   --o-trimmed-sequences reads_qza/reads_trimmed.qza

1.7 可视化修剪后的序列

接下来我们可以用 demux summarize 命令可视化修剪之后测序数据，生成的报告内容主要包括每个样品的 reads 数量以及这些 reads 的质量分布。与 FASTQC + MultiQC 生成的报告相比，这虽有些简陋，但也足够。

qiime demux summarize \
   --i-data reads_qza/reads_trimmed.qza \
   --o-visualization reads_qza/reads_trimmed_summary.qzv

.qzv 文件是一种特殊的 artifact 文件，主要用来可视化分析的 summary。我们可以直接将它拖到 QIIME2 view 网页[2] 进行可视化。也可直接解压，点 index.html查看报告。

2. 用 deblur 将 reads 降噪为 ASV （amplicon sequence variants）

这一步我们可使用两种 pipeline，deblur 或 DADA2，任选其一即可。下面的步骤以 deblur 流程为例。

2.1 合并双端序列

目前 deblur 尚不支持未合并的双端序列，所以要先用 VSEARCH 合并双端序列。但该命令比 PEAR 等方法更为严格，若这一步输出的 reads 过少，我们也可以可导出后用 PEAR 进行合并，然后重新导入数据。

qiime vsearch join-pairs \
   --i-demultiplexed-seqs reads_qza/reads_trimmed.qza \
   --o-joined-sequences reads_qza/reads_trimmed_joined.qza

2.2 过滤低质量的 reads

基于默认选项过滤掉低质量的 reads：

qiime quality-filter q-score-joined \
   --i-demux reads_qza/reads_trimmed_joined.qza \
   --o-filter-stats filt_stats.qza \
   --o-filtered-sequences reads_qza/reads_trimmed_joined_filt.qza

对过滤后的数据进行统计和可视化，看看发生了哪些改变：

qiime demux summarize \
   --i-data reads_qza/reads_trimmed_joined_filt.qza \
   --o-visualization reads_qza/reads_trimmed_joined_filt_summary.qzv

2.3 运行 deblur

运行 deblur 流程得到 ASVs。下面的命令将保留 singletons，除非我们设置参数 --p-min-reads 1 。

注意，如果这步出现报错，需使用 --p-trim-length 参数将所有序列修剪到相同的长度。建议根据上面生成的 .qzv 网页报告确定正确的修剪长度。最后，在 deblur 选项中输入该长度 --p-trim-length，然后重新运行命令。

qiime deblur denoise-16S \
   --i-demultiplexed-seqs reads_qza/reads_trimmed_joined_filt.qza \
   --p-trim-length -1 \
   --p-sample-stats \
   --p-jobs-to-start $NCORES \
   --p-min-reads 1 \
   --output-dir deblur_output

2.4 对 deblur 流程的结果统计可视化

qiime feature-table summarize \
   --i-table deblur_output/table.qza \
   --o-visualization deblur_output/deblur_table_summary.qzv

在这个报告中，我们需要留意运行 deblur 后是否还保留足够多的 reads，稍后我们也将根据这个报告来确定过滤的 cut-off 值。该降噪工具会过滤与已知噪声匹配或与预期扩增子区域的相似性不匹配的 reads。若运行 deblur 后，样品的深度非常低（与输入深度相比），这表明，你可能哪里设置错了，或数据中可能包含大量噪音亦或是 deblur 并不适合于你的数据集。

3. 对 ASV 进行物种注释

在这一步中我们将使用 SILVA 数据库训练 Naive-Bayes 分类器来对 ASV 进行物种注释。

3.1 建立或获取分类器

这种方法要求事先基于参考数据库训练分类器。QIIME 2 团队建议为不同的引物组合建立专门的分类器。对于一些大家常用的引物组合，可直接在仓库中下载（ http://kronos.pharmacology.dal.ca/public_files/taxa_classifiers/qiime2-2020.2_classifiers/ ），没有的话则需要自己手动建立分类器：

•16S V4/V5 region (classifier_silva_132_99_16S_V4.V5_515F_926R.qza)•16S V3/V4 region (classifier_silva_132_99_16S_V3.V4_341F_805R.qza)•16S V6/V8 region (classifier_silva_132_99_16S_V6.V8_B969F_BA1406R.qza)•16S V6/V8 region targeting archaea (classifier_silva_132_99_16S_V6.V8_A956F_A1401R.qza)•16S V3/V4 region targeting cyanobacteria (classifier_silva_132_99_16S_V3.V4_CYA359F_CYA781R.qza)•18S V4 region (classifier_silva_132_99_18S_V4_E572F_E1009R.qza)•Full ITS - fungi only (classifier_sh_refs_qiime_ver8_99_s_02.02.2019_ITS.qza)•Full ITS - all eukaryotes (classifier_sh_refs_qiime_ver8_99_s_all_02.02.2019_ITS.qza)

仓库中的 SILVA 版本为 132，若想使用最新的 SILVA 138 版本可参考：https://github.com/mikerobeson/make_SILVA_db 关于如何手动根据自己的引物建立分类器请参阅：https://github.com/LangilleLab/microbiome_helper/wiki/Creating-QIIME-2-Taxonomic-Classifiers

此外，在使用这些自定义分类器时，我们应仔细检查它们在数据集上是否正确执行，手动检查分类器对 ASV 的分类尤为重要。理论上，使用特定于引物的分类器，可以改进物种注释的效果，但仍建议你在首次运行自定义 16S 分类器时同时运行全长 16S 分类器进行比较。

3.2 运行分类器

这一步是 SOP 中运行时间最长、内存用量最大的命令之一。

qiime feature-classifier classify-sklearn \
   --i-reads deblur_output/representative_sequences.qza \
   --i-classifier /home/shared/taxa_classifiers/qiime2-2020.2_classifiers/classifier_silva_132_99_16S.qza \
   --p-n-jobs $NCORES \
   --output-dir taxa

导出输出文件以查看分类和置信度得分：

qiime tools export \
   --input-path taxa/classification.qza --output-path taxa

3.3 用 BLAST 评估物种注释

我们可以把一些 ASV 的注释结果与其 BLASTn 的结果进行比较，来大致评估分类器的性能：

qiime feature-table tabulate-seqs --i-data deblur_output/representative_sequences.qza \
                                  --o-visualization deblur_output/representative_sequences.qzv

用 QIIME2 view 打开 representative_sequences.qzv，可以看到 ASV 的序列。可在不同门类中选择大约 5 个 ASV 用 BLAST 进行验证。

4. 过滤结果

对去噪后的结果进行筛选是微生物组数据分析的一个重要步骤。你可以在 QIIME 2过滤教程[3] 中看到关于这个过程的更多细节。

4.1 过滤掉罕见的 ASVs

我们可以根据步骤 2.4 中创建的摘要可视化，选择一个 cut-off 值。比如我们可以移除所有频率小于平均采样深度 0.1% 的 ASVs。要计算这个 cut-off 值，你需要查看步骤 2.4 生成的报告，确定平均样本深度，乘以0.001，舍入最近的整数即可。

一旦确定了 cut-off 值，可以用下面的命令进行过滤 :

qiime feature-table filter-features \
   --i-table deblur_output/table.qza \
   --p-min-frequency X \
   --p-min-samples 1 \
   --o-filtered-table deblur_output/deblur_table_filt.qza

4.2 过滤污染物及未分类的 ASVs

既然我们已经注释了 ASV，那就可以用这些信息来过滤一些污染物及未分类的 ASV。在16S 测序数据中，两种常见的污染物是线粒体和叶绿体的16S 序列，这些序列可以通过过滤他们的分类标签来去除。另外，我们还需要过滤在门类级别未分类的 ASV，因为这些序列更有可能是噪音（例如可能为嵌合体序列）。注意，如果你的数据不是基于 SILVA 数据库进行分类的，那么就需要将命令中的 D_1__ 更改为特定字符串，以便能够识别 phylum-level 的注释，或者干脆省略该行。当然，如果你的目的是研究环境中的新物种，也可以省略这一行。

qiime taxa filter-table \
   --i-table deblur_output/deblur_table_filt.qza \
   --i-taxonomy taxa/classification.qza \
   --p-include D_1__ \
   --p-exclude mitochondria,chloroplast \
   --o-filtered-table deblur_output/deblur_table_filt_contam.qza

4.3 排除低深度的样本

经过这些过滤步骤后，通常某些样品的深度会很低，那么这些样品也可以排除在下游分析之外，因为它们会在很大程度上增加噪音。我们通常使用的最小 cut-off 在 1000 到 4000 reads 范围内。如果你希望保留所有样本，除了那些完全失败的样本( depth < 50 reads )之外，那就直接用一个较低的 cut-off 值即可。

首先，我们需要绘制稀释曲线，以确定在什么测序深度下样本的丰富度达到瓶颈，然后选择一个尽可能接近这个值的 cut-off 值。

绘制稀疏曲线前需要先可视化之前一步生成的数据：

qiime feature-table summarize \
   --i-table deblur_output/deblur_table_filt_contam.qza \
   --o-visualization deblur_output/deblur_table_filt_contam_summary.qzv

从这个表中，我们可以确定样品的最大深度。然后用这个值设置 --p-max-depth 参数，绘制稀疏曲线：

qiime diversity alpha-rarefaction \
   --i-table deblur_output/deblur_table_filt_contam.qza \
   --p-max-depth X \
   --p-steps 20 \
   --p-metrics 'observed_otus' \
   --o-visualization rarefaction_curves_test.qzv

观察这些曲线，以确定保留样本的最小深度，设置 --p-min-frequency 进行过滤:

qiime feature-table filter-samples \
   --i-table deblur_output/deblur_table_filt_contam.qza \
   --p-min-frequency SET_CUTOFF \
   --o-filtered-table deblur_output/deblur_table_final.qza

如果你不想排除任何样本，则可使用上一步的输出 deblur_table_filt_contam.qza 作为最终结果文件（cp deblur_output/deblur_table_filt_contam.qza deblur_output/deblur_table_final.qza）。

4.4 取子集和汇总过滤后的结果

一旦我们有了最终筛选后的数据，我们将需要生成 ASV 序列的 QZA 格式子集。用如下命令从序列文件中排除低频率的 ASV：

qiime feature-table filter-seqs \
   --i-data deblur_output/representative_sequences.qza \
   --i-table deblur_output/deblur_table_final.qza \
   --o-filtered-data deblur_output/rep_seqs_final.qza

进行可视化：

qiime feature-table summarize \
   --i-table deblur_output/deblur_table_final.qza \
   --o-visualization deblur_output/deblur_table_final_summary.qzv

5. 使用 SEPP QIIME 2 插件对 ASV 建树

SEPP 是将短序列放入参考系统发生树的一种方法。使用 SEPP QIIME 2 插件对 ASV 建树：

qiime fragment-insertion sepp \
   --i-representative-sequences deblur_output/rep_seqs_final.qza \
   --i-reference-database /home/shared/rRNA_db/16S/sepp-refs-gg-13-8.qza \
   --o-tree asvs-tree.qza \
   --o-placements insertion-placements.qza \
   --p-threads $NCORES

sepp-refs-gg-13-8.qza 下载地址：https://library.qiime2.org/plugins/q2-fragment-insertion/16/ 。

6. 生成稀疏曲线

绘制所有样本的稀疏曲线：

qiime diversity alpha-rarefaction \
   --i-table deblur_output/deblur_table_final.qza \
   --p-max-depth X \
   --p-steps 20 \
   --i-phylogeny asvs-tree.qza \
   --m-metadata-file $METADATA \
   --o-visualization rarefaction_curves.qzv

7. 生成相对丰度堆叠柱状图

用以下命令可输出交互式的堆叠条形图：

qiime taxa barplot \
   --i-table deblur_output/deblur_table_final.qza \
   --i-taxonomy taxa/classification.qza \
   --m-metadata-file $METADATA \
   --o-visualization taxa/taxa_barplot.qzv

QIIME 2 默认是绘制每个样本的条形图。我们也可以根据 metadata 中的分类信息绘制堆叠柱状图：

qiime feature-table group \
   --i-table deblur_output/deblur_table_final.qza \
   --p-axis sample \
   --p-mode sum \
   --m-metadata-file $METADATA \
   --m-metadata-column CATEGORY \
   --o-grouped-table deblur_output/deblur_table_final_CATEGORY.qza

qiime taxa barplot \
   --i-table deblur_output/deblur_table_final_CATEGORY.qza \
   --i-taxonomy taxa/classification.qza \
   --m-metadata-file CATEGORY_METADATA \
   --o-visualization taxa/taxa_barplot_CATEGORY.qzv

8. 计算多样性及可视化

在计算这些指标之前，该命令还将稀释所有样本到样本排序深度(X 为最低合理样本深度; 将排除深度低于 cut-off 的样本)。

qiime diversity core-metrics-phylogenetic \
   --i-table deblur_output/deblur_table_final.qza \
   --i-phylogeny asvs-tree.qza \
   --p-sampling-depth X \
   --m-metadata-file $METADATA \
   --p-n-jobs $NCORES \
   --output-dir diversity

接下来我们可以比较不同组的 alpha 多样性并绘制箱线图。例如，生成箱线图比较 Shannon alpha 多样性：

qiime diversity alpha-group-significance \
   --i-alpha-diversity diversity/shannon_vector.qza \
   --m-metadata-file $METADATA \
   --o-visualization diversity/shannon_compare_groups.qzv

9. 用 ANCOM 识别差异特征

ANCOM[4] 是一种测试样本集合之间特征的相对f丰度差异的方法。但它要求所有的特征为非零值，所以首先需要添加一个假计数（一般选择1） :

qiime composition add-pseudocount \
   --i-table deblur_output/deblur_table_final.qza \
   --p-pseudocount 1 \
   --o-composition-table deblur_output/deblur_table_final_pseudocount.qza

用下面的命令运行 ANCOM，CATEGORY 为需要比较的分类。

qiime composition ancom \
   --i-table deblur_output/deblur_table_final_pseudocount.qza \
   --m-metadata-file $METADATA \
   --m-metadata-column CATEGORY \
   --output-dir ancom_output

10. 导出文件

输出 ASV 的 FASTA 文件:

qiime tools export \
   --input-path deblur_output/rep_seqs_final.qza \
   --output-path deblur_output_exported

导出 BIOM 表，并加入将物种分类注释信息：

sed -i -e '1 s/Feature/#Feature/' -e '1 s/Taxon/taxonomy/' taxa/taxonomy.tsv

qiime tools export \
   --input-path deblur_output/deblur_table_final.qza \
   --output-path deblur_output_exported

biom add-metadata \
   -i deblur_output_exported/feature-table.biom \
   -o deblur_output_exported/feature-table_w_tax.biom \
   --observation-metadata-fp taxa/taxonomy.tsv \
   --sc-separated taxonomy

biom convert \
   -i deblur_output_exported/feature-table_w_tax.biom \
   -o deblur_output_exported/feature-table_w_tax.txt \
   --to-tsv \
   --header-key taxonomy

11. 下游分析

虽然 QIIME 2 也可以进行一些数据可视化，但对于更多样化且专业的分析还是得用 R 来完成。

11.1 将 QIIME2 生成的 `.qza` 文件导入 R

安装 qiime2R ：

devtools::install_github("jbisanz/qiime2R")

导入 .qza 文件为 phyloseq 对象：

library("qiime2R")

physeq<-qza_to_phyloseq(
  features="./deblur_output/deblur_table_final.qza",
  tree="./asvs-tree.qza",
  taxonomy="./taxa/classification.qza",
  metadata = "./metadata.txt"
)

11.2 用 Phyloseq 包可视化微生物组数据

此部分内容参见我之前写的 Phyloseq 教程。引用链接

[1] QIIME2 “Moving Pictures” tutorial: https://docs.qiime2.org/2020.2/tutorials/moving-pictures/ [2] QIIME2 view 网页: https://view.qiime2.org/ [3] QIIME 2过滤教程: https://docs.qiime2.org/2020.2/tutorials/filtering/ [4] ANCOM: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4450248/

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2020-06-14，如有侵权请联系 cloudcommunity@tencent.com 删除

https

网络安全

html

本文分享自生信菜鸟团微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

https

网络安全

html

登录后参与评论

0 条评论

热度