基于語義特征向量的DNA與轉(zhuǎn)錄因子結(jié)合特異性預(yù)測(cè)

孫曉雨 權(quán)麗君 梅 杰 黃立群 呂 強(qiáng)*

1(蘇州大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院 江蘇 蘇州 215006) 2(江蘇省計(jì)算機(jī)信息處理技術(shù)重點(diǎn)實(shí)驗(yàn)室 江蘇 蘇州 215006)

0 引 言

表達(dá)是結(jié)構(gòu)基因在生物體內(nèi)的轉(zhuǎn)錄、翻譯及加工的過程，而大量TF與基因組的作用能夠影響基因表達(dá)[1]，其中DNA與特定TF蛋白的結(jié)合，即DNA序列特異性，對(duì)于轉(zhuǎn)錄和選擇性剪切這類基因的調(diào)控過程研究具有特別重要的意義。由于數(shù)據(jù)量擴(kuò)大和特征維度增加，通過生物實(shí)驗(yàn)手段獲得DNA序列特異性費(fèi)時(shí)費(fèi)力，因此，目前只有很少一部分的TF和DNA結(jié)合特異性得到明確的表征。近年來，隨著高通量基因組測(cè)序數(shù)據(jù)的發(fā)展和計(jì)算性能的提高，類似于ChIP-Seq[2-3]等實(shí)驗(yàn)產(chǎn)生的高通量數(shù)據(jù)對(duì)于使用計(jì)算手段預(yù)測(cè)DNA序列特異性、繼而尋找發(fā)現(xiàn)新的生物學(xué)機(jī)制起到了至關(guān)重要的作用。

目前，多種基于序列的計(jì)算方法已經(jīng)被提出來預(yù)測(cè)TF與DNA的結(jié)合特異性，包括以權(quán)重矩陣為代表的傳統(tǒng)計(jì)算方法和以機(jī)器學(xué)習(xí)為代表的新型學(xué)習(xí)方法等。位置權(quán)重矩陣(Position Weight Matrix，PWM)[4]和二核苷酸權(quán)重矩陣(Dinucleotide Weight Matrix，DWM)[5]描述了每一個(gè)堿基在每一個(gè)位置上出現(xiàn)的頻率。它們的不足在于假定每個(gè)或每兩個(gè)核苷酸對(duì)TF與DNA序列的結(jié)合貢獻(xiàn)是獨(dú)立的，這種粗暴的可加性假設(shè)給預(yù)測(cè)準(zhǔn)確度帶來了較大的偏差。傳統(tǒng)機(jī)器學(xué)習(xí)方法gkm-SVM[6]是一種基于支持向量機(jī)的分類模型，輸入DNA序列，輸出該序列的打分。但是gkm-SVM只能用于小于5 000的樣本數(shù)，而LS-gkm[7]在數(shù)據(jù)量方面進(jìn)行了升級(jí)，改善了這一問題。最近，深度學(xué)習(xí)模型已廣泛應(yīng)用于該領(lǐng)域。例如，DeepBind[8]成功應(yīng)用卷積神經(jīng)網(wǎng)絡(luò)(Convolution Neural Network，CNN)[9]來捕獲主題特征，其性能優(yōu)于以前的方法。DanQ[10]成功地將遞歸神經(jīng)網(wǎng)絡(luò)(Recurrent Neural Network，RNN)[11-12]附加到基于CNN的模型，以便可以捕獲核苷酸堿基的長程相互作用[13]。但是，目前所有現(xiàn)有的基于機(jī)器學(xué)習(xí)的方法都在數(shù)據(jù)集上分別訓(xùn)練多個(gè)模型對(duì)應(yīng)不同的TF和實(shí)驗(yàn)環(huán)境CELL，也就是不包括上下文限制。此外，這些方法也缺少對(duì)自身提取的特征值進(jìn)行生物意義上的解析，進(jìn)而影響對(duì)預(yù)測(cè)問題的深入研究。

本文提出一種結(jié)合深度學(xué)習(xí)和傳統(tǒng)機(jī)器學(xué)習(xí)方法的模型semanticSVC，預(yù)測(cè)TF與DNA序列的結(jié)合特異性。首先，對(duì)來自ENCODE項(xiàng)目[14]的ChIP-seq數(shù)據(jù)集進(jìn)行整合，使用以EMB-CNN-RNN-MLP為框架的深度學(xué)習(xí)模型進(jìn)行訓(xùn)練，獲取具有多TF全局語義的序列特征向量，本文稱作semanticFeature。然后以semantic-Feature作為DNA序列特征，通過支持向量分類機(jī)(Support Vector Classify，SVC)進(jìn)行目標(biāo)訓(xùn)練，其最終預(yù)測(cè)結(jié)果AUC(Area Under roc Curve)分布要優(yōu)于其他先進(jìn)的方法。此外，通過對(duì)semanticFeature的降維可視化處理，不僅解析了其具有的生化意義，而且證明了本文方法應(yīng)用的廣泛性。

