基于卷積神經(jīng)網(wǎng)絡(luò)的細(xì)菌轉(zhuǎn)錄終止子預(yù)測

金 冬，張 萌，賈藏芝

(大連海事大學(xué) 理學(xué)院,遼寧 大連 116026)

在遺傳學(xué)中，轉(zhuǎn)錄終止子通常位于poly(A)位點(diǎn)的下游，提供終止轉(zhuǎn)錄的信號。它通過對新合成的轉(zhuǎn)錄本RNA提供信號來介導(dǎo)轉(zhuǎn)錄終止[1-2]。一般來說，原核生物的轉(zhuǎn)錄終止子可分為兩類：一類依賴于ρ(Rho)因子才能實(shí)現(xiàn)終止作用，記作Rho-dependent；另一類則不依賴ρ因子便能實(shí)現(xiàn)終止作用，記為Rho-independent。ρ因子是一種解旋酶，可以破壞mRNA-DNA-RNA聚合酶的轉(zhuǎn)錄復(fù)合體。通常在細(xì)菌和噬菌體中發(fā)現(xiàn)Rho-dependent轉(zhuǎn)錄終止子[3-4]。Rho-independent 的終止位點(diǎn)位于翻譯終止密碼子的下游，由mRNA上一個(gè)無結(jié)構(gòu)的、富含GC堿基對的序列組成[5]。Rho-independent[6-7]因子包含7-20個(gè)GC-rich區(qū)域，后跟一個(gè)短poly-T或T-stretch，在延長轉(zhuǎn)錄本上形成自退火發(fā)夾結(jié)構(gòu)，轉(zhuǎn)錄中的RNA聚合酶遇到發(fā)夾結(jié)構(gòu)將會暫停前進(jìn)。它通過破壞mRNA-DNA-RNA聚合酶三元復(fù)合物，最終使轉(zhuǎn)錄終止。

在傳統(tǒng)的實(shí)驗(yàn)中，轉(zhuǎn)錄終止子是否存在通常是通過測定mRNA的長度來確定的，而這種方法往往無法精確的識別終止位點(diǎn)[8]。因此，許多預(yù)測轉(zhuǎn)錄終止子的方法被開發(fā)出來。近年來，Yada等[9]利用隱馬爾可夫模型預(yù)測了大腸桿菌基因的轉(zhuǎn)錄終止子。Ermolaeva[10]和Unniraman等[11]分別使用TransTerm算法和GeSTer算法預(yù)測了轉(zhuǎn)錄終止子。2001年，Lesnik等[12]提出了一種基于熱力學(xué)評分系統(tǒng)預(yù)測大腸桿菌K-12基因組終止子的方法。隨著機(jī)器學(xué)習(xí)技術(shù)的發(fā)展，許多分類任務(wù)得到了解決。Feng等[8]提取偽k-tuple核苷酸組成特征，并通過二項(xiàng)分布進(jìn)行特征選擇。隨后將選擇的特征與支持向量機(jī)(SVM)相結(jié)合，構(gòu)建了一個(gè)名為iTerm-PseKNC[8]的計(jì)算方法來預(yù)測轉(zhuǎn)錄終止子。最近，F(xiàn)an等[13]采用k核苷酸的位置信息、核苷酸的含量、核苷酸的47種物理化學(xué)性質(zhì)作為特征向量，并結(jié)合XGBoost分類算法構(gòu)建了預(yù)測模型iterb-PPse,取得了相當(dāng)不錯(cuò)的效果。值得注意的是，原有的預(yù)測方法都是采用傳統(tǒng)的機(jī)器學(xué)習(xí)方法作為分類算法。近幾年，深度學(xué)習(xí)種的卷積神經(jīng)網(wǎng)絡(luò)框架在生物信息領(lǐng)域得到了廣泛應(yīng)用，并且取得了令人滿意的分類性能[14-15]。因此，我們嘗試將卷積神經(jīng)網(wǎng)絡(luò)應(yīng)用于細(xì)菌終止子的預(yù)測。

