磷酸化修飾位點(diǎn)預(yù)測分析研究

鄧曉政，徐瑞杰，陳宇

（青島大學(xué) 數(shù)據(jù)科學(xué)與軟件工程學(xué)院，山東青島 266071）

蛋白質(zhì)磷酸化修飾是生物體內(nèi)一種普遍存在的翻譯后修飾類型[1]，至少有三分之一的細(xì)胞蛋白中存在磷酸化修飾[2]。該過程具有可逆性,與許多細(xì)胞內(nèi)生化過程相關(guān)，如蛋白-蛋白相互作用[3]、免疫反應(yīng)[4]、有絲分裂/細(xì)胞周期[5]等。磷酸化修飾的失調(diào)可能會導(dǎo)致多種疾病發(fā)生[6]。研究發(fā)現(xiàn)，大多癌癥以及肌肉疾病與磷酸化位點(diǎn)及其協(xié)同位點(diǎn)的異常修飾相關(guān)[7]。因此，準(zhǔn)確定位磷酸化修飾位點(diǎn)對人類疾病的研究具有重大意義。

近年來，磷酸化修飾位點(diǎn)的定位更多地依賴于機(jī)器學(xué)習(xí)算法，如支持向量機(jī)、人工神經(jīng)網(wǎng)絡(luò)、隱形馬爾可夫模型和隨機(jī)森林等。2014年，Dou等人[8]通過8種不同的序列特征打分以及支持向量機(jī)算法構(gòu)建預(yù)測模型。Hamid等人[9]在2016年開發(fā)了一種基于序列結(jié)構(gòu)屬性以及隨機(jī)森林算法的磷酸化位點(diǎn)預(yù)測工具。2017年，Tan等人[10]以Dou的磷酸化位點(diǎn)數(shù)據(jù)為模板開發(fā)了基于位置的卡方表特征和偽位置特定評分矩陣的支持向量機(jī)分類器。雖然結(jié)合機(jī)器學(xué)習(xí)算法對磷酸化修飾位點(diǎn)的預(yù)測取得了一定進(jìn)展，但是許多模型不具備跨物種普適預(yù)測能力，對不同物種的預(yù)測存在偏差。因此針對于單個(gè)物種利用機(jī)器學(xué)習(xí)算法開發(fā)精度更高的蛋白質(zhì)翻譯后修飾位點(diǎn)預(yù)測器十分必要。

本文提出了一種專門用于定位人類磷酸化修飾位點(diǎn)的預(yù)測方法，該方法采用動態(tài)篩選機(jī)制針對數(shù)據(jù)集選取最優(yōu)序列氨基酸理化性質(zhì)和氨基酸組成等特征，并基于支持向量機(jī)算法進(jìn)行磷酸化修飾位點(diǎn)的準(zhǔn)確定位，相關(guān)算法以及預(yù)測工具命名為HPSP（Human Phosphorylation Site Predictor）。

1 材料與方法

1.1 材料

從UniprotKB/Swiss-Prot中下載經(jīng)過驗(yàn)證的人類蛋白質(zhì)序列數(shù)據(jù)共計(jì)20 402條，篩選出關(guān)于絲氨酸（S）、蘇氨酸（T）和酪氨酸（Y）磷酸化位點(diǎn)的蛋白質(zhì)序列。使用CD-HIT程序去除同一性閾值為30%的冗余序列，最后獲得絲氨酸、蘇氨酸、酪氨酸各4 917條、2 002條以及609條蛋白質(zhì)序列片段。

