基于機器學習的蛋白質編碼區識別

包曉娜,何黎黎,崔景安

(北京建筑大學 理學院 北京102616)

大多數真核生物的編碼區是不連續的,編碼蛋白質的序列在基因序列中被非編碼序列隔開(見圖1)。編碼的序列又稱為外顯子(Exon),攜帶著遺傳信息,能夠決定和指導生物的性狀;非編碼序列又稱為內含子(Intron)[1]。如果一個基因有n個內含子,一般總是把基因的外顯子分隔成n+1個部分。且內含子的核苷酸數量比外顯子多許多倍[1-2]。因此,外顯子和內含子的準確識別是一個具有挑戰性的研究。外顯子和內含子區分也有助于研究基因功能、基因表達、基因注釋、基因轉錄調控,對于內含子功能的研究也具有一定的輔助作用[3-4],故外顯子和內含子的分類具有重要的意義。

多年來,學者們已經提出了基因編碼區(外顯子)預測的多種方法。一般可以分為基于同源比對的方法和不依賴同源比對的方法?；谛蛄型葱缘姆椒ㄊ且袁F有的基因數據庫為標準,對待檢測DNA序列進行相似性識別,從而根據已有經驗判斷未知序列的外顯子和內含子區域。BLAST[5]、MUSCLE[6]是常見的比對工具,近年來也有諸如GeMoMa[7]的基因預測程序被提出?；谛蛄型葱缘姆椒蚀_率較高,但測序成本高、比對效率等因素制約了該項技術的發展?；诖?許多的學者將研究重點轉向不依賴比對技術的模型。數字信號處理技術在該領域發揮著關鍵的作用[8]。且數字信號處理前通常需對DNA序列進行數值映射[9]。VOSS[10]是一種廣泛使用的固定映射技術,它將DNA序列轉化為4個二進制指示符序列XA[n],XC[n],XG[n],XT[n]。核苷酸在特定堿基位置出現用1表示,未出現用0表示。Z曲線理論[11]是基于物理化學性質的映射方式。利用傳統四面體的對稱性開發,它將DNA或RNA序列映射到折疊曲線中。Z曲線表示出DNA序列攜帶的所有信息[8],可用于基因鑒定和DNA或RNA序列分析[12]、識別細菌和古細菌基因組中蛋白質編碼基因[13]等。此外,在眾多序列編碼方法中,k字符相對頻率技術(k-mer)[14]是較常見和簡便的方法。圖2展示了當k為4步長為1時的短序列的k-mer生成過程。機器學習的迅猛發展也為蛋白質編碼區的識別帶來了許多新的算法。如CNN-MGP[15]、GeneMark EP+[16]、DBN[17]。CNN-MGP[15]是用于宏基因組學基因預測的卷積神經網絡,能夠提取編碼區和非編碼區的特征。GeneMark EP+[16]是用于真核基因預測的算法和工具。深度置信網絡DBN[17]通過多層玻爾茲曼機對DNA序列進行數值轉換,用深度置信網絡模型對外顯子和內含子分類判別。盡管已經有許多的外顯子與內含子分類方法被提出,但是準確率、敏感度、特異度、AUC值等評價參數還有待提升。

圖1 真核生物外顯子與內含子交替示意圖Fig.1 Schematic diagram of exon intron alternation in eukaryotic coding region

圖2 k字符相對頻率技術提取k-mer示意圖(k=4)Fig.2 Schematic diagram of k-mer extraction by k-character relative frequency technology(k=4)

將數值映射和機器學習分類器相結合,提出了一個組合算法(具體流程見圖3)。首先,給定一個外顯子或內含子,將其通過密碼子與氨基酸的對應轉換為特定的氨基酸序列,此處的轉換不同于標準的翻譯過程。然后,利用經典的k-mer技術獲取序列的特征向量。最后,將外顯子與內含子的特征向量輸入邏輯回歸分類器中,訓練模型并識別蛋白質編碼區(外顯子)。利用真核生物基準數據集HMR195和BG570對模型進行了五折交叉驗證,AUC值分別達到了0.981 3和0.987 4。將兩個數據集合并計算時,敏感度和特異度分別為0.954 1、0.942 8。通過對比發現,新算法的識別結果明顯優于VOSSDFT、傳統的貝葉斯判別等方法。新算法識別HMR195和BG570數據集的時間為1.46 s、3.58 s,表明組合模型能夠高效又準確地鑒定真核生物的外顯子和內含子。

