呂成偉， 樊永顯

(桂林電子科技大學(xué) 計(jì)算機(jī)與信息安全學(xué)院，廣西 桂林 541004)

自經(jīng)過(guò)修飾的RNA核糖核酸第一次發(fā)現(xiàn)以來(lái)，被人類(lèi)已知的RNA修飾類(lèi)型已經(jīng)達(dá)到了約150種[1]。研究表明，RNA修飾是基因調(diào)控的關(guān)鍵組成部分[2]，其參與了轉(zhuǎn)錄后的各種生物過(guò)程，如蛋白質(zhì)翻譯和定位、mRNA剪接等，并發(fā)揮著重要作用[3]。但是，RNA修飾在其他方面的功能對(duì)人們來(lái)說(shuō)仍然是未知的。因此，預(yù)測(cè)RNA修飾位點(diǎn)對(duì)于理解它們的分子機(jī)制和功能起著至關(guān)重要的作用。

新一代測(cè)序技術(shù)的出現(xiàn)為在全基因組范圍內(nèi)研究RNA修飾提供了契機(jī)。如N1-甲基甘氨酸(m1A)、N6-甲基甘氨酸(m6A)和5-甲基胞嘧啶(m5C)圖譜可用于人類(lèi)轉(zhuǎn)錄組。雖然這些基于生物實(shí)驗(yàn)技術(shù)的高通量測(cè)序方法在理解生物功能和RNA修飾方面起到了積極推動(dòng)的作用，但這些方法有很大局限性，其實(shí)驗(yàn)成本高、耗時(shí)長(zhǎng)。為了解決該問(wèn)題，一些用于識(shí)別RNA修飾位點(diǎn)的基于高分辨率實(shí)驗(yàn)數(shù)據(jù)的計(jì)算方法被提出。針對(duì)m6A修飾位點(diǎn)的識(shí)別問(wèn)題，Chen等[4]提出了基于序列的iRNAMethy方法，使用了偽二核苷酸組分(pseudo dinucleotide composition，簡(jiǎn)稱(chēng)PseDNC)編碼方式，在特征提取方式上取得了一些突破。Chen等[5]在編碼方式上進(jìn)行了創(chuàng)新，提出了m6Apred方法，在原有序列信息的基礎(chǔ)上計(jì)算出核苷酸的頻率信息，且加入了其化學(xué)分類(lèi)特征，進(jìn)一步提高了m6A修飾位點(diǎn)的識(shí)別準(zhǔn)確率。Zhang等[6]提出了一種新的提取特征的方法m6A-HPCS，該方法的思想與Chen等[5]提出的m6Apred基本一致，從23種核苷酸物理化學(xué)性質(zhì)中出尋找一個(gè)最優(yōu)子集，結(jié)合自協(xié)方差和互協(xié)方差變換提取序列特征。

上述幾種方法為基于序列的RNA修飾位點(diǎn)的預(yù)測(cè)開(kāi)辟了道路，并取得了一定的成功。但這幾種方法對(duì)m6A和m5C修飾位點(diǎn)的識(shí)別準(zhǔn)確率不夠理想，仍有較大的提升空間。針對(duì)上述問(wèn)題，在采用PseKNC的編碼方式對(duì)樣本序列進(jìn)行編碼的基礎(chǔ)上，添加了位置特異性單核苷酸及二核苷酸偏好特征，并基于XGBoost集成算法構(gòu)建了預(yù)測(cè)模型。實(shí)驗(yàn)證明，該模型的魯棒性好，且針對(duì)上述提出的3種RNA修飾位點(diǎn)的預(yù)測(cè)都取得了較高的準(zhǔn)確率。

1 材料與方法

為了開(kāi)發(fā)一種新的預(yù)測(cè)方法，文獻(xiàn)[7-11]遵循了Chou[12]提出的一些原則，并明確以下4個(gè)步驟：

