Miseq 16S amplicon V3V4 PE300測序是目前菌群結構譜研究最為常用的測序手段。本文將以此類測序的下機數據為例展示“如何從Miseq測序數據中快速提取出可以用來統(tǒng)計分析的菌屬相對豐度表”的工作流程。該豐度表是做菌群研究最為基本的數據，要想發(fā)文章還必須做大量的統(tǒng)計分析。在以后的文章中我們會繼續(xù)推出一些統(tǒng)計分析方法，敬請期待！

軟件準備

Linux環(huán)境：

1.       FastQC

質檢下機數據

軟件地址：https://www.bioinformatics./projects/fastqc/

2.       Cutadapt

切除測序引物

軟件地址：https://cutadapt./en/stable/

3.       QIIME2

序列拼接、質控、比對、注釋

軟件版本：QIIME2 2018.4或QIIME2 2018.8

軟件地址：https:///

Windows環(huán)境：

1.       Filezilla

下載Linux環(huán)境中的數據或上傳數據到Linux環(huán)境

軟件地址：https:///

2.       Excel

數據處理

3.  QIIME2 view

查看QIIME2輸出的以.qzv為后綴的文件

網頁地址：https://view./

數據準備

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./2019.1/data-resources/

下載得到gg_13_8_otus.tar.gz（最新版，大小為305M）后將其解壓得到99_otus.fasta（序列）和99_otu_taxonomy.txt（分類）兩個文件，文件獲取如下：

流程概覽

工作流程

1. 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下載結果文件并打開

2. 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來處理

3. 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。

4. 用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是質檢統(tǒng)計，如下表；

table.qza是細菌特征豐度表；

rep-seqs.qza是細菌特征代表性序列

4.2. DADA2統(tǒng)計結果可視化

qiime metadata tabulate \

  --m-input-file denoising-stats.qza \

  --o-visualization denoising-stats.qzv

最后一列是Clean data，它將被用于下游分析

5. 引物特異性菌群比對數據庫

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. 細菌分類學注釋

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. 豐度統(tǒng)計

將細菌特征豐度表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. 數據處理

7.1. 獲取屬水平菌豐度表

QIIME2 view網頁中打開taxa-bar-plots.qzv，下載level6的CSV文件，如下：

7.2. 標準化菌屬豐度表

把CSV文件導入到Excel中進行標準化，即每個菌屬的原始豐度除以該菌所在樣本的總菌屬豐度得到標準相對菌屬相對豐度

處理方法如下：

以上就是我個人總結的“從Miseq測序數據中快速提取出可以用來統(tǒng)計分析的菌屬相對豐度表”全部工作流程，如有問題可以留言交流，以后會繼續(xù)推出“如何利用該菌屬相對豐度表進行統(tǒng)計分析”的文章，如有興趣可以關注。

轉自生信草堂公眾號，已授權

數據請在公眾號獲取~