圖3 本文算法的框圖Fig.3 Block diagram of the algorithm in this article

1 數據

1.1 數據的獲取

本文對真核生物的DNA序列進行編碼區判別分析,實驗中用到2個基準數據集,分別是HMR195[18]和BG570[19]數據。HMR195數據由195個哺乳動物DNA序列組成,包括人類、小鼠和大鼠,共948個外顯子。BG570是指由570個脊椎動物序列組成的基因組測試數據集,共2 649個外顯子。兩個數據集可從網址http：//www.imtech.res.in/raghava/genebench中獲取?；鶞蕯祿拈L度范圍、外顯子數目和內含子數目如表1所示。為了保證對外顯子和內含子分類的全面性,將短的(長度低于20 bp)外顯子和內含子序列也加入了實驗中。

表1 基準數據的外顯子和內含子分布表Table 1 Exon and intron distribution table ofbenchmark data

1.2 數據的預處理

1.2.1 DNA序列的數值轉化

在實現外顯子和內含子的精準分類與預測前,通常需要對DNA序列進行數值映射,即將DNA序列轉化為一個數值形式的表示[17]。本文提出了一個全新的DNA序列數值化映射方法,結合k-mer技術[14],將DNA序列中的外顯子和內含子分別轉化為一個特征向量。下面介紹特征向量的提取過程：

給定一個外顯子ACAGCGACC：第1步,從第一個核苷酸A處開始,通過每次僅移動一個核苷酸,將外顯子轉化為一段特定氨基酸序列,具體為,‘ACA’對應氨基酸T,‘CAG’對應氨基酸Q,‘AGC’對應氨基酸S,‘GCG’對應氨基酸A,‘CGA’對應氨基酸R,‘GAC’對應氨基酸D,‘ACC’對應氨基酸T,由此得到一段特定氨基酸序列為TQSARDT;第2步,結合經典的k字符相對頻率技術,規定k值從1至5變化,將TQSARDT轉化為特征向量。假設k=2,則2-mer種類包括TQ、QS、SA、AR、RD、DT。特征向量由2-mer頻數構成,即(fTQ,fQS,fSA,fAR,fRD,fDT)=(1,1,1,1,1,1),其中fTQ表示TQ的頻數。

1.2.2 DNA序列的特征提取

DNA序列特征提取源于特定氨基酸序列k-mer的種類和數值頻率。具體來說,特征向量的元素(即所有的k-mer種類)是DNA序列的特征,即1.2.1節提到的TQ、QS等。通常來講,一段氨基酸序列中的k-mer種類數為20k。但是,由于特定氨基酸序列的轉化不同于生物學中標準的翻譯過程,且存在不同的密碼子對應同一種氨基酸,所以本算法的k-mer種類數遠遠少于20k,這大大節約了計算時的內存消耗。以脯氨酸P為例,如圖4,它由4個密碼子編譯CCT、CCC、CCA、CCG,由本文的轉化過程,P后的下一個氨基酸共5種,分別為亮氨酸L、組氨酸H、谷氨酰胺Q、脯氨酸P、精氨酸R,遠少于20種。因此,每個氨基酸后可能出現的氨基酸種類少于20種。最終k-mer組合種類數隨之大大減少,也就是說本文的轉化過程大大降低了特征向量的維度。表2列出了k=2時的全部95種特征向量的元素特征。

圖4 脯氨酸P后面會出現的氨基酸種類示意圖Fig.4 Schematic diagram of amino acid types that will appear after P

表2 k=2時,特征向量的95種元素Table 2 When k = 2, 95 elements of eigenvector

表3列出了部分短外顯子或內含子的特征向量(以k=2為例)詳細求解過程。規定外顯子類別為1,內含子類別為0。

