Windows下16S rRNA基因擴(kuò)增子測(cè)序數(shù)據(jù)分析的簡(jiǎn)易流程

方梅梅，王禹煊，王明月，鄭和龍，張 璐，向沙沙，張國(guó)慶，李余動(dòng)

(浙江工商大學(xué) 食品與生物工程學(xué)院,杭州 310018)

地球上的各種環(huán)境中，如土壤，海洋或人體腸道等，都生存著成千上萬(wàn)的微生物，包含細(xì)菌、真菌、病毒和古菌等。微生物組(Microbiome)指某一特定生態(tài)環(huán)境中全部微生物及其遺傳物質(zhì)的總和[1-2]。人類(lèi)微生物組的組成與功能具有非常豐富的多樣性，并會(huì)隨著環(huán)境的變化而動(dòng)態(tài)變化。如在人類(lèi)口腔中定殖的微生物就有超過(guò)600種，口腔微生物菌群失衡會(huì)引起齲齒、口臭等口腔健康問(wèn)題[3]。隨著高通量測(cè)序技術(shù)的快速發(fā)展，宏基因組方法可以系統(tǒng)研究不同環(huán)境中微生物群落的組成與多樣性變化[4]。大量研究表明，微生物組結(jié)構(gòu)變化與人類(lèi)疾病，如糖尿病、癌癥、消化不良等密切相關(guān)[5]。微生物組學(xué)已成為最熱門(mén)的研究方向之一，目前有多項(xiàng)國(guó)際合作大的研究項(xiàng)目正在開(kāi)展，如人類(lèi)微生物組計(jì)劃(HMP)、地球微生物組計(jì)劃(EMP)和中國(guó)微生物組計(jì)劃等[2，6]。

16S rRNA基因全長(zhǎng)約1 542 bp，因在結(jié)構(gòu)和功能上都具有高度的保守性，存在于所有的細(xì)菌中，因此常被用于微生物分類(lèi)鑒定的標(biāo)志物。16S rRNA基因有9個(gè)高變區(qū)(V1～V9)中，V4區(qū)是最準(zhǔn)確且能識(shí)別到屬水平的可變區(qū)，而V3～V4區(qū)測(cè)序長(zhǎng)度(約465 bp)比V4區(qū)更長(zhǎng)(約290 bp)，相較單獨(dú)V4區(qū)測(cè)序更加準(zhǔn)確[7]。新一代高通量測(cè)序平臺(tái)Illumina MiSeq克服了傳統(tǒng)高通量測(cè)序方法的讀長(zhǎng)較短與誤差較大的缺陷，因而用于16S rRNA的V3～V4區(qū)測(cè)序可獲得更準(zhǔn)確的微生物物種分類(lèi)信息。

高通量測(cè)序項(xiàng)目一般會(huì)產(chǎn)生海量的測(cè)序原始數(shù)據(jù)，對(duì)于原始數(shù)據(jù)的過(guò)濾篩選，加工以及統(tǒng)計(jì)分析，對(duì)計(jì)算機(jī)的軟硬件要求較高，使沒(méi)有計(jì)算機(jī)背景的研究人員難以分析數(shù)據(jù)和解釋研究結(jié)果[8]。尤其是，常用的微生物組分析軟件(如QIIME、Mothur等)需要在Linux系統(tǒng)下運(yùn)行，而普通的生物研究人員完全掌握Linux系統(tǒng)的使用非常困難[9]。因此，如何使大多數(shù)科研人員能進(jìn)行生物信息分析是生物大數(shù)據(jù)時(shí)代的一個(gè)挑戰(zhàn)[10]。盡管現(xiàn)在已經(jīng)開(kāi)發(fā)有一些基于網(wǎng)頁(yè)的數(shù)據(jù)分析平臺(tái)，如MG-RAST[11]，SEED2[12]等，但是這些平臺(tái)往往有許多限制，如數(shù)據(jù)不能太大，網(wǎng)絡(luò)傳輸速度慢，分析任務(wù)等待時(shí)間長(zhǎng)等。尤其是網(wǎng)絡(luò)分析平臺(tái)不適合那些需要對(duì)樣本數(shù)據(jù)嚴(yán)格保密的項(xiàng)目。

