基于組合模型的轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)構(gòu)建算法研究*

劉曉燕，張誠誠，郭茂祖，邢林林

1.哈爾濱工業(yè)大學(xué) 計算機科學(xué)與技術(shù)學(xué)院，哈爾濱 150001

2.北京建筑大學(xué) 電氣與信息工程學(xué)院，北京 100044

1 引言

轉(zhuǎn)錄調(diào)控是基因調(diào)控的主要過程，即轉(zhuǎn)錄因子通過結(jié)合位點進(jìn)而控制目標(biāo)基因的表達(dá)[1]。構(gòu)建高精度的基因轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)一直以來都是研究熱點，其主要研究利用實驗數(shù)據(jù)重構(gòu)網(wǎng)絡(luò)[2]。隨著第二代大規(guī)模基因測序技術(shù)和高通量基因表達(dá)分析技術(shù)的發(fā)展，使得方便地獲取基因表達(dá)數(shù)據(jù)、基因序列數(shù)據(jù)成為可能。此外，隨著近年來對機器學(xué)習(xí)的不斷研究，對構(gòu)建基因轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)的研究進(jìn)入了一個新的階段，構(gòu)建高精度基因轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)也成為可能。

轉(zhuǎn)錄調(diào)控的本質(zhì)就是轉(zhuǎn)錄因子通過結(jié)合啟動子位點，進(jìn)而控制目標(biāo)基因的轉(zhuǎn)錄水平來完成相應(yīng)的功能。基因表達(dá)數(shù)據(jù)反映的是直接或間接測量得到的基因轉(zhuǎn)錄產(chǎn)物mRNA在樣本中的豐度，這些數(shù)據(jù)可以用于分析哪些基因的表達(dá)發(fā)生了改變，基因之間有何相關(guān)性，在不同條件下基因的活動是如何受影響的。啟動子是基因的一個組成部分，控制基因表達(dá)（轉(zhuǎn)錄）的起始時間和表達(dá)的程度。啟動子本身并不控制基因活動，而是通過與轉(zhuǎn)錄因子結(jié)合，從而控制基因活動。

目前，構(gòu)建基因轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)的算法主要有如下幾種。Bornholdt在2008年使用布爾網(wǎng)絡(luò)模型構(gòu)建調(diào)控網(wǎng)絡(luò)，在網(wǎng)絡(luò)中每個基因有開、關(guān)兩個狀態(tài)，狀態(tài)“開”表示一個基因轉(zhuǎn)錄表達(dá)，形成基因產(chǎn)物，而狀態(tài)“關(guān)”則代表一個基因未轉(zhuǎn)錄[3]。布爾網(wǎng)絡(luò)雖然簡單，但是建模能力有所欠缺。Werhli等人在2007年使用貝葉斯網(wǎng)絡(luò)模型構(gòu)建轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)[4]。貝葉斯網(wǎng)絡(luò)雖然可以建模比較復(fù)雜的網(wǎng)絡(luò)，但是計算復(fù)雜度過高。Huynh-Thu等人在2011年使用基于樹的模型 GENIE3（gene network inference with ensemble of trees）[5]，將構(gòu)建轉(zhuǎn)錄調(diào)控問題轉(zhuǎn)換為機器學(xué)習(xí)里面的回歸問題，在樹節(jié)點分裂過程中得到基因與基因之間調(diào)控重要性的排名。2015年，Huynh-Thu等人提出了新的算法模型Jump3[6]，對GENIE3算法進(jìn)行了改進(jìn)，之前的算法只用了基因表達(dá)數(shù)據(jù)，在新模型中Huynh-Thu引入了啟動子狀態(tài)這一變量，并且對樹分裂節(jié)點標(biāo)準(zhǔn)進(jìn)行了修改，從而帶來了精度上的提升，但是計算時間復(fù)雜度過高，程序運行時間非常長，訓(xùn)練n個基因的轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)Jump3算法至少需要n個子模型，即隨機森林。Gillani[7]和Maetschke[8]等人在2014年將基因轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)構(gòu)建問題轉(zhuǎn)為二分類問題，使用支持向量機進(jìn)行訓(xùn)練，并比較了不同核函數(shù)的性能，實驗中并未比較其他算法且只用了基因表達(dá)數(shù)據(jù)。Robiolo等人在2015年利用神經(jīng)網(wǎng)絡(luò)對轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)進(jìn)行建模[9]，利用迭代次數(shù)、誤差構(gòu)建評分體系，進(jìn)而計算基因之間的調(diào)控可能性，此方法同樣存在計算復(fù)雜度過高的問題，得到n個基因的轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)至少需要(n×(n-1))/2個神經(jīng)網(wǎng)絡(luò)模型。

