导读
Miseq 16S amplicon V3V4 PE300测序是目前菌群结构谱研究最为常用的测序手段。本文将以此类测序的下机数据为例展示“如何从Miseq测序数据中快速提取出可以用来统计分析的菌属相对丰度表”的工作流程。该丰度表是做菌群研究最为基本的数据，要想发文章还必须做大量的统计分析。以后推出更多方法。
一、软件准备
Linux环境：
1. FastQC
用途：质检下机数据
软件地址：https://www.bioinformatics.babraham.ac.uk/projects/fastqc/
2. Cutadapt
用途：切除测序引物
软件地址：https://cutadapt.readthedocs.io/en/stable/
3. QIIME2
序列拼接、质控、比对、注释
软件地址：https://qiime2.org/
Windows环境：
1. Filezilla
用途：下载Linux环境中的数据或上传数据到Linux环境
软件地址：https://filezilla-project.org/
2. Excel
用途：数据处理
3. QIIME2 view
用途查看QIIME2输出的以.qzv为后缀的文件
网页地址：https://view.qiime2.org/
二、数据准备
1. Miseq 16S amplicon V3V4测序下机数据
R1.fastq，R2.fastq
2. 测序引物
p1 -> CCTACGGGNGGCWGCAG
p2 -> GACTACHVGGGTATCTAATCC
3. 表型文件metadata.txt
准备存放样本信息的表型文件，以tab键为分隔符。可以先在Excel中做表，然后保存为tsv文件。
格式如下：
4. Greengenes细菌16S全长序列数据库
下载地址：https://docs.qiime2.org/2019.1/data-resources/
下载得到gg_13_8_otus.tar.gz（最新版，大小为305M）后将其解压得到99_otus.fasta（序列）和99_otu_taxonomy.txt（分类）两个文件，文件获取如下：
三、流程概览
四、工作流程
- FastQC质检
1.1. 合并R1、R2
cat *R1.fastq > merge.R1.fastq
cat *R2.fastq > merge.R2.fastq
1.2. FastQC质检
可以使用-t：指定线程数和-o：指定输出位置
fastq -t 8 merge.R1.fastq
fastq -t 8 merge.R2.fastq
1.3. 使用Filezilla下载结果文件并打开，如下：
- Cutadapt切引物
2.1. 检查引物的位置和序列
位置：5’端
序列：p1 -> CCTACGGGNGGCWGCAG; p2 -> GACTACHVGGGTATCTAATCC
2.2. 切引物
cutadapt -g CCTACGGGNGGCWGCAG -G GACTACHVGGGTATCTAATCC -o *R1.fastq -p *R2.fastq *r1.fastq *r2.fastq --core=2
由下图可见99%以上的Reads都包含引物
2.3.质检
Reads起始端质量明显提高，末端的低质量碱基可利用下面的DADA2来处理
- QIIME2数据导入
3.1. 制作manifest.txt列表文件，存放每一个样本的信息，格式如下：
sample-id,absolute-filepath,direction
sample-1,/filepath/sample1_r1.fastq,forward
sample-1,/filepath/sample1_r2.fastq,reverse
sample-2,/filepath/sample2_r1.fastq,forward
sample-2,/filepath/sample2_r2.fastq,reverse
注意不能有空格、换行符、制表符等
3.2. 格式转换
qiime tools import \
--type 'SampleData[PairedEndSequencesWithQuality]' \
--input-path manifest.txt \
--output-path manifest.qza \
--source-format PairedEndFastqManifestPhred33
3.3. 可视化检查
qiime demux summarize \
--i-data manifest.qza \
--o-visualization manifest.qzv
manifest.qzv文件需要从Linux中下载后再拖拽到qiime2 view网页中才能打开。此处可以得到质检矢量图，通过放大观察可以清楚的判断碱基质量明显下降的位置，从而辅助确定下一步中的reads1_cutpoint和reads2_cutpoint。
- 用DADA2进行切割、去嵌合体、拼接等
4.1. 使用10个线程运行DADA2
为保证碱基质量这里再次要对Reads进行切割。Reads起始端质量很高时N1 N2设为0即可；观察manifest.qzv确定reads1_cutpoint和reads2_cutpoint，这里我将其分别设为275和250。
qiime dada2 denoise-paired \
--i-demultiplexed-seqs manifest.qza \
--p-trim-left-f N1 \
--p-trim-left-r N2 \
--p-trunc-len-f reads1_cutpoint \
--p-trunc-len-r reads2_cutpoint \
--o-table table.qza \
--o-representative-sequences rep-seqs.qza \
--o-denoising-stats denoising-stats.qza \
--p-n-threads 10
此步骤会生成三个新文件：
denoising-stats.qza是质检统计，如下表；
table.qza是细菌特征丰度表；
rep-seqs.qza是细菌特征代表性序列
4.2. DADA2统计结果可视化
qiime metadata tabulate \
--m-input-file denoising-stats.qza \
--o-visualization denoising-stats.qzv
最后一列是Clean data，它将被用于下游分析。
- 引物特异性菌群比对数据库
5.1. 转换格式
将99_otus.fasta（序列）和99_otu_taxonomy.txt（分类）两个文件的格式转换QIIME2能识别和利用的格式
qiime tools import \
--type 'FeatureData[Sequence]' \
--input-path 99_otus.fasta \
--output-path 99_otus.qza
qiime tools import \
--type 'FeatureData[Taxonomy]' \
--input-format HeaderlessTSVTaxonomyFormat \
--input-path 99_otu_taxonomy.txt \
--output-path 99-taxonomy.qza
5.2. 抽提V3V4模板序列
用测序引物序列从Greengenes数据库中的16S全长序列99_otus.qza中抽提出引物特异性的细菌参考序列，就会得到本研究特异性的参考序列
qiime feature-classifier extract-reads \
--i-sequences 99_otus.qza \
--p-f-primer CCTACGGGNGGCWGCAG \
--p-r-primer GACTACHVGGGTATCTAATCC \
--o-reads 99-v3v4-seqs.qza
5.3. 训练V3V4分类器
qiime feature-classifier fit-classifier-naive-bayes \
--i-reference-reads 99-v3v4-seqs.qza \
--i-reference-taxonomy 99-taxonomy.qza \
--o-classifier gg-13-8-99-v3v4-classifier.qza
- 细菌分类学注释
6.1. 比对
把DADA2分析得到的细菌特征代表性序列rep-seqs.qza和训练好的分类器gg-13-8-99-v3v4-classifier.qza进行比对，获得具体的细菌分类信息taxonomy.qza
qiime feature-classifier classify-sklearn \
--i-classifier gg-13-8-99-v3v4-classifier.qza \
--i-reads rep-seqs.qza \
--o-classification taxonomy.qza
6.2. 丰度统计
将细菌特征丰度表table.qza和细菌分类信息taxonomy.qza进行整合获得完整的细菌分类丰度表，包含界、门、纲、目、科、属、种多水平的细菌丰度信息
qiime taxa barplot \
--i-table table.qza \
--i-taxonomy taxonomy.qza \
--m-metadata-file metadata.tsv \
--o-visualization taxa-bar-plots.qzv
- 数据处理
7.1. 获取属水平菌丰度表
QIIME2 view网页中打开taxa-bar-plots.qzv，下载level6的CSV文件，如下：
7.2. 标准化菌属丰度表
把CSV文件导入到Excel中进行标准化，即每个菌属的原始丰度除以该菌所在样本的总菌属丰度得到标准相对菌属相对丰度
处理方法如下：
以上就是我个人总结的“从Miseq测序数据中快速提取出可以用来统计分析的菌属相对丰度表”全部工作流程，如有问题可以留言交流，以后会继续推出“如何利用该菌属相对丰度表进行统计分析”的文章，如有兴趣可以关注。
相关阅读：
16S测序分析（一）菌属丰度表获取
16S测序分析（二）菌群多样性分析
16S测序分析（三）用LEfSe寻找组间差异细菌
16S测序分析（四）用MaAsLin寻找组间差异细菌
16S测序分析（五）用RandomForest寻找关键细菌
16S测序分析（六）用PICRUSt预测菌群KEGG代谢通路
