基于卷積神經網絡的細菌轉錄終止子預測

金 冬，張 萌，賈藏芝

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

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

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

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

1 材料和方法

1.1 數據收集和預處理

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

表1 不同物種的數據集

1.2 特征提取

特征提取在開發基于機器學習算法的計算模型中起著非常重要的作用。本文從序列中提取了兩類特征：一個是k-mer，另一個是EIIP[17]。

1.2.1k-mer核苷酸組成

給出一個DNA序列D，它的直觀表達式是[18]：

D=R1R2R3R4R5R6R7…RL

(1)

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

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

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

(2)

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

1.2.2 EIIP

EIIP作為特征已被廣泛的用來預測基因序列[17]?；贓IIP的識別方法廣泛應用于基因結構識別的關鍵部分，如F56F11.4基因的預測[24]、囊性纖維化基因的預測和和增強子的識別[25]等。

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

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

(3)

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

1.3 特征選擇

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

1.3.1 MRMD

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

(4)

其中，

(5)

(6)

(7)

(8)

(9)

從而，第i個特征的最大相關值MR定義為：

(10)

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

(11)

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

maxMDi=EDi(1≤i≤M)

(12)

根據以上結果，第i個特征MRMD值定義為max(MRi+MDi)，根據此值的大小，對特征進行排序。數值越大，表明此特征與目標標簽的相關性越強[26]。

1.3.2F-score

第j個特征的F-score定義為：

(13)

1.3.3 二項分布

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

(14)

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

(15)

(16)

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

最后，根據以下公式計算訓練數據集中的第j個k-mer核苷酸的概率:

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

(17)

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

2 卷積神經網絡

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

圖1 神經網絡模型的架構

3 模型訓練和性能評估

3.1 參數選擇

采用貝葉斯對卷積神經網絡的神經元個數(a)、批次(b)、dropout(c)、學習率(d)、激活函數(e)以及全連接層數(f)這六種參數進行優化。除上述參數外，所涉及到的其它參數均按照scikit-learn庫的默認值。其中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]中取值。根據貝葉斯優化方法對參數組合進行五十次尋優，耗時1小時11分鐘，優化過程以及最佳參數結果(見圖2)。貝葉斯優化器建立了搜索空間的替代模型，并在此維度內進行搜索，而不是在實際搜索空間內進行搜索。優化參數的二維圖，最終選取損失值最小(0.122 5)的參數組合。

圖2 參數選擇的結果

3.2 性能評估

為了評估轉錄終止預測模型的性能，使用準確性(Acc)、靈敏度(Sn)、特異性(Sp)和馬修相關系數(MCC)作為五折交叉驗證和獨立數據集測試的評估標準。

(18)

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

4 結果分析

4.1 選擇最優的特征子集

為尋找使CNN分類器達到最優性能的k-mer(k=4,5,6,7)特征，利用五折交叉驗證對每類特征進行測試。如圖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的特征維數為4096，高維度特征可能包含冗余信息，導致過擬合。因此，使用了MRMD、F-Score和二項分布這三種常用的特征選擇方法來尋找最優特征子集。

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

第1步，根據F-score值、MRMD值和二項分布值對6-mer向量中的4 096個元素進行排序；

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

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

4.2 模型對比

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

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

4.3 獨立測試集表現

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

圖4 在獨立測試中模型與iTerm-PseKNC的準確率比較

4.4 特征可視化

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

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

5 結 論

1)在這項研究中，提出了一種新的計算模型TermCNN可以快速準確地識別轉錄終止子；

2)將代表性的6-mer特征子集和EIIP作為輸入參數，利用CNN對模型進行訓練和優化；

3)五折交叉驗證和多個獨立測試結果證明了模型的競爭力，其性能結果明顯優于其他算法和現有計算工具iTerm-PseKNC，但是在靈敏度方面比iterb-PPse稍低。