表3 當k=2時,部分外顯子和內含子序列的特征向量Table 3 When k=2, eigenvectors of some exon and intron sequences

2 模型的構建

2.1 二分類算法的選擇

在完成DNA序列的數值轉換后,為了找到最適合特征向量的二分類模型,本文對五種分類器進行了嘗試和驗證,分別是隨機森林(Random forest)[20]、邏輯回歸(Logistic regression)[21]、高斯樸素貝葉斯(Gaussian naive bayes)[22]、支持向量機(SVM)[23]、k最鄰近分類算法(KNN)[24]。計算時,采用五折交叉驗證[25]。五折交叉驗證是判斷分類器性能的一種統計分析方法。它將原始數據分為5組,不重復地抽取其中4組作為訓練集,剩余的1組作為測試集,共得到5種測試結果,最終取用平均數。

為了對5種不同的算法進行有效的對比和測度,此處使用三個評價指標ROC(Receiver operating characteristic)曲線、AUC值和近似相關系數AC值。ROC曲線[26]是以假陽率(False positive rate)作為橫軸線(成本),以真陽率(True positive rate)作為縱軸線(收益),來說明在各種閾值條件下的假陽率和真陽率的關系曲線。ROC曲線與對角線的距離愈接近,表明試驗中識別編碼區與非編碼區的能力愈弱,亦即該方法的分類預測能力愈弱。為了更準確地描述算法的判別能力,通常將ROC曲線下的區域面積用AUC[26]進行定量和比較,AUC數值愈接近1,說明分類的有效性越好。近似相關系數AC[26]是一種得到普遍認可的綜合評估指標,TP(True positive)為外顯子被正確預測的個數,FP(False positive)為預測為外顯子但實際為內含子的個數,TN(True negative)為內含子被正確預測的個數,FN(False negative)為預測為內含子但實際為外顯子的個數。此外,為了檢驗結果的統計學顯著性,采用Delong檢驗[27]對ROC-AUC進行成對比較,p<0.05被認為具有統計學意義。

(1)

具體實驗結果如圖5、圖6和表4所示。圖5中,k=2時,在HMR195數據集,邏輯回歸的AUC平均數分別為0.981 3,明顯高于其他模型的結果。如圖6,BG570數據集也得到類似的結果,邏輯回歸算法在所有k值優于其他算法。

圖5 HMR195數據集中,五種算法AUC值的熱圖Fig.5 In HMR195 data set, heatmaps of AUC values of five algorithms

AC值對比結果如表4所示,k=2時,邏輯回歸取得了最大的AC值。AC值兼顧了TP、TN、FP、FN四個參數的值。AC值越大,表明分類效果越好。同時可以發現,當k取其他值時,邏輯回歸算法相較其余四種方法也具有明顯的優勢。因此,由特征向量與邏輯回歸組合的分類模型較準確。

2.2 組合模型的確定

最終,組合模型確定為特征向量與邏輯回歸分類器的結合。首先,將DNA序列轉化為特定的氨基酸序列;其次,由特定氨基酸序列得到特征向量。最后,將特征向量放入邏輯回歸分類器中,獲得外顯子和內含子的預測結果。

如圖7,選取五折交叉驗證中的一次實驗結果,畫出ROC曲線圖(k=2)?？梢悦黠@看出,組合模型最貼近面積為1的四邊形線,分類效果較好。并且,HMR195的結果具有統計學顯著性(邏輯回歸VS隨機森林：p=5.07×10-8;邏輯回歸VS樸素貝葉斯：p=4.99×10-16;邏輯回歸VS支持向量機：p=7.74×10-10;邏輯回歸VS k最鄰近算法：p=8.91×10-7)。BG570數據的試驗結果也顯著(邏輯回歸VS隨機森林：p=8.05×10-16;邏輯回歸VS樸素貝葉斯：p=3.70×10-54;邏輯回歸VS支持向量機：p=4.67×10-9;邏輯回歸VS k最鄰近算法：p=1.24×10-7)。