微軟在2016年推出Windows的Linux子系統(tǒng)-Windows Subsystem for Linux (WSL)，原生支持在windows10上運(yùn)行Linux程序，包括大多數(shù)的Linux命令與工具，如bash, grep, ssh等。這樣大量在Linux平臺(tái)運(yùn)行的生物信息軟件都可以直接安裝在WSL下運(yùn)行。WSL使Linux軟件安裝與使用變得容易，只需要在Windows系統(tǒng)下，學(xué)習(xí)Linux系統(tǒng)的基本命令，不需要專(zhuān)門(mén)在電腦上或用虛擬機(jī)(如VirtualBox)安裝Linux系統(tǒng)。

為方便研究人員分析微生物組數(shù)據(jù)，我們開(kāi)發(fā)了一套在WSL-Ubuntu下分析樣品16S rRNA基因擴(kuò)增子高通量測(cè)序的簡(jiǎn)易流程。本流程基于Win10的Linux子系統(tǒng)，避免了安裝Linux的復(fù)雜過(guò)程，并簡(jiǎn)化了QIIME[13]、VSEARCH[14]及相關(guān)軟件的安裝過(guò)程，只要運(yùn)行一個(gè)安裝腳本就自動(dòng)安裝并配置運(yùn)行所需的環(huán)境變量，方便沒(méi)有計(jì)算機(jī)背景的研究人員進(jìn)行相關(guān)微生物組數(shù)據(jù)分析。

1 材料與方法

1.1 唾液樣本采集與DNA測(cè)序數(shù)據(jù)

本實(shí)驗(yàn)數(shù)據(jù)來(lái)源于一項(xiàng)針對(duì)在校大學(xué)生的口腔微生物研究項(xiàng)目。實(shí)驗(yàn)簡(jiǎn)要說(shuō)明如下：首先在實(shí)驗(yàn)的前期階段，我們進(jìn)行樣本信息的搜集，篩選和整理，后送實(shí)驗(yàn)室提取DNA和質(zhì)量檢驗(yàn)，其中，唾液樣本的采集參考陳瑜等人的方法[15]，DNA的提取使用QIAGEN DNA Mini KIT試劑盒提取，經(jīng)PCR擴(kuò)增16S DNA的V3～V4區(qū)后，進(jìn)行瓊脂糖凝膠電泳和核酸蛋白定量檢測(cè)來(lái)檢測(cè)提取的DNA的純度及濃度。對(duì)于檢測(cè)合格的樣本DNA送測(cè)序公司采用Illumina HiSeq測(cè)序平臺(tái)進(jìn)行16S rDNA雙端(Paired-End)測(cè)序。公司返回的測(cè)序結(jié)果為FASTQ格式文件的壓縮包，包含所有的樣本正反兩個(gè)方向的序列信息。