1)構(gòu)造一個(gè)高質(zhì)量的基準(zhǔn)數(shù)據(jù)集用來(lái)訓(xùn)練和測(cè)試預(yù)測(cè)模型；

2)對(duì)生物序列進(jìn)行編碼，使其能夠被預(yù)測(cè)模型識(shí)別；

3)選擇或者自主研發(fā)一種魯棒的算法來(lái)建立預(yù)測(cè)模型；

4)進(jìn)行交叉驗(yàn)證，并客觀地評(píng)價(jià)預(yù)測(cè)模型的好壞。

1.1 基準(zhǔn)數(shù)據(jù)集

本研究所使用的數(shù)據(jù)集包括m1A、m6A、m5C 3種RNA序列[13]，這3種RNA序列分別包含6 366、1 130和120個(gè)正樣本。為了平衡正負(fù)樣本，從對(duì)應(yīng)的負(fù)樣本中隨機(jī)選取了6 366、1 130和120個(gè)樣本分別作為m1A、m6A和m5C的負(fù)樣本。

1.2 構(gòu)建RNA序列樣本

生物信息學(xué)面臨的一個(gè)極具挑戰(zhàn)的問(wèn)題是：如何對(duì)生物序列進(jìn)行有效編碼，使得其編碼后的序列盡可能包含序列模式特征，并能被現(xiàn)有的機(jī)器學(xué)習(xí)模型直接識(shí)別。幾乎所有的機(jī)器學(xué)習(xí)算法都是為了處理向量而開(kāi)發(fā)的，并不能直接識(shí)別序列樣本。常用的獨(dú)熱編碼雖然能解決不能被機(jī)器學(xué)習(xí)模型識(shí)別的問(wèn)題，但這種編碼方式會(huì)丟失許多重要的序列模式特征，從而導(dǎo)致最終的預(yù)測(cè)模型預(yù)測(cè)能力低下。為了解決蛋白質(zhì)序列的問(wèn)題，研究人員開(kāi)發(fā)了偽氨基酸組成物(pseudo amino acid composition，簡(jiǎn)稱(chēng)PseAAC)[14-18]。自PseAAC引入以來(lái)，它的概念幾乎滲透到計(jì)算蛋白質(zhì)組學(xué)的每個(gè)領(lǐng)域。在PseAAC概念的啟發(fā)和巨大成功的鼓舞下，PseKNC[28]被提出來(lái)，并應(yīng)用到基因組分析的各個(gè)領(lǐng)域。

假設(shè)RNA序列由L個(gè)核苷酸殘基組成，

S={R1,R2,…,RL}，

(1)

其中R1、R2分別為RNA序列的第1、2個(gè)核苷酸殘基。

本研究用PseKNC對(duì)RNA序列進(jìn)行編碼，最終得到具有4k個(gè)分量的向量，即

(2)

對(duì)于式(2)，當(dāng)k=2時(shí)，表示RNA序列由二核苷酸配置而成，于是有

D=[f(AA)f(AC)f(AG)…f(UU)]T=

(3)

D=[f(AAA)f(AAC)f(AAG)f(UUU)]Τ=

(4)

1.3 XGBoost

XGBoost是一種基于GBDT(gradient boosting decision tree，簡(jiǎn)稱(chēng)GBDT)梯度下降框架的集成學(xué)習(xí)算法。GBDT是將梯度下降和決策樹(shù)相結(jié)合，基于前一個(gè)分類(lèi)器殘差減少的方向上，構(gòu)造新的分類(lèi)器，通過(guò)多次迭代構(gòu)造一組弱分類(lèi)器，對(duì)弱分類(lèi)器輸出結(jié)果進(jìn)行加權(quán)累加，累加結(jié)果作為強(qiáng)分類(lèi)器輸出[21]。XGBoost與GBDT相比，其優(yōu)點(diǎn)在于改變了GBDT基于Boosting串行序列化求解問(wèn)題的方式，利用CPU多線程分布式并行計(jì)算，并通過(guò)對(duì)殘差進(jìn)行泰勒二次展開(kāi)進(jìn)行求解，從而打破了現(xiàn)有庫(kù)的計(jì)算速度和精度，使得數(shù)據(jù)處理和運(yùn)算的速度得到了提升。

