王皆恒，李 校

(四川大學(xué)生命科學(xué)學(xué)院 四川省分子生物與生物技術(shù)重點(diǎn)實(shí)驗(yàn)室，成都610064)

1 引 言

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

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

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

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

2 材料與方法

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

我們使用了PDNA_62和PDNA_224數(shù)據(jù)集來(lái)評(píng)估我們的預(yù)測(cè)能力，兩個(gè)數(shù)據(jù)集的來(lái)源如下：

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

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

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

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

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

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

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

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

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

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

Score(α,β)=

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

αi,j=eMi,jλupi

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

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

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

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

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

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

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

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

2.4 對(duì)預(yù)測(cè)結(jié)果進(jìn)行歸一化整合

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

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

3 結(jié)果分析

PdDNA，即本論文所研究得出的蛋白質(zhì)-DNA相互作用位點(diǎn)位點(diǎn)預(yù)測(cè)方法，主要整合了SVM支持向量機(jī)打分和序列匹配度打分.在SVM模型訓(xùn)練過(guò)程中，我們使用PDNA_62標(biāo)準(zhǔn)實(shí)驗(yàn)數(shù)據(jù)集進(jìn)行測(cè)試，下表為五交叉檢驗(yàn)結(jié)果隨讀碼窗口長(zhǎng)度變化所得到的預(yù)測(cè)信息.我們的評(píng)估指標(biāo)包括敏感度(Sensitivity，SN)、特異性(Specificity，SP)、強(qiáng)度(Strength)、準(zhǔn)確度(ACC)和馬太相關(guān)系數(shù)(MCC)和，其具體計(jì)算方法如下.最終確定長(zhǎng)度11為最佳窗口長(zhǎng)度(表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ù)測(cè)精度隨SVM窗口長(zhǎng)度變化情況Tab.1 The prediction result of the PDNA_62 dataset with different window sizes of the SVM

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

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

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

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

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

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

4 討 論

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

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

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