導(dǎo)讀
Miseq 16S amplicon V3V4 PE300測(cè)序是目前菌群結(jié)構(gòu)譜研究最為常用的測(cè)序手段疮方。本文將以此類測(cè)序的下機(jī)數(shù)據(jù)為例展示“如何從Miseq測(cè)序數(shù)據(jù)中快速提取出可以用來統(tǒng)計(jì)分析的菌屬相對(duì)豐度表”的工作流程。該豐度表是做菌群研究最為基本的數(shù)據(jù)，要想發(fā)文章還必須做大量的統(tǒng)計(jì)分析。以后推出更多方法岳服。
一、軟件準(zhǔn)備
Linux環(huán)境：
1. FastQC
用途：質(zhì)檢下機(jī)數(shù)據(jù)
軟件地址：https://www.bioinformatics.babraham.ac.uk/projects/fastqc/
2. Cutadapt
用途：切除測(cè)序引物
軟件地址：https://cutadapt.readthedocs.io/en/stable/
3. QIIME2
序列拼接惫恼、質(zhì)控鲸沮、比對(duì)畅形、注釋
軟件地址：https://qiime2.org/
Windows環(huán)境：
1. Filezilla
用途：下載Linux環(huán)境中的數(shù)據(jù)或上傳數(shù)據(jù)到Linux環(huán)境
軟件地址：https://filezilla-project.org/
2. Excel
用途：數(shù)據(jù)處理
3. QIIME2 view
用途查看QIIME2輸出的以.qzv為后綴的文件
網(wǎng)頁地址：https://view.qiime2.org/
二、數(shù)據(jù)準(zhǔn)備
1. Miseq 16S amplicon V3V4測(cè)序下機(jī)數(shù)據(jù)
R1.fastq诉探，R2.fastq
2. 測(cè)序引物
p1 -> CCTACGGGNGGCWGCAG
p2 -> GACTACHVGGGTATCTAATCC
3. 表型文件metadata.txt
準(zhǔn)備存放樣本信息的表型文件日熬，以tab鍵為分隔符∩隹瑁可以先在Excel中做表竖席，然后保存為tsv文件。
格式如下：
4. Greengenes細(xì)菌16S全長(zhǎng)序列數(shù)據(jù)庫
下載地址：https://docs.qiime2.org/2019.1/data-resources/
下載得到gg_13_8_otus.tar.gz（最新版敬肚，大小為305M）后將其解壓得到99_otus.fasta（序列）和99_otu_taxonomy.txt（分類）兩個(gè)文件毕荐，文件獲取如下：
三、流程概覽
四艳馒、工作流程
- FastQC質(zhì)檢
1.1. 合并R1憎亚、R2
cat *R1.fastq > merge.R1.fastq
cat *R2.fastq > merge.R2.fastq
1.2. FastQC質(zhì)檢
可以使用-t：指定線程數(shù)和-o：指定輸出位置
fastq -t 8 merge.R1.fastq
fastq -t 8 merge.R2.fastq
1.3. 使用Filezilla下載結(jié)果文件并打開，如下：
- 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.質(zhì)檢
Reads起始端質(zhì)量明顯提高弄慰，末端的低質(zhì)量堿基可利用下面的DADA2來處理
- QIIME2數(shù)據(jù)導(dǎo)入
3.1. 制作manifest.txt列表文件第美，存放每一個(gè)樣本的信息，格式如下：
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. 格式轉(zhuǎn)換
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網(wǎng)頁中才能打開。此處可以得到質(zhì)檢矢量圖慌闭，通過放大觀察可以清楚的判斷堿基質(zhì)量明顯下降的位置别威，從而輔助確定下一步中的reads1_cutpoint和reads2_cutpoint。
- 用DADA2進(jìn)行切割驴剔、去嵌合體省古、拼接等
4.1. 使用10個(gè)線程運(yùn)行DADA2
為保證堿基質(zhì)量這里再次要對(duì)Reads進(jìn)行切割。Reads起始端質(zhì)量很高時(shí)N1 N2設(shè)為0即可丧失；觀察manifest.qzv確定reads1_cutpoint和reads2_cutpoint豺妓，這里我將其分別設(shè)為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
此步驟會(huì)生成三個(gè)新文件：
denoising-stats.qza是質(zhì)檢統(tǒng)計(jì)利花，如下表科侈；
table.qza是細(xì)菌特征豐度表；
rep-seqs.qza是細(xì)菌特征代表性序列
4.2. DADA2統(tǒng)計(jì)結(jié)果可視化
qiime metadata tabulate \
--m-input-file denoising-stats.qza \
--o-visualization denoising-stats.qzv
最后一列是Clean data炒事，它將被用于下游分析。
- 引物特異性菌群比對(duì)數(shù)據(jù)庫
5.1. 轉(zhuǎn)換格式
將99_otus.fasta（序列）和99_otu_taxonomy.txt（分類）兩個(gè)文件的格式轉(zhuǎn)換QIIME2能識(shí)別和利用的格式
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模板序列
用測(cè)序引物序列從Greengenes數(shù)據(jù)庫中的16S全長(zhǎng)序列99_otus.qza中抽提出引物特異性的細(xì)菌參考序列蔫慧，就會(huì)得到本研究特異性的參考序列
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. 訓(xùn)練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
- 細(xì)菌分類學(xué)注釋
6.1. 比對(duì)
把DADA2分析得到的細(xì)菌特征代表性序列rep-seqs.qza和訓(xùn)練好的分類器gg-13-8-99-v3v4-classifier.qza進(jìn)行比對(duì)挠乳，獲得具體的細(xì)菌分類信息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)計(jì)
將細(xì)菌特征豐度表table.qza和細(xì)菌分類信息taxonomy.qza進(jìn)行整合獲得完整的細(xì)菌分類豐度表，包含界、門睡扬、綱盟蚣、目、科卖怜、屬屎开、種多水平的細(xì)菌豐度信息
qiime taxa barplot \
--i-table table.qza \
--i-taxonomy taxonomy.qza \
--m-metadata-file metadata.tsv \
--o-visualization taxa-bar-plots.qzv
- 數(shù)據(jù)處理
7.1. 獲取屬水平菌豐度表
QIIME2 view網(wǎng)頁中打開taxa-bar-plots.qzv，下載level6的CSV文件马靠，如下：
7.2. 標(biāo)準(zhǔn)化菌屬豐度表
把CSV文件導(dǎo)入到Excel中進(jìn)行標(biāo)準(zhǔn)化奄抽，即每個(gè)菌屬的原始豐度除以該菌所在樣本的總菌屬豐度得到標(biāo)準(zhǔn)相對(duì)菌屬相對(duì)豐度
處理方法如下：
以上就是我個(gè)人總結(jié)的“從Miseq測(cè)序數(shù)據(jù)中快速提取出可以用來統(tǒng)計(jì)分析的菌屬相對(duì)豐度表”全部工作流程，如有問題可以留言交流甩鳄，以后會(huì)繼續(xù)推出“如何利用該菌屬相對(duì)豐度表進(jìn)行統(tǒng)計(jì)分析”的文章逞度，如有興趣可以關(guān)注。
相關(guān)閱讀：
16S測(cè)序分析（一）菌屬豐度表獲取
16S測(cè)序分析（二）菌群多樣性分析
16S測(cè)序分析（三）用LEfSe尋找組間差異細(xì)菌
16S測(cè)序分析（四）用MaAsLin尋找組間差異細(xì)菌
16S測(cè)序分析（五）用RandomForest尋找關(guān)鍵細(xì)菌
16S測(cè)序分析（六）用PICRUSt預(yù)測(cè)菌群KEGG代謝通路
