|
Miseq 16S amplicon V3V4 PE300測(cè)序是目前菌群結(jié)構(gòu)譜研究最為常用的測(cè)序手段。本文將以此類測(cè)序的下機(jī)數(shù)據(jù)為例展示“如何從Miseq測(cè)序數(shù)據(jù)中快速提取出可以用來(lái)統(tǒng)計(jì)分析的菌屬相對(duì)豐度表”的工作流程。該豐度表是做菌群研究最為基本的數(shù)據(jù)，要想發(fā)文章還必須做大量的統(tǒng)計(jì)分析。在以后的文章中我們會(huì)繼續(xù)推出一些統(tǒng)計(jì)分析方法，敬請(qǐng)期待！ 軟件準(zhǔn)備 Linux環(huán)境： 1. FastQC 質(zhì)檢下機(jī)數(shù)據(jù) 軟件地址：https://www.bioinformatics./projects/fastqc/ 2. Cutadapt 切除測(cè)序引物 軟件地址：https://cutadapt./en/stable/ 3. QIIME2 序列拼接、質(zhì)控、比對(duì)、注釋 軟件版本：QIIME2 2018.4或QIIME2 2018.8 軟件地址：https:/// Windows環(huán)境： 1. Filezilla 下載Linux環(huán)境中的數(shù)據(jù)或上傳數(shù)據(jù)到Linux環(huán)境 軟件地址：https:/// 2. Excel 數(shù)據(jù)處理 3. QIIME2 view 查看QIIME2輸出的以.qzv為后綴的文件 網(wǎng)頁(yè)地址：https://view./ 數(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鍵為分隔符?？梢韵仍贓xcel中做表，然后保存為tsv文件。 格式如下：  4. Greengenes細(xì)菌16S全長(zhǎng)序列數(shù)據(jù)庫(kù) 下載地址：https://docs./2019.1/data-resources/ 下載得到gg_13_8_otus.tar.gz（最新版，大小為305M）后將其解壓得到99_otus.fasta（序列）和99_otu_taxonomy.txt（分類）兩個(gè)文件，文件獲取如下：  流程概覽  工作流程 1. 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é)果文件并打開(kāi)  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 由下圖可見(jiàn)99%以上的Reads都包含引物  2.3. 質(zhì)檢 Reads起始端質(zhì)量明顯提高，末端的低質(zhì)量堿基可利用下面的DADA2來(lái)處理  3. 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)頁(yè)中才能打開(kāi)。此處可以得到質(zhì)檢矢量圖，通過(guò)放大觀察可以清楚的判斷堿基質(zhì)量明顯下降的位置，從而輔助確定下一步中的reads1_cutpoint和reads2_cutpoint。  4. 用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，它將被用于下游分析  5. 引物特異性菌群比對(duì)數(shù)據(jù)庫(kù) 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ù)庫(kù)中的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 6. 細(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 7. 數(shù)據(jù)處理 7.1. 獲取屬水平菌豐度表 QIIME2 view網(wǎng)頁(yè)中打開(kāi)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ù)中快速提取出可以用來(lái)統(tǒng)計(jì)分析的菌屬相對(duì)豐度表”全部工作流程，如有問(wèn)題可以留言交流，以后會(huì)繼續(xù)推出“如何利用該菌屬相對(duì)豐度表進(jìn)行統(tǒng)計(jì)分析”的文章，如有興趣可以關(guān)注。 轉(zhuǎn)自生信草堂公眾號(hào)，已授權(quán) 數(shù)據(jù)請(qǐng)?jiān)诠娞?hào)獲取~ |
|
|
