劉曉燕,張誠誠,郭茂祖,邢林林
1.哈爾濱工業大學 計算機科學與技術學院,哈爾濱 150001
2.北京建筑大學 電氣與信息工程學院,北京 100044
轉錄調控是基因調控的主要過程,即轉錄因子通過結合位點進而控制目標基因的表達[1]。構建高精度的基因轉錄調控網絡一直以來都是研究熱點,其主要研究利用實驗數據重構網絡[2]。隨著第二代大規模基因測序技術和高通量基因表達分析技術的發展,使得方便地獲取基因表達數據、基因序列數據成為可能。此外,隨著近年來對機器學習的不斷研究,對構建基因轉錄調控網絡的研究進入了一個新的階段,構建高精度基因轉錄調控網絡也成為可能。
轉錄調控的本質就是轉錄因子通過結合啟動子位點,進而控制目標基因的轉錄水平來完成相應的功能。基因表達數據反映的是直接或間接測量得到的基因轉錄產物mRNA在樣本中的豐度,這些數據可以用于分析哪些基因的表達發生了改變,基因之間有何相關性,在不同條件下基因的活動是如何受影響的。啟動子是基因的一個組成部分,控制基因表達(轉錄)的起始時間和表達的程度。啟動子本身并不控制基因活動,而是通過與轉錄因子結合,從而控制基因活動。
目前,構建基因轉錄調控網絡的算法主要有如下幾種。Bornholdt在2008年使用布爾網絡模型構建調控網絡,在網絡中每個基因有開、關兩個狀態,狀態“開”表示一個基因轉錄表達,形成基因產物,而狀態“關”則代表一個基因未轉錄[3]。布爾網絡雖然簡單,但是建模能力有所欠缺。Werhli等人在2007年使用貝葉斯網絡模型構建轉錄調控網絡[4]。貝葉斯網絡雖然可以建模比較復雜的網絡,但是計算復雜度過高。Huynh-Thu等人在2011年使用基于樹的模型 GENIE3(gene network inference with ensemble of trees)[5],將構建轉錄調控問題轉換為機器學習里面的回歸問題,在樹節點分裂過程中得到基因與基因之間調控重要性的排名。2015年,Huynh-Thu等人提出了新的算法模型Jump3[6],對GENIE3算法進行了改進,之前的算法只用了基因表達數據,在新模型中Huynh-Thu引入了啟動子狀態這一變量,并且對樹分裂節點標準進行了修改,從而帶來了精度上的提升,但是計算時間復雜度過高,程序運行時間非常長,訓練n個基因的轉錄調控網絡Jump3算法至少需要n個子模型,即隨機森林。Gillani[7]和Maetschke[8]等人在2014年將基因轉錄調控網絡構建問題轉為二分類問題,使用支持向量機進行訓練,并比較了不同核函數的性能,實驗中并未比較其他算法且只用了基因表達數據。Robiolo等人在2015年利用神經網絡對轉錄調控網絡進行建模[9],利用迭代次數、誤差構建評分體系,進而計算基因之間的調控可能性,此方法同樣存在計算復雜度過高的問題,得到n個基因的轉錄調控網絡至少需要(n×(n-1))/2個神經網絡模型。
本文主要研究的是轉錄調控網絡的構建算法。一方面,目前研究方法只用了基因表達數據;另一方面,目前對于基因功能研究方面已經形成了比較完善的GO注釋。GeneOntology(基因本體)是一種用于功能注釋的樹型結構數據[10],具有相同GO注釋的基因功能往往相似。于是本文利用基因表達數據、基因序列數據、基因注釋數據這3個數據來源。此外,綜合考慮精度和計算復雜度,本文將轉錄調控網絡構建問題轉化為二分類問題,提出了基于深度自編碼器和組合模型構建基因轉錄調控網絡。
本文組織結構如下:第2章介紹了數據獲取、數據預處理以及負例集的構建。第3章詳細闡述了本文提出的基于深度自編碼器的XGBoost和邏輯回歸組合模型DAXL(combined model with XGBoost and logistic regression based on deep AutoEncoder)。第 4章以擬南芥調控網絡的構建為例,對DAXL方法進行了驗證,并與GENIE方法進行了對比,給出了實驗結果。第5章對全文進行總結。
轉錄因子活性會明顯影響目標基因的表達水平,因此兩者在轉錄水平上的表達譜具有相關性[11]。本文從GEO(gene expression omnibus)上獲取編號為GSE41212[12]的基因表達數據,該數據集從時間和空間兩個角度詳細分析了擬南芥種子發芽的過程,記錄了擬南芥種子的4個組織在9個時間點的基因表達值。基因序列數據主要包括兩部分:啟動子序列數據、蛋白質的氨基酸序列數據。其中啟動子的序列數據是從TAIR上下載,氨基酸序列數據是從http://www.arabidopsis.org/上下載。本文的GO數據可以從http://geneontology.org/page/download-annotations上下載。正例集是已經被實驗驗證的,確定有基因轉錄調控關系的轉錄因子、目標基因對,也是后續進行模型訓練必要的數據。AGRIS收錄了19 013對具有直接調控關系的轉錄因子和目標基因。本文選取該數據集作為正例集,可從http://arabidopsis.med.ohiostate.edu/上下載。
本文收集的生物數據有兩個主要問題:一是基因序列數據、基因注釋數據都是字符串,無法直接進行模型訓練;二是把基因表達數據、基因序列數據、基因注釋數據以基因ID作為主鍵進行鏈接時存在缺失數據。
本文獲得基因表達數據、基因序列數據、基因注釋數據的具體處理如下。
從GEO上獲取的基因表達數據是數值型的,在與基因序列數據、基因注釋數據進行以基因ID為主鍵的鏈接時,存在部分缺失的情況。因為模型中包括集成樹模型,對缺失值不敏感,所以本文直接用-1填充缺失值,從而得到116維特征:

