基于雙向門控循環(huán)單元的5-甲基胞嘧啶位點預測

黃修威，方中純，李海榮

（1 內(nèi)蒙古科技大學 信息工程學院，包頭 014010；2 內(nèi)蒙古科技大學 工程訓練中心（創(chuàng)新創(chuàng)業(yè)教育學院），包頭 014010）

迄今為止，已有170 多種化學修飾在RNA 中被發(fā)現(xiàn)［1-3］.5-甲基胞嘧啶（m5C）是最常見的轉錄后修飾之一，它是通過RNA 甲基轉移酶催化將甲基轉移到胞嘧啶環(huán)的第5 位碳原子形成的.m5C 位點參與多種生物過程，包括調(diào)節(jié)蛋白質(zhì)翻譯和RNA-蛋白質(zhì)相互作用，調(diào)節(jié)核mRNA 輸出和RNA 可變剪接，增強RNA 穩(wěn)定性等［4-7］.此外，有研究還證明m5C 修飾與許多疾病有關，如乳腺癌［8］、常染色體隱性智力殘疾［9］、肌萎縮側索硬化［10］和帕金森氏?。?1］.因此，準確鑒定m5C 位點在RNA 中的位置是開展相應疾病和生物學功能研究的首要和關鍵任務.

目前已經(jīng)開發(fā)出了數(shù)種方法來預測RNA 中的m5C位點.在實驗中，已經(jīng)開發(fā)了幾種傳統(tǒng)的高通量測序技術，如亞硫酸氫鹽測序［12］，m5C-RIP［13］，mi-CLIP［14］和RBS-seq［15］等實驗方法.雖然這些技術在鑒定不同物種的RNA 中的m5C 位點方面取得了一定的成功，然而，這些技術非常耗時和昂貴，難以跟上后基因組時代RNA 序列數(shù)量的急劇增加.因此又提出了多種基于機器學習算法來預測RNA 序列中m5C 位點的模型.FENG 等［16］首先在人類中通過偽二核苷酸組成（PseDNC）的三個理化特征設計了基于支持向量機（Support Vector Machine，SVM）的m5CPseDNC 模型.后來，QIU 等［17］開發(fā)了一個稱為iRNAm5C-PseDNC 的模型，該模型建立在不平衡和冗余的數(shù)據(jù)集上并通過使用隨機森林來預測人類的m5C位點.之后，F(xiàn)ANG 等［18］開發(fā)了一個基于SVM的預測模型RNAm5CPred，來預測人類的m5C位點，并且使用獨立的測試集來評估不同的方法.對于擬南芥中的m5C 位點，SONG 等［19］率先開發(fā)了預測模型PEA-M5C，在該模型中，結合了三種特征編碼技術以提高綜合性能.LI 等［20］從GEO 數(shù)據(jù)庫中收集數(shù)據(jù)，基于隨機森林算法開發(fā)了一個名為RNAm5 Cfinder 的模型，可用于預測小鼠和人類的m5C 位點.2020 年CHEN 等［21］通過引入位置特異性傾向相關特征，建立了一個新的模型m5CPred-SVM，來預測三種不同物種的RNA m5C 位點.同年DOU 等［22］提出了一種新的基于支持向量機（SVM）的工具，稱為iRNAm5C_SVM，其通過組合多個序列特征來識別擬南芥中的m5C 位點.最近，CHAI 等［23］通過采用特征融合策略來利用信息序列，并堆疊集成學習框架，結合了五種流行的機器學習算法提出了Staem5模型，用于小鼠和擬南芥m5C位點的識別.

傳統(tǒng)機器學習方法的性能高度依賴于特征提取的有效性，然而特定RNA 修飾最相關的特征組合很難判斷.基于深度學習的計算架構能夠在沒有任何人工干預的情況下從序列中提取最重要的特征，并進行端到端預測.通過多個神經(jīng)元層和激活函數(shù)，它將獲得從原始輸入到潛在表征的映射能力，而潛在表征是由標記數(shù)據(jù)訓練的.在這樣的數(shù)據(jù)驅動模型下，與RNA 修飾位點語義信息相關的RNA序列深層特征就會顯現(xiàn)出來.

為了彌補現(xiàn)有計算模型在性能和計算成本方面的缺點，本文提出了一種簡單有效的基于雙向門控循環(huán)單元（BiGRU）的預測模型pm5C-BGRU，用于鑒定RNA 序列中的m5C 位點.輸入RNA 序列由獨熱編碼（One-hot encoding）和核苷酸化學性質(zhì)（NCP）的組合來表示.核苷酸的化學性質(zhì)是核苷酸在官能團、氫鍵和環(huán)結構方面最基本的表征.BiGRU能夠從RNA 序列中學習序列前后特征間的長期依賴關系，這使得pm5C-BGRU 能夠更準確地識別m5C 位點.最后使用10 倍交叉驗證和獨立數(shù)據(jù)集測試的方法來評估和比較pm5C-BGRU的性能.