本文主要研究的是轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)的構(gòu)建算法。一方面，目前研究方法只用了基因表達(dá)數(shù)據(jù)；另一方面，目前對于基因功能研究方面已經(jīng)形成了比較完善的GO注釋。GeneOntology（基因本體）是一種用于功能注釋的樹型結(jié)構(gòu)數(shù)據(jù)[10]，具有相同GO注釋的基因功能往往相似。于是本文利用基因表達(dá)數(shù)據(jù)、基因序列數(shù)據(jù)、基因注釋數(shù)據(jù)這3個數(shù)據(jù)來源。此外，綜合考慮精度和計算復(fù)雜度，本文將轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)構(gòu)建問題轉(zhuǎn)化為二分類問題，提出了基于深度自編碼器和組合模型構(gòu)建基因轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)。

本文組織結(jié)構(gòu)如下：第2章介紹了數(shù)據(jù)獲取、數(shù)據(jù)預(yù)處理以及負(fù)例集的構(gòu)建。第3章詳細(xì)闡述了本文提出的基于深度自編碼器的XGBoost和邏輯回歸組合模型DAXL（combined model with XGBoost and logistic regression based on deep AutoEncoder）。第 4章以擬南芥調(diào)控網(wǎng)絡(luò)的構(gòu)建為例，對DAXL方法進(jìn)行了驗證，并與GENIE方法進(jìn)行了對比，給出了實驗結(jié)果。第5章對全文進(jìn)行總結(jié)。

2 數(shù)據(jù)獲取、預(yù)處理及負(fù)例集的構(gòu)建

2.1 數(shù)據(jù)獲取

轉(zhuǎn)錄因子活性會明顯影響目標(biāo)基因的表達(dá)水平，因此兩者在轉(zhuǎn)錄水平上的表達(dá)譜具有相關(guān)性[11]。本文從GEO（gene expression omnibus）上獲取編號為GSE41212[12]的基因表達(dá)數(shù)據(jù)，該數(shù)據(jù)集從時間和空間兩個角度詳細(xì)分析了擬南芥種子發(fā)芽的過程，記錄了擬南芥種子的4個組織在9個時間點的基因表達(dá)值。基因序列數(shù)據(jù)主要包括兩部分：啟動子序列數(shù)據(jù)、蛋白質(zhì)的氨基酸序列數(shù)據(jù)。其中啟動子的序列數(shù)據(jù)是從TAIR上下載，氨基酸序列數(shù)據(jù)是從http://www.arabidopsis.org/上下載。本文的GO數(shù)據(jù)可以從http://geneontology.org/page/download-annotations上下載。正例集是已經(jīng)被實驗驗證的，確定有基因轉(zhuǎn)錄調(diào)控關(guān)系的轉(zhuǎn)錄因子、目標(biāo)基因?qū)Γ彩呛罄m(xù)進(jìn)行模型訓(xùn)練必要的數(shù)據(jù)。AGRIS收錄了19 013對具有直接調(diào)控關(guān)系的轉(zhuǎn)錄因子和目標(biāo)基因。本文選取該數(shù)據(jù)集作為正例集，可從http://arabidopsis.med.ohiostate.edu/上下載。

2.2 數(shù)據(jù)預(yù)處理

本文收集的生物數(shù)據(jù)有兩個主要問題：一是基因序列數(shù)據(jù)、基因注釋數(shù)據(jù)都是字符串，無法直接進(jìn)行模型訓(xùn)練；二是把基因表達(dá)數(shù)據(jù)、基因序列數(shù)據(jù)、基因注釋數(shù)據(jù)以基因ID作為主鍵進(jìn)行鏈接時存在缺失數(shù)據(jù)。

本文獲得基因表達(dá)數(shù)據(jù)、基因序列數(shù)據(jù)、基因注釋數(shù)據(jù)的具體處理如下。