在生物學中,相鄰的3個核苷酸被稱為密碼子,每個密碼子編碼一個氨基酸,因為對于啟動子序列數據而言,每一條啟動子序列,都可以得到全部密碼子的頻次,對于每一個轉錄因子都可以得到如下64維特征:

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

轉錄因子通過調控目標基因參與同一生物過程,因此兩者往往具有相同的功能。轉錄因子在基因本體(GO)數據庫中,每個GO ID是一個唯一的標識符,對應某一術語名稱。因此,可以利用GO術語建立特征向量[13]。本文利用擬南芥的GO數據,最終大約取得3 112維特征。具體步驟如下:
(1)根據選定的數據集中的基因ID,提取對應每一個ID的所有GO術語。
(2)將提取出的GO術語依次編號。
(3)假設共有n個GO術語,若基因含有i(i=1,2,…,n)號GO術語,則為1,否則為0,則每個基因功能信息可以表示成一個n維向量:

通過數據預處理,對于每一對轉錄因子、目標基因,可以得到如下一條向量:

將轉錄調控網絡問題轉換為機器學習二分類問題,需要構建負樣本。本文從生物背景出發,構建負樣本采用Yu等人在文獻[14]中提出的方法。構建負樣本的基本思想是:假設基因為G,轉錄因子為T,如果基因G不包含轉錄因子T的任意一個結合位點,則認為(T,G)這一對基因為負樣本。
首先使用深度自編碼器訓練原始高維基因注釋數據,得到新的低維可以表征基因注釋數據的向量;然后把基因表達數據、基因序列數據、新的向量交給XGBoost進行訓練,對于訓練好的XGBoost模型,可以得到每一條樣本最終落在每一個樹上葉子節點的信息,從而進行01編碼,即樣本落在某一個葉子節點上為1,否則為0;最后將01編碼向量交給邏輯回歸訓練進行分類,最終得到轉錄調控網絡。整體算法流程如圖1所示。

Fig.1 Flow chart of DAXL圖1 DAXL整體流程圖
經過數據預處理的特征有6 540維,其中基因注釋數據特征特別稀疏,不利于基于集成樹的模型進行訓練。為了解決這個問題,本文利用深度自編碼器[15-16]進行Embedding,即先訓練深度自編碼器模型,然后對于訓練好的模型直接取其中隱層作為新的輸入樣本。
深度自動編碼器即多層自編碼器,其輸入輸出一致,內部經過多層的編碼、解碼,具有很好的學習能力。本文使用深度自編碼器完成高維特征的embedding目標。利用深度自編碼器隱層單元數量可以自行定義的特性,將隱層節點設置為從高到低,再從低到高。然后提取最中間隱層學習到的向量作為原始輸入,如圖2所示,進而解決基因注釋數據高維且稀疏的問題。

Fig.2 DeepAutoEncoder圖2 深度自編碼器
如圖2所示,本文使用五層的深度自編碼器,第一層(L1)和第五層(L5)是相同的向量,中間有3個隱藏層。以第一層與第五層向量差的平方作為學習目標,使用反向傳播算法學習參數,通過多輪迭代后,得到最終的網絡結構。本文為了獲得可以表征原始輸入向量的低維向量,用第三層(L3)的輸出表示原始輸入向量。
XGBoost模型最早在2014年由Chen等人提出[17],它是boosting的一個擴展變體。XGBoost在以決策樹為基礎的基學習器構建boosting集成的基礎上,使用泰勒展開式近似目標函數,利用一階、二階導數,從而提高精度。從圖3中可以看出,XGBoost由多個弱分類器組成,且每一分類器的學習都相互獨立,比較容易進行數據并行、特征并行,因此并行粒度大。此外,樹模型是一個典型的if-else結構,可以處理非線性特征問題。

