一種基于宏基因組模擬數據的生物標志物篩選方法

王曉君++滕琳

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

關鍵詞：宏基因組；生物標志物；mRMR；自助抽樣法

中圖分類號： Q789文獻標志碼： A文章編號：1002-1302（2016）05-0056-04

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

1材料與方法

1.1模擬數據集的產生

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

在真標志物中，一個小方格是一個25（樣本）×10（變量）的矩陣。矩陣每一列的值都是由正態分布函數或者伽馬分布函數（利用R語言中rnorm或者rgamma函數實現）產生的。表格中填充淺灰色的格子表示由伽馬分布函數產生，填充深灰色的格子表示由正態分布函數產生。而假標志物組中，每一個都是一個25（樣本）×900（變量）的矩陣，其數值由正態分布函數產生。

每個包含在真標志物中小方格都是一個20（樣本）×10（變量）的矩陣。矩陣每列的值都由伽馬分布函數（利用R語言中rgamma函數實現）產生。但對于假標志物組，每個格子都是一個20（樣本）×300（變量）的矩陣，其數值也是由伽馬分布函數產生。

1.2分析流程

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

2結果與討論

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

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

冗余率=冗余的標志物數目選擇標志物總數目×100%；（1）

非冗余率=特異的真生物標記物數選擇標志物總數目×100%。（2）

2.2去冗余性分析

由圖4可知，對于數據集S2、S3，本研究的新方法得到了最好的分析結果（表3），同時在數據集S1中，也得到了很好的區分效果。此外，本研究基于自助抽樣的新方法較其他方法得到了更多的非冗余真標志物。宏基因組數據量龐大，各種各樣的微生物基因片段都包含其中，表征微生物種屬特性及其功能的特異性標準是研究生物標志物的意義所在。在復雜的數據庫中尋找特異的生物標志物來重構菌群的復雜性，因此其選擇的冗余性不可避免。本試驗基于自助抽樣方法很好地解決了冗余性這個難題，對于后續宏基因組工作有重要的應用價值。

2.3穩健性分析

基于3個模擬數據集，筆者分析比較了本方法與其他已經在宏基因組研究中應用的方法（如LEfSe、Metastats、Wilcoxon）在穩健性方面存在的差異。對于每種方法，選擇100個生

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

2.4分類準確性分析

分類準確性是生物標志物選擇方法是否具有競爭力的重要指標。分類準確率計算公式如下：

分類準確率=準確分類的樣本數目測試樣本中樣本總數×100%。（3）

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

3結論

目前宏基因組數據缺乏生物標志物的背景信息，使得利用各種方法預測宏基因組生物標志物變得困難[7]。筆者提出了將基于自助抽樣的方法用于宏基因組生物標志物的鑒定，它是一個自上而下的方法，結合了mRMR方法和自助重抽樣技術?；谀M數據集，通過其與2種自上而下的方法（Metastats、LEfSe）以及自下而上的方法（Wilcoxon秩和檢驗）進行對比，表明本方法可以在較高準確率的基礎上更加穩健地選擇更多的非冗余生物標志物。但本方法在鑒定功能性的生物標志物方面不是非常理想，還需進一步完善。

參考文獻：

[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]高岳. 應用宏基因組技術從微生物中獲得活性物質的研究進展[J]. 江蘇農業科學，2014，42（1）：5-8.趙劍波，郭繼英，姜全，等. 桃抗重茬砧木GF677組培快繁技術[J]. 江蘇農業科學，2016，44（5）：60-61，68.