1 實(shí)驗(yàn)方法

本文的semanticSVC將DNA序列作為輸入，并根據(jù)TF、實(shí)驗(yàn)環(huán)境CELL信息，即TF所在的細(xì)胞環(huán)境，選擇相應(yīng)的SVC模型進(jìn)行特異性預(yù)測(cè)。結(jié)構(gòu)如圖1所示。

圖1 semanticSVC結(jié)構(gòu)

semanticSVC主要包括兩部分，第一部分是基于CNN-RNN深度框架的多TF統(tǒng)一模型，名為semanticDeep，該模型可以有效地對(duì)DNA序列進(jìn)行建模，不僅可以獲取序列的長程依賴關(guān)系，而且可以解讀與多個(gè)TF相關(guān)的特異性信息，這種特異性信息稱作“語義特征向量”，即semanticFeature。第二部分根據(jù)不同TF，分別對(duì)基于semanticFeature的SVC模型進(jìn)行構(gòu)建，最終預(yù)測(cè)出DNA與TF的結(jié)合性。我們將semanticFeature進(jìn)行可視化，可以直觀地幫助讀者進(jìn)行語義分析。

1.1 數(shù)據(jù)來源

本文數(shù)據(jù)均來自ENCODE(ENCyclopedia Of DNA Elements)項(xiàng)目[14](https://www.encodeproject.org)的608個(gè)不同的測(cè)定實(shí)驗(yàn)獲取的數(shù)據(jù)，測(cè)定實(shí)驗(yàn)主要根據(jù)TF和實(shí)驗(yàn)環(huán)境CELL進(jìn)行分類。為了驗(yàn)證semantic-Feature的優(yōu)越性，ENCODE項(xiàng)目數(shù)據(jù)被分為兩部分：506個(gè)用于semanticDeep模型的構(gòu)建，取名為D506；102個(gè)用于SVC模型的構(gòu)建，取名為D102。這兩部分的實(shí)驗(yàn)數(shù)據(jù)都有各自的訓(xùn)練集和測(cè)試集。

每個(gè)數(shù)據(jù)集包含n條DNA序列和其對(duì)應(yīng)的標(biāo)簽，標(biāo)簽為1表示該序列與TF可以結(jié)合，也就是正樣本；標(biāo)簽是0表示該序列與TF不能結(jié)合，即負(fù)樣本。序列長度為101，其中正樣本是通過ChIP-seq數(shù)據(jù)的峰值測(cè)量獲得，而負(fù)樣本序列通過采用dinucleotide shuffling經(jīng)典方法[8]模擬生成。

1.2 semanticSVC模型設(shè)計(jì)

本文的semanticDeep根據(jù)數(shù)據(jù)集D506進(jìn)行構(gòu)建。如圖2所示，模型的輸入是DNA序列的詞向量和TF的one-hot編碼。如果把一條DNA序列看作一個(gè)句子，那么DNA序列與TF的結(jié)合位點(diǎn)，即motif，可以看作句子的一個(gè)單詞，一個(gè)motif的長度一般是8，所以這里選取8作為詞向量的k-mer[15]。輸出是DNA與TF特異性的兩種分類結(jié)果:0或1。對(duì)于一個(gè)序列s、TF和CELL，semanticDeep可以抽象為四個(gè)階段，該過程表示為：

f(s,tf,cell)=

MLP(RNN(CNN(EMB(s))),(tf,cell))

(1)

圖2 semanticDeep結(jié)構(gòu)

模型第一層是EMB層(Embedding)。經(jīng)過EMB層，把每個(gè)序列替換為向量的索引，在訓(xùn)練過程中，向量會(huì)得到更新，通過索引尋找就更加快捷方便。

第二層是CNNs層。其由兩層1維CNN組成，第一層CNN有32個(gè)卷積核，第二層CNN有64個(gè)卷積核。如果把堿基比作字母，那么一條DNA序列就是一條沒有分詞的句子，CNN可以有效地將句子分成詞。每層卷積神經(jīng)網(wǎng)絡(luò)后又加了一層BatchNormalization[16]，提高網(wǎng)絡(luò)的泛化能力，避免梯度下降和梯度爆炸。

第三層是BRNN層。該層由三層RNN包括一層雙向LSTM和兩層雙向門循環(huán)單元(Gated Recurrent Unit，GRU)構(gòu)成。

LSTM中設(shè)計(jì)了輸入門(input gate)、忘記門(forget gate)和輸出門(output gate)對(duì)信息進(jìn)行長期處理，以獲得信息之間的長程依賴關(guān)系。LSTM模型的計(jì)算公式如下：

it=σ(Wixt+Uiht-1)

(2)

ft=σ(Wfxt+Ufht-1)

(3)

ot=σ(Woxt+Uoht-1)

(4)

(5)

(6)

ht=ot×tanhct

(7)

GRU的設(shè)計(jì)是將LSTM中的input gate和forget gate更改為更新門(update gate)，將output gate更改為重置門(reset gate)。GRU模型的計(jì)算公式如下：

zt=σ(Wzxt+Uzht-1)

(8)

rt=σ(Wrxt+Urht-1)

(9)

(10)

ht=(1-zt)×ht+zt×ht-1

(11)

式中：z表示update gate；r表示reset gate。

由于雙向GRU比雙向LSTM少一個(gè)門，所以參數(shù)更少、需要的存儲(chǔ)空間更少、訓(xùn)練時(shí)間更少；而雙向LSTM在數(shù)據(jù)集較大時(shí)的性能更好。在雙向LSTM后面加一層雙向GRU是為了減少模型參數(shù)、降低模型的過擬合。最后，經(jīng)過綜合考慮訓(xùn)練性能和速度，semanticDeep采用了雙向LSTM和雙向GRU混用的方式。DNA序列上的一個(gè)堿基與其他堿基是息息相關(guān)的，也就是這個(gè)基本單元不僅與它前面的基本單元相關(guān)，也與它后面基本單元相關(guān)。因此使用BRNN可以建模序列的雙向關(guān)系，即獲取CNNs層導(dǎo)出的分詞之間的相關(guān)性。

第四層由多層感知機(jī)(Multi-Layer Perception，MLP)構(gòu)成，映射一組輸入向量到一組輸出向量。每一層全連接到下一層，除了輸入節(jié)點(diǎn)，每個(gè)節(jié)點(diǎn)都帶有非線性激活函數(shù)。為了避免出現(xiàn)過擬合現(xiàn)象，我們還在MLP之間加入了Dropout[17]層，最后使用Sigmoid函數(shù)將輸出變成一個(gè)0至1之間的值。

對(duì)于D102數(shù)據(jù)集，通過semanticDeep的第三層BRNN層將DNA序列轉(zhuǎn)化為具有語義的64維向量semanticFeature，最后將semanticFeature和其對(duì)應(yīng)的序列標(biāo)簽作為本文SVC模型的輸入進(jìn)行訓(xùn)練。經(jīng)訓(xùn)練，每個(gè)TF和CELL對(duì)應(yīng)一個(gè)SVC模型，每個(gè)模型可以預(yù)測(cè)輸出基于semanticFeature的DNA與該TF的結(jié)合特異性。通過SVC，可以獲得數(shù)據(jù)集中的支持向量(Support Vector，SV)，支持向量意味著這些DNA樣例的預(yù)測(cè)可靠性不高。

對(duì)于semanticDeep模型訓(xùn)練，本文使用Keras作為框架，TensorFlow[18]作為后臺(tái)，NVIDIA GeForce GTX 1080Ti的GPU作為計(jì)算單元，SVC模型使用sklearn.svm[19]包實(shí)現(xiàn)。

1.3 semanticFeature可視化

為了更好更直觀地探索semanticFeature的可解釋性和生物意義, 通過降維可視化的方法，對(duì)它們的分布和位置變化進(jìn)行語義分析。因?yàn)?4維的向量無法直觀地體現(xiàn)semanticFeature的語義信息，為了更好地對(duì)semanticFeature在2維平面可視化，本文使用t-SNE降維工具[20]將64維的semanticFeature降至2維，其中每個(gè)2維向量代表一條DNA序列。然后將這些2維向量以坐標(biāo)的形式映射到2-D圖上，因此每個(gè)TF都可以形成一個(gè)散點(diǎn)圖，圖中的每個(gè)點(diǎn)都代表了對(duì)應(yīng)序列的2維向量。根據(jù)序列的標(biāo)簽將非支持向量的點(diǎn)分成兩種形狀，其中：三角形的點(diǎn)代表不結(jié)合的序列；Ⅹ型的點(diǎn)代表結(jié)合的序列；支持向量的點(diǎn)則用普通的圓點(diǎn)表示。該圖可以直觀地展現(xiàn)本文方法提取出的semanticFeature的語義信息，即可以將序列分成是否結(jié)合特定TF的相反語義信息類；此外，DNA序列之間存在的相關(guān)性信息，特別是當(dāng)序列發(fā)生突變時(shí)所引起的軌跡變化也可以直觀展現(xiàn)。