Fig.3 Schematic diagram of XGBoost圖3 XGBoost示意圖
邏輯回歸模型(logistic regression)是最大熵模型的一個特例[18],它的雛形是線性回歸,是一個線性模型,不可以直接處理非線性問題。通過sigmoid函數進行轉換,將預測值映射到[0,1],從而轉變為分類模型,比較適合01編碼訓練數據。可以通過極大似然估計、梯度下降得到相對較優解。從圖4中可以看出,邏輯回歸可以通過行并行、列并行的方式計算梯度。
為了充分發揮非線性、線性模型的優勢,借鑒Facebook 2014年在計算廣告領域提出的新模型[19],本文構建了組合模型,如圖5所示。
樣本X先通過XGBoost進行訓練,XGBoost訓練完成后,對于每一條輸入向量,可以根據葉子節點的信息得到新的一維向量。舉例來說,在圖5中,訓練好的XGBoost模型有兩棵樹,第一棵樹有3個葉子節點,第二棵樹有2個葉子節點。如果一條樣本在XGBoost模型的第一棵樹上落在第一個葉子節點,在XGBoost模型的第二棵樹上落在第一個葉子節點,在0-1 Code層可以得到(1,0,0,1,0)這個新向量,然后這個向量再通過邏輯回歸模型進行訓練。
用TSNE(t-distributed stochastic neighbor embedding)[20]將高維向量映射到二維空間,如圖6、圖7所示。圖6是原始基因表達數據,圖7是經過01編碼之后得到的新的數據。其中紅色的點是正例,藍色的點是負例,從而可以看出經過01編碼之后,數據集正負例集相對分開,更加利于模型進行分類。
原始基因注釋數據經過離散化后有3 112維,深度自編碼器在一定程度上可以降維且可以學習到一個更好的表達形式,從而提高后續組合模型的精度。圖8展示了深度自編碼器中間層節點數在10、50、100、200、500時對XGBoost、XGBoost+LR最終預測結果F1值的影響,其中節點數為0時表示沒有使用深度自編碼器對基因注釋數據進行處理。一方面說明組合模型DAXL精度比XGBoost高;另一方面說明深度自編碼器可以對離散化的基因注釋進行充分學習,從而得到更好的表達形式。
由于構建負例集時存在假陽性問題,即實際是正例可能被認為是負例進行訓練,本文使用F1作為評價指標。F1評價指標綜合考慮準確率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 深度自編碼器中間層節點數對XGBoost、XGBoost+LR的影響

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

Table 1 Result of experiments inArabidopsis dataset表1 在擬南芥數據集上的實驗結果
從表2中可以看出,DAXL方法比只使用基因表達數據的GENIE3方法好。同時,在相同的數據集上,DAXL模型也比支持向量機方法好,因為支持向量機對缺失值敏感,且核函數的選擇對精度影響很大。

Table 2 Result of experiments inArabidopsis dataset表2 在擬南芥數據集上的實驗結果
為了進一步驗證本文算法的可靠性以及本文工作的意義,從http://arabidopsis.med.ohio-sta te.edu/得到最新的轉錄調控關系,其中有69對轉錄調控關系在原來的訓練驗證過程中沒有出現過,本文使用DAXL模型對這69對調控關系進行預測。如圖9所示,發現有62對的預測概率在0.5以上,有47對的預測概率在0.8以上。

Fig.9 Predictive probability of novel transcriptional regulatory relationships圖9 新的轉錄調控關系預測概率
本文分析了目前經典的基因轉錄調控網絡構建算法,針對存在的問題,提出了基于深度自編碼器和組合模型的方法。本文方法充分利用了基因表達數據、基因序列數據、基因注釋數據,用深度自編碼器巧妙地解決了基因注釋數據高維稀疏的問題;使用XGBoost與邏輯回歸的組合模型,結合生物背景知識,把基因轉錄調控網絡構建問題轉為二分類問題,從而構建了基因轉錄調控網絡。本文在真實的擬南芥數據集上對提出的方法進行了驗證,并與現有方法進行了比較。結果說明本文方法比單一采用LR、XGBoost更加有效,較對比方法GENIE3提高了15個百分點。