一種基于宏基因組模擬數(shù)據(jù)的生物標(biāo)志物篩選方法

王曉君++滕琳

摘要：鑒于生物圈中微生物資源的巨大開(kāi)發(fā)潛力以及測(cè)序技術(shù)不斷發(fā)展，宏基因組學(xué)研究的不斷深入，微生物群落已經(jīng)被看作一個(gè)整體來(lái)進(jìn)行分析并且已經(jīng)得到廣泛應(yīng)用。然而由于微生物的多樣性以及微生物菌群的復(fù)雜性，使得精確確定和定量宏基因組數(shù)據(jù)中的分類單元成為宏基因組數(shù)據(jù)分析的難點(diǎn)。已有的宏基因組數(shù)據(jù)標(biāo)記分析工具無(wú)法解決微生物群落預(yù)測(cè)結(jié)果重現(xiàn)的穩(wěn)健性、準(zhǔn)確性以及處理非冗余標(biāo)記物方面遇到的問(wèn)題。筆者提出了一個(gè)新的基于宏基因組自助抽樣（metagenomic bootstrap）的生物標(biāo)志物選擇方法，它結(jié)合了mRMR（minimal redundancy maximal relevance）和自助抽樣方法（bootstrapping），可以更加穩(wěn)健、準(zhǔn)確而有效地通過(guò)對(duì)宏基因組數(shù)據(jù)的挖掘?qū)崿F(xiàn)非冗余標(biāo)記物的篩選?；谀M數(shù)據(jù)集，通過(guò)其與2種自上而下的方法（Metastats、LEfSe）以及自下而上的方法（Wilcoxon秩和檢驗(yàn)）進(jìn)行對(duì)比，表明本方法可以在較高準(zhǔn)確率的基礎(chǔ)上更加穩(wěn)健地選擇更多的非冗余生物標(biāo)志物。

關(guān)鍵詞：宏基因組；生物標(biāo)志物；mRMR；自助抽樣法

中圖分類號(hào)： Q789文獻(xiàn)標(biāo)志碼： A文章編號(hào)：1002-1302（2016）05-0056-04

微生物一直被人們視為巨大的生物資源，尤其是其龐大的基因組數(shù)據(jù)包含有大量不為人知的新功能基因，將對(duì)人類的生產(chǎn)、生活做出卓越貢獻(xiàn)[1]。然而，微生物資源中九成以上的微生物是不可培養(yǎng)的，也就意味著在新基因探索的道路上，人類面臨著不小的困難。新一代測(cè)序技術(shù)的出現(xiàn)將幫助人們揭示不可（或難）培養(yǎng)微生物的基因組信息，從而發(fā)現(xiàn)新的微生物或新的功能基因。隨著微生物基因組數(shù)據(jù)庫(kù)的不斷壯大，人們普遍意識(shí)到宏基因組數(shù)據(jù)分析的難點(diǎn)，宏基因組數(shù)據(jù)中生物標(biāo)志物的鑒定以及應(yīng)用非常重要。但宏基因組數(shù)據(jù)分析并不簡(jiǎn)單，研究顯示，微生物群落展現(xiàn)出了非同一般的主體間可變性，更不可思議的是，此可變性竟然出現(xiàn)在人類和環(huán)境菌群中[2-3]。目前，人們已知的宏基因組生物標(biāo)志物的鑒定方法有2種：一種是自下而上的方法，主要包括Wilcoxon秩和檢驗(yàn)[4]，測(cè)試每個(gè)分類單元，選擇群體間具有差異的元素作為標(biāo)志物；另外一種是自上而下的方法，主要包括Metastats、LEfSe。雖然這2種方法都可以用來(lái)統(tǒng)計(jì)評(píng)估宏基因組數(shù)據(jù)的差異，對(duì)生物標(biāo)志物進(jìn)行鑒定，但這些方法很難解決數(shù)據(jù)分析結(jié)果重現(xiàn)的穩(wěn)健性、冗余性等問(wèn)題。筆者提出一個(gè)自上而下的結(jié)合mRMR[5]和自助抽樣法從微生物宏基因組樣本中篩選生物標(biāo)志物的方法，此方法首先分析微生物群落的整體分布，然后進(jìn)行生物標(biāo)志物篩選，不同于傳統(tǒng)生物標(biāo)志物篩選的是，它結(jié)合了mRMR，能更為有效地避免了生物冗余標(biāo)志物這一難題。