2 實(shí)驗(yàn)與結(jié)果分析

2.1 評(píng)價(jià)指標(biāo)

本文性能分析使用的評(píng)價(jià)指標(biāo)是一種二分類模型的評(píng)價(jià)指標(biāo)：AUC指標(biāo)。

(12)

式中：positives表示正樣本集;M是正樣本的數(shù)目；N是負(fù)樣本的數(shù)目；將樣本的score值按照從小到大的順序排序，最小的score對(duì)應(yīng)了樣本的rank1，以此類推得到ranki[2]。

2.2 性能分析

如圖3所示是semanticSVC和Deepbind、DanQ的AUC分布。

圖3 三種方法的AUC分布

統(tǒng)計(jì)D102數(shù)據(jù)集上SVC模型測(cè)試數(shù)據(jù)集的AUC分布，發(fā)現(xiàn)Deepbind方法的AUC分布在0.69～0.99之間，平均分?jǐn)?shù)為0.94，DanQ方法的AUC分布在0.71～0.99之間，平均分?jǐn)?shù)為0.90，而本文的semanticSVC方法的AUC分布在0.82～0.99之間，平均分?jǐn)?shù)達(dá)到0.96，平均準(zhǔn)確率達(dá)到0.92。盡管三個(gè)模型的AUC最高都接近于1.00，但是從圖3可以看出，本文semanticSVC方法的AUC分布比其他兩個(gè)模型更加集中，這說明本文模型的訓(xùn)練結(jié)果不僅優(yōu)于其他兩個(gè)模型，而且預(yù)測(cè)準(zhǔn)確度更加穩(wěn)定。此外，也驗(yàn)證了基于深度學(xué)習(xí)方法的semanticDeep抽取出的semanticFeature的優(yōu)越性。