1 數(shù)據(jù)和方法

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

高質(zhì)量的數(shù)據(jù)集對于訓練和評估預測模型都極其重要.從最近發(fā)表的文獻中收集了三個物種的m5C 數(shù)據(jù).擬南芥的m5C 數(shù)據(jù)從NCBI 基因表達綜合（GEO）數(shù)據(jù)庫獲得，登錄號為GSE94065［24］.因為擬南芥的m5C 樣本足夠大，所以使用CD-HIT 程序［25］通過將序列同一性的閾值設置為80%來去除冗余樣本.最終，在擬南芥基因組中獲得了6289 個陽性樣本.擬南芥的陰性樣品從它們的基因組中獲得，為了平衡陽性和陰性樣本，隨機挑選相同數(shù)量的擬南芥陰性樣本并使用CD-HIT 程序與陽性樣本進行相同操作.小鼠和人類數(shù)據(jù)集來自于CHEN等［21］的研究工作，其中分別包含小鼠和人類的5563和269 個陽性樣本及相同數(shù)量的陰性樣本.所用樣本均已使用CD-HIT程序去除數(shù)據(jù)集中的相似序列.

基準數(shù)據(jù)集通常分為兩部分，一個是訓練數(shù)據(jù)集，另一個是獨立測試集.訓練數(shù)據(jù)集用于模型構建、交叉驗證和確定學習算法的超參數(shù).獨立測試數(shù)據(jù)集用于測試性能和模型泛化能力.三個物種的訓練和獨立測試數(shù)據(jù)集數(shù)量如表1所示.

表1 不同物種的m5C數(shù)據(jù)集Tab.1 m5C data sets of different species

樣本中每個RNA 序列是長度為41 bp 的RNA片段，這是現(xiàn)有工具中最常用的長度.m5C數(shù)據(jù)集中的每個RNA片段可以表示如下：

其中中心X 是靶向位點，表示m5C 的C（胞嘧啶核苷酸）.N1到N20代表朝向靶位點的上游側翼核苷酸，N22到N41而表示其下游側翼核苷酸.

1.2 RNA片段的編碼

pm5C-BGRU 采用RNA 序列作為輸入.RNA 序列由獨熱編碼和核苷酸化學性質(zhì)（NCP）兩種常用編碼技術的組合來表示.獨熱編碼是序列分析中最常見和最有效的編碼方式之一，它將RNA 中的每種核糖核苷酸序列編碼為一個獨熱向量，RNA 片段的A（腺嘌呤核苷酸）表示為［1，0，0，0］，C（胞嘧啶核苷酸）表示為［0，1，0，0］，G（鳥嘌呤核苷酸）表示為［0，0，1，0］，U（尿嘧啶核苷酸）表示為［0，0，0，1］.而NCP 是RNA 序列中每個核苷酸在三維笛卡爾坐標系中基于它們的三個化學基團的表示.RNA 序列中的四個核苷酸每一個都有不同的化學特性，如表2所示.考慮到環(huán)結構，A 和G 是由兩個環(huán)組成的嘌呤，而C和U 是由一個環(huán)組成的嘧啶；就二級結構形成而言，A 和U 之間的鍵是弱氫鍵，而C 和G 之間的鍵是強氫鍵；關于化學官能團，A和C屬于氨基，而G和U屬于酮基.根據(jù)這三種化學性質(zhì)，這四種核苷酸可以通過賦值1或0分為三個不同的組，即A表示為［1，1，1］，C 表示為［0，1，0］，G 表示為［1，0，0］，U 表示為［0，0，1］.通過將獨熱編碼和NCP 產(chǎn)生的向量連接在一起，RNA 序列中每個核苷酸均被編碼成一個7 維的特征向量，如圖1 所示，即A=［1，0，0，0，1，1，1］，C=［0，1，0，0，0，1，0］，G=［0，0，1，0，1，0，0］，U=［0，0，0，1，0，0，1］，每個RNA 序列由41 個這樣的列表組成，即一條RNA 片段可以編碼成一個7×41的矩陣.

圖1 RNA序列編碼Fig.1 RNA sequence coding

表2 核苷酸化學性質(zhì)Tab.2 Chemical properties of nucleotides

1.3 門控循環(huán)單元