為方便學(xué)習(xí)本分析流程，本文只使用部分唾液樣本(10個(gè))的微生物組數(shù)據(jù)，測(cè)序數(shù)據(jù)可從本流程網(wǎng)站下載(http://www.ligene.cn/win16s/)。整個(gè)分析過(guò)程只需要約2 h。

1.2 電腦系統(tǒng)與軟件安裝

1.2.1 裝有最新版Windows 10系統(tǒng)的筆記本電腦

為了運(yùn)行本分析流程，首先需要一臺(tái)Windows 10系統(tǒng)的電腦，并且將該系統(tǒng)更新到最新版本(版本1 709及以上)。升級(jí)過(guò)程可以選擇電腦自動(dòng)升級(jí)，也可以在微軟網(wǎng)站下載升級(jí)程序，進(jìn)行手動(dòng)升級(jí)(Windows易升軟件: https://www.microsoft.com/en-us/software-download/windows10)。在系統(tǒng)升級(jí)完畢后，就可開(kāi)啟系統(tǒng)的Windows Subsystem for Linux(WSL)功能：在控制面板-“應(yīng)用”-“程序與功能”-“啟用或關(guān)閉Windows功能”-“Windows Subsystem for Linux”中開(kāi)啟“Windows Subsystem for Linux”,從而開(kāi)啟WSL的功能。

1.2.2 WSL+Ubuntu安裝

完成WSL的開(kāi)啟之后，在Windows10自帶的微軟應(yīng)用商店“Microsoft Store”中搜索并下載Ubuntu。當(dāng)Ubuntu app下載完畢之后，點(diǎn)擊“Launch”啟動(dòng)安裝。然后窗口會(huì)提示你創(chuàng)建一個(gè)Linux的賬戶(hù)，以后在Linux環(huán)境下進(jìn)行數(shù)據(jù)分析都是在這個(gè)用戶(hù)名下。這個(gè)用戶(hù)名和密碼與電腦微軟系統(tǒng)的賬戶(hù)和密碼沒(méi)有任何關(guān)系，因此不必設(shè)置為相同的。安裝完成之后，即可以通過(guò)開(kāi)始菜單的Ubuntu程序圖標(biāo)來(lái)運(yùn)行Ubuntu Linux系統(tǒng)。

1.2.3 Xming安裝

目前WSL不支持圖形界面，為了實(shí)現(xiàn)數(shù)據(jù)分析結(jié)果的可視化，需要安裝一個(gè)Xming軟件，通過(guò)Xming官網(wǎng)(http://sourceforge.net/projects/xming/)下載安裝包(最新版本6.9.0.31)，進(jìn)行安裝即可。

1.2.4 安裝win16s流程

此流程是基于VSEARCH (v2.7.1)與QIIME (v1.9.1)等開(kāi)源軟件，因此要在WSL-Ubuntu的命令行終端中輸入不同軟件的安裝指令，所有命令集見(jiàn)附件1(流程安裝)。為方便初學(xué)者安裝，我們提供一個(gè)bash腳本(win16s_install.sh)，下載此文件后，在ubuntu bash終端使用 “bash win16s_install.sh”命令運(yùn)行此腳本，即可一次安裝所有流程運(yùn)行依賴(lài)軟件，并能自動(dòng)配置合適的環(huán)境變量。

Win16s所有軟件與相關(guān)腳本可以從網(wǎng)址下載：http://www.ligene.cn/win16s/.

1.3 win16s分析流程

本分析流程的基本步驟如圖1所示。主要步驟為：樣本雙端測(cè)序運(yùn)用VSEARCH軟件進(jìn)行序列拼接， OTU聚類(lèi)，后通過(guò)比對(duì)RDP分類(lèi)數(shù)據(jù)庫(kù)注釋物種，最后用QIIME相關(guān)腳本進(jìn)行α/β多樣性分析(如PCoA等)。

圖1 16S rRNA基因擴(kuò)增子測(cè)序數(shù)據(jù)分析流程圖Fig.1 Flowchart of 16S rRNA gene amplicon data analysis pipeline

2 結(jié)果與討論

本結(jié)果是基于對(duì)10個(gè)唾液樣本的16S rRNA基因測(cè)序數(shù)據(jù)分析的教學(xué)過(guò)程。首先，每個(gè)學(xué)生按流程分析樣本數(shù)據(jù)，可以得到樣本的OTU聚類(lèi)與物種分類(lèi)信息；在查看樣品的細(xì)菌種類(lèi)后，統(tǒng)計(jì)每個(gè)樣品的差異菌種。然后，根據(jù)問(wèn)卷調(diào)查結(jié)果對(duì)樣本進(jìn)行分組：葷組與素組(見(jiàn)表1)，通過(guò)微生物組多樣性相關(guān)分析命令，統(tǒng)計(jì)不同組樣品間的菌群結(jié)構(gòu)差異。

2.1 樣本信息

為了分析不同飲食習(xí)慣大學(xué)生的唾液微生物的種群差異性，我們首選通過(guò)調(diào)查問(wèn)卷來(lái)篩選出合適的志愿者，然后再統(tǒng)一組織開(kāi)展取樣工作。取樣人群的基本信息，包括樣本編號(hào)，性別，年齡，BMI值，取樣方法，所屬組別(見(jiàn)表1)。為了避免長(zhǎng)時(shí)間的儲(chǔ)存對(duì)樣本產(chǎn)生的影響，對(duì)于取得的唾液樣本將直接用于后續(xù)研究(具體參見(jiàn)材料與方法1.1部分)。

表1 樣本信息表Table 1 Basic information of the samples

2.2 數(shù)據(jù)質(zhì)控與預(yù)處理

為保證數(shù)據(jù)質(zhì)量和避免錯(cuò)誤結(jié)果，16s rDNA擴(kuò)增子測(cè)序數(shù)據(jù)需要質(zhì)控及預(yù)處理后，再用于下游的分析。首先常用FastQC軟件對(duì)測(cè)序數(shù)據(jù)進(jìn)行質(zhì)量評(píng)估，由于一般公司的測(cè)序報(bào)告文件會(huì)附帶評(píng)估報(bào)告，質(zhì)量太差會(huì)重測(cè)，此步非本實(shí)驗(yàn)必須步聚。然后，利用常用的NGS質(zhì)控工具(如Trimmomatic[16])進(jìn)行預(yù)處理，如過(guò)濾低質(zhì)量reads(Q<20)，切除read末端低質(zhì)量堿基與未識(shí)別堿基N、切除測(cè)序接頭(Adapters)、去除長(zhǎng)度太短的序列等。

本實(shí)驗(yàn)數(shù)據(jù)為Illumina HiSeq2500平臺(tái)Paired-end測(cè)序(PE250)，并已去除引物與Barcodes后的原始序列，每個(gè)樣品有兩個(gè)正反向測(cè)序reads數(shù)據(jù)文件(注意文件名書(shū)寫(xiě)格式，如SAA_R1.fastq與SAA_R2.fastq)。正向Reads(R1)和反向Reads(R2)兩段序列有部分重疊可進(jìn)行拼接。一般默認(rèn)用fastq-join方法拼接，所得序列再過(guò)濾低質(zhì)量reads與切除末端低質(zhì)量堿基。由于16S rRNA的V3～V4區(qū)長(zhǎng)度為465 bp，本實(shí)驗(yàn)過(guò)濾長(zhǎng)度小于380 bp的序列片段。原始序列經(jīng)過(guò)濾后，最終每個(gè)樣本得到約15 000條有效reads(見(jiàn)表2)用于后續(xù)分析，好的質(zhì)量控制可以得到準(zhǔn)確的分類(lèi)結(jié)果。

2.3 OTU聚類(lèi)

OTU(Operational Taxonomic Units)是在系統(tǒng)發(fā)生學(xué)研究中，為了便于進(jìn)行分析，人為給某一個(gè)分類(lèi)單元(品系，種，屬，分組等)設(shè)置的同一標(biāo)志。要了解一個(gè)樣品測(cè)序結(jié)果中的菌種、菌屬等數(shù)目信息，就需要對(duì)序列進(jìn)行聚類(lèi)操作(Clustering)。通過(guò)聚類(lèi)將序列按照彼此的相似性97%以上歸為一個(gè)OTU。一般建議采用開(kāi)放參考(Open reference)聚類(lèi)方式，但是PICRUSt預(yù)測(cè)需要以默認(rèn)Greengene數(shù)據(jù)庫(kù)的封閉參考(Closed reference)聚類(lèi)分析結(jié)果。

本實(shí)驗(yàn)對(duì)樣品數(shù)據(jù)進(jìn)行質(zhì)量控制、序列拼接及去除嵌合體后，選擇序列相似度值97%進(jìn)行OTU聚類(lèi)。細(xì)菌的16S rRNA分類(lèi)一般使用97%相似度，而真核生物的18S rRNA分類(lèi)一般使用96%相似度，使用不同的算法或相似度閥值可能得到完全不同的OTU數(shù)[17]。得到OTU后，采用RDP算法與Greengenes rDNA數(shù)據(jù)庫(kù)比對(duì)，注釋各OTU的分類(lèi)單元。大多數(shù)OTU序列能注釋到屬水平，一部分OTU可以鑒定到種水平。

本實(shí)驗(yàn)共得到了212個(gè)OTU，在門(mén)和屬兩個(gè)水平上對(duì)同類(lèi)的OTU進(jìn)行分類(lèi)鑒定，歸屬于12個(gè)門(mén)，94個(gè)屬 (見(jiàn)表2)。通過(guò)OTU 聚類(lèi)及分類(lèi)鑒定得到的結(jié)果，可以作為后續(xù)菌群統(tǒng)計(jì)分析的基礎(chǔ)。

2.4 細(xì)菌種群的豐富度

不同人體的口腔微生物種群存在著一定的差異，通過(guò)分析可以得到不同分類(lèi)水平的細(xì)菌種群結(jié)構(gòu)及其豐富度，可以清晰地得到不同樣本在同一分類(lèi)水平的差異。分析流程在進(jìn)行可視化操作后，在taxa_summary文件夾中生成子文件夾taxa_summary_plots,打開(kāi)其中的文件bar_charts.html,即可觀察到各樣品自門(mén)水平(L2)到屬水平(L6)的分類(lèi)及其組成比例。圖2是在門(mén)和屬水平對(duì)菌群的差異性的結(jié)果。具體菌群分類(lèi)與豐度信息可以在多樣性分析輸出目錄(“taxa_summary/”)，查看在不同分類(lèi)水平下的細(xì)菌種群的豐富度。由結(jié)果可知，部分菌種在所有樣本中都存在；而少數(shù)菌種則只在部分樣本中存在。

通過(guò)對(duì)不同樣本分組的種群豐富度分析，本實(shí)驗(yàn)共檢出優(yōu)勢(shì)菌涉及12個(gè)門(mén)(見(jiàn)圖2(a))，其中Firmicutes、Proteobacteria、Bacteroidetes、Actinobacteria不存在顯著性差異，而Fusobacteria存在顯著性組間差異(wilcox.test,P< 0.05)。本實(shí)驗(yàn)共檢出優(yōu)勢(shì)菌涉及94個(gè)屬(見(jiàn)圖2(b))， 葷食組中Streptococcus比素食組顯著增加，而Fusobacterium比素食組顯著減少(wilcox.test,P< 0.05)。

表2 唾液細(xì)菌測(cè)序的基本信息Table 2 Basic sequencing information of bacteria in saliva

圖2 實(shí)驗(yàn)樣品菌群分類(lèi)與相對(duì)豐度Fig.2 Taxonomy and relative abundance of the bacterial populations across samples

2.5 其它統(tǒng)計(jì)分析

最后得到的biom文件(otus_table_tax.biom)可用于后續(xù)統(tǒng)計(jì)分析。通?？梢杂?jì)算兩類(lèi)多樣性指數(shù)：Alpha多樣性指數(shù)(α-diversity)與Beta多樣性指數(shù)(-diversity)。α多樣性指數(shù)度量的是群落內(nèi)部的物種多樣性，多樣性指數(shù)度量的是群落間的物種多樣性。按照生態(tài)學(xué)的觀點(diǎn)，多樣性高的生態(tài)群落，抵抗力和穩(wěn)定性都可能更強(qiáng)。這里以QIIME流程中的core_diversity_analyses.py腳本分析，并通過(guò)最小測(cè)序量樣本的reads(Counts)數(shù)作為-e參數(shù)(Rarefaction抽樣數(shù))。core_diversity_analyses.py分析運(yùn)行一系列標(biāo)準(zhǔn)的alpha與beta多樣性分析，可以查看分析輸出文件夾(core_diversity/)下的文件查看相關(guān)結(jié)果。最后，可以將OTU表的biom格式文件轉(zhuǎn)換成文本文件，用于其它統(tǒng)計(jì)分析軟件的分析，如STAMP[18]或R軟件包phyloseq, vena等。

2.5.1 樣本菌群Alpha多樣性分析

本實(shí)驗(yàn)統(tǒng)計(jì)了樣本α多樣性分析中的不同多樣性指數(shù)。圖3(a)是系統(tǒng)發(fā)育多樣性指數(shù)(PD)，兩組的多樣性非常相似，雖然兩條曲線(xiàn)(PD值)區(qū)分明顯，但誤差棒有重疊，即沒(méi)有顯著統(tǒng)計(jì)差異(P> 0.05)。另外，葷食組和素食組Simpson多樣性指數(shù)分別為(0.94±0.04)與(0.95±0.02)，Shannon多樣性指數(shù)分別為(5.54±0.63)與(5.78±0.30)，可見(jiàn)兩組樣本有著相似的多樣性指數(shù)。

2.5.2 樣本菌群Beta多樣性分析

本實(shí)驗(yàn)利用QIIME的腳本對(duì)得到的OTU表做差異距離矩陣，后用主坐標(biāo)分析(Principal coordinates analysis, PCoA)將距離矩陣可視化。PCoA結(jié)果如圖3(b)所示,其中一個(gè)點(diǎn)代表一個(gè)樣本，相同的顏色表示是同一組的樣本，兩個(gè)點(diǎn)之間的距離越近，說(shuō)明兩個(gè)樣品的微生物群落差異越小。由PCoA分析結(jié)果(見(jiàn)圖3(b))可知兩組樣本均比較分散，主成分區(qū)分度不是特別明顯，但在各自象限均占主要優(yōu)勢(shì)，說(shuō)明兩組樣本的主成分有一定差異，但沒(méi)有顯著的差異。這可能與本研究的采樣對(duì)象樣本數(shù)量太少有關(guān)，還需要擴(kuò)大取樣人群范圍進(jìn)一步研究。

圖3 樣本分組的菌群差異Fig.3 Bacterial population difference between sample groups

生活習(xí)慣和飲食習(xí)慣的不同，形成具有差異的口腔微生物群落結(jié)構(gòu)，存在不同的優(yōu)勢(shì)菌和劣勢(shì)菌。它們影響著口腔的環(huán)境，對(duì)人體造成不同的影響。通過(guò)分類(lèi)統(tǒng)計(jì)分析研究不同組別的菌群組成差異，進(jìn)一步研究菌群差異與樣本之間存在的關(guān)系，是微生物組研究的意義和目的。

3 結(jié) 論

通過(guò)整合開(kāi)源軟件VSEARCH與QIIME等分析流程，我們開(kāi)發(fā)了在Windows系統(tǒng)下分析16S rRNA基因擴(kuò)增子測(cè)序數(shù)據(jù)的簡(jiǎn)易流程。本流程不僅適用于16S rRNA數(shù)據(jù)分析，還可以用于18S rRNA或ITS等標(biāo)記基因的分析。以后隨著研究的深入，我們將提供更多元化的分析腳本，不斷更新現(xiàn)有分析流程，并同時(shí)在網(wǎng)上公布，供初學(xué)者學(xué)習(xí)微生物組分析方法[19]。

WSL是微軟在Windows 10中推出的Linux子系統(tǒng)功能,允許Linux ELF64二進(jìn)制文件運(yùn)行在Windows系統(tǒng)上。隨著微軟對(duì)WSL的改進(jìn)，越來(lái)越多的生物信息軟件將支持在WSL運(yùn)行，從而方便廣大生物學(xué)背景的人員在Windows系統(tǒng)下進(jìn)行生物信息數(shù)據(jù)分析[20]。然而WSL也有一定局限，如它只支持64位程序，一些老的32位程序?qū)⒉荒茉赪SL下運(yùn)行，如QIIME默認(rèn)的uclust程序不能在WSL下編譯成功。好在目前大多數(shù)32位程序都有相應(yīng)的64位替代程序，如uclust、USEARCH[20]可以用免費(fèi)的64位程序VSEARCH代替。另外，在WSL下，程序?qū)τ脖P(pán)的讀寫(xiě)效率還比較低，程序運(yùn)行的效率還是比在原生Linux系統(tǒng)下運(yùn)行的速度慢。但WSL適合作為數(shù)據(jù)分析入門(mén)的學(xué)習(xí)環(huán)境，當(dāng)研究人員習(xí)慣了Linux命令操作，可以很方便遷移到原生Linux計(jì)算機(jī)上進(jìn)行實(shí)際項(xiàng)目數(shù)據(jù)的分析研究工作。