造成XGBoost模型精度高、運(yùn)行速度快的另一個(gè)因素是選用分類(lèi)回歸樹(shù)(classification and regression tree，簡(jiǎn)稱(chēng)CART)作為決策樹(shù)。由于CART樹(shù)的葉子節(jié)點(diǎn)對(duì)應(yīng)的值是一個(gè)實(shí)際的分?jǐn)?shù)，而非一個(gè)確定的類(lèi)別，這使得優(yōu)化算法的實(shí)現(xiàn)變得更加高效。將XGBoost模型表示為如下數(shù)學(xué)形式：

(5)

其中：F為所有可能的CART樹(shù)的集合；f為一棵具體的CART樹(shù)；K為樹(shù)的棵數(shù)。該XGBoost模型由k棵CART樹(shù)組成。

1.4 性能評(píng)估

通過(guò)如下步驟判斷一個(gè)模型的優(yōu)劣：

1)采取交叉驗(yàn)證的方法測(cè)試模型。如文獻(xiàn)[22]采用交叉驗(yàn)證的方法對(duì)模型進(jìn)行測(cè)試，得到了廣泛的認(rèn)可和使用。交叉驗(yàn)證對(duì)于本實(shí)驗(yàn)的模型測(cè)試也同樣不失為一種好的方法。為了減少隨意性，并使得實(shí)驗(yàn)結(jié)果更準(zhǔn)確，采用10次十折交叉驗(yàn)證。

2)采用Chou[12]在研究信號(hào)肽預(yù)測(cè)中使用的4個(gè)度量參數(shù)評(píng)價(jià)模型。根據(jù)文獻(xiàn)[8,23-24]中的定義，敏感性N、特異性P、準(zhǔn)確率A和馬修斯相關(guān)系數(shù)M分別表示為：

(6)

其中：NTP為含有修飾位點(diǎn)的樣本被正確預(yù)測(cè)為含有該修飾位點(diǎn)數(shù)；NTN為不含有修飾位點(diǎn)的樣本被正確預(yù)測(cè)為不含有該修飾位點(diǎn)數(shù)；NFP為不含有修飾位點(diǎn)的樣本被錯(cuò)誤預(yù)測(cè)為含有該修飾位點(diǎn)數(shù)；NFN為含有修飾位點(diǎn)的樣本被錯(cuò)誤預(yù)測(cè)為不含有該修飾位點(diǎn)數(shù)。

1.5 位置特異性核苷酸偏好特征

位置特異性偏好思想在生物信息學(xué)得到了廣泛應(yīng)用，在功能位點(diǎn)的識(shí)別及啟動(dòng)子位點(diǎn)的識(shí)別方面都取得了非常不錯(cuò)的效果[25-27]。其原理是統(tǒng)計(jì)生物序列中某些關(guān)鍵位置或某種核苷酸出現(xiàn)的概率[26]，將得到的概率值作為位點(diǎn)識(shí)別的特征。受此啟發(fā)，將位置特異性單核苷酸和雙核苷酸偏好特征特征應(yīng)用到RNA修飾位點(diǎn)中。

1.5.1 位置特異性單核苷酸偏好特征

由式(1)可知，每個(gè)樣本由L個(gè)核苷酸組成，對(duì)于一個(gè)基準(zhǔn)數(shù)據(jù)集的所有樣本，可分別計(jì)算出第j(j=1,2,…,L)個(gè)位置上4種核苷酸出現(xiàn)的概率，并用一個(gè)長(zhǎng)度為4的向量表示：

(7)