循環(huán)神經(jīng)網(wǎng)絡（Recurrent Neural Network，RNN）是一種能夠記憶上下文信息的網(wǎng)絡架構，非常適合生物序列分析［26］.門控循環(huán)單元（Gated Recurrent Unit，GRU）為RNN 的精簡版，不但能記住上一時刻的輸入信息，而且能夠抑制梯度爆炸或梯度消失的現(xiàn)象.門控循環(huán)單元結構簡單，只有兩個門控單元，分別是重置門和更新門.重置門決定著上一時刻隱含層被遺忘信息的大小，更新門衡量過去信息的重要程度，決定著傳遞到當前隱含層信息的大小.GRU網(wǎng)絡的基本結構如圖2所示.

圖2 GRU網(wǎng)絡基本結構圖Fig.2 Basic Structure of GRU Network

圖2 中x1為模型輸入，ht-1為上一時刻隱含狀態(tài)，ht為當前時刻隱含狀態(tài)，為候選隱含狀態(tài)，rt代表重置門，zt代表更新門.具體計算過程如下：

式中，σ、Tanh為Sigmoid激活函數(shù)和Tanh激活函數(shù)；Wz、Wr、W為遺忘門、更新門、隱藏層對應的權值.

對于序列特征不明顯的RNA 修飾位點預測問題，要考慮的不僅僅是在之前的信息，之后的信息也同樣重要.雙向門控循環(huán)單元（BiGRU）能夠同時在發(fā)生位點前后考慮輸入數(shù)據(jù)的特征，并集成它們，其結構如圖3所示.

圖3 BiGRU網(wǎng)絡Fig.3 BiGRU network

BiGRU 網(wǎng)絡多了一層隱藏層，信息可以通過正向和反向輸入，第x個的輸出為：

本研究中通過將RNA 序列信息輸入到BiGRU網(wǎng)絡，來找尋修飾位點周圍核糖核苷酸排列的潛在規(guī)律.本模型pm5C-BGRU 搭建了兩個尺寸為64 個神經(jīng)元的BiGRU 網(wǎng)絡層，在BiGRU 層之后，設置了一個尺寸為64 個神經(jīng)元的全連接（fully connection）層.兩層之間用一個 flatten 層來連接，同時為防止過擬合，在三個網(wǎng)絡層的后面引入了一個參數(shù)為0.3的 dropout 層.所有網(wǎng)絡層的激活函數(shù)均為ReLU 函數(shù)，它能夠產(chǎn)生稀疏輸出并加速收斂.優(yōu)化器使用Adam 優(yōu)化器，并且將學習率設置為2e-5.為了避免過擬合，設置early-stopping函數(shù)的規(guī)則是：當損失函數(shù)在300 次訓練過程中都沒有減小時，提前結束訓練.網(wǎng)絡結構具體參數(shù)如表3 所示，模型框架如圖4所示.

圖4 模型框架示意圖Fig.4 Model framework

表3 BiGRU網(wǎng)絡詳細超參數(shù)Tab.3 BiGRU network structure and hyper-parameters

1.4 評價指標

為了測量模型的性能，本文計算了敏感性（Sn）、特異性（Sp）、準確性（Acc）和馬修斯相關系數(shù)（MCC）.這些性能指標用公式表示如下：

其中TP為真陽性，TN 為真陰性，F(xiàn)P為假陽性，F(xiàn)N 為假陰性.Sn 反映了被正確識別的真陽性類別的比例，Sp反映正確識別的真陰性類別的比例，Acc根據(jù)正確檢測的正樣本和負樣本來呈現(xiàn)模型的總體性能.MCC 是一個比較均衡的指標，描述預測結果與實際結果之間的相關系數(shù).此外，受試者工作特征曲線（ROC）的曲線下面積AUC 也是分類模型常見的一種性能度量，它表示不同類之間可分離的程度，可評估預測的整體性能.

1.5 實驗環(huán)境

操作系統(tǒng)Windows10，CPU 型號為Intel?Core?i7-10700 CPU@2.90 GHz，運行內(nèi)存為32 GB.該預測模型是基于Pytorch 深度學習框架，編程語言為Python.

2 結果與討論

2.1 序列組成分析

