車丹丹 郭 順 姜青山
1(中國科學院深圳先進技術研究院 深圳 518055)
2(中國科學院大學深圳先進技術學院 深圳 518055)
解析基因調控網(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)計方法以提高模型準確率。
調控網(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
1989 年,一種基于統(tǒng)計分析的重構基因調控網(wǎng)絡方法被提出[9],該方法使用相關系數(shù)來定義基因之間的相似性度量并以此推導基因網(wǎng)絡。這種方法的主要缺點是相關系數(shù)很難識別基因之間更加復雜的依賴關系。
為解決基于統(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ù)的評分來估計基因之間的依賴關系。以上僅基于信息論的方法均存在一個缺點——識別率有限,對于大量冗余和不相關特征無法做有效識別。
自 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ù)的時候。
自 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)計學方法處理。
(2)學習過程:為了識別所有潛在的調控關系,學習過程將每個子集定義為統(tǒng)計上的特征選擇問題。對于小樣本的數(shù)據(jù),直接通過特征選擇方法來求解效果并不理想,本文選擇 XGBoost 模型,通過改變訓練樣本的權重,同時學習多個基模型,并將多個基模型線性組合以提高性能。
(3)融合過程:該過程將學習過程得到的每個子集的調控關系進行融合,最終將所有的調控關系根據(jù)權重大小進行排序。對于不同子問題分別建模得到的權重,存在量綱不一致問題,本文在初始排序基礎上增加歸一化和統(tǒng)計方法以提高模型準確率。
圖 2 調控網(wǎng)絡推斷方法流程圖Fig. 2 Flow chart of regulation network inference method
圖 3 基因表達數(shù)據(jù)處理流程圖Fig. 3 Gene expression data processing flow chart
本文選擇 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)構建的。
由于每個子問題的建模過程都是獨立的,所以不能簡單地使用從每個子問題評估的監(jiān)管關系的可信度進行全局排名。因此,本文采用基于 L2 范數(shù)的規(guī)范化方法來解決這個問題,并且每個子問題的權重 wi,j被規(guī)范化為:
其中,p 為每個子問題中基因 j 的候選調控因子的數(shù)目。
此外,本文還使用統(tǒng)計方法,通過更新全局權重 wi,j來進一步細化推斷 GRN。這種改進是基于這樣一個假設:如果一個候選調控因子 i 調控多個靶基因,那么它將是一個重要的調控因子,并且所有調控關系的可信度都應該提高。根據(jù)這一點,本文將權重 wi,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ù)量。
前人的研究表明,隨機化和正則化在重建 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ù)“基于樹的”方法的默認設置相同。
在 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 則順序地學習這些弱學習器(每個基礎模型都依賴于前面的模型),并按照某種確定性的策略將它們組合起來。
現(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
本文提出一種基于 Boosting 集成模型的特征選擇框架:針對基因靜態(tài)數(shù)據(jù)進行調控網(wǎng)絡推斷,建立特有的集成框架,同時通過在目標函數(shù)中加入 2 種不同的懲罰項,以及限制樹結構和剪枝過程,盡可能地避免過擬合,較大幅度地提高了預測準確率。另外,對于不同子問題分別建模得到的權重存在量綱不一致的問題,本文在初始排序基礎上增加歸一化和統(tǒng)計方法以提高模型準確率。使用來自 DREAM5 挑戰(zhàn)的基準數(shù)據(jù)集測試結果表明,本文所提出的 XGBoost 比現(xiàn)有其他方法獲得更好的性能。在未來的研究中,將進一步研究利用先驗信息進行特征選擇,充分結合多方數(shù)據(jù),并增強該方法在其他生物物種基因調控關系應用中的可移植性。