其中，MA,j、MC,j、MG,j和MU,j分別為A、C、G和U這4種核苷酸在第j個(gè)位置出現(xiàn)的概率。

將j從1取值到L得到的位置特異性向量Mj組合在一起，構(gòu)成一個(gè)4×L維的位置特異性單核苷酸偏好矩陣M：

(8)

1.5.2 位置特異性雙核苷酸偏好特征

取2個(gè)相鄰的核苷酸為一個(gè)單元，則式(1)樣本序列可表示為

S′={N1,N2,…,NL-1}，

(9)

其中，Nj={Rj,Rj+1}，j=1,2,…,L-1表示第j個(gè)位置雙核苷酸的類(lèi)型，而雙核苷酸的種類(lèi)共有16種，即Nj∈{AA,AC,AG,AU,CA,…,UU}。

與M的計(jì)算過(guò)程類(lèi)似，可計(jì)算出一個(gè)維度為16×(L-1)的位置特異性雙核苷酸偏好矩陣：

(10)

1.6 特征選擇

雖然用PseKNC對(duì)RNA序列進(jìn)行編碼是一個(gè)不錯(cuò)的選擇，但這種編碼方式有一個(gè)缺點(diǎn)，即數(shù)據(jù)樣本維度會(huì)出現(xiàn)爆炸式的增長(zhǎng)，從而會(huì)出現(xiàn)以下問(wèn)題：1)容易造成過(guò)擬合導(dǎo)致模型的泛化能力偏低；2)信息冗余和噪聲會(huì)導(dǎo)致模型準(zhǔn)確率低下，達(dá)不到預(yù)期效果；3)高維度樣本使得計(jì)算機(jī)的運(yùn)算負(fù)荷加重，從而使模型的運(yùn)行時(shí)間大大增加。

為了解決上述問(wèn)題，對(duì)數(shù)據(jù)樣本進(jìn)行特征選擇，進(jìn)而降低樣本維度。本研究采用F-score[28]特征選擇方法，定義如下：

(11)

2 結(jié)果

2.1 PseKNC編碼方式中最優(yōu)k值的確定

對(duì)m1A、m6A、m5C這3種RNA序列采用PseKNC的編碼方式進(jìn)行編碼，但k取不同的值時(shí)，結(jié)果也不同，為了找到最優(yōu)的k值，對(duì)k=2,3,4,5,6,7時(shí)分別進(jìn)行編碼，其十折交叉驗(yàn)證的結(jié)果如圖1所示。

圖1 PseKNC編碼中不同的k值對(duì)m1A, m6A 和 m5C修飾位點(diǎn)的識(shí)別準(zhǔn)確率的影響

從圖1可看出，k取不同值時(shí)，m1A、m6A、m5C這3種RNA修飾位點(diǎn)的識(shí)別準(zhǔn)確率也隨之波動(dòng)，其中m1A對(duì)于不同的k值波動(dòng)較小，而k值的不同對(duì)m6A和m5C的影響較大。當(dāng)k=5時(shí)，m1A修飾位點(diǎn)的識(shí)別準(zhǔn)確率最高；k=4時(shí)，m6A修飾位點(diǎn)的識(shí)別準(zhǔn)確率最高；k=3時(shí)m5C修飾位點(diǎn)的識(shí)別準(zhǔn)確率最高。在進(jìn)行特征選擇之前，m1A、m6A、m5C的PseKNC編碼中k的最優(yōu)值分別為5、4、3。

為了進(jìn)一步提高識(shí)別準(zhǔn)確率，采用F-score特征選擇方法得到k最終的最優(yōu)解。

在m1A、m6A、m5C取不同k值(k=2,3,4,5,6,7)編碼后的序列基礎(chǔ)上，融合位置特異性單核苷酸偏好特征和位置特異性雙核苷酸偏好特征，再對(duì)其進(jìn)行特征選擇。十折交叉驗(yàn)證的結(jié)果如表1所示。

