基于XGBoost的RNA修飾位點的識別

呂成偉， 樊永顯

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

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

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

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

1 材料與方法

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

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

2)對生物序列進行編碼，使其能夠被預(yù)測模型識別；

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

4)進行交叉驗證，并客觀地評價預(yù)測模型的好壞。

1.1 基準數(shù)據(jù)集

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

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

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

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

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

(1)

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

本研究用PseKNC對RNA序列進行編碼，最終得到具有4k個分量的向量，即

(2)

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

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

(5)

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

1.4 性能評估

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

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

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

(6)

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

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

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

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

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

(7)

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

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

(8)

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

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

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

(9)

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

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

(10)

1.6 特征選擇

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

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

(11)

2 結(jié)果

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

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

圖1 PseKNC編碼中不同的k值對m1A, m6A 和 m5C修飾位點的識別準確率的影響

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

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

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

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

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

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

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

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

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

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

圖2 XGBoost模型在m1A,m6A和m5C修飾位點上的識別性能

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

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

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

3 結(jié)束語

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