1材料與方法

1.1模擬數(shù)據(jù)集的產(chǎn)生

S1模擬數(shù)據(jù)集：根據(jù)文獻(xiàn)，微生物群落的分類分布都遵循正態(tài)分布，故而基于正態(tài)分布，產(chǎn)生模擬數(shù)據(jù)集S1（S1未列出，僅說(shuō)明特性，其具體的結(jié)構(gòu)類似于下面即將產(chǎn)生的數(shù)據(jù)集S3，只是在生成數(shù)據(jù)時(shí)產(chǎn)生的是正態(tài)分布的數(shù)據(jù)，不同分類之間的差異指的是均值差異。）。S1中共有1 000個(gè)變量和120個(gè)樣本，包含2個(gè)分類（每個(gè)分類包含3個(gè)亞類，每個(gè)亞類包含20個(gè)樣本）。對(duì)于每一個(gè)樣本來(lái)說(shuō)，都包含10個(gè)真標(biāo)志物組（10個(gè)變量/組）和1個(gè)假標(biāo)志物組（900個(gè)變量/假標(biāo)志物）。數(shù)據(jù)集S1的特性是真標(biāo)志物中的2個(gè)分類組均值差異較大，在每個(gè)分類內(nèi)部，亞類之間的差異很?。ㄔ诿總€(gè)標(biāo)志物組內(nèi)，雖然生成數(shù)據(jù)時(shí)沒(méi)有差異，但是由于隨機(jī)函數(shù)的緣故，差異在所難免）。S2模擬數(shù)據(jù)集：筆者分析以前本實(shí)驗(yàn)室口腔微生物宏基因組數(shù)據(jù)[6]發(fā)現(xiàn)，微生物群落的宏基因組數(shù)據(jù)的分布不單是正態(tài)分布這么簡(jiǎn)單，往往會(huì)有10%的變量符合正態(tài)分布和伽瑪分布2種混合分布模式，因此基于正態(tài)和伽瑪混合分布產(chǎn)生模擬數(shù)據(jù)集S2（表1）。數(shù)據(jù)集S2有2個(gè)重要特性：第一，對(duì)于真標(biāo)志物，2個(gè)分類組參數(shù)shape（伽瑪分布中的1個(gè)重要參數(shù)）或者均值差異較大，每個(gè)分類內(nèi)部亞類之間的差異較小；第二，對(duì)于假標(biāo)志物，它們?cè)诜诸?、亞類之間均值沒(méi)有差異（每個(gè)標(biāo)志物組內(nèi)隨機(jī)差異如S1所述）。此外，處在相同標(biāo)志物組內(nèi)的變量被認(rèn)為是冗余的變量。S3模擬數(shù)據(jù)集：根據(jù)之前口腔樣本數(shù)據(jù)發(fā)現(xiàn)，超過(guò)40%的變量?jī)H符合伽瑪分布，因此基于伽瑪分布產(chǎn)生模擬數(shù)據(jù)集S3（表2）。數(shù)據(jù)集S3區(qū)別于S2的特性在于真標(biāo)志物中S3數(shù)據(jù)集的2個(gè)分類組在參數(shù)shape上差異較大，在每個(gè)分類內(nèi)，亞類之間的差異較小。

