王皆恒，李 校

(四川大學生命科學學院 四川省分子生物與生物技術(shù)重點實驗室，成都610064)

1 引 言

對生命活動起到關(guān)鍵作用的兩類生物大分子就是蛋白質(zhì)和DNA.而這二者之間的相互作用則是大量生理生化反應(yīng)的核心，它在基因表達調(diào)控[1]、組蛋白修飾[2]，DNA復制、修復和重組[3]等細胞過程中發(fā)揮著極其重要的作用.蛋白質(zhì)中有一類序列被稱為“功能殘基”，它們只占龐大的功能蛋白質(zhì)序列中的一小部分，但卻真正的參與了與其他生物大分子的相互作用以及各類生理生化反應(yīng).也正因如此解析這類“功能殘基”的真正功能和作用位點變得極為關(guān)鍵.例如：如果想要深入地了解轉(zhuǎn)錄過程的機理，不可避免地需要首先了解轉(zhuǎn)錄過程中的DNA相互作用位點.關(guān)于蛋白質(zhì)功能殘基位點的尋找主要有兩種方法：第一種方法即為傳統(tǒng)的通過位點突變來確定蛋白質(zhì)功能殘基的實驗方法，實驗方法的準確率較高但需要大量人力物力的投入，同時實驗周期較長，完全無法滿足目前生物學數(shù)據(jù)的增長速度.由于蛋白質(zhì)的序列結(jié)構(gòu)和功能之間有著緊密的聯(lián)系，若兩種蛋白質(zhì)的功能殘基序列或結(jié)構(gòu)相似，那么從生物意義上來講這二者很有可能存在相似的生物學功能，因此僅僅通過單純的數(shù)據(jù)計算也可以預(yù)測得到有意義的功能預(yù)測結(jié)果.

隨著基因組/轉(zhuǎn)錄組測序技術(shù)的不斷發(fā)展和成本不斷降低，大量的生物學數(shù)據(jù)也需要這樣一種高效的方法快速篩選需要的生物學數(shù)據(jù)，進而快速而精準的預(yù)測蛋白質(zhì)上的功能殘基位點.

在蛋白質(zhì)轉(zhuǎn)移至靶細胞形式功能的過程中，蛋白質(zhì)的部分基團承載著前往靶細胞的定位信息，通過這些基團內(nèi)的信息尋找靶細胞相關(guān)相互作用位點的過程也被稱之為亞細胞定位[4].在亞細胞定位信息領(lǐng)域，通過支持向量機(SVM)進行亞細胞定位分析已經(jīng)取得顯著成果.同時，支持向量機也可以對蛋白質(zhì)中簡單的四種超二級結(jié)構(gòu)進行預(yù)測[5].經(jīng)過改進的Weighted SVM[6]對于蛋白質(zhì)磷酸化位點也有著不錯的預(yù)測效果.而多重支持向量機首尾相連進行大量特種歸類也被稱為神經(jīng)網(wǎng)絡(luò)算法[7].生物信息學分析著眼于序列比對、結(jié)構(gòu)預(yù)測、分子進化等領(lǐng)域[8]，隨著數(shù)據(jù)量的增加神經(jīng)網(wǎng)絡(luò)算法勢必能夠為其他更深入的科學工作提供更可靠的實驗數(shù)據(jù)分析.

日前，在PDB數(shù)據(jù)庫[9]中儲存的蛋白質(zhì)-DNA復合體三維結(jié)構(gòu)數(shù)據(jù)日益增加，這些結(jié)構(gòu)數(shù)據(jù)也為功能殘基位點預(yù)測提供了大量的數(shù)據(jù)來源.在生物信息學中我們可以在這些數(shù)據(jù)庫中挖掘提煉已有數(shù)據(jù)，對其進行不同特點和功能的分類和統(tǒng)計分析，從而得到潛在的結(jié)構(gòu)或者序列共同點，并以此為基礎(chǔ)整理成自動化算法對未知蛋白的功能殘基位點進行預(yù)測.在本論文當中我們共設(shè)計了兩種預(yù)測蛋白質(zhì)功能殘基的計算方法：基于相似序列的一級結(jié)構(gòu)預(yù)測、基于機器學習算法的支持向量機(SVM)預(yù)測.并對其結(jié)果進行歸一化整合，得到最佳的預(yù)測結(jié)果.

2 材料與方法

2.1 數(shù)據(jù)集構(gòu)建