表4 k從1至5,5種算法的AC平均值Table 4 K from 1 to 5, mean AC value of 5 algorithms

圖7 五個算法模型的ROC曲線圖Fig.7 ROC curves of 5 algorithm models

3 實驗結果

3.1 單獨數據集對比分析

為了說明本文新方法與其余方法的優劣,將其與經典的VOSSDFT[10,28]、EIIPDFT[28-29]、SPDFT[28,30]和Code13-Marple[28]進行了比較。VOSSDFT、EIIPDFT、SPDFT均是基于離散傅里葉變換的技術(Discrete Fourier Transform, DFT)來區分真核生物外顯子和內含子[10,29-30]。Code13-Marple是一種基于自回歸譜分析和小波變換的集成算法。由表5,以HMR195為例,新方法(k=2)的AUC值達到了0.981 3,比其余四種方法分別高出了0.418 7、0.470 0、0.385 1、0.263 4;在BG570數據集上,AUC和AC值也遠遠超過其余四種模型中的最大值。新算法明顯優于其他三種傳統的基于DFT的方法和Code13-Marple。

表5 組合模型與其他方法的比較Table 5 Comparison of eigenvector method with other methods

3.2 合并數據集對比分析

為驗證算法在較大數據集上的分類效果,將HMR195和BG570兩組數據合并得到合并數據集,共3 597個外顯子、4 354個內含子。此外,為了更加全面的評估組合模型的性能,增加了準確率、敏感度、特異度以及運行時間這四個對比維度,并與經典的貝葉斯判別法[31]進行比較。貝葉斯判別法是進行判別分析的一種多元統計分析方法。合并數據集后,k值取3時本文算法得到最好的預測結果。

(2)

(3)

(4)

表6是兩種方法的對比分析表,其中準確率acc[26]為全部序列中被正確預測的序列的比例;敏感度Sn[26]為所有實際外顯子中被正確預測為外顯子的比例;特異度Sp[26]為所有真實的內含子被正確預測為內含子的比例。在合并后的較大數據集上,組合模型的敏感度Sn為0.954 1遠遠大于貝葉斯判別法的0.787 2。在運行時間方面,組合模型只需要8.91 s,而貝葉斯判別法需要27.28 s。因此,本文方法不僅適用于小數據集,在較大數據集上同樣表現優異,并且運行速度快于貝葉斯判別法。本文組合模型以及貝葉斯判別法的計算基于處理器為Intel(R) Core(TM) i7-8550U CPU@1.80 GHz和16.0 GB RAM的設備,使用Python3.8編程獲得。

表6 二種模型的比對結果分析表Table 6 Analysis table of comparison results of two models

4 結論及展望

本研究提出了一個基于特征向量的數值映射方法,之后結合邏輯回歸算法,對基因外顯子和內含子實現了精確的分類。將組合模型運用到編碼區識別,給出了一個全新的研究視角。為了證明組合模型的可行性,利用HMR195和BG570兩個真核生物數據集,將其與現有的成熟方法進行了對比(見表5和表6),均證明了它的有效性。此外,為證實模型在更大數據集上的效果,本文新收集了462條人類DNA序列[32]進行試驗,共包含2 843個外顯子,2 381個內含子。全部數據可從網址https：//www.fruitfly.org/sequence/human-datasets.html獲取。當全部數據共同訓練時,共6 440個外顯子,6 735個內含子。本文方法實驗結果：acc、Sn、Sp、AC、AUC的值分別為0.957 7、0.966 6、0.949 0、0.915 5、0.989 4(k=2)。當擴大數據集后,組合模型對于外顯子和內含子依然能起到很好的識別效果。其次,1.2.2節中特征向量的提取過程充分利用了密碼子的簡并性,降低了特征向量的維度。然而本文還未將外顯子和內含子的結構信息作為特征的重要因素,之后的研究中會考慮加入結構信息,從而進一步提升模型的性能。并且,本文后續研究仍將擴大樣本量,嘗試更多更全面物種的蛋白質編碼區分類,爭取構建快速便捷的外顯子與內含子識別工具。