博文
QIIME 2教程之生成特征表和特征序列(案例一)
|
案例一（地球微生组计划(EMP)单端数据）
数据来源研究文章《Moving pictures of the human microbiome》在2011年发表于Genome Biology。（QIIME 2在五个时间点对来自两个人（subject-1，2）四个身体部位（gut (feces), mouth(tongue), and skin (left and right palms)）的微生物组样本进行分析，只在第1个时间点使用抗生素。）
本教程中使用的数据基于Illumina HiSeq产出，使用地球微生物组计划扩增16S rRNA基因高变区4（V4）测序的方法。
（一） 样本元数据
开始分析之前，阅读样本元数据，熟悉本研究中使用的样本信息。元数据以制表符分隔的文本格式。（元数据之详情见：https://mp.weixin.qq.com/s/o9FLzogXGMurh9EcBGaQHA）
下载元数据：
wget -c http://210.75.224.110/github/QIIME2ChineseManual/2020.11/moving-pictures/sample-metadata.tsv
元数据可视化：
qiime metadata tabulate --m-input-file sample-metadata.tsv --o-visualization sample-metadata.qzv
（二） 下载和导入数据
1. 下载数据
本教程中，我们将处理完整的序列数据的一小部分，以便命令能够快速运行(减少等待时间)。
mkdir -p emp-single-end-sequences
# 3.6M
wget \
-O "emp-single-end-sequences/barcodes.fastq.gz" \
"https://data.qiime2.org/2020.11/tutorials/moving-pictures/emp-single-end-sequences/barcodes.fastq.gz"
# 24M
wget \
-O "emp-single-end-sequences/sequences.fastq.gz" \
"https://data.qiime2.org/2020.11/tutorials/moving-pictures/emp-single-end-sequences/sequences.fastq.gz"
2. 导入数据（导入到QIIME 2对象中）
EMPSingleEndSequences是包含多样本混合的序列文件，这意味着序列尚未分配给样本（因此包括sequences.fastq.gz和barcode.fastq.gz文件，其中barcode.fastq.gz包含与sequences.fastq.gz（302581条reads）中的每个序列相关联的条形码。3个文件关联如下：
time qiime tools import \
--type EMPSingleEndSequences \
--input-path emp-single-end-sequences \
--output-path emp-single-end-sequences.qza
输出对象：
emp-single-end-sequences.qza（包含barcode和序列的信息）
不同数据类型的导入，详见：https://blog.sciencenet.cn/home.php?mod=space&uid=994715&do=blog&id=1362408
（三） 拆分样本
EMP数据为混合样本序列，为了混合序列进行样本拆分，我们需要知道哪个条形码序列与每个样本相关联。此信息包含在样品元数据文件中。
time qiime demux emp-single \ --i-seqs emp-single-end-sequences.qza \ --m-barcodes-file sample-metadata.tsv \ --m-barcodes-column barcode-sequence \ --o-per-sample-sequences demux.qza \ --o-error-correction-details demux-details.qza 输出结果： demux-details.qza demux.qza | demux emp-single命令指的是这些序列是根据地球微生物组计划标准方法添加的条形码，并且是单端序列 --m-barcodes-column是元数据的barcodes列的表头名
demux.qza包含样本拆分后的序列 demux-details.qza包括Golay标签错误校正的详细，在本教程中不作讨论 |
* 拆分样本后，以下命令获得对应样本序列文件(sample.fastq.gz)以及清单文件(MANIFEST)
qiime tools export --input-path demux.qza --output-path exported-seqs
* 统计结果可视化
qime qiime demux summarize \
--i-data demux.qza \
--o-visualization demux.qzv
| 【Overview】 | 【Interactive Qaulity Plot】 |
|   |
主要分为三部分：上部为摘要；中部为样本不同数据量分布频率柱状图，可下载PDF，下部为每个样本的测序量。 | 263931个序列中10000个序列的随机抽样。 同样为三部分：上部为每个位置碱基的质量分布交互式箱线图，鼠标悬停在上面，即可在下面(中部)文字和表格中显示鼠标所在位置碱基质量的详细信息；下部为拆分样本的长度摘要(一般等长测序无差别)。 |
（四） 序列质控和生成特征表
QIIME 2插件多种质量控制方法可选，包括DADA2、Deblur和基于基本质量分数的过滤。
在本教程中，我们使用DADA2和Deblur两种方法分别介绍这个步骤。这两种方法的结果为2个：一个QIIME 2特征表FeatureTable[Frequency]和一个代表性序列FeatureData[Sequence]对象。1. Frequency对象包含数据集中每个样本中每个唯一序列的计数（频率）；2. Sequence对象将FeatureTable中的特征ID与序列对应。
【备注】此步主要有DADA2和Deblur两种方法可选，推荐使用DADA2(相较USEARCH的UPARSE算法，目前DADA2方法仅去噪去嵌合，不再按相似度聚类，结果与真实物种的序列更接近)。
1. 方法1-DADA2
time 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 生成三个输出文件 stats-dada2.qza dada2计算统计结果 table-dada2.qza 特征表 rep-seqs-dada2.qza 代表序列 特征表统计进行进行可视化 qiime metadata tabulate \ --m-input-file stats-dada2.qza \ --o-visualization stats-dada2.qzv | -p-trim-left m 去除每个序列的前m个碱基(如引物、标签序列barcode)。截取左端低质量序列，我们看上图中箱线图，左端质量都很高，无低质量区，设置为0 --p-trunc-len n 序列截取长度n，也是为了去除右端低质量序列，我们看到大于120以后，质量下降极大，甚至中位数都下降至20以下，需要全部去除，综合考虑决定设置为120 |
统计结果解释：
结果为每个样本，输入、过滤、去噪和非嵌合的统计，并支持按列排序，检索和功能，用于样本异常筛选，特征表抽平标准化等。
【补充】下游分析，将继续使用dada2的结果，需要将它们统一重命名，以便后续分析：
mv rep-seqs-dada2.qza rep-seqs.qza
mv table-dada2.qza table.qza
2. 方法2-Deblur
Deblur使用序列错误配置文件将错误的序列与从其来源的真实生物序列相关联，从而得到高质量的序列变异数据，主要为两个步骤。
* 第一步，应用基于质量分数的初始质量过滤过程，是Bokulich等人2013年发表的质量过滤方法。
time qiime quality-filter q-score \ --i-demux demux.qza \ --o-filtered-sequences demux-filtered.qza \ --o-filter-stats demux-filter-stats.qza 输出对象: demux-filtered.qza: 序列质量过滤后结果。 demux-filter-stats.qza: 序列质量过滤后结果统计。 统计结果可视化 qiime metadata tabulate \ --m-input-file demux-filter-stats.qza \ --o-visualization demux-filter-stats.qzv 输出结果 demux-filter-stats.qzv质量过程统计表，同上面提到的stats-dada2.qzv统计表类似。 【备注】在Deblur的论文中，作者使用了当时推荐的过滤参数。而这里使用的参数基于最新的经验，效果更好。 |
统计结果解释：
示例如上：包括6列，第1列为样本名称，2-6列分别为总输入读长、总保留高读长、截断的读长、截断后太短的读长和超过最大模糊碱基的读长的数量统计。我们通常只关注2，3列数量即可，其它列常用于异常的输助判断。
* 第二步，使qiime deblur denoise-16S方法应用于Deblur工作流程。deblur去噪16S过程，输入文件为质控后的序列，设置截取长度参数，生成结果文件有代表序列、特征表、样本统计。
time 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 输出结果: deblur-stats.qza 过程统计 table-deblur.qza特征表 rep-seqs-deblur.qza代表序列 结果可视化： qiime deblur visualize-stats \ --i-deblur-stats deblur-stats.qza \ --o-visualization deblur-stats.qzv 输出结果 demux-filter-stats.qzv 分析统计表，有分析中每个步骤的统计表 | 此方法需要一个用于质量过滤的参数，即截断位置n长度的序列的--p-trim-length n。通常，Deblur开发人员建议将该值设置为质量分数中位数开始下降至低质量区时的长度。在本次数据上，质量图（在质量过滤之前,第3部分箱线图）表明合理的选择是在115至130序列位置范围内。这是一个主观的评估。你可能不采用该建议的一种原因是存在多个批次测序的元分析。在这种情况的元分析中，比较所有批次的序列长度是否相同，以避免人为引入特定的偏差，全局考虑这些是非常重要的。由于我们已经使用修剪长度为120 bp用于qiime dada2 denoise-single分析，并且由于120 bp是基于质量图的结果，这里我们将使用--p-trim-length 120参数。下一个命令可能需要10分钟才能运行完成。 |
统计结果解释：
以上为deblur去噪和鉴定ASV处理过程统计结果，最终结果reads-hit-reference列（the number of Deblur reads which recruited to the postitive reference）。
【备注】deblur最大缺点就是慢，本次只分析了33个样品，共177,092条序列。而实际研究中大项目会有成千上万的样本，1亿-10亿条序列，此步分析可能需要几个月甚至根本无法完成，不推荐。
【补充】下游分析，将继续使用deblur的结果，需要将它们统一改名，以便继续分析：
mv rep-seqs-deblur.qza rep-seqs.qza
mv table-deblur.qza table.qza
【备注】
以上两种方法只选择一种即可。推荐dada2速度更快一些，步骤也少一些。有精力的条件下，可以两种方法都试试，比较一下两种方法哪个结果更适合自己。其实每种方法都有存在的意义，而且也有适用的范围，要在具体的项目中，结合背景知识分析哪种方法结果更好时才知道。
（五） 特征表和特征序列汇总
在质量筛选步骤完成之后，你将希望探索数据结果。可以使用以下两个命令进行此操作，这两个命令将创建数据的可视化摘要。特性表汇总命令（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 | 输出结果 table.qzv: 特征表统计。 rep-seqs.qzv: 代表序列统计，可点击序列跳转NCBI blast查看相近序列的信息 |
【特性表汇总命令（feature-table summarize）】
提供关于与每个样品和每个特性相关联的序列数量、这些分布的直方图以及一些相关的汇总统计数据的信息。
【特征表序列（feature-table tabulate-seqs）】
提供特征ID到序列的映射，并提供链接以针对NCBI nt数据库轻松BLAST每个序列。
【参考】
特征(feature):“特征"是指一个观测单位，比如一个OTU、一个序列变异(ASV)、一个基因、一个代谢物等。QIIME 2使用这种通用的术语（指”特征“）是由于QIIME 2支持许多类型的"特征”。
https://docs.qiime2.org/2020.11/tutorials/moving-pictures/
QIIME 2教程. 04人体微生物组分析Moving Pictures(2020.11)
https://blog.sciencenet.cn/blog-994715-1362530.html
上一篇：QIIME 2教程之物种组成分析
下一篇：QIIME 2教程之生成特征表和特征序列(案例二)