在真標(biāo)志物中，一個(gè)小方格是一個(gè)25（樣本）×10（變量）的矩陣。矩陣每一列的值都是由正態(tài)分布函數(shù)或者伽馬分布函數(shù)（利用R語(yǔ)言中rnorm或者rgamma函數(shù)實(shí)現(xiàn)）產(chǎn)生的。表格中填充淺灰色的格子表示由伽馬分布函數(shù)產(chǎn)生，填充深灰色的格子表示由正態(tài)分布函數(shù)產(chǎn)生。而假標(biāo)志物組中，每一個(gè)都是一個(gè)25（樣本）×900（變量）的矩陣，其數(shù)值由正態(tài)分布函數(shù)產(chǎn)生。

每個(gè)包含在真標(biāo)志物中小方格都是一個(gè)20（樣本）×10（變量）的矩陣。矩陣每列的值都由伽馬分布函數(shù)（利用R語(yǔ)言中rgamma函數(shù)實(shí)現(xiàn)）產(chǎn)生。但對(duì)于假標(biāo)志物組，每個(gè)格子都是一個(gè)20（樣本）×300（變量）的矩陣，其數(shù)值也是由伽馬分布函數(shù)產(chǎn)生。

1.2分析流程

歸一化：為了減少原始數(shù)據(jù)的噪聲，增強(qiáng)mRMR方法選擇具有識(shí)別能力的變量，模擬數(shù)據(jù)集需要進(jìn)行離散化，即用原始數(shù)據(jù)的均值（μ）和標(biāo)準(zhǔn)差（σ）對(duì)數(shù)據(jù)進(jìn)行離散化。任何數(shù)據(jù)大于μ+σ/2轉(zhuǎn)換為1，小于μ-σ/2轉(zhuǎn)換為-1，其他數(shù)據(jù)轉(zhuǎn)換為0。同時(shí)，原始的讀長(zhǎng)數(shù)目需要進(jìn)行歸一化，轉(zhuǎn)換為相對(duì)豐度，即每個(gè)變量的讀長(zhǎng)數(shù)除以所有樣本在該變量中的讀長(zhǎng)總數(shù)，每個(gè)變量的總和為1（變量中80%都是0將被忽略）。主要分析流程：歸一化后的數(shù)據(jù)采用變量篩選和自助重抽樣2個(gè)步驟進(jìn)行去冗余，具體流程見(jiàn)圖1。第一步的參數(shù)為1～M，其中M為第一次變量篩選時(shí)被mRMR篩選出的候選變量，用于區(qū)分不同樣本（可能含有冗余變量）；第二步為自助重抽樣，參數(shù)為2～B；第三步為變量排序，參數(shù)是3～M′，這些變量是上一步抽樣中被mRMR選出的，當(dāng)所有的自助重抽樣與變量選取完成后，按照變量出現(xiàn)次數(shù)進(jìn)行排序，選取最終M′個(gè)變量作為最終用戶需要的變量（M>M′）。

2結(jié)果與討論

2.1基于宏基因組的自助抽樣方法的參數(shù)選擇

此方法過(guò)程主要包括3步：變量篩選步驟、自助重抽樣和變量篩選過(guò)程以及變量排序，整個(gè)過(guò)程包含3個(gè)主要參數(shù)，分別是M、M′、B，它們對(duì)于選擇生物標(biāo)記物的質(zhì)量有重大影響。對(duì)于模擬數(shù)據(jù)集S1來(lái)說(shuō)，參數(shù)M設(shè)置為50。當(dāng)M等于50時(shí)，幾乎全部的非冗余變量都會(huì)被mRMR從1 000個(gè)變量中選出，考慮到計(jì)算的效率，50已經(jīng)足夠，因此沒(méi)有選擇更高的標(biāo)準(zhǔn)。對(duì)于參數(shù)B的選擇，筆者設(shè)置了一系列自助重抽樣次數(shù)的梯度，結(jié)果顯示，當(dāng)B超過(guò)40時(shí)，被選擇出來(lái)的真標(biāo)記物s不再增加（由于原始數(shù)據(jù)s的不固定性，因此選擇多個(gè)s來(lái)表征數(shù)據(jù)的變化趨勢(shì)）（圖2）。B值設(shè)為40。同樣的道理，對(duì)于數(shù)據(jù)集S2、S3中M′的選擇，結(jié)果與S1具有一致性（圖3）。由于S1只包含10個(gè)真標(biāo)記物組，因此參數(shù)M′設(shè)置為10（最為理想的結(jié)果是每個(gè)標(biāo)志物組中含有1個(gè)變量M′）。因此，將整體數(shù)據(jù)集參數(shù)M、B、M′分別設(shè)置為50、40、10。對(duì)本研究中基于自助抽樣的生物標(biāo)志物選擇方法進(jìn)行了去冗余性和準(zhǔn)確性分析，來(lái)考察本方法是否更適合于宏基因組數(shù)據(jù)分析。冗余率、非冗余率計(jì)算公式如下：