分別以磷酸化位點(diǎn)為中心，截取前后各13個(gè)氨基酸殘基的對稱窗口作為編碼特征的磷酸化序列片段，窗口內(nèi)殘基總數(shù)為27個(gè)。如果磷酸化位點(diǎn)在N端或C端附近，使用大寫字母“O”代表缺失的氨基酸。經(jīng)過上述處理后分別得到3種序列片段各22 763條、3 805條和1 122條作為正樣本集。假設(shè)在同一蛋白質(zhì)上沒有任何磷酸化信息標(biāo)記的絲氨酸/蘇氨酸/酪氨酸殘基都是非磷酸化位點(diǎn)，本文將同一蛋白質(zhì)中沒有任何磷酸化信息標(biāo)記的絲氨酸/蘇氨酸/酪氨酸殘基作為負(fù)樣本集。從各自蛋白質(zhì)序列片段中截取未被磷酸化的絲氨酸、蘇氨酸、酪氨酸前后各13個(gè)殘基的對稱窗口作為負(fù)樣本集，分別為319 176條、91 936 條、11 628 條。

為避免正負(fù)樣本集數(shù)目極不平衡導(dǎo)致模型過擬合的情況，對負(fù)樣本集進(jìn)行隨機(jī)抽樣，使正負(fù)樣本集保持相同的數(shù)量規(guī)模。

1.2 方法

1.2.1 物理化學(xué)屬性特征

氨基酸的物理化學(xué)屬性（PCP）如疏水性、分子量、可及表面等，在蛋白質(zhì)的結(jié)構(gòu)和功能研究中具有重要作用。蛋白質(zhì)序列的結(jié)構(gòu)和功能在某種程度上與組成蛋白質(zhì)的每一個(gè)氨基酸的物理化學(xué)屬性都有很大關(guān)聯(lián)。而且PCP已經(jīng)成功應(yīng)用于蛋白質(zhì)磷酸化修飾位點(diǎn)預(yù)測工作中[11]。從AAindex數(shù)據(jù)庫[12]中獲取544種物理化學(xué)性質(zhì)，利用F值檢驗(yàn)進(jìn)行物理化學(xué)屬性的篩選，具體公式如式1所示。

式中，μi+、μi-分別表示正樣本和負(fù)樣本中第i個(gè)物理化學(xué)屬性的平均值，σi+，σi-分別表示正樣本和負(fù)樣本中第i個(gè)物理化學(xué)屬性的方差。

某一個(gè)物理化學(xué)屬性的F值越高，則認(rèn)為該物理化學(xué)屬性特征能很好地區(qū)分磷酸化位點(diǎn)和非磷酸化位點(diǎn)。對544個(gè)物理化學(xué)屬性分別計(jì)算F值，取F值最大的前20個(gè)物理化學(xué)屬性作為特征。由于AAindex中氨基酸指數(shù)范圍較廣，可能對訓(xùn)練模型的預(yù)測造成影響，因此對氨基酸指數(shù)進(jìn)行歸一化處理，具體的處理過程如式2所示。

式中，Pi表示某一物理化學(xué)屬性在AAindex中的值，Pinorm表示經(jīng)過歸一化處理后的結(jié)果。Pmax、Pmin分別表示某一物理化學(xué)屬性在AAindex中的最大值和最小值。

最后每一個(gè)序列片段可以得到20維物理化學(xué)屬性特征 PCP(P)=(p1,p2,p3…p20)。

1.2.2 KNN距離

KNN（K最鄰近分類算法）距離從可能的磷酸化位點(diǎn)周圍局部序列中提取序列相似性信息，這些序列相似性信息反映了蛋白質(zhì)序列結(jié)構(gòu)上的差異。序列相似性越高，序列之間存在的功能就越相似。利用KNN距離提取相似性信息，首先計(jì)算未知蛋白片段到正負(fù)樣本集的距離。距離D(S1,S2)定義為式3。

式中，p表示序列片段中心兩側(cè)氨基酸殘基的數(shù)目，i表示序列片段中氨基酸的位置。Blo是基于BLOSUM62矩陣[13]得到的氨基酸相似性打分矩陣。具體的轉(zhuǎn)換過程定義為式4。