在本研究中，根據(jù)Feng等[8]的工作，引入了一種新的轉(zhuǎn)錄終止子預(yù)測模型,稱為TermCNN。首先從大腸桿菌DNA序列中提取k-mer(k= 4,5,6,7)核苷酸組成特征作為CNN的輸入向量。在五折交叉驗(yàn)證中，挑選出準(zhǔn)確率最高的6-mer核苷酸組成特征。然后采用最大相關(guān)-最大距離(MRMD)、二項(xiàng)分布和F-score這三種特征選擇方法來尋找6-mer特征的最優(yōu)特征子集，以減少無用信息和節(jié)省運(yùn)行時(shí)間。最后將選擇出的最優(yōu)特征子集與電子-離子相互作用偽電位(EIIP)特征相結(jié)合，輸入到CNN進(jìn)行訓(xùn)練，構(gòu)建高精度模型。五折交叉驗(yàn)證以及五個(gè)獨(dú)立測試數(shù)據(jù)集的實(shí)驗(yàn)結(jié)果一致顯示了本文提出的預(yù)測模型TermCNN的有效性，特別是用于區(qū)分不同種類的終止子。

1 材料和方法

1.1 數(shù)據(jù)收集和預(yù)處理

一個(gè)客觀的基準(zhǔn)數(shù)據(jù)集是建立終止子預(yù)測模型的基礎(chǔ)。從RegulonDB[16]中收集大腸桿菌的終止子，去冗余后得到286個(gè)Rho-independent 終止子和19個(gè)Rho-dependent終止子[8]。與之前的數(shù)據(jù)集相比，RegulonDB新增了25個(gè)轉(zhuǎn)錄終止子。將新發(fā)現(xiàn)的25個(gè)轉(zhuǎn)錄終止子視為一個(gè)獨(dú)立的測試集，命名為E_Ter_25。對于訓(xùn)練數(shù)據(jù)集，采用了與Feng[8]相同的數(shù)據(jù)集，包含280個(gè)終止子和560個(gè)非終止子，便于評估和比較不同預(yù)測器的性能。對于獨(dú)立測試，F(xiàn)eng[8]使用了兩個(gè)終止子獨(dú)立測試集，分別是E_Ter_147和B_Ter_425。從Fan等[13]的工作中選取兩個(gè)均為負(fù)樣本構(gòu)成的獨(dú)立測試集，樣本是分別從大腸桿菌和枯草芽孢桿菌的上游截取的，記為E_Nonter_159和B_Nonter_122。在縮寫中，E表示來自大腸桿菌的序列，B表示來自枯草芽孢桿菌的序列，數(shù)字表示每個(gè)數(shù)據(jù)集中的樣本數(shù)量(見表1)。

表1 不同物種的數(shù)據(jù)集

1.2 特征提取

特征提取在開發(fā)基于機(jī)器學(xué)習(xí)算法的計(jì)算模型中起著非常重要的作用。本文從序列中提取了兩類特征：一個(gè)是k-mer，另一個(gè)是EIIP[17]。

1.2.1k-mer核苷酸組成

給出一個(gè)DNA序列D，它的直觀表達(dá)式是[18]：

D=R1R2R3R4R5R6R7…RL

(1)

其中Ri表示在DNA序列中第i個(gè)位置的核苷酸。

k-mer核苷酸是將DNA序列轉(zhuǎn)化為數(shù)字向量的一種簡單而常用的方法，這一方法具有重要的生物學(xué)意義，在DNA調(diào)控元件識別中已得到了廣泛的應(yīng)用[19-23]。k-mer可以將任何DNA序列表示為4k維的向量如下：

R=[φ1φ2…φu…φ4k]

(2)

其中φu(u=1,2…,4k)為沿著序列第u個(gè)k-mer的頻率。在本工作中，k=1、2、3、4、5、6、7，并與EIIP相結(jié)合進(jìn)行測試，尋找最優(yōu)的特征集。

1.2.2 EIIP

EIIP作為特征已被廣泛的用來預(yù)測基因序列[17]?；贓IIP的識別方法廣泛應(yīng)用于基因結(jié)構(gòu)識別的關(guān)鍵部分，如F56F11.4基因的預(yù)測[24]、囊性纖維化基因的預(yù)測和和增強(qiáng)子的識別[25]等。

四個(gè)核苷酸的EIIP值分別為，A: 0.126，G:0.0806，C: 0.134，T:0.133。計(jì)算每條序列中A、T、G、C的平均EIIP值，構(gòu)造特征向量為[25]:

P=[EIIPAAA·fAAA,EIIPAAC·fAAC,…,EIIPTTT·fTTT]

(3)

其中fXYZ為任意三核苷酸XYZ的頻率，EIIPXYZ是三核苷酸XYZ的EIIP值之和，X,Y,Z∈{A,C,G,T}。

1.3 特征選擇