冗余率=冗余的標(biāo)志物數(shù)目選擇標(biāo)志物總數(shù)目×100%；（1）

非冗余率=特異的真生物標(biāo)記物數(shù)選擇標(biāo)志物總數(shù)目×100%。（2）

2.2去冗余性分析

由圖4可知，對(duì)于數(shù)據(jù)集S2、S3，本研究的新方法得到了最好的分析結(jié)果（表3），同時(shí)在數(shù)據(jù)集S1中，也得到了很好的區(qū)分效果。此外，本研究基于自助抽樣的新方法較其他方法得到了更多的非冗余真標(biāo)志物。宏基因組數(shù)據(jù)量龐大，各種各樣的微生物基因片段都包含其中，表征微生物種屬特性及其功能的特異性標(biāo)準(zhǔn)是研究生物標(biāo)志物的意義所在。在復(fù)雜的數(shù)據(jù)庫(kù)中尋找特異的生物標(biāo)志物來(lái)重構(gòu)菌群的復(fù)雜性，因此其選擇的冗余性不可避免。本試驗(yàn)基于自助抽樣方法很好地解決了冗余性這個(gè)難題，對(duì)于后續(xù)宏基因組工作有重要的應(yīng)用價(jià)值。

2.3穩(wěn)健性分析

基于3個(gè)模擬數(shù)據(jù)集，筆者分析比較了本方法與其他已經(jīng)在宏基因組研究中應(yīng)用的方法（如LEfSe、Metastats、Wilcoxon）在穩(wěn)健性方面存在的差異。對(duì)于每種方法，選擇100個(gè)生

物標(biāo)志物（等于每個(gè)數(shù)據(jù)集中真生物標(biāo)志物數(shù)目）計(jì)算100個(gè)生物標(biāo)志物的百分率，結(jié)果見(jiàn)圖5、表4。在已有的研究方法中，Wilcoxon在3個(gè)模擬數(shù)據(jù)集中的穩(wěn)健性是最高的，本方法與Wilcoxon方法在3個(gè)數(shù)據(jù)集上相當(dāng)，甚至表現(xiàn)更好?；诤昊蚪M數(shù)據(jù)生物標(biāo)志物選擇的方法，選擇出的生物標(biāo)志物具有較少的冗余固然重要，但是能夠選擇出在不同分組樣本中有差異的生物標(biāo)志物是前提。本方法的穩(wěn)健性能夠保證選出的生物標(biāo)志物能夠代表或者區(qū)分不同的樣本，只有這樣的生物標(biāo)志物才有生物學(xué)意義。

2.4分類準(zhǔn)確性分析

分類準(zhǔn)確性是生物標(biāo)志物選擇方法是否具有競(jìng)爭(zhēng)力的重要指標(biāo)。分類準(zhǔn)確率計(jì)算公式如下：

分類準(zhǔn)確率=準(zhǔn)確分類的樣本數(shù)目測(cè)試樣本中樣本總數(shù)×100%。（3）