接下來從數(shù)據(jù)集D102中隨機(jī)選擇5個(gè)測(cè)試集的AUC結(jié)果進(jìn)行比較。結(jié)果如表1所示。

表1 5個(gè)測(cè)試數(shù)據(jù)集性能對(duì)比

各個(gè)測(cè)試數(shù)據(jù)集的信息描述見表2，測(cè)試集所在環(huán)境指的是測(cè)試集對(duì)應(yīng)的TF和CELL。

表2 5個(gè)測(cè)試數(shù)據(jù)集信息所在環(huán)境

在測(cè)試集1上DeepBind方法的AUC在0.85以下，DanQ方法的AUC在0.9以下，而本文方法AUC達(dá)到了0.96；在測(cè)試集3和測(cè)試集5上，DeepBind方法和DanQ方法均沒有達(dá)到0.95，而本文方法達(dá)到了0.97；在其余兩個(gè)測(cè)試集上，三個(gè)方法的AUC差距不大，均接近于1，可見本文方法在一定程度上優(yōu)于其他兩個(gè)方法。

2.3 semanticFeature的可視化分析

圖4列舉了ENCODE項(xiàng)目中的ChIP-seq數(shù)據(jù)對(duì)應(yīng)的TF SP1 CELL K562和TF E2F1 CELL HeLa-S3的semanticFeature的散點(diǎn)圖，其中顏色最淺的Ⅹ型點(diǎn)代表序列標(biāo)簽是1的不屬于支持向量的序列，顏色較深的三角形點(diǎn)代表序列標(biāo)簽是0的不屬于支持向量的序列，顏色最深的中間的圓點(diǎn)代表屬于支持向量的序列。可以清晰地看到，屬于SV的序列幾乎把其他所有的序列都正確地分成兩部分，左半邊代表可以與TF結(jié)合的序列，右半部分代表不與TF結(jié)合的序列。已知在二分類問題的SVC算法中，SV可以決定出分割超平面的位置，也就是分類的決策邊界。而屬于SV的DNA序列的semanticFeature也恰好將正負(fù)樣本分割在它們的兩側(cè)。同時(shí)，其他100個(gè)TF對(duì)應(yīng)的ChIP-seq數(shù)據(jù)的semanticFeature的散點(diǎn)圖也具有同樣規(guī)律。

(a) TF SP1 CELL K562

(b) TF E2F1 CELL HeLa-S3圖4 semanticFeature可視化