特征選擇方法可以降低特征向量的維數(shù)，為訓(xùn)練分類器找到最優(yōu)特征子集。近幾年，最大相關(guān)-最大距離(MRMD)[26]、F-score[27]和二項(xiàng)分布(BD)[28]，方法在改善預(yù)測器性能上具有顯著成效，已廣泛應(yīng)用于生物信息學(xué)領(lǐng)域。

1.3.1 MRMD

MRMD利用皮爾遜相關(guān)系數(shù)計(jì)算特征子集與目標(biāo)類的相關(guān)性，并使用歐氏距離函數(shù)計(jì)算特征子集的冗余度，相關(guān)性與距離的和最大的特征被選擇到最終的特征子集中。首先定義兩個(gè)向量的相關(guān)系數(shù)如下：

(4)

其中，

(5)

(6)

(7)

(8)

(9)

從而，第i個(gè)特征的最大相關(guān)值MR定義為：

(10)

本文中，兩個(gè)特征的距離采用歐式距離，定義如下：

(11)

最大距離就是取所有歐氏距離中的最大值，記為：

maxMDi=EDi(1≤i≤M)

(12)

根據(jù)以上結(jié)果，第i個(gè)特征MRMD值定義為max(MRi+MDi)，根據(jù)此值的大小，對特征進(jìn)行排序。數(shù)值越大，表明此特征與目標(biāo)標(biāo)簽的相關(guān)性越強(qiáng)[26]。

1.3.2F-score

第j個(gè)特征的F-score定義為：

(13)

1.3.3 二項(xiàng)分布

Feng[8]和Su[29]采用基于二項(xiàng)分布的技術(shù)，通過SVM分類器進(jìn)行五折交叉驗(yàn)證的性能結(jié)果對特征進(jìn)行選擇。這里，先驗(yàn)概率qi定義為k-mer核苷酸的頻率，如下所示:

(14)

其中mi(i=1,-1)分別表示正、負(fù)訓(xùn)練數(shù)據(jù)集(即終止子和非終止子數(shù)據(jù)集)中的k-mer片段總數(shù)。M為全部訓(xùn)練數(shù)據(jù)集中k-mer片段的總數(shù)。因此，第j個(gè)k-mer核苷酸(j=1,2,…, 4k)在正樣本和負(fù)樣本中的概率可以定義為：

(15)

(16)

其中Nj表示終止子和非終止子訓(xùn)練數(shù)據(jù)集中第j個(gè)k-mer核苷酸的總數(shù)。n1,j和n-1,j分別表示正、負(fù)訓(xùn)練數(shù)據(jù)集中第j個(gè)k-mer核苷酸的總數(shù)。

最后，根據(jù)以下公式計(jì)算訓(xùn)練數(shù)據(jù)集中的第j個(gè)k-mer核苷酸的概率:

Pj=min(p(n1,j),p(n-1,j)

(17)

所有的k-mer核苷酸可以根據(jù)概率的大小進(jìn)行排序，也就是說，Pj越小，相應(yīng)的k-mer核苷酸對分類效果越有效。

2 卷積神經(jīng)網(wǎng)絡(luò)

CNN已被廣泛應(yīng)用于各種分類任務(wù)中，其在圖像識別、圖像檢測、語音識別等方面表現(xiàn)出良好的性能。隨著深度學(xué)習(xí)的深入研究[30]，CNN還用于預(yù)測啟動(dòng)子[31]、蛋白質(zhì)泛素化位點(diǎn)[32]、蛋白質(zhì)翻譯后修飾位點(diǎn)的capsule網(wǎng)絡(luò)[33]、RNA假尿苷位點(diǎn)[34]。在本研究中，借助Keras工具，使用CNN模型識別轉(zhuǎn)錄終止子。TermCNN由兩個(gè)卷積層、兩個(gè)池化層和連接層組成(見圖1)。轉(zhuǎn)錄終止子包含更多的GC堿基對，因此使用了一個(gè)平均池化層，池化大小為3×3，這適合于獲取序列的GC含量。還使用dropout來防止模型的過擬合。對于隨機(jī)梯度下降法，選擇了Adam優(yōu)化算法。整個(gè)程序在Python 3.6中使用，實(shí)驗(yàn)環(huán)境為：主機(jī)CPU型號為AMD Ryzen 74 800 H with Radeon Graphics，主頻為2.90 GHz，物理內(nèi)存為16 GB，操作系統(tǒng)為64位Windows10，深度學(xué)習(xí)框架為TensorFlow 2.0.0。

圖1 神經(jīng)網(wǎng)絡(luò)模型的架構(gòu)

3 模型訓(xùn)練和性能評估

3.1 參數(shù)選擇

采用貝葉斯對卷積神經(jīng)網(wǎng)絡(luò)的神經(jīng)元個(gè)數(shù)(a)、批次(b)、dropout(c)、學(xué)習(xí)率(d)、激活函數(shù)(e)以及全連接層數(shù)(f)這六種參數(shù)進(jìn)行優(yōu)化。除上述參數(shù)外，所涉及到的其它參數(shù)均按照scikit-learn庫的默認(rèn)值。其中a在[8,64]中取值，b在集合[8,128]中取值，c在集合{0.1,0.3,0.5,0.7}中取值，d在集合[0.000 1,0.01]中取值，e的選取有兩種情況relu和sigmoid，f在集合[1,10]中取值。根據(jù)貝葉斯優(yōu)化方法對參數(shù)組合進(jìn)行五十次尋優(yōu)，耗時(shí)1小時(shí)11分鐘，優(yōu)化過程以及最佳參數(shù)結(jié)果(見圖2)。貝葉斯優(yōu)化器建立了搜索空間的替代模型，并在此維度內(nèi)進(jìn)行搜索，而不是在實(shí)際搜索空間內(nèi)進(jìn)行搜索。優(yōu)化參數(shù)的二維圖，最終選取損失值最小(0.122 5)的參數(shù)組合。

圖2 參數(shù)選擇的結(jié)果

3.2 性能評估

為了評估轉(zhuǎn)錄終止預(yù)測模型的性能，使用準(zhǔn)確性(Acc)、靈敏度(Sn)、特異性(Sp)和馬修相關(guān)系數(shù)(MCC)作為五折交叉驗(yàn)證和獨(dú)立數(shù)據(jù)集測試的評估標(biāo)準(zhǔn)。

(18)

其中TP、TN、FP和FN分別代表真陽性、真陰性、假陽性和假陰性的數(shù)量。

4 結(jié)果分析

4.1 選擇最優(yōu)的特征子集

為尋找使CNN分類器達(dá)到最優(yōu)性能的k-mer(k=4,5,6,7)特征，利用五折交叉驗(yàn)證對每類特征進(jìn)行測試。如圖3所示，6-mer與CNN的整合得到了最好的MCC和Acc。6-mer的MCC值為0.942，比4-mer高0.101，比5-mer高0.004，比7-mer高0.035?？紤]到6-mer的特征維數(shù)為4096，高維度特征可能包含冗余信息，導(dǎo)致過擬合。因此，使用了MRMD、F-Score和二項(xiàng)分布這三種常用的特征選擇方法來尋找最優(yōu)特征子集。

圖3 四種具有不同數(shù)量特征的模型性能

第1步，根據(jù)F-score值、MRMD值和二項(xiàng)分布值對6-mer向量中的4 096個(gè)元素進(jìn)行排序；

第2步，設(shè)k=30作為初值。需要指出的是，為特征子集選擇的維數(shù)是某個(gè)數(shù)字k的平方，特征向量可以轉(zhuǎn)換成一個(gè)方陣作為輸入。因此，選取排名靠前的k2-64元素與EIIP特征結(jié)合形成長度為k的方0陣，然后將一維特征向量轉(zhuǎn)換為二維方陣作為CNN的輸入。

以步長為5，在特征方陣長度為k+5，通過五折交叉驗(yàn)證尋找準(zhǔn)確率最高的特征子集。最后在精度最高的特征維數(shù)周圍使用步長為1篩選出最優(yōu)特征子集。并比較最優(yōu)特征集和無特征選擇的特征集的結(jié)果。對于F-score、二項(xiàng)分布和MRMD，k值分別為63、64和51時(shí)的準(zhǔn)確率最高。相對時(shí)間成本和準(zhǔn)確性，將MRMD方法選擇的6-mer特征向量中前2 537個(gè)元素與64個(gè)EIIP特征相結(jié)合，Acc為97.62%，Sn為92.86%，Sp為100%，MCC為0.947。這表明所建立的模型TermCNN具有良好的識別轉(zhuǎn)錄終止子的能力。

4.2 模型對比

