鄧曉政，徐瑞杰，陳宇

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

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

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

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

1 材料與方法

1.1 材料

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

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

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

1.2 方法

1.2.1 物理化學屬性特征

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

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

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

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

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

1.2.2 KNN距離

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

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

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

1.2.3 信息熵與熵密度

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

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

1.2.4 氨基酸組成

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

2 結果與分析

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

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

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

表1 訓練集上的預測結果

3 討論

3.1 特征分析

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

表2 基于不同特征的酪氨酸磷酸化修飾位點預測性能

3.2 工具對比

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

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

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

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

4 結論

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