車丹丹 郭 順 姜青山

1(中國科學院深圳先進技術研究院 深圳 518055)

2(中國科學院大學深圳先進技術學院 深圳 518055)

1 引 言

解析基因調控網(wǎng)絡(Gene Regulatory Networks，GRNs)的結構對生物信息學至關重要，因為它為生物有機體的發(fā)展機理及功能等提供了一個新的研究視角。隨著微陣列技術的發(fā)展，全基因組范圍上的基因表達都可被觀測到，這為從基因表達數(shù)據(jù)上推導調控網(wǎng)絡拓撲結構提供了契機?；蛘{控建模是根據(jù)基因表達數(shù)據(jù)所蘊含的信息而建立的反映基因與基因之間調控關系的網(wǎng)絡，而重構基因調控網(wǎng)絡是后基因組時代非常重要的研究領域[1]?；蛘{控網(wǎng)絡對所有生物種類和系統(tǒng)的作用是顯而易見的，因為其在維持生物有機體的功能方面發(fā)揮著重要作用[2]。因此，重構基因調控網(wǎng)絡有著廣泛的應用前景，它為藥物設計和醫(yī)療相關領域等提供了重要信息。

基因是遺傳信息的基本載體。雖然一個有機體中的所有基因都是相同的，但它們可以根據(jù)基因間的相互作用及網(wǎng)絡在不同組織中準確表達并執(zhí)行特定的功能[3]。在研究和認識基因之間及相應網(wǎng)絡的相互作用工作中，重構的基因調控網(wǎng)絡可作為一種工作模型，為研究者在實驗設計上提供輔助和形成新的假說。例如，Hecker 等[4]將重構基因調控網(wǎng)絡應用于卵巢癌，產(chǎn)生了一系列可檢驗的假說，并發(fā)現(xiàn)了一個潛在的藥物靶點。Camacho 等[5]提出機器學習和網(wǎng)絡生物學相結合的交叉學科有望在疾病生物學、藥物發(fā)現(xiàn)、微生物研究和合成生物學等領域取得重大突破。

Boosting 方法是一種強大且常用的統(tǒng)計學習方法，主要貢獻在于可以將弱學習算法提升為強學習算法，具有許多傳統(tǒng)方法所沒有的優(yōu)點，故其在基因表達數(shù)據(jù)上的應用十分廣泛。在分類學習中，Boosting 方法通過反復修改訓練數(shù)據(jù)的權值分布，構建一系列基本分類器，并將這些分類器進行線性組合，構成一個強分類器。對于回歸問題，Boosting 方法通過多次對訓練樣本做重抽樣建模，學習多個回歸器，并將這些回歸器進行組合，可較大幅度提高回歸模型的性能[6]。

本文針對靜態(tài)基因表達數(shù)據(jù)，致力于研究基因調控網(wǎng)絡重構方法，目標在于構建具有更好可靠性及準確率的基因調控網(wǎng)絡。具體地，針對真實的靜態(tài)基因表達數(shù)據(jù)，建立特征選擇集成框架，選擇 Boosting 模型計算基因調控關系的初始權重，并在初始排序基礎上增加歸一化和統(tǒng)計方法以提高模型準確率。

2 基因調控網(wǎng)絡推斷方法研究現(xiàn)狀

調控網(wǎng)絡已廣泛應用于各種生物系統(tǒng)的結構建模，且已開發(fā)出一系列方法來構建可靠的生物網(wǎng)絡[7]。Lee 和 Tzou[8]通過分析從基因表達數(shù)據(jù)上推導調控網(wǎng)絡的不同計算方法，表明現(xiàn)有方法均具有不同程度的準確性和復雜性，雖然實用性已有所提高，但準確率仍待提高。由于是從基因表達數(shù)據(jù)來推導調控網(wǎng)絡拓撲結構，因此這些方法被稱為基因調控網(wǎng)絡的逆向工程(Reverse-Engineering)，也可稱為重構基因調控網(wǎng)絡(Reconstructing GRNs)。

按照時間順序，國內外學者從基因表達數(shù)據(jù)上推導調控網(wǎng)絡拓撲結構的問題研究工作主要有 4 個階段(見圖 1):基于統(tǒng)計分析的方法[9]、基于信息論的方法[10-12]、基于概率圖模型的方法[13-17]和基于機器學習的方法[18-26]。

圖 1 現(xiàn)有研究方法發(fā)展歷程Fig. 1 Development of existing research methods

2.1 基于統(tǒng)計分析的方法

1989 年，一種基于統(tǒng)計分析的重構基因調控網(wǎng)絡方法被提出[9]，該方法使用相關系數(shù)來定義基因之間的相似性度量并以此推導基因網(wǎng)絡。這種方法的主要缺點是相關系數(shù)很難識別基因之間更加復雜的依賴關系。

2.2 基于信息論的方法

為解決基于統(tǒng)計分析方法的局限性，研究者提出了一些基于信息論的方法——使用基因之間的互信息(Mutual Information)作為度量并以此推導調控網(wǎng)絡[10-12]。由于在相關網(wǎng)絡中，一般存在間接調控關系，故需采用一些方法來消除間接調控關系的影響。例如，CMI2NI[10]通過計算包含和排除兩個基因之間邊緣的假設分布之間的 Kullback-Leibler 差異，來量化給定兩個基因之間的相互信息；BC3NET(Bagging C3NET)[11]使用基于對互信息值連同最大步長(Maximization Step)的估計來進一步提高所重構基因調控網(wǎng)絡的識別率；ANOVerence[12]引入元信息(Meta-Information)來推導調控網(wǎng)絡并使用基于相關系數(shù)的評分來估計基因之間的依賴關系。以上僅基于信息論的方法均存在一個缺點——識別率有限，對于大量冗余和不相關特征無法做有效識別。

2.3 基于概率圖模型的方法

自 2000 年起，許多基于概率圖模型的方法(如貝葉斯網(wǎng)絡方法)被廣泛應用于重構基因調控網(wǎng)絡。Friedman 等[13]首次提出了基于貝葉斯網(wǎng)絡來表示統(tǒng)計依賴關系的框架——利用貝葉斯網(wǎng)絡學習工具從微陣列數(shù)據(jù)中恢復基因相互作用；Liu 等[14]提出了局部貝葉斯網(wǎng)絡方法，利用網(wǎng)絡分解策略和偽正邊緣消除方法從基因表達數(shù)據(jù)中推斷 GRNs。DBNCS(Dynamic Bayesian Network Comprehensive Score)算法將綜合評分與動態(tài)貝葉斯模型相結合，首次構造出具有多重時延的 GRNs[15]；Xing 等[16]提出了一種基于互信息和斷點檢測的候選自動選擇算法(CAS)來限制搜索空間，以加速貝葉斯網(wǎng)絡的學習過程，但效率仍低于同期部分其他方法。De Campos 等[17]利用結構約束的貝葉斯網(wǎng)絡學習算法來對基因表達數(shù)據(jù)構建 GRNs。

然而，靜態(tài)貝葉斯網(wǎng)絡方法只能構建無環(huán)網(wǎng)絡，并需將數(shù)據(jù)進行離散化處理從而造成信息丟失，而動態(tài)貝葉斯網(wǎng)絡方法則通常局限于在時間序列數(shù)據(jù)上的應用。另外，無論是從理論上還是從計算效率上來看，貝葉斯網(wǎng)絡的結構學習都是一項很大的挑戰(zhàn)，尤其是當該類方法應用于維度特別高的基因表達數(shù)據(jù)的時候。

2.4 基于機器學習的方法

自 2009 年起，出現(xiàn)了許多對重構基因調控網(wǎng)絡方法比較評估方面的研究。其中，DREAM(Dialogue for Reverse Engineering Assessments and Methods)挑戰(zhàn)[18]為相關領域的研究者們提供了基準數(shù)據(jù)集來驗證評估他們的工作。其中，基于特征選擇框架的集成方法在 DREAM 挑戰(zhàn)基準數(shù)據(jù)集上有較為突出的識別率。例如，GENIE3(GEne Network Inference with Ensemble of Trees)[19]在學習過程中使用隨機森林(Random Forests)的特征選擇方法，但該方法在理論上不具有較好的可理解性。為此，TIGRESS[20]在學習過程中使用最小角回歸方法(Least Angle Regression)的特征選擇方法并結合穩(wěn)定性選擇(Stability Selection)來解決重構基因調控網(wǎng)絡問題。NIMEFI(Network Inference Using Multiple Ensemble Feature Importance Algorithms)[21]考慮了不同基于特征選擇框架的集成方法的互補性，將 GENIE3、E-SVR(Ensemble Support Vector Regression)以及 E-EL(Ensemble Elasitc Net)等方法置于統(tǒng)一的框架下進行學習，并以此解決重構基因調控網(wǎng)絡問題。然而，NIMEFI 最大的問題在于其所需的參數(shù)遠大于其他方法，這給模型的參數(shù)選擇帶來很大的挑戰(zhàn)，并且極大地增加了該模型的不確定性。Guo 等[22]利用基于偏最小二乘的線性方法對 GRNs 進行建模，但在真實實驗的數(shù)據(jù)集上該線性模型仍存在明顯局限性，模型準確率低；Chi 和 Liu[23]利用基于模糊邏輯和神經(jīng)網(wǎng)絡的認知模糊影響圖(FCMs)對 GRNs 進行建模，但預測結果準確率不高；Deng 等[24]利用微分方程進行網(wǎng)絡推斷，并引入了一個具有自適應數(shù)值微分的線性微分方程模型，該模型可擴展到非常大的調節(jié)網(wǎng)絡；Petralia 等[25]提出了一個靈活統(tǒng)一的集成框架，允許將來自異類數(shù)據(jù)的信息共同考慮用于 GRNs 推斷；Zheng 等[26]提出了一種結合互信息的集成框架，首先對候選調控基因進行預加權，然后利用 MARS(Multivariate Adaptive Regression Splines)檢測非線性調控鏈，但該方法存在模型過擬合問題，無法準確提取基因調控關系。

由于基因表達數(shù)據(jù)通常含有大量冗余和不相關特征，而現(xiàn)有基于機器學習的模型或結構復雜參數(shù)過多導致過擬合嚴重，模型不穩(wěn)定；或使用線性模型等簡單模型進行擬合，效果差、準確率低。本文提出一種基于 Boosting 集成模型的方法[6](XGBoost)，應用隨機化和正則化來解決模型過擬合問題，同時針對不同子問題建模所得權重不一致問題，對初始權重增加歸一化和統(tǒng)計學方法處理。

3 基于 XGBoost 的基因靜態(tài)數(shù)據(jù)調控網(wǎng)絡推斷方法

(2)學習過程:為了識別所有潛在的調控關系，學習過程將每個子集定義為統(tǒng)計上的特征選擇問題。對于小樣本的數(shù)據(jù)，直接通過特征選擇方法來求解效果并不理想，本文選擇 XGBoost 模型，通過改變訓練樣本的權重，同時學習多個基模型，并將多個基模型線性組合以提高性能。

(3)融合過程:該過程將學習過程得到的每個子集的調控關系進行融合，最終將所有的調控關系根據(jù)權重大小進行排序。對于不同子問題分別建模得到的權重，存在量綱不一致問題，本文在初始排序基礎上增加歸一化和統(tǒng)計方法以提高模型準確率。

3.1 分解過程

圖 2 調控網(wǎng)絡推斷方法流程圖Fig. 2 Flow chart of regulation network inference method

圖 3 基因表達數(shù)據(jù)處理流程圖Fig. 3 Gene expression data processing flow chart

3.2 學習過程

本文選擇 Boosting 模型評估基因調控關系的重要性，多次對訓練樣本做重抽樣并建模，學習多個回歸器，并將這些回歸器進行組合，提高回歸性能。由于傳統(tǒng)的 Boosting 集成學習方法需要學習多個弱學習器，訓練時間相對較長，而 XGBoost 模型[6]使用二階導數(shù)的信息來幫助迭代訓練，損失函數(shù)值將能更快地下降，提高訓練速度，獲得高性能模型。XGBoost 的目標函數(shù)可以表述為:

其中，T 為樹結構中的樹的葉子節(jié)點個數(shù)； 和 λ 為控制收縮的參數(shù)；w 為葉子權重。正則化項可以對最終學習得到的權重進行平滑，即 XGBoost 模型通過 LASSO(L1)和 Ridge(L2)正則化來懲罰更復雜的模型，從而避免過擬合問題。

本文應用特征變量 Gi的個數(shù) Ni來分割所有樹結構中的目標變量，作為 Gi的重要性。分割標準和其他細節(jié)可參考文獻[6]。針對靜態(tài)表達式數(shù)據(jù)，本文選擇 XGBoost 模型作為公式(1)中的 f 來解決問題，其中模型是通過公式(2)構建的。

3.3 融合過程

由于每個子問題的建模過程都是獨立的，所以不能簡單地使用從每個子問題評估的監(jiān)管關系的可信度進行全局排名。因此，本文采用基于 L2 范數(shù)的規(guī)范化方法來解決這個問題，并且每個子問題的權重 wi,j被規(guī)范化為:

其中，p 為每個子問題中基因 j 的候選調控因子的數(shù)目。

此外，本文還使用統(tǒng)計方法，通過更新全局權重 wi,j來進一步細化推斷 GRN。這種改進是基于這樣一個假設:如果一個候選調控因子 i 調控多個靶基因，那么它將是一個重要的調控因子，并且所有調控關系的可信度都應該提高。根據(jù)這一點，本文將權重 wi,j更新表述為:

4 結果分析與評估

4.1 數(shù)據(jù)集與評價指標

本課題的數(shù)據(jù)集及評價指標均來自于 DREAM 挑戰(zhàn)平臺(基因調控網(wǎng)絡挑戰(zhàn))[18]。該平臺為生物學和醫(yī)學的研究者們提出了眾多挑戰(zhàn)，并提供了基準數(shù)據(jù)集來驗證評估他們的工作。其中，DREAM5 是第一個利用大規(guī)模真實數(shù)據(jù)集構建 GRNs 的挑戰(zhàn)，其中靶基因達到 O(103)數(shù)量級、調控基因則為 O(102)。

DREAM5 數(shù)據(jù)集共包含 4 個網(wǎng)絡:網(wǎng)絡 1 是通過 in-silico 模擬導出的，另外 3 個網(wǎng)絡則是從不同實驗中獲得的。其中，網(wǎng)絡 2 來源于金葡萄桿菌(S.Aureus)相關實驗，網(wǎng)絡 3 來源于原核生物(E.coli)實驗，網(wǎng)絡 4 則來源于真核生物(S.cerevisiae)實驗。網(wǎng)絡 3 和 4 調控關系的金標準來自于 RegulonDB 和 Gene Ontology(GO)兩個數(shù)據(jù)庫，而網(wǎng)絡 2 沒有對應的金標準，故本文暫不考慮。最終選取 DREAM5 挑戰(zhàn)的 1、3 和 4 網(wǎng)絡數(shù)據(jù)作為實驗數(shù)據(jù)。

本文選取來自 DREAM5 挑戰(zhàn)的靜態(tài)基因表達數(shù)據(jù)集，基本信息如表 1 所示。其中，insilico 網(wǎng)絡包括 1 643 個靶基因、195 個轉錄因子(Transcription Factors)，金標準中共有 4 012 條調控關系；E.coli 網(wǎng)絡包括 4 511 個靶基因、334 個轉錄因子，金標準中共有 2 066 條調控關系；S.cerevisiae 網(wǎng)絡包括 5 950 個靶基因、333 個轉錄因子，金標準中共有 3 940 條調控關系。

表 1 靜態(tài)基因表達數(shù)據(jù)集Table 1 Steady gene expression data set

為評價該方法的性能，本文考慮兩種常用的評價指標:接受者操作特征曲線面積(Area Under Receiver Operating Characteristic，AUROC)和正確率-召回率曲線面積(Area Under Precision-Recall Curves，AUPR)。其中，AUROC 是基于真陽率(True Positive Rate，TPR)與假陽率(False Positive Rate，F(xiàn)PR)的接受者操作特征(ROC)范圍，AUPR 是根據(jù)精確度(Precision)與召回率(Recall)得出的領域。真陽率(TPR)為檢測出來的真陽性樣本數(shù)除以所有真實陽性樣本數(shù)；假陽率(FPR)為檢測出來的假陽性樣本數(shù)除以所有真實陰性樣本數(shù)；召回率為真陽性樣本在所有檢測正確樣本總數(shù)中的占比。

其中，TP(True Positive)為真陽性的數(shù)量；TN(True Negative)為真陰性的數(shù)量；FP(False Positive)為假陽性的數(shù)量；FN(False Negative)為假陰性的數(shù)量。

4.2 隨機化和正則化

前人的研究表明，隨機化和正則化在重建 GRNs 中是有效的。其中，隨機化包括樣本抽樣和特征選擇，如引導程序和子特征。正則化是通過將懲罰項添加到目標函數(shù)，進而控制模型的復雜性，回歸模型中最常用的正則化技術是 LASSO(L1)和 Ridge(L2)。

本文方法 XGBNet 基于 XGBoost。其中，XGBoost python 軟件包提供了用于實現(xiàn)的各種參數(shù)，本文選擇決策樹作為基學習器。參數(shù) max_depth 和 min_child_weight 與模型中每棵樹的結構相關，且都設置為 4；將控制每棵樹中訓練樣本比率的參數(shù) subsample 設置為 0.7；參數(shù)colsample_bytree 控制每棵樹中特征(候選調節(jié)器)的比率，并在此處設置為 0.9；學習率 eta 設置為 0.000 8，樹的數(shù)量設置為 1 000，與大多數(shù)“基于樹的”方法的默認設置相同。

4.3 實驗結果

在 DREAM5 數(shù)據(jù)集中采用 in-silico 模擬數(shù)據(jù)、E.coli 和 S.cerevisiae 實驗數(shù)據(jù)來對所提出模型的性能進行評估。更進一步地，選擇了幾種最新的 GRNs 推斷方法，包括 iRafNet[25]、HiDi[24]、PLSNET[22]和 DREAM 挑戰(zhàn)賽[18]的獲勝者與本文結果進行對比分析，具體結果如表 2 所示。

表 2 為不同方法在 DREAM5 數(shù)據(jù)集中 3 種網(wǎng)絡上的 AUPR 和 AUROC。其中，KO(Knock Out Data)為剔除某個基因后的基因表達數(shù)據(jù)；SS(Steady-State Expression Data)為平穩(wěn)基因表達數(shù)據(jù)。從表 2 可以看出，iRafNet 和 HiDi 都集成了穩(wěn)態(tài)表達數(shù)據(jù)和剔除數(shù)據(jù)，而 DREAM5 挑戰(zhàn)賽的冠軍僅使用穩(wěn)態(tài)表達數(shù)據(jù)，本方法也僅需用到穩(wěn)態(tài)表達數(shù)據(jù)，但就網(wǎng)絡 1 的 AUPR 和 AUROC 以及網(wǎng)絡 3 和 4 的 AUROC 而言，XGBNet 仍然比其他方法具有更好的性能。一個主要原因可能是，由于不完整的 KO 數(shù)據(jù)所提供的信息很少，對于推斷出 GRNs 所起的作用很小。此外，由于生物實驗數(shù)據(jù)的采集誤差大、獲取途徑不一致等問題，導致 5 種方法在網(wǎng)絡 3 和 4 的 AUPR 值均較低。

從表 2 數(shù)值可看出，本文方法在 in-silico 生成的模擬數(shù)據(jù)集中，AUPR 和 AUROC 兩個評估指標均顯著優(yōu)于現(xiàn)有方法；在 E.coli 和 S.cerevisiae 兩種生物的真實實驗數(shù)據(jù)中，AUROC 指標均高于現(xiàn)有最優(yōu)方法。推測原因在于，本文所提出方法建立特有的集成框架，通過在目標函數(shù)中加入 2 種不同的懲罰項，以及限制樹結構和剪枝過程，盡可能地避免過擬合，較大幅度地提高了預測準確率；同時用于基因調控網(wǎng)絡的多為 bagging 集成方法(以隨機森林為典型代表)，而以 XGBoost 和 AdaBoost 為主的 Boosting 集成方法性能則通常優(yōu)于 bagging 集成。兩種方法的主要區(qū)別在于 bagging 的每個弱學習器都是獨立并行學習的，而 Boosting 則順序地學習這些弱學習器(每個基礎模型都依賴于前面的模型)，并按照某種確定性的策略將它們組合起來。

4.4 討論與分析

現(xiàn)階段對基因靜態(tài)數(shù)據(jù)調控網(wǎng)絡推斷的研究，主要集中在基于信息論、概率圖和機器學習的方法中。Zhang 等[10]、Ricardo 等[11]和 Küffner 等[12]使用基因之間的互信息作為度量并以此推導調控網(wǎng)絡，來進一步消除間接調控關系的影響，然而僅使用互信息仍不能足夠多地提取基因調控的有效信息；Liu 等[14]、Yu 等[15]、Xing 等[16]和 de Campos 等[17]利用葉斯網(wǎng)絡學習工具從靜態(tài)基因數(shù)據(jù)中計算基因相互作用，從最初的靜態(tài)貝葉斯到改良的動態(tài)貝葉斯，做出了很多的突破，然而效率仍低于現(xiàn)有部分其他方法，同時模型準確率也有待提高；Huynh-Thu 等[19]、Haury 等[20]、Ruyssinck 等[21]和 Guo 等[22]構建基于特征選擇框架的集成方法，且用于基因調控網(wǎng)絡的多為 bagging 集成方法(以隨機森林為典型代表)。與上述基因靜態(tài)數(shù)據(jù)調控網(wǎng)絡推斷研究不同的是，本文選擇基于 XGBoost 的集成方法，每個弱學習器都依賴于前一個模型，并按照某種確定性的策略將其組合起來，且著重解決模型過擬合問題，并對集成模型的結果增加歸一化和統(tǒng)計學方法處理。當然，本文還存在不足之處，未來工作的重點是進一步研究利用先驗信息進行特征選擇，并增強方法的可移植性。

表 2 不同方法結果比較Table 2 Results of different methods

5 結 論

本文提出一種基于 Boosting 集成模型的特征選擇框架:針對基因靜態(tài)數(shù)據(jù)進行調控網(wǎng)絡推斷，建立特有的集成框架，同時通過在目標函數(shù)中加入 2 種不同的懲罰項，以及限制樹結構和剪枝過程，盡可能地避免過擬合，較大幅度地提高了預測準確率。另外，對于不同子問題分別建模得到的權重存在量綱不一致的問題，本文在初始排序基礎上增加歸一化和統(tǒng)計方法以提高模型準確率。使用來自 DREAM5 挑戰(zhàn)的基準數(shù)據(jù)集測試結果表明，本文所提出的 XGBoost 比現(xiàn)有其他方法獲得更好的性能。在未來的研究中，將進一步研究利用先驗信息進行特征選擇，充分結合多方數(shù)據(jù)，并增強該方法在其他生物物種基因調控關系應用中的可移植性。