表1 特征選擇后，不同的k值對(duì)m1A、m6A、m5C修飾位點(diǎn)的識(shí)別結(jié)果

從表1可看出，融合了位置特異性單核苷酸偏好特征和位置特異性雙核苷酸偏好特征并經(jīng)過(guò)特征選擇后，m1A、m6A、m5C這3種RNA修飾位點(diǎn)的識(shí)別準(zhǔn)確率均有較大提升。其中：m1A在k=7時(shí)，識(shí)別準(zhǔn)確率達(dá)到最高，為99.9%；m6A在k=5時(shí)，識(shí)別準(zhǔn)確率達(dá)到最高，為88.1%；m5C在k=3時(shí)，識(shí)別準(zhǔn)確率達(dá)到最高，為88.2%。而在此之前，m1A、m6A、m5C這3種RNA修飾位點(diǎn)的識(shí)別準(zhǔn)確率最高時(shí)其PseKNC編碼對(duì)應(yīng)的k值分別為5、4、3，顯然，除了m5C的PseKNC編碼中最優(yōu)k值未發(fā)生改變，其他2個(gè)均發(fā)生了改變。

2.2 基于網(wǎng)格搜索的XGBoost模型參數(shù)尋優(yōu)

在確定了PseKNC編碼的最優(yōu)k值后，繼續(xù)對(duì)XGBoost預(yù)測(cè)模型的參數(shù)進(jìn)行尋優(yōu)。要想完全發(fā)揮XGBoost的強(qiáng)大性能，對(duì)其進(jìn)行調(diào)參是必不可少的一項(xiàng)工作。這里選用網(wǎng)格搜索的方法對(duì)其進(jìn)行調(diào)參。網(wǎng)格搜索的原理是：在所有候選的參數(shù)中，通過(guò)循環(huán)遍歷，嘗試每種可能性，交叉驗(yàn)證后，表現(xiàn)最好的參數(shù)組合就是最終結(jié)果。該方法的優(yōu)點(diǎn)是結(jié)果準(zhǔn)確，但缺點(diǎn)是當(dāng)參數(shù)數(shù)量過(guò)多時(shí)，參數(shù)尋優(yōu)的計(jì)算過(guò)程非常耗時(shí)。XGBoost中的參數(shù)主要分為通用參數(shù)、學(xué)習(xí)任務(wù)參數(shù)和命令行參數(shù)3大類(lèi)。其中通用參數(shù)有20多個(gè)，學(xué)習(xí)任務(wù)參數(shù)有4個(gè)，命令行參數(shù)有十多個(gè)，若對(duì)這些參數(shù)進(jìn)行網(wǎng)格搜索尋優(yōu)，將會(huì)非常耗時(shí)，使調(diào)參的工作面臨巨大挑戰(zhàn)。

為了應(yīng)對(duì)這一挑戰(zhàn)，針對(duì)性地選取一些核心參數(shù)，將對(duì)模型性能影響不大的參數(shù)剔除。據(jù)此，在通用參數(shù)中選取booster、learning_rate、max_depth、min_child_weight、subsample、colsample_bytree、gamma、lambda和alpha這8個(gè)參數(shù)，在學(xué)習(xí)任務(wù)參數(shù)中選取objective參數(shù)，在命令行參數(shù)中選取num_round參數(shù)。其中：max_depth和min_child_weight對(duì)決策樹(shù)的構(gòu)建起約束作用；subsample和colsample_bytree均是關(guān)于采樣的參數(shù)；lambda和alpha均是正則項(xiàng)。為了進(jìn)一步減少運(yùn)算量，采用組合分批網(wǎng)格搜索的策略，即將具有相同作用的參數(shù)組合在一起，對(duì)組合進(jìn)行網(wǎng)格搜索，將最優(yōu)組合與其他參數(shù)組合在一起，再進(jìn)行網(wǎng)格搜索。調(diào)參后的XGBoost模型在jackknife測(cè)試下的識(shí)別結(jié)果如表2所示。從表2可看出，m6A修飾位點(diǎn)的識(shí)別準(zhǔn)確率從88.1%提升到了92.6%，m5C修飾位點(diǎn)的識(shí)別準(zhǔn)確率從88.2%提升到了89.6%，雖然m1A修飾位點(diǎn)的識(shí)別準(zhǔn)確率并未得到提升，但在此之前已經(jīng)達(dá)到了99.9%這樣一個(gè)非常理想的水平。