其中，M表示BLOSUM62矩陣，min和max分別為BLOSUM62矩陣中的最小值和最大值。計(jì)算未知蛋白片段到正負(fù)樣本集的KNN距離并排序后，選擇前k個(gè)樣本并統(tǒng)計(jì)這k個(gè)樣本中正樣本所占的百分比，即為最終的KNN特征值。KNN特征提取中，k值的選擇對于分類效果具有很大影響，k值過大會使算法時(shí)間復(fù)雜度越高，而k值過小則可能丟棄真正有意義的序列。因此本文選取5個(gè)不同的k值，每條蛋白質(zhì)序列可得到 5 維 KNN 特征 KNN(P)=(k1,k2,k3,k4,k5)。

1.2.3 信息熵與熵密度

1948年，Shannon[14]提出信息熵以度量給定系統(tǒng)信息含量。越是混亂的系統(tǒng)，信息熵越高；而越是有序的系統(tǒng)，信息熵越低。對于蛋白質(zhì)序列片段而言，保守信息隨著位置變化而變化，而磷酸化修飾位點(diǎn)附近的某些殘基對磷酸化位點(diǎn)的識別有重要影響。因此信息熵是對蛋白質(zhì)序列片段中各個(gè)位置殘基保守性信息量的度量。將信息熵與熵密度作為衡量磷酸化可能性高低的標(biāo)準(zhǔn)，具體計(jì)算方法如式5和式6所示。

其中，X表示蛋白質(zhì)序列片段，fi(X)表示在該蛋白質(zhì)片段中第i個(gè)氨基酸殘基的出現(xiàn)頻率。最后每一條蛋白質(zhì)序列片段得到1維的信息熵H(P)=x和20維的熵密度特征S(P)=(s1,s2,s2…s20)。

1.2.4 氨基酸組成

氨基酸組成（AAC）是最經(jīng)典的蛋白質(zhì)特征編碼方法之一。Lee等[15]利用修飾位點(diǎn)周圍的氨基酸組成信息作為預(yù)測蛋白質(zhì)翻譯后修飾位點(diǎn)的一個(gè)重要特征。AAC與每一個(gè)氨基酸的物理化學(xué)屬性能夠反映蛋白質(zhì)序列的生物化學(xué)環(huán)境，與磷酸化序列片段具有相似的生物化學(xué)環(huán)境的片段很有可能存在潛在的磷酸化位點(diǎn)。每個(gè)蛋白質(zhì)片段可被編碼為一個(gè)20維的向量 AAC(P)=(a1,a2,a3…a20)。

2 結(jié)果與分析

本文使用LIBSVM進(jìn)行預(yù)測模型的搭建。首先從已獲得數(shù)據(jù)集合中選取部分?jǐn)?shù)據(jù)來進(jìn)行訓(xùn)練，分別從訓(xùn)練集中隨機(jī)抽取十次正負(fù)樣本，選取比例為1∶ 1。

在此基礎(chǔ)上，對支持向量機(jī)徑向基（RBF）核函數(shù)和c、g參數(shù)進(jìn)行優(yōu)化。在利用LIBSVM進(jìn)行預(yù)測之前，利用grid.py得到最優(yōu)cost值和gamma值，使用交叉驗(yàn)證對比選取RBF核類型和C-SVC類型來創(chuàng)建模型。RBF核函數(shù)如式7所示。

采用十倍交叉驗(yàn)證法對預(yù)測方法的性能進(jìn)行了評價(jià)。利用精度（Pr）、靈敏度（Sn）、特異性（Sp）、準(zhǔn)確度（ACC）和馬太相關(guān)系數(shù)（MCC）對該預(yù)測系統(tǒng)的性能進(jìn)行評價(jià)。MCC是反映正負(fù)樣本成功預(yù)測的綜合指標(biāo)，其值范圍為-1～1，數(shù)值越趨近于1表示預(yù)測性能越好，數(shù)值越趨近于-1表示預(yù)測性能越差。該模型在訓(xùn)練集上的預(yù)測結(jié)果如表1所示。

表1 訓(xùn)練集上的預(yù)測結(jié)果