從GEO上獲取的基因表達(dá)數(shù)據(jù)是數(shù)值型的，在與基因序列數(shù)據(jù)、基因注釋數(shù)據(jù)進(jìn)行以基因ID為主鍵的鏈接時，存在部分缺失的情況。因為模型中包括集成樹模型，對缺失值不敏感，所以本文直接用-1填充缺失值，從而得到116維特征:

在生物學(xué)中，相鄰的3個核苷酸被稱為密碼子，每個密碼子編碼一個氨基酸，因為對于啟動子序列數(shù)據(jù)而言，每一條啟動子序列，都可以得到全部密碼子的頻次，對于每一個轉(zhuǎn)錄因子都可以得到如下64維特征：

蛋白質(zhì)序列可以表示成由20個氨基酸字符組成的字符串，其字符集合為{A,V,L,I,F,P,M,S,T,C,W,Y,N,Q,D,E,K,R,H,G}，對于每一條蛋白質(zhì)序列數(shù)據(jù)，統(tǒng)計這20個字符的出現(xiàn)頻次，從而得到如下的20維特征：

轉(zhuǎn)錄因子通過調(diào)控目標(biāo)基因參與同一生物過程，因此兩者往往具有相同的功能。轉(zhuǎn)錄因子在基因本體（GO）數(shù)據(jù)庫中，每個GO ID是一個唯一的標(biāo)識符，對應(yīng)某一術(shù)語名稱。因此，可以利用GO術(shù)語建立特征向量[13]。本文利用擬南芥的GO數(shù)據(jù)，最終大約取得3 112維特征。具體步驟如下：

（1）根據(jù)選定的數(shù)據(jù)集中的基因ID，提取對應(yīng)每一個ID的所有GO術(shù)語。

（2）將提取出的GO術(shù)語依次編號。

（3）假設(shè)共有n個GO術(shù)語，若基因含有i(i=1,2,…,n)號GO術(shù)語，則為1，否則為0，則每個基因功能信息可以表示成一個n維向量：

通過數(shù)據(jù)預(yù)處理,對于每一對轉(zhuǎn)錄因子、目標(biāo)基因，可以得到如下一條向量：

2.3 構(gòu)建負(fù)例集

將轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)問題轉(zhuǎn)換為機器學(xué)習(xí)二分類問題，需要構(gòu)建負(fù)樣本。本文從生物背景出發(fā)，構(gòu)建負(fù)樣本采用Yu等人在文獻(xiàn)[14]中提出的方法。構(gòu)建負(fù)樣本的基本思想是：假設(shè)基因為G，轉(zhuǎn)錄因子為T，如果基因G不包含轉(zhuǎn)錄因子T的任意一個結(jié)合位點，則認(rèn)為(T,G)這一對基因為負(fù)樣本。

3 基于深度自編碼器的XGBoost、邏輯回歸組合模型DAXL

首先使用深度自編碼器訓(xùn)練原始高維基因注釋數(shù)據(jù)，得到新的低維可以表征基因注釋數(shù)據(jù)的向量；然后把基因表達(dá)數(shù)據(jù)、基因序列數(shù)據(jù)、新的向量交給XGBoost進(jìn)行訓(xùn)練，對于訓(xùn)練好的XGBoost模型，可以得到每一條樣本最終落在每一個樹上葉子節(jié)點的信息，從而進(jìn)行01編碼，即樣本落在某一個葉子節(jié)點上為1，否則為0；最后將01編碼向量交給邏輯回歸訓(xùn)練進(jìn)行分類，最終得到轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)。整體算法流程如圖1所示。

Fig.1 Flow chart of DAXL圖1 DAXL整體流程圖

3.1 基于深度自編碼器的基因注釋數(shù)據(jù)降維方法

經(jīng)過數(shù)據(jù)預(yù)處理的特征有6 540維，其中基因注釋數(shù)據(jù)特征特別稀疏，不利于基于集成樹的模型進(jìn)行訓(xùn)練。為了解決這個問題，本文利用深度自編碼器[15-16]進(jìn)行Embedding，即先訓(xùn)練深度自編碼器模型，然后對于訓(xùn)練好的模型直接取其中隱層作為新的輸入樣本。