我們使用了PDNA_62和PDNA_224數(shù)據(jù)集來評估我們的預(yù)測能力，兩個數(shù)據(jù)集的來源如下：

PDNA_62：PDNA_62是一個由Ahma and Sarai建立的經(jīng)典非冗余蛋白質(zhì)-DNA復合物數(shù)據(jù)集[10]，主要集中了Protein Data Bank (PDB)[9]數(shù)據(jù)庫中分辨率高于3.5 ?的蛋白質(zhì)-DNA復合物，并去除序列相似度高于25%的冗余序列.PDNA_62數(shù)據(jù)集中共包含1 215個DNA相互作用位點，以及6 948個非DNA相互作用位點.

PDNA_224：根據(jù)最新版本的PDB數(shù)據(jù)庫(2019年7月4日)，檢索DNA結(jié)合蛋白，并將X晶體衍射分辨率設(shè)置為3.0 ?，這樣從數(shù)據(jù)庫里面得到978個蛋白質(zhì)-DNA復合物.利用PISCES software過濾[11]對候選復合物序列進行過濾，并剔除相似度高于25%的冗余序列，再去除與PDNA_62中的同源序列，最終得到224條非冗余蛋白質(zhì)序列.其中共包含3 778個DNA相互作用位點和53 570個非DNA相互作用位點.

2.2 基于SVM的預(yù)測器

由于蛋白質(zhì)中每個殘基的特征可以用一個讀碼窗來描述，讀碼窗由中心殘基(central residue)及其n個相鄰殘基的位置組合而成，若相鄰殘基不存在，則用0來表示.假設(shè)讀碼窗的長度取L，每個中心殘基擁有P=24個特征數(shù)，那么預(yù)測蛋白質(zhì)-DNA相互作用位點的特征總數(shù)是L×P=L×24，最后將這些特征量格式化為支持向量機的輸入數(shù)據(jù)，得到預(yù)測結(jié)果.我們所利用的是公用的LibSVM軟件，該軟件從http://www.csie.nut.edu.tw/～cjlin/libsvm[12]下載得到.LibSVM將其輸出結(jié)果轉(zhuǎn)化成條件概率通過使sigmoid函數(shù)，sigmoid函數(shù)的定義如下：

其中x是支持向量機的輸入特征，sdv(x)代表輸入特征x的閾值，P(Y=1|x)是條件概率結(jié)果，A和B分別是sigmoid函數(shù)的斜率和偏移參數(shù).官方文檔推薦A=-2.0，B=-0.5.此外，由于我們的數(shù)據(jù)集中正負樣本差異過大，因此我們使用了隨機過抽樣(over-sampling)分析，并使用LibSVM的RBF(radial basis function)內(nèi)核進行超平面分類器進行建模，以期盡可能消除樣本數(shù)量對SVM帶來的偏差.

為了是其預(yù)測結(jié)果能夠更好地與后續(xù)序列匹配預(yù)測器結(jié)果整合，我們將上述方程進行反函數(shù)處理，使最終結(jié)果為SVM預(yù)測閾值.

2.3 基于序列匹配的預(yù)測器

PDNA_62和PDNA_224中的序列文件保存為FASTA格式，但其并不能很好地反映序列之間的替換關(guān)系.因此利用PSI-BLAST，根據(jù)BLOSUM62替代矩陣對數(shù)據(jù)庫進行初始化打分，生成PSSM蛋白質(zhì)特征矩陣(圖1).

圖1 蛋白質(zhì)1TC3中C段肽鏈的PSSM特征矩陣Fig.1 PSSM matrix of C-chain in protein 1TC3

使用初始化后的PSSM數(shù)據(jù)文件進行計算.序列的提取與SVM中的信息窗口類似，使用一個長度為n的讀碼框(n=3、5、7、9、11……)，若兩邊數(shù)據(jù)不足則用0補齊.讀碼框在已知(B)和未知(A)序列中同時滑動，并在相應(yīng)位點進行計算.相應(yīng)的計算公式我們使用了鋅指結(jié)構(gòu)預(yù)測中CHDEs氨基酸得分的計算公式[13].

Score(α,β)=