為了更好地預測人類、小鼠和擬南芥的m5C 位點，使用Two Sample Logos［27］比較三個物種的m5C位點和非m5C 位點周圍核苷酸的位置保守性，研究m5C 位點側翼核苷酸的分布和偏好.三個物種序列位置保守的標志如圖5所示.顯然，就核苷酸分布而言，含m5C 的序列與非m5C 樣品有著顯著不同.小鼠序列中G 和C 在含有m5C 位點的序列中顯著富集，而A和U在不含m5C的序列中顯著富集.對于擬南芥，在m5C位點的上游序列有大量的C，相反嘌呤在非m5C上游序列中富集.另外，與擬南芥中C的高豐度相比，其他兩個物種在含有m5C 位點的序列中有核苷酸G 的富集.圖5表明，在m5C位點周圍存在特定的核苷酸偏好，因此通過使用序列信息構建預測模型是合理的.此外，不同物種之間m5C 位點上下游核苷酸的差異說明應當構建物種特異性預測模型.

圖5 m5C和非m5C位點周圍的核苷酸分布圖Fig.5 Nucleotide distribution around m5C and non m5C sites

2.2 不同編碼方法的結果比較

在三個物種訓練數(shù)據(jù)集上使用10 折交叉驗證對所提出的模型進行評估.為了研究One-hot encoding 和NCP 整合后的有效性，進行了三組實驗.第一組實驗只使用了One-hot encoding 表示RNA 序列，第二組實驗僅使用NCP 來表示RNA 序列，而第三組實驗整合了這兩種編碼方法.這些實驗的結果如表4所示.可以看出，與單獨使用One-hot encoding或NCP 相比，將兩種編碼方法進行整合后的模型，在敏感性、特異性、準確性、MCC 方面都均有提升，可以在m5C位點的鑒定中產(chǎn)生更好的結果.

表4 不同編碼方式下模型的性能比較Tab.4 Performance comparison of models under different coding methods

2.3 與現(xiàn)有預測工具對比

本節(jié)在獨立測試集上將本文提出的pm5CBGRU 模型與近期已發(fā)表的預測模型進行了比較.因為不同預測方法使用了不同的基準數(shù)據(jù)集，所以使用獨立測試集可以更客觀.對人類用pm5C-BGRU與RNAm5CPred［18］和m5CPred-SVM［21］比較；小鼠和擬南芥則用pm5C-BGRU 則與m5CPred-SVM 和Staem5［23］進行比較.表5 顯示了不同模型對相應物種的預測結果.

表5 pm5C-BGRU與現(xiàn)有方法的性能比較結果Tab.5 Performance comparison results of pm5C-BGRU and existing methods

從表5 可以看出，pm5C-BGRU 在人類的獨立測試集上幾乎在所有評估指標上都是最佳的.在小鼠的獨立測試集上，Staem5 的Sp 要高于pm5C-BGRU，但其他評估指標還是pm5C-BGRU 最優(yōu).然而，在擬南芥的獨立測試集上，m5CPred-SVM 與pm5CBGRU 取得了相近的預測性能，Acc 僅僅相差0.2.盡管兩個模型在擬南芥獨立測試集上的評價指標相似，但與m5CPred-SVM 相比，pm5C-BGRU 沒有手工特征提取的過程，使得模型更加簡便.所有這些結果表明，pm5C-BGRU 模型在預測m5C 位點方面比其他現(xiàn)有方法表現(xiàn)得更好.

2.4 跨物種預測

本研究分別為人類、小家鼠和擬南芥三個物種建立了三個模型.為了評估這些模型的物種特異性和可移植性，本節(jié)通過使用其中任一模型來測試該模型對其他物種中的m5C 位點的識別性能，測試結果如圖6所示，其中橫坐標為相應模型.

圖6 跨物種預測性能Fig.6 Cross species prediction performance

由圖6 可知，基于自身數(shù)據(jù)構建的模型在自己的獨立測試集上表現(xiàn)良好，而在其他物種上的表現(xiàn)不理想.另外，人類與小鼠相互間的測試結果要優(yōu)于人類與擬南芥、小鼠與擬南芥相互間的測試結果，這可能是因為人類和小鼠均為哺乳動物，親緣關系較近，而擬南芥為植物的緣故.

3 結論與分析

本研究提出了基于深度學習框架的模型pm5CBGRU，用于準確地識別RNA 序列中的m5C 位點.pm5C-BGRU 通過利用獨熱編碼和用于表示RNA 序列的核苷酸化學性質(zhì)組合來提取重要的特征，這種組合有助于實現(xiàn)高效的m5C 位點鑒定結果.并且，與現(xiàn)有的機器學習預測模型相比，pm5C-BGRU表現(xiàn)出了更好的性能.該模型通過精確預測m5C 位點，將對m5C 的生物學功能和相關疾病的研究起到幫助.在后續(xù)的工作中，可以考慮將與特定疾病相關的m5C 位點使用深度學習方法進行預測，但擴展相關疾病的m5C數(shù)據(jù)集也是一項挑戰(zhàn).