此部分只采用S2及S3作為驗(yàn)證分類準(zhǔn)確率與否的數(shù)據(jù)集，由于S1數(shù)據(jù)集內(nèi)部區(qū)分非常明顯，對(duì)于任何一種區(qū)分方法都能實(shí)現(xiàn)很好的分類結(jié)果，因此在后2個(gè)數(shù)據(jù)集中分析比較這幾種方法的優(yōu)劣更有意義。分類時(shí)，使用這4種方法選擇的10個(gè)標(biāo)志物來(lái)建模。其中，每個(gè)數(shù)據(jù)集都有2類，每類含有60個(gè)樣本，采用50個(gè)樣本作為訓(xùn)練數(shù)據(jù)集，10個(gè)樣本作為檢驗(yàn)數(shù)據(jù)集，結(jié)果顯示，在2個(gè)數(shù)據(jù)集準(zhǔn)確性的分析中，基于自助抽樣的方法較其他3種方法具有更高的分類準(zhǔn)確性以及最小的區(qū)分結(jié)果變異性，即最小的s（標(biāo)準(zhǔn)方差）值（圖6）。分類準(zhǔn)確性是筆者選擇方法的一個(gè)重要指標(biāo)，基于自助抽樣方法與其他生物標(biāo)志物選擇方法相比，在分類準(zhǔn)確性方面具有非常明顯的優(yōu)勢(shì)，在今后對(duì)于宏基因組研究中，本方法可以很好地實(shí)現(xiàn)對(duì)于生物標(biāo)志物的選擇。

3結(jié)論

目前宏基因組數(shù)據(jù)缺乏生物標(biāo)志物的背景信息，使得利用各種方法預(yù)測(cè)宏基因組生物標(biāo)志物變得困難[7]。筆者提出了將基于自助抽樣的方法用于宏基因組生物標(biāo)志物的鑒定，它是一個(gè)自上而下的方法，結(jié)合了mRMR方法和自助重抽樣技術(shù)。基于模擬數(shù)據(jù)集，通過(guò)其與2種自上而下的方法（Metastats、LEfSe）以及自下而上的方法（Wilcoxon秩和檢驗(yàn)）進(jìn)行對(duì)比，表明本方法可以在較高準(zhǔn)確率的基礎(chǔ)上更加穩(wěn)健地選擇更多的非冗余生物標(biāo)志物。但本方法在鑒定功能性的生物標(biāo)志物方面不是非常理想，還需進(jìn)一步完善。

參考文獻(xiàn)：

[1]Ndimba B K，Ndimba R J，Johnson T S，et al.Biofuels as a sustainable energy source：An update of the applications of proteomics in bioenergy crops and algae[J]. Journal of Proteomics，2013，93：234-244.

[2]Pedros-Alio C. Marine microbial diversity：can it be determined？ [J]. Trends in Microbiology，2006，14（6）：257-263.

[3]Liao，L，Xu X W，Jiang X W，et al. Microbial diversity in deep-sea sediment from the cobalt-rich crust deposit region in the Pacific Ocean[J]. Microbiology Ecology，2011，78（3）：565-585.

[4]Bauer D F. Constructing confidence sets using rank statistics[J]. Journal of the American Statistical Association，1972，67（339）：687-690.

[5]Ding C，Peng H C. Minimum redundancy feature selection from microarray gene expression data[C]. Proceedings of the 2003 IEEE Bioinformatics Conference，2003：523-528.

[6]Huang S，Li R，Zeng X W，et al. Predictive modeling of gingivitis severity and susceptibility via oral microbiota[J]. The ISME Journal，2014，8（9）：1768-1780.

[7]高岳. 應(yīng)用宏基因組技術(shù)從微生物中獲得活性物質(zhì)的研究進(jìn)展[J]. 江蘇農(nóng)業(yè)科學(xué)，2014，42（1）：5-8.趙劍波，郭繼英，姜全，等. 桃抗重茬砧木GF677組培快繁技術(shù)[J]. 江蘇農(nóng)業(yè)科學(xué)，2016，44（5）：60-61，68.