深度自動編碼器即多層自編碼器，其輸入輸出一致，內(nèi)部經(jīng)過多層的編碼、解碼，具有很好的學(xué)習(xí)能力。本文使用深度自編碼器完成高維特征的embedding目標(biāo)。利用深度自編碼器隱層單元數(shù)量可以自行定義的特性，將隱層節(jié)點設(shè)置為從高到低，再從低到高。然后提取最中間隱層學(xué)習(xí)到的向量作為原始輸入，如圖2所示，進(jìn)而解決基因注釋數(shù)據(jù)高維且稀疏的問題。

Fig.2 DeepAutoEncoder圖2 深度自編碼器

如圖2所示，本文使用五層的深度自編碼器，第一層（L1）和第五層（L5）是相同的向量，中間有3個隱藏層。以第一層與第五層向量差的平方作為學(xué)習(xí)目標(biāo)，使用反向傳播算法學(xué)習(xí)參數(shù)，通過多輪迭代后，得到最終的網(wǎng)絡(luò)結(jié)構(gòu)。本文為了獲得可以表征原始輸入向量的低維向量，用第三層（L3）的輸出表示原始輸入向量。

3.2 基于XGBoost的01編碼器

XGBoost模型最早在2014年由Chen等人提出[17]，它是boosting的一個擴(kuò)展變體。XGBoost在以決策樹為基礎(chǔ)的基學(xué)習(xí)器構(gòu)建boosting集成的基礎(chǔ)上，使用泰勒展開式近似目標(biāo)函數(shù)，利用一階、二階導(dǎo)數(shù)，從而提高精度。從圖3中可以看出，XGBoost由多個弱分類器組成，且每一分類器的學(xué)習(xí)都相互獨立，比較容易進(jìn)行數(shù)據(jù)并行、特征并行，因此并行粒度大。此外，樹模型是一個典型的if-else結(jié)構(gòu)，可以處理非線性特征問題。

Fig.3 Schematic diagram of XGBoost圖3 XGBoost示意圖

3.3 基于邏輯回歸的轉(zhuǎn)錄調(diào)控關(guān)系預(yù)測方法

邏輯回歸模型（logistic regression）是最大熵模型的一個特例[18]，它的雛形是線性回歸，是一個線性模型，不可以直接處理非線性問題。通過sigmoid函數(shù)進(jìn)行轉(zhuǎn)換，將預(yù)測值映射到[0,1]，從而轉(zhuǎn)變?yōu)榉诸惸Ｐ停容^適合01編碼訓(xùn)練數(shù)據(jù)。可以通過極大似然估計、梯度下降得到相對較優(yōu)解。從圖4中可以看出，邏輯回歸可以通過行并行、列并行的方式計算梯度。

3.4 DAXL算法

為了充分發(fā)揮非線性、線性模型的優(yōu)勢，借鑒Facebook 2014年在計算廣告領(lǐng)域提出的新模型[19]，本文構(gòu)建了組合模型，如圖5所示。

樣本X先通過XGBoost進(jìn)行訓(xùn)練，XGBoost訓(xùn)練完成后，對于每一條輸入向量，可以根據(jù)葉子節(jié)點的信息得到新的一維向量。舉例來說，在圖5中，訓(xùn)練好的XGBoost模型有兩棵樹，第一棵樹有3個葉子節(jié)點，第二棵樹有2個葉子節(jié)點。如果一條樣本在XGBoost模型的第一棵樹上落在第一個葉子節(jié)點，在XGBoost模型的第二棵樹上落在第一個葉子節(jié)點，在0-1 Code層可以得到（1,0,0,1,0）這個新向量，然后這個向量再通過邏輯回歸模型進(jìn)行訓(xùn)練。

4 實驗

4.1 01編碼有效性驗證

用TSNE（t-distributed stochastic neighbor embedding）[20]將高維向量映射到二維空間，如圖6、圖7所示。圖6是原始基因表達(dá)數(shù)據(jù)，圖7是經(jīng)過01編碼之后得到的新的數(shù)據(jù)。其中紅色的點是正例，藍(lán)色的點是負(fù)例，從而可以看出經(jīng)過01編碼之后，數(shù)據(jù)集正負(fù)例集相對分開，更加利于模型進(jìn)行分類。

4.2 深度自編碼器中間層參數(shù)對結(jié)果的影響