為了證明使用深度學(xué)習(xí)識別轉(zhuǎn)錄終止子的優(yōu)越性，將CNN與決策樹、多層感知器、邏輯回歸、樸素貝葉斯、基于SVM的iTerm-PseKNC、iterb-PPSE和CNN+LSTM進(jìn)行了比較。結(jié)果如表2所示?？梢姡跍\層機(jī)器學(xué)習(xí)中，基于SVM的iTerm-PseKNC達(dá)到了最好的Acc(95.71%)，MCC(0.888)，CNN實(shí)現(xiàn)了較好的性能，達(dá)到97.98%的Acc和0.955的MCC，但是基于XGBoost的iterb-PPse給出了最好的結(jié)果，Acc為99.88%，MCC為0.999。TermCNN比iterb-PPse稍微遜色的原因有兩個(gè)：1)提取的特征過于單一。本文僅僅考慮的終止子序列的6-mer特征，沒有考慮位置及核苷酸的物理化學(xué)性質(zhì)；2)數(shù)量的樣本數(shù)量較少，不能夠體現(xiàn)CNN的優(yōu)越性。隨著越來越多終止子序列的發(fā)現(xiàn)，也將會繼續(xù)優(yōu)化我們的模型。

表2 不同分類器在五折交叉驗(yàn)證中識別終止子的比較

4.3 獨(dú)立測試集表現(xiàn)

為了更好地評價(jià)模型的泛化能力，進(jìn)一步測試了五個(gè)獨(dú)立的數(shù)據(jù)集E_Ter_147、B_Ter_425，E_Ter_25，E_Nonter_159和B_Nonter_122。對于E_Ter_147，TermCNN正確預(yù)測了147個(gè)終止子，iTerm-PseKNC也正確預(yù)測了147個(gè)終止子。對于B_Ter_425，TermCNN型正確預(yù)測了417個(gè)終止子(98.12%)，而iTerm-PseKNC僅正確預(yù)測了372個(gè)終止子(87.53%)。對于新的獨(dú)立測試集，TermCNN正確預(yù)測了所有25個(gè)終止子(100%)，而iTerm-PseKNC正確預(yù)測了24個(gè)終止子(96%)，如圖4所示。為了多方面檢驗(yàn)所建立模型的有效性，從iterb-PPse中選取兩個(gè)負(fù)樣本數(shù)據(jù)集E_Nonter_159和B_Nonter_122。對于E_Nonter_159，TermCNN預(yù)測了158個(gè)非轉(zhuǎn)錄終止子(99.37%)。對于B_Nonter_122，TermCNN預(yù)測了121個(gè)非轉(zhuǎn)錄終止子(99.18%)。相比于iterb-PPse, TermCNN預(yù)測對的數(shù)目少一個(gè)。比較遺憾的是，由于iTerm-PseKNC提供的網(wǎng)絡(luò)服務(wù)器不能正常使用，因此無法和它進(jìn)行比較。

圖4 在獨(dú)立測試中模型與iTerm-PseKNC的準(zhǔn)確率比較

4.4 特征可視化

為了更加直觀的可以看到特征的有效性，通過采用t分布隨機(jī)鄰居嵌入(t-SNE)進(jìn)行特征可視化。圖5中每個(gè)點(diǎn)代表一個(gè)樣本，藍(lán)色點(diǎn)表示轉(zhuǎn)錄終止子位點(diǎn)，紅色點(diǎn)表示非轉(zhuǎn)錄終止子位點(diǎn)。一開始可以清晰的看到只用原始特征表示的兩類點(diǎn)很難分開，后經(jīng)過神經(jīng)網(wǎng)絡(luò)層層訓(xùn)練，在全連接層的輸出向量可以比較明顯的劃分兩類。因此，顯示CNN處理轉(zhuǎn)錄終止子數(shù)據(jù)很有效。

圖5 t-SNE可視化特征表示

5 結(jié) 論

1)在這項(xiàng)研究中，提出了一種新的計(jì)算模型TermCNN可以快速準(zhǔn)確地識別轉(zhuǎn)錄終止子；

2)將代表性的6-mer特征子集和EIIP作為輸入?yún)?shù)，利用CNN對模型進(jìn)行訓(xùn)練和優(yōu)化；

3)五折交叉驗(yàn)證和多個(gè)獨(dú)立測試結(jié)果證明了模型的競爭力，其性能結(jié)果明顯優(yōu)于其他算法和現(xiàn)有計(jì)算工具iTerm-PseKNC，但是在靈敏度方面比iterb-PPse稍低。