此公示中α,β表示兩條序列，j表示長度為n的讀碼框中第j位點的匹配的分，i則表示每個位點的氨基酸與20個標準氨基酸進行匹配，αni、βni則表示這些氨基酸的匹配的分，關(guān)于單個位點與20個氨基酸得分對應(yīng)關(guān)系，依舊沿用BLOSUM62氨基酸得分矩陣，具體的計算方法如下.Pi則代表氨基酸i出現(xiàn)的背景頻率.最終將所有讀碼框的得分相加，總得分即能夠代表功能未知序列(A)和功能已知序列(B)之間的相似關(guān)系.

αi,j=eMi,jλupi

其中Mi,j即為BLOSUM62矩陣中兩氨基酸對應(yīng)得分，λμ即為PSSM格式化后的α文件中standard ungapped Lambda value.

由于此方法得到大量得分結(jié)果，最理想的情況中，我們希望這些整理后的坐標呈強烈的線性關(guān)系，這樣才能夠進一步說明兩個片段之間擁有強烈的相關(guān)性.因此我們以單個蛋白質(zhì)肽鏈前x%、整個訓練集得分的第y%作為閾值，由于皮爾遜系數(shù)可以用來表示坐標點的相關(guān)性，因此我們使用皮爾遜相關(guān)性約束作為參數(shù)，同時調(diào)整窗口長度(n=7,9,11,……)，建立模型進行訓練，并使用訓練結(jié)果最優(yōu)的模型進行得分篩選，得到初步的待處理匹配位置，并將匹配位置整理成坐標形式.篩選模型可分為四類：

type1：用單個蛋白質(zhì)肽鏈前x%最為閾值，并進行皮爾遜相關(guān)系數(shù)的約束；

type2：用整體蛋白質(zhì)肽鏈前y%作為閾值，并進行皮爾遜相關(guān)系數(shù)的約束；

type3：用單個蛋白質(zhì)肽鏈前x%最為閾值，不進行皮爾遜相關(guān)系數(shù)的約束；

type4：用整體蛋白質(zhì)肽鏈前y%作為閾值，不進行皮爾遜相關(guān)系數(shù)的約束.

最終我們將確定模型篩選出來的位置坐標進行線性匹配，以求找到更多的坐標點位于同一條直線.那么這些位于同一直線的坐標點所代表的序列位點，即可視為擁有強烈的序列相似性和功能相關(guān)性.

圖2 不同坐標點的皮爾遜系數(shù)[14]Fig.2 Pearson coefficient of different dataset[14]

2.4 對預(yù)測結(jié)果進行歸一化整合

由于預(yù)測結(jié)果各異，且有著各自的優(yōu)缺點，因此我們使用數(shù)學期望對兩種預(yù)測器進行數(shù)據(jù)整合，根據(jù)數(shù)學期望公式，我們可整合支持向量機打分sdv(x)與序列匹配打分Score(α,β)，具體函數(shù)如下：

最終，對長度為L的α序列，其于β序列的結(jié)合得分為E(x).

3 結(jié)果分析

PdDNA，即本論文所研究得出的蛋白質(zhì)-DNA相互作用位點位點預(yù)測方法，主要整合了SVM支持向量機打分和序列匹配度打分.在SVM模型訓練過程中，我們使用PDNA_62標準實驗數(shù)據(jù)集進行測試，下表為五交叉檢驗結(jié)果隨讀碼窗口長度變化所得到的預(yù)測信息.我們的評估指標包括敏感度(Sensitivity，SN)、特異性(Specificity，SP)、強度(Strength)、準確度(ACC)和馬太相關(guān)系數(shù)(MCC)和，其具體計算方法如下.最終確定長度11為最佳窗口長度(表1).

SN=TP/(TP+FN)

SP=TN/(TN+FP)

Strength=(SN+SP)/2

Acc=(TP+TN)/(TP+FP+TN+FN)

MCC=(TP×TN-FP×FN)/

表1 PDNA_62數(shù)據(jù)集預(yù)測精度隨SVM窗口長度變化情況Tab.1 The prediction result of the PDNA_62 dataset with different window sizes of the SVM

對于數(shù)據(jù)集PDNA_62，在基于PdDNA運算模型的測試中，其準確度(ACC)為85.15%，馬太相關(guān)系數(shù)(MCC)0.55為，強度(Strength)為85.89%，正集預(yù)測成功率為84.91%，負集預(yù)測成功率為86.87%.與其他主流預(yù)測方法對比結(jié)果如下(表2).同時，最終預(yù)測結(jié)果的ROC曲線中(圖3)也可明顯看出，結(jié)合SVM與序列匹配的預(yù)測準確率比單獨的SVM預(yù)測高出5%～10%不等.