原始基因注釋數(shù)據(jù)經(jīng)過離散化后有3 112維，深度自編碼器在一定程度上可以降維且可以學(xué)習(xí)到一個更好的表達(dá)形式，從而提高后續(xù)組合模型的精度。圖8展示了深度自編碼器中間層節(jié)點數(shù)在10、50、100、200、500時對XGBoost、XGBoost+LR最終預(yù)測結(jié)果F1值的影響，其中節(jié)點數(shù)為0時表示沒有使用深度自編碼器對基因注釋數(shù)據(jù)進(jìn)行處理。一方面說明組合模型DAXL精度比XGBoost高；另一方面說明深度自編碼器可以對離散化的基因注釋進(jìn)行充分學(xué)習(xí)，從而得到更好的表達(dá)形式。

4.3 轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)構(gòu)建結(jié)果

由于構(gòu)建負(fù)例集時存在假陽性問題，即實際是正例可能被認(rèn)為是負(fù)例進(jìn)行訓(xùn)練，本文使用F1作為評價指標(biāo)。F1評價指標(biāo)綜合考慮準(zhǔn)確率P和召回率R[21]。

Fig.4 Logistic regression圖4 邏輯回歸

Fig.5 Combined model圖5 組合模型

Fig.6 Initial image based on TSNE圖6 原始TSNE圖

Fig.7 0-1 code image based on TSNE圖7 0-1編碼TSNE圖

Fig.8 Influence of the number of nodes in inter layer of deepAutoEncoder on XGBoost and XGBoost+LR圖8 深度自編碼器中間層節(jié)點數(shù)對XGBoost、XGBoost+LR的影響

從表1中可以看出，XGBoost模型比邏輯回歸模型精度高，因為XGBoost具有非線性擬合能力。DAXL模型的精度比XGBoost高，因為采用了01編碼，同時綜合了XGBoost、邏輯回歸模型的優(yōu)點。

Table 1 Result of experiments inArabidopsis dataset表1 在擬南芥數(shù)據(jù)集上的實驗結(jié)果

從表2中可以看出，DAXL方法比只使用基因表達(dá)數(shù)據(jù)的GENIE3方法好。同時，在相同的數(shù)據(jù)集上，DAXL模型也比支持向量機方法好，因為支持向量機對缺失值敏感，且核函數(shù)的選擇對精度影響很大。

Table 2 Result of experiments inArabidopsis dataset表2 在擬南芥數(shù)據(jù)集上的實驗結(jié)果

4.4 預(yù)測新的轉(zhuǎn)錄調(diào)控關(guān)系

為了進(jìn)一步驗證本文算法的可靠性以及本文工作的意義，從http://arabidopsis.med.ohio-sta te.edu/得到最新的轉(zhuǎn)錄調(diào)控關(guān)系，其中有69對轉(zhuǎn)錄調(diào)控關(guān)系在原來的訓(xùn)練驗證過程中沒有出現(xiàn)過，本文使用DAXL模型對這69對調(diào)控關(guān)系進(jìn)行預(yù)測。如圖9所示，發(fā)現(xiàn)有62對的預(yù)測概率在0.5以上，有47對的預(yù)測概率在0.8以上。

Fig.9 Predictive probability of novel transcriptional regulatory relationships圖9 新的轉(zhuǎn)錄調(diào)控關(guān)系預(yù)測概率

5 總結(jié)

本文分析了目前經(jīng)典的基因轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)構(gòu)建算法，針對存在的問題，提出了基于深度自編碼器和組合模型的方法。本文方法充分利用了基因表達(dá)數(shù)據(jù)、基因序列數(shù)據(jù)、基因注釋數(shù)據(jù)，用深度自編碼器巧妙地解決了基因注釋數(shù)據(jù)高維稀疏的問題；使用XGBoost與邏輯回歸的組合模型，結(jié)合生物背景知識，把基因轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)構(gòu)建問題轉(zhuǎn)為二分類問題，從而構(gòu)建了基因轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)。本文在真實的擬南芥數(shù)據(jù)集上對提出的方法進(jìn)行了驗證，并與現(xiàn)有方法進(jìn)行了比較。結(jié)果說明本文方法比單一采用LR、XGBoost更加有效，較對比方法GENIE3提高了15個百分點。