3 討論

3.1 特征分析

為了獲得磷酸化修飾位點(diǎn)預(yù)測的最優(yōu)特征，以酪氨酸正負(fù)樣本集為例，進(jìn)一步驗(yàn)證四類特征單獨(dú)預(yù)測以及全融合特征預(yù)測在訓(xùn)練集上的性能，具體表現(xiàn)如表2所示。對于單個(gè)特征預(yù)測而言，KNN距離特征的預(yù)測性能最優(yōu)，MCC為0.54，比其他三類特征預(yù)測的MCC值高0.2左右，這也印證了前述特征分析中正負(fù)樣本KNN打分值差異性較大的結(jié)論。雖然全融合特征的Sp值相比于KNN距離特征沒有明顯的提高，但是在其他四個(gè)評價(jià)指標(biāo)上存在顯著提高。由此可見，四類特征全融合表現(xiàn)最優(yōu)，最終選取四類特征全融合來搭建預(yù)測模型。

表2 基于不同特征的酪氨酸磷酸化修飾位點(diǎn)預(yù)測性能

3.2 工具對比

磷酸化位點(diǎn)預(yù)測工具有很多，為了進(jìn)一步說明本算法在已知磷酸化位點(diǎn)預(yù)測上的優(yōu)勢，將本文提出的預(yù)測模型HPSP與RF-Phos方法在獨(dú)立測試集上進(jìn)行比較。從文獻(xiàn)[16]中收集獨(dú)立測試數(shù)據(jù)集，從400個(gè)磷酸化修飾位點(diǎn)中除去訓(xùn)練集中存在的數(shù)據(jù)，剩下的作為獨(dú)立測試集。正樣本數(shù)量分別為188條、82條、100條，負(fù)樣本數(shù)量分別為364條、246條、128條。

為了更加客觀地比較，使用四種評價(jià)指標(biāo)對預(yù)測性能進(jìn)行評估。在獨(dú)立測試集上，雖然在酪氨酸上的MCC、ACC值與RF-Phos相當(dāng)，但是在特異性Sp上有較大的提高。在絲氨酸與蘇氨酸上，靈敏度和特異性都高于RF-Phos。三種磷酸化修飾的MCC與ACC值均高于RF-Phos，其中在蘇氨酸磷酸化位點(diǎn)上的性能表現(xiàn)最優(yōu)，準(zhǔn)確率ACC提高了0.3%，MCC從0.56提高到0.63，是性能評估中一個(gè)全面的綜合指標(biāo)。表3所示的測試結(jié)果表明本模型HPSP整體優(yōu)于RF-Phos的預(yù)測算法。

表3 模型HPSP與RF-Phos工具性能比較結(jié)果

該結(jié)果充分說明模型HPSP對磷酸化修飾位點(diǎn)的預(yù)測具有良好的性能。今后，可以探索其他的方法來進(jìn)一步提高模型的預(yù)測能力，如增加其他具有代表意義的特征、改變滑動窗口長度大小或者配合其他機(jī)器學(xué)習(xí)分類預(yù)測算法等方法；同時(shí)也可以將激酶信息整合到預(yù)測模型中以識別激酶特異性磷酸化修飾位點(diǎn)。

4 結(jié)論

本文基于支持向量機(jī)開發(fā)了一種針對人類蛋白質(zhì)序列數(shù)據(jù)的非特異性磷酸化位點(diǎn)預(yù)測方法HPSP，該方法整合了信息熵和密度熵等四類特征，利用F值檢驗(yàn)方法篩選出能夠顯著區(qū)分磷酸化位點(diǎn)和非磷酸化位點(diǎn)的特征，利用支持向量機(jī)算法進(jìn)行預(yù)測模型的搭建。在獨(dú)立測試集，對絲氨酸、蘇氨酸、酪氨酸磷酸化修飾位點(diǎn)的預(yù)測準(zhǔn)確值分別達(dá)76.33%、81.10%、71.50%。