表2 對于PDNA_62數(shù)據(jù)集，PdDNA預(yù)測結(jié)果Tab.2 The prediction result of the PDNA_62 dataset of PdDNA algorithm

圖3 結(jié)合SVM預(yù)測器和基于序列匹配預(yù)測器，通過ROC曲線展示數(shù)據(jù)集PDNA_62的最終的預(yù)測結(jié)果Fig.3 ROC curves for the DNA-binding sites prediction in PDNA_62 dataset by combining SVM predictor with sequence-based predictor

同時，對于自建庫PDNA_224也擁有更好的預(yù)測情況(表3).從ROC曲線(圖4)中可以看出，其SVM與序列匹配整合結(jié)果也比單獨的SVM預(yù)測略有提高，但相比PDNA_62數(shù)據(jù)集提升幅度較小，這可能是由于我們的正負集樣本數(shù)量差距進一步擴大，進而造成了SVM分類器出現(xiàn)偏移造成的.

表3 對于PDNA_224數(shù)據(jù)集，PdDNA預(yù)測結(jié)果Tab.3 The prediction result of the PDNA_224 dataset of PdDNA algorithm

圖4 結(jié)合SVM預(yù)測器和基于序列匹配預(yù)測器，通過ROC曲線展示數(shù)據(jù)集PDNA_224的最終的預(yù)測結(jié)果Fig.4 ROC curves for the DNA-binding sites prediction in PDNA_224 dataset by combining SVM predictor with sequence-based predictor

4 討 論

在本論文中，我們通過SVM支持向量機方法、序列匹配算法的融合，得到了一種全新的用來預(yù)測蛋白質(zhì)-DNA相互作用位點的實驗方法—PdDNA.并根據(jù)這個方法對標準實驗數(shù)據(jù)庫PDNA_62進行了數(shù)據(jù)分析，結(jié)果表明我們的PdDNA程序所得到的預(yù)測結(jié)果為：ACC，85.15%；MCC，0.55；Strength，85.89%；SP，84.91%；SN：86.87%.其預(yù)測準確率高于目前蛋白質(zhì)預(yù)測領(lǐng)域任何其他算法.其高達86%以上的預(yù)測準確度也證明了我們的預(yù)測方法具有高度的可用性和實用性，能夠給蛋白質(zhì)功能組學、分子生物學等其他領(lǐng)域的科學研究帶來重要的預(yù)測參考，進一步提高科研工作的效率

同時，我們還建立了具有參考意義的PDNA_224蛋白質(zhì)-DNA相互作用位點數(shù)據(jù)庫.對于自建數(shù)據(jù)庫的預(yù)測效果也令人滿意，其最好的ACC為83.07，MCC為0.42，Str為83.05.以期對他人研究和算法設(shè)計提供有價值的數(shù)據(jù)來源.

但我們的數(shù)據(jù)集中仍存在正負樣本不平衡的問題，由于客觀實驗限制我們很難在數(shù)據(jù)庫中找到足量的正集數(shù)據(jù).正集樣本過少的情況下，其樣本無法廣泛分布，因此SVM預(yù)測的分類器確會朝向正集偏移，導致其敏感度下降.根據(jù)LibSVM的官方文檔[18]中的說明，在特征數(shù)遠小于樣本數(shù)的情況下建議使用RBF內(nèi)核.在RBF內(nèi)核中LibSVM會啟用超平面分類器，并啟用了隨機過抽樣方法，以此對抗樣本數(shù)量不平衡的問題，這也是本文中采用的解決方法.但這種方法對計算量需求極大.例如在PDNA_224數(shù)據(jù)集中，負集樣本為正集樣本的10倍，此時負集樣本被隨機抽樣為10份并分別與正集樣本進行建模，最終對模型進行擬合.雖精度較高但計算時間會陡增十倍.若在計算資源不足的情況下，也可以考慮手動設(shè)置懲罰因子C，對分類結(jié)果進行加權(quán).仍以PDNA_224數(shù)據(jù)集舉例，正負數(shù)據(jù)集樣本比例為1∶10，則懲罰因子C比例可定位10∶1，在LibSVM中的參數(shù)則應(yīng)設(shè)置為-w1 10;-w-1 1，則可在節(jié)省計算資源的情況下得到大致相同的分類結(jié)果.缺點則是會偏離原始數(shù)據(jù)的概率分布，因而本文中沒有采用.