表2 參數(shù)調(diào)整后的XGBoost模型的識(shí)別結(jié)果

ROC(receiver operating characteristic)曲線是一個(gè)能直觀展現(xiàn)模型性能的另一個(gè)重要指標(biāo)[29]。ROC曲線下側(cè)包含的面積(AUC)越大，模型的性能越好[30]。圖2為XGBoost預(yù)測(cè)模型分別對(duì)m1A、m6A、m5C這3種RNA修飾位點(diǎn)進(jìn)行識(shí)別后生成的ROC曲線。從圖2可看出，m1A、m6A、m5C所對(duì)應(yīng)的AUC值分別為0.998 6、0.931 2和0.955 8，表明XGBoost預(yù)測(cè)模型的魯棒性非常好。

圖2 XGBoost模型在m1A,m6A和m5C修飾位點(diǎn)上的識(shí)別性能

2.3 不同方法的識(shí)別結(jié)果對(duì)比

將XGBoost預(yù)測(cè)模型的識(shí)別結(jié)果與文獻(xiàn)[13]使用的SVM預(yù)測(cè)模型的識(shí)別結(jié)果進(jìn)行對(duì)比，2個(gè)模型經(jīng)過(guò)jackknife測(cè)試后的結(jié)果如表3所示。從表3可看出，XGBoost預(yù)測(cè)模型和SVM預(yù)測(cè)模型在m1A修飾位點(diǎn)的識(shí)別上均取得了較好的結(jié)果，準(zhǔn)確率分別達(dá)到了99.9%、99.1%；在m6A修飾位點(diǎn)的識(shí)別上，XGBoost預(yù)測(cè)模型的準(zhǔn)確率為92.6%，SVM預(yù)測(cè)模型的準(zhǔn)確率為90.4%，提升了2.2%；在m5C修飾位點(diǎn)的識(shí)別上，XGBoost預(yù)測(cè)模型的準(zhǔn)確率達(dá)到了89.6%，遠(yuǎn)高于SVM預(yù)測(cè)模型的77.5%，提升了12.1%。

表3 XGBoost與SVM的識(shí)別結(jié)果比較

3 結(jié)束語(yǔ)

為了更快速、準(zhǔn)確地識(shí)別RNA序列中的修飾位點(diǎn)m1A、m6A、m5C，提出了一種融合位置特異性單核苷酸及雙核苷酸偏好特征的PseKNC編碼方式，并構(gòu)建了一個(gè)基于XGBoost的RNA修飾位點(diǎn)的預(yù)測(cè)模型。與現(xiàn)有的SVM預(yù)測(cè)模型相比，其準(zhǔn)確率和馬修斯相關(guān)系數(shù)均取得了明顯提升，其中，對(duì)于在現(xiàn)有SVM預(yù)測(cè)模型上識(shí)別效果相對(duì)較差的m5C修飾位點(diǎn)，在XGBoost預(yù)測(cè)模型上取得了較大的突破，識(shí)別準(zhǔn)確率從77.5%提高到了89.6%，馬修斯相關(guān)系數(shù)從0.552提高到了0.792，此外，敏感性和特異性也分別從0.758和0.792提高到了0.913和0.880。XGBoost預(yù)測(cè)模型的提出為RNA修飾位點(diǎn)的識(shí)別提供了高效、可靠的方法。