基于機(jī)器學(xué)習(xí)的蛋白質(zhì)編碼區(qū)識(shí)別

包曉娜,何黎黎,崔景安

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

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

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

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

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

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

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

1 數(shù)據(jù)

1.1 數(shù)據(jù)的獲取

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

表1 基準(zhǔn)數(shù)據(jù)的外顯子和內(nèi)含子分布表Table 1 Exon and intron distribution table ofbenchmark data

1.2 數(shù)據(jù)的預(yù)處理

1.2.1 DNA序列的數(shù)值轉(zhuǎn)化

在實(shí)現(xiàn)外顯子和內(nèi)含子的精準(zhǔn)分類與預(yù)測(cè)前,通常需要對(duì)DNA序列進(jìn)行數(shù)值映射,即將DNA序列轉(zhuǎn)化為一個(gè)數(shù)值形式的表示[17]。本文提出了一個(gè)全新的DNA序列數(shù)值化映射方法,結(jié)合k-mer技術(shù)[14],將DNA序列中的外顯子和內(nèi)含子分別轉(zhuǎn)化為一個(gè)特征向量。下面介紹特征向量的提取過(guò)程：

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

1.2.2 DNA序列的特征提取

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

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

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

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

表3 當(dāng)k=2時(shí),部分外顯子和內(nèi)含子序列的特征向量Table 3 When k=2, eigenvectors of some exon and intron sequences

2 模型的構(gòu)建

2.1 二分類算法的選擇

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

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

(1)

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

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

AC值對(duì)比結(jié)果如表4所示,k=2時(shí),邏輯回歸取得了最大的AC值。AC值兼顧了TP、TN、FP、FN四個(gè)參數(shù)的值。AC值越大,表明分類效果越好。同時(shí)可以發(fā)現(xiàn),當(dāng)k取其他值時(shí),邏輯回歸算法相較其余四種方法也具有明顯的優(yōu)勢(shì)。因此,由特征向量與邏輯回歸組合的分類模型較準(zhǔn)確。

2.2 組合模型的確定

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

如圖7,選取五折交叉驗(yàn)證中的一次實(shí)驗(yàn)結(jié)果,畫(huà)出ROC曲線圖(k=2)?？梢悦黠@看出,組合模型最貼近面積為1的四邊形線,分類效果較好。并且,HMR195的結(jié)果具有統(tǒng)計(jì)學(xué)顯著性(邏輯回歸VS隨機(jī)森林：p=5.07×10-8;邏輯回歸VS樸素貝葉斯：p=4.99×10-16;邏輯回歸VS支持向量機(jī)：p=7.74×10-10;邏輯回歸VS k最鄰近算法：p=8.91×10-7)。BG570數(shù)據(jù)的試驗(yàn)結(jié)果也顯著(邏輯回歸VS隨機(jī)森林：p=8.05×10-16;邏輯回歸VS樸素貝葉斯：p=3.70×10-54;邏輯回歸VS支持向量機(jī)：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 五個(gè)算法模型的ROC曲線圖Fig.7 ROC curves of 5 algorithm models

3 實(shí)驗(yàn)結(jié)果

3.1 單獨(dú)數(shù)據(jù)集對(duì)比分析

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

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

3.2 合并數(shù)據(jù)集對(duì)比分析

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

(2)

(3)

(4)

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

表6 二種模型的比對(duì)結(jié)果分析表Table 6 Analysis table of comparison results of two models

4 結(jié)論及展望

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