因此，可以得到這樣的結(jié)論，通過semanticDeep獲得的semanticFeature可以抽取到與DNA特異性相關(guān)卻與TF無關(guān)的重要元級(jí)屬性。從而也解釋了為什么semanticSVC在僅僅使用了淺層學(xué)習(xí)方法的前提下，其預(yù)測(cè)結(jié)果卻優(yōu)于其他基于深度學(xué)習(xí)框架的預(yù)測(cè)器。

2.4 單點(diǎn)突變應(yīng)用分析

為了驗(yàn)證本文方法對(duì)于單點(diǎn)突變的敏感性，選取來自HGMD(The Human Gene Mutation Database)[21]的兩個(gè)基因HNF4A和BRCA1的四個(gè)突變序列對(duì)(Accession Number：CR1211032，CR133305，CR103902和CR090197)作為實(shí)驗(yàn)數(shù)據(jù)，DNA序列長度均為61，突變點(diǎn)位于序列中間。如圖4所示，“★”代表原始DNA序列， “◆”代表突變DNA序列，矩形框里包含的是CR1211032(圖4(a))和CR090197(圖4(b))，圓框中包含的CR133305，三角框中包含的CR103902。其中：CR1211032、CR133305、CR103902對(duì)其基因表達(dá)的減少和功能水平的降低具有影響[22-24]，而CR090197增加了啟動(dòng)因子的活性從而增強(qiáng)了與TF的結(jié)合性[25]。根據(jù)TRANSFAC database[26]存儲(chǔ)的實(shí)驗(yàn)數(shù)據(jù)，分別選擇與這三個(gè)結(jié)合位點(diǎn)有相互作用的TF SP1 CELL K562和TF E2F1 CELL HeLa-S3所對(duì)應(yīng)的模型進(jìn)行DNA特異性預(yù)測(cè)。根據(jù)semanticSVC預(yù)測(cè)可知，突變CR1211032(結(jié)合概率由0.99變?yōu)?.14)、CR133305(結(jié)合概率由0.55變?yōu)?.10)、CR103902(結(jié)合概率由0.57變?yōu)?.18)發(fā)生之后，其DNA序列與TF SP1的結(jié)合變低了，而突變CR076712(結(jié)合概率由0.88變?yōu)?.99)發(fā)生之后，該DNA與TF E2F1的結(jié)合變強(qiáng)了，以上預(yù)測(cè)結(jié)果都印證文獻(xiàn)[22-25]中的實(shí)驗(yàn)結(jié)果。除此之外，通過本文方法的semanticDeep抽取出突變序列對(duì)的semanticFeature，也正確地將突變前后的變化通過二維圖展現(xiàn)出來，例如，圖4(a)中，基因HNF4A的CR1211032突變對(duì)中的原始DNA序列所在的位置處于結(jié)合區(qū)域，而突變后的DNA序列位于非結(jié)合區(qū)域，因此，我們可以推斷出本文方法得到的語義特征向量semanticFeature，雖與特定的單個(gè)TF無關(guān)，但蘊(yùn)涵了所有參與訓(xùn)練的多個(gè)TF的有用信息，從而可以直觀地反映序列突變導(dǎo)致特異性變化的潛在因素。

3 結(jié) 語

基于生物序列建模是計(jì)算生物學(xué)的一個(gè)基礎(chǔ)性難題，本文針對(duì)前人工作的不足之處，提出一個(gè)深度學(xué)習(xí)與淺層學(xué)習(xí)模型結(jié)合的方法：基于深度模型的語義特征提取和基于語義特征的SVC模型構(gòu)建。通過可視化分析可得：對(duì)給定的DNA序列生成的語義特征向量，與特定單個(gè)TF無關(guān)，但蘊(yùn)涵了所有參與訓(xùn)練的TF的有用信息。與現(xiàn)有方法相比，基于語義特征的SVC模型在預(yù)測(cè)ENCODE項(xiàng)目中的ChIP-seq數(shù)據(jù)的結(jié)合特異性方面取得了最優(yōu)的結(jié)果。實(shí)驗(yàn)結(jié)果不僅證明了本文方法的有效性，而且展示了淺層機(jī)器學(xué)習(xí)模型可以通過semanticFeature快速訓(xùn)練其他任務(wù)。在未來的工作中我們將在序列模型中引入無監(jiān)督訓(xùn)練機(jī)制，以使模型達(dá)到更好的泛化能力。此外我們也將嘗試引入并行以及超參數(shù)優(yōu)化技術(shù)來加速模型具體超參數(shù)的搜索。我們工作的最終目的是希望設(shè)計(jì)出更加符合生物意義的序列模型，從而促進(jìn)計(jì)算生物學(xué)相關(guān)問題的解決。