PCA方法在蛋白質(zhì)亞細(xì)胞定位中應(yīng)用

馬 軍 偉， 史 舵， 顧 宏， 張 杰

(1.大連理工大學(xué) 控制科學(xué)與工程學(xué)院，遼寧 大連 116024；2.山西省電力公司 電力通信中心，山西 太原 030001；3.安徽工業(yè)大學(xué) 數(shù)理學(xué)院，安徽 馬鞍山 243002)

0 引 言

蛋白質(zhì)功能的研究是蛋白質(zhì)工程的重要環(huán)節(jié)，它對(duì)人們生活生產(chǎn)有著重要的意義，廣泛應(yīng)用于食品安全、新藥研制、工業(yè)生產(chǎn)等領(lǐng)域，并在科學(xué)進(jìn)步和社會(huì)發(fā)展中扮演推進(jìn)劑的角色.如何更好地分析研究蛋白質(zhì)的功能成為當(dāng)前的主要問題.研究發(fā)現(xiàn)，蛋白質(zhì)的功能與它的亞細(xì)胞定位密切相關(guān)，可以根據(jù)蛋白質(zhì)的亞細(xì)胞位置來推斷蛋白質(zhì)的功能.傳統(tǒng)蛋白質(zhì)亞細(xì)胞定位預(yù)測(cè)方法主要是生物實(shí)驗(yàn)方法，如細(xì)胞分餾法、熒光顯微法等.但是此類方法一般費(fèi)用較高，而且比較費(fèi)時(shí)[1].隨著生物數(shù)據(jù)庫中蛋白質(zhì)序列數(shù)量的急劇膨脹，通過實(shí)驗(yàn)的方法獲得蛋白質(zhì)的亞細(xì)胞位置信息越來越不現(xiàn)實(shí).因此，急需發(fā)展計(jì)算方法對(duì)蛋白質(zhì)進(jìn)行亞細(xì)胞定位預(yù)測(cè).

序列編碼技術(shù)是蛋白質(zhì)亞細(xì)胞定位預(yù)測(cè)技術(shù)的基礎(chǔ)，一般來說，按照基于氨基酸組成(amino acid composition，AAC)的編碼方法，蛋白質(zhì)將被表示為一個(gè)20維的特征向量，每一維對(duì)應(yīng)一種氨基酸，以這種氨基酸在蛋白質(zhì)序列中出現(xiàn)的頻率為該維元素的值.然而，若用此種方法表示蛋白質(zhì)，氨基酸之間的相對(duì)位置及序列長(zhǎng)度都將遺失，從而造成表示上的固有局限性.為了解決此問題，Chou[2]提出了一個(gè)更高級(jí)的蛋白質(zhì)表示模型——偽 氨 基 酸 組 成 (pseudo－amino acid composition，PseAAC)，PseAAC 的應(yīng)用顯著提高了預(yù)測(cè)精度.PseAAC包含20＋λ個(gè)成分，其中前20個(gè)成分按照AAC蛋白質(zhì)序列方法編碼，后λ個(gè)成分代表λ個(gè)不同級(jí)別的序列次序相關(guān)因子，這些離散的數(shù)列能近似地表示蛋白質(zhì)的序列次序效應(yīng)，從而提高了預(yù)測(cè)精度.近年來，PseAAC被廣泛應(yīng)用于生物預(yù)測(cè)模型[3、4].

一般來說，λ的值越大，包含的序列效益越大，但λ的值不能超過蛋白質(zhì)的長(zhǎng)度.若λ的值過大，則會(huì)造成如冗余和溢出之類的統(tǒng)計(jì)預(yù)測(cè)錯(cuò)誤.因此，對(duì)于不同的訓(xùn)練數(shù)據(jù)集，應(yīng)有不同的最優(yōu)λ值.若數(shù)據(jù)集中最短蛋白質(zhì)鏈的氨基酸殘基數(shù)比較大，那么選擇合適的λ值會(huì)比較耗時(shí).本文在此情況下用主成分分析法提取主特征來解決這個(gè)問題.

1 理論分析

1.1 PseAAC

蛋白質(zhì)的AAC特征表示由20個(gè)元素組成，分別代表20種不同氨基酸在蛋白質(zhì)序列中出現(xiàn)的頻率.偽氨基酸組成不但包含這些信息，還包含一些其他成分，通過這些成分近似反映蛋白質(zhì)的序列順序效應(yīng).

假設(shè)蛋白質(zhì)鏈有L個(gè)氨基酸殘基:

通過一系列序列順序相關(guān)因子可以近似地反映序列次序效應(yīng)，相關(guān)因子的定義如下:

其中λ＜L.θ1稱為第一級(jí)相關(guān)因子，反映蛋白質(zhì)序列中相鄰氨基酸相關(guān)性；θ2稱為第二級(jí)相關(guān)因子，反映蛋白質(zhì)序列中所有每間隔一個(gè)氨基酸的相關(guān)性；θλ稱為第λ級(jí)相關(guān)因子，反映蛋白質(zhì)序列中所有每間隔λ－1個(gè)氨基酸的相關(guān)性(圖1).其相關(guān)函數(shù)Ci，j定義為

其中H1(Rj)、H2(Rj)和M(Rj)分別是Rj的疏水性值、親水性值和側(cè)鏈氨基酸質(zhì)量.然后取代普通的基于氨基酸編碼的方法，蛋白質(zhì)X序列可以通過下式表示:

式中

其中fi表示蛋白質(zhì)X中氨基酸的出現(xiàn)頻率，θj表示蛋白質(zhì)X第j級(jí)序列次序效應(yīng)相關(guān)因子；w是序列次序效應(yīng)的權(quán)重因子，在這里取w＝0.05；λ表示此模型中采用的相關(guān)因子類型數(shù)量，且0≤λ＜L.當(dāng)λ＝0時(shí)，說明模型中沒有了能反映序列次序效應(yīng)的相關(guān)因子，PseAAC便退化為普通的AAC模型；當(dāng)0＜λ＜L，從式(3)、(4)中可以看出，前20個(gè)分量反映氨基酸組成效應(yīng)，后λ個(gè)分量反映序列次序效應(yīng).由于λ個(gè)相關(guān)因子同序列次序效應(yīng)緊密相關(guān)，也就是說λ越大，所包含的序列信息就越多，但是λ也有上限，必須小于蛋白質(zhì)的氨基酸殘基數(shù)目.此外，λ過大也有可能降低蛋白質(zhì)的聚類性能，從而影響分類正確率.因此，對(duì)于一個(gè)給定的數(shù)據(jù)集，必存在一個(gè)最優(yōu)的λ值.如果通過試驗(yàn)法獲得，需要多次反復(fù)地試驗(yàn)才能找到最優(yōu)值，將會(huì)特別費(fèi)時(shí)費(fèi)力.本文采用一種新方法來解決此問題.

圖1 PseAAC示意圖Fig.1 Schematic drawing of PseAAC

1.2 PCA

在PseAAC中前20個(gè)分量反映了蛋白質(zhì)的組成信息，當(dāng)其中一個(gè)量較大時(shí)同時(shí)會(huì)影響其他量的大小，后λ個(gè)分量反映了序列長(zhǎng)度及次序效應(yīng)，所以具有一定的相關(guān)性，從而說明數(shù)據(jù)在一定程度上有信息的重疊.主成分分析(PCA)采用一種降維的方式，找出幾個(gè)綜合因子來代表原來眾多的特征，使這些綜合因子盡可能地反映原來變量的信息，而且彼此之間互不相關(guān)，只研究樣本特征中少數(shù)幾個(gè)能最大程度保留原始特征變化方面信息的特征組合[5].這樣也避免了求取最優(yōu)λ值所帶來的低效性.

假設(shè)所討論問題有n個(gè)指標(biāo)，可以看成n個(gè)隨機(jī)變量，記為(X1X2…Xn)，主成分分析的實(shí)質(zhì)是線性組合這n個(gè)指標(biāo)構(gòu)成新指標(biāo)

同時(shí)滿足

(1)每一個(gè)主成分系數(shù)平方和為1，

(2)主成分之間相互獨(dú)立，

(3)主成分的方差依次遞減，即重要性依次遞減，

在求取時(shí)，可設(shè)X＝ (X1X2…Xn)T，新指標(biāo)Y＝ (Y1Y2…Yn)T.則要在保持主成分之間相互獨(dú)立且方差依次遞減的原則下，找到變換矩陣U，有

可設(shè)X的協(xié)方差矩陣為根據(jù)方差的性質(zhì)可知，ΣX是非負(fù)對(duì)稱矩陣，所以必 存 在 正 交 矩 陣V， 使 得VTΣXV＝中(λ1λ2…λn)是ΣX的特征向量，可以假設(shè)λ1≥λ2≥…≥λn，V正好是特征值對(duì)應(yīng)的特征向量矩陣.可以證明V便是所要求的U，且各個(gè)主成分的方差var(Yi)＝λi(i＝1，2，…，n)[5].

所以，第i個(gè)主分量的方差貢獻(xiàn)率可以定義為反映了樣本數(shù)據(jù)的信息變化情況.前m(m＜n)個(gè)主分量的累積方差貢獻(xiàn)率定義為，可以選取前m個(gè)主分量使其累積方差貢獻(xiàn)率達(dá)到一定的要求(如85%～95%)，這樣便可達(dá)到降低原始數(shù)據(jù)維數(shù)的目的.

2 實(shí)驗(yàn)準(zhǔn)備

2.1 分類器

為了說明本文方法的有效性，選用兩個(gè)比較常用的分類器:k近鄰(k－NN)算法及反向傳播網(wǎng)絡(luò)(back－propagation network，BP network).

k－NN算法是一種比較簡(jiǎn)單成熟的分類算法，但其良好的分類性能并不弱于其他較復(fù)雜算法，已廣泛應(yīng)用于蛋白質(zhì)二級(jí)結(jié)構(gòu)預(yù)測(cè)及亞細(xì)胞定位.這里k取3，它被認(rèn)為比1－NN有更好的可解釋性[6]，·2用來測(cè)量樣本之間的相似度.

BP網(wǎng)絡(luò)學(xué)習(xí)算法對(duì)于逼近實(shí)數(shù)值、離散值和向量值的目標(biāo)函數(shù)有很強(qiáng)的健壯性，已被成功應(yīng)用于人臉識(shí)別、視覺導(dǎo)航和生物信息學(xué)等領(lǐng)域.本文選用具有8個(gè)神經(jīng)元的單隱藏層網(wǎng)絡(luò)結(jié)構(gòu)，基于變學(xué)習(xí)率的后向傳播算法用來更新權(quán)值.最大訓(xùn)練次數(shù)設(shè)為300，同時(shí)為了防止過擬合，在每一次訓(xùn)練時(shí)都會(huì)估計(jì)泛化誤差，若連續(xù)5次訓(xùn)練誤差沒有變化則訓(xùn)練會(huì)被提前終止[7].

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

本文選用2個(gè)不同的數(shù)據(jù)集來驗(yàn)證算法性能，分別進(jìn)行PCA優(yōu)化并觀察優(yōu)化之后的結(jié)果相比優(yōu)化之前結(jié)果的改進(jìn)程度.

第一個(gè)數(shù)據(jù)集由Chen等[8]構(gòu)造，含有315個(gè)蛋白質(zhì)序列(去掉兩個(gè)已不用的蛋白質(zhì)序列)，具有6類亞細(xì)胞:細(xì)胞質(zhì)蛋白質(zhì)序列110個(gè)，細(xì)胞質(zhì)膜蛋白質(zhì)序列55個(gè)，線粒體蛋白質(zhì)序列34個(gè)，分泌蛋白質(zhì)序列17個(gè)，細(xì)胞核蛋白質(zhì)序列52個(gè)，內(nèi)質(zhì)網(wǎng)蛋白質(zhì)序列47個(gè).第二個(gè)數(shù)據(jù)集由Gardy等[9]構(gòu)造，已廣泛應(yīng)用于蛋白質(zhì)亞細(xì)胞定位中.其含有541個(gè)蛋白質(zhì)序列，具有4類亞細(xì)胞:細(xì)胞質(zhì)蛋白質(zhì)序列194個(gè)，細(xì)胞質(zhì)膜蛋白質(zhì)序列103個(gè)，細(xì)胞壁蛋白質(zhì)序列61個(gè)，細(xì)胞外蛋白質(zhì)序列183個(gè).為方便使用，這兩個(gè)數(shù)據(jù)庫分別記作CH315和GA541.

2.3 精度評(píng)測(cè)

采用蛋白質(zhì)亞細(xì)胞定位研究中常用的5折交叉驗(yàn)證法.具體做法如下:將數(shù)據(jù)集均分成5等份，選擇其中4個(gè)子集作為訓(xùn)練集，然后選擇第5個(gè)子集作為測(cè)試集，重復(fù)5次保證每一個(gè)子集都擔(dān)任過測(cè)試集.實(shí)際上在神經(jīng)網(wǎng)絡(luò)實(shí)驗(yàn)中，用了5次5折交叉驗(yàn)證法來盡可能保證結(jié)果的合理性，這是因?yàn)樘荻认陆邓惴〞?huì)收斂到相對(duì)于網(wǎng)絡(luò)權(quán)值的局部極小值，而不是全局最小值.同時(shí)選用總體預(yù)測(cè)精度At作為評(píng)價(jià)指標(biāo)，它被看作判斷一個(gè)分類器好壞最重要的評(píng)價(jià)標(biāo)準(zhǔn).定義如下:

其中P(i)表示第i類蛋白質(zhì)序列被正確識(shí)別的數(shù)量，N是蛋白質(zhì)序列的總數(shù)量，k是類別數(shù)量.

3 實(shí)驗(yàn)結(jié)果

不同的數(shù)據(jù)集中最短蛋白質(zhì)序列長(zhǎng)度是不同的，根據(jù)1.1節(jié)分析知λ必須小于最短序列長(zhǎng)度，而且λ越大所包含的序列順序效應(yīng)越大.在數(shù)據(jù)集CH315中蛋白質(zhì)序列長(zhǎng)度最短為87，在GA541中最短為40.所以在使用PCA對(duì)數(shù)據(jù)集進(jìn)行優(yōu)化時(shí)，λ取最大值，分別為86和39.相應(yīng)地，根據(jù)式(3)，其輸入向量維數(shù)分別是20＋86＝106和20＋39＝59.在實(shí)驗(yàn)中，累積方差貢獻(xiàn)率設(shè)置為90%以提取關(guān)鍵主成分.

圖2顯示的是k－NN算法在CH315與GA541數(shù)據(jù)集中的預(yù)測(cè)精度與λ之間的關(guān)系.為了增加對(duì)比效應(yīng)，PCA作用后的效果圖已畫出(虛線部分).可以看出預(yù)測(cè)精度隨著λ取值不同在不斷變化.在CH315數(shù)據(jù)集中，有3個(gè)最優(yōu)λ值，分別是73、74和75，對(duì)應(yīng)精度約為73.80%.在GA541數(shù)據(jù)集中，最優(yōu)λ值是2，對(duì)應(yīng)精度約為82.25%.經(jīng)過PCA優(yōu)化后，預(yù)測(cè)精度分別達(dá)到74.15%和82.46%，顯然要高于最優(yōu)λ值對(duì)應(yīng)的精度.更重要的是，免去了求解最優(yōu)值所帶來的低效率.

圖2 k－NN算法在兩個(gè)數(shù)據(jù)集上的總體預(yù)測(cè)精度Fig.2 Total prediction accuracy using k－NN algorithm on the two datasets

圖3 是BP神經(jīng)網(wǎng)絡(luò)在CH315與GA541數(shù)據(jù)集上的預(yù)測(cè)精度與λ關(guān)系圖(虛線對(duì)應(yīng)著PCA作用后的情況).不僅給出了每種情況下5次交叉驗(yàn)證的均值，同時(shí)畫出了方差效果圖.可以看出PCA作用后均值明顯提升，同時(shí)方差變化范圍明顯減小，尤其在CH315數(shù)據(jù)集上表現(xiàn)更為突出，顯示了PCA的優(yōu)勢(shì).

圖3 BP網(wǎng)絡(luò)在兩個(gè)數(shù)據(jù)集上的總體預(yù)測(cè)精度Fig.3 Total prediction accuracy using BP network on the two datasets

4 結(jié)果討論

當(dāng)λ＝0時(shí)，PseAAC退化為AAC，從上述實(shí)驗(yàn)可以看出PseAAC確實(shí)要優(yōu)于AAC.但是，對(duì)于不同的λ此模型會(huì)產(chǎn)生不同的預(yù)測(cè)精度，為了得到最優(yōu)解，可以采用試驗(yàn)的方法逐一探求λ值，但是這種方法既繁瑣又低效.實(shí)際上，為了解決這一問題，Chou等[10]提出了集成方法，通過將不同λ值的PseAAC融合成一個(gè)整體，集成了210個(gè)獨(dú)立分類器，但是并沒有指出基底分類器的數(shù)目對(duì)輸出結(jié)果的影響.此外，集成如此大量的分類器難免產(chǎn)生較大的計(jì)算復(fù)雜度.本文采用主成分分析法，通過提取關(guān)鍵主特征來解決這一問題，試驗(yàn)結(jié)果顯示此方法確實(shí)可以提高蛋白質(zhì)亞細(xì)胞定位預(yù)測(cè)的準(zhǔn)確度.

進(jìn)一步，也嘗試PCA與其他分類器的結(jié)合.作者構(gòu)造了神經(jīng)網(wǎng)絡(luò)集成器定位亞細(xì)胞.在文獻(xiàn)[11]中，601個(gè)蛋白質(zhì)序列從革蘭陰性桿菌庫中提出，其中包括140個(gè)細(xì)胞質(zhì)序列，74個(gè)細(xì)胞外序列，280個(gè)內(nèi)膜序列，38個(gè)外膜序列，69個(gè)周質(zhì)序列.最短序列長(zhǎng)度為29，所以在PCA過程中λ設(shè)為28，90%的累積方差貢獻(xiàn)率依然被用于提取關(guān)鍵主成分.表1顯示了分類器的性能比較，同時(shí)，為了增加對(duì)比度，性能強(qiáng)大的PSORTb分類器預(yù)測(cè)結(jié)果也已給出.從表中可以看出，經(jīng)過PCA優(yōu)化后，各個(gè)類的精度都是最高的，總體預(yù)測(cè)精度達(dá)到82.0%，比PSORTb分類器和集成學(xué)習(xí)器各高16.6%和5.0%，說明PCA確實(shí)抓住了蛋白質(zhì)樣本的主特征.

表1 3種不同方法的預(yù)測(cè)結(jié)果Tab.1 Predicted results by the three different methods

5 結(jié) 語

本文試驗(yàn)中用了90%的累積方差貢獻(xiàn)率，當(dāng)然這不是最優(yōu)選擇，在其他情況下也可以選擇其他值，一般來說，大于等于85%的累積方差貢獻(xiàn)率是合理的.另外，PCA只是一種線性變換，當(dāng)數(shù)據(jù)高度復(fù)雜的時(shí)候效果并不明顯，有時(shí)還有可能降低預(yù)測(cè)性能，下一步準(zhǔn)備嘗試kernel PCA(KPCA)、kernel independent component analysis(KICA)等非線性變換來替換主成分分析，對(duì)提升分類器的預(yù)測(cè)性能有一定幫助.

[1]SHEN Hong－bin， CHOU Kuo－chen. Predicting protein subnuclear location with optimized evidencetheoreticK－nearest classifier and pseudo amino acid composition [J]. Biochemical and Biophysical Research Communications，2005，337(3):752－756

[2]CHOU Kuo－chen. Prediction of protein cellular attributes using pseudo－amino acid composition [J].Proteins:Structure，F(xiàn)unction，and Bioinformatics，2001，43(3):246－255

[3]DING Y S，ZHANG T L.Using Chou′s pseudo amino acid composition to predict subcellular localization of apoptosis proteins:an approach with immune genetic algorithm－based ensemble classifier[J].Pattern Recognition Letters，2008，29(13):1887－1892

[4]ZENG Y，GUO Y，XIAO R，etal.Using the augmented Chou′s pseudo amino acid composition for predicting protein submitochondria locations based on auto covariance approach [J].Journal of Theoretical Biology，2009，259(2):366－372

[5]李弼程，邵美珍，黃 潔.模式識(shí)別原理與應(yīng)用[M].西安:西安電子科技大學(xué)出版社，2008

[6]VEENMAN C，REINDERS M.The nearest subclass classifier:A compromise between the nearest mean and nearest neighbor classifier [J]. IEEE Transactions on Pattern Analysis and Machine Intelligence，2005，27(9):1417－1429

[7]ZHOU Z H，WU J，TANG W.Ensembling neural networks:Many could be better than all [J].Artificial Intelligence，2002，137(1－2):239－263

[8]CHEN Y L，LI Q Z.Prediction of the subcellular location of apoptosis proteins [J].Journal of Theoretical Biology，2007，245(4):775－783

[9]GARDY J，LAIRD M，CHEN F，etal.PSORTb v.2.0:expanded prediction of bacterial protein subcellular localization and insights gained from comparative proteome analysis [J].Bioinformatics，2005，21(5):617－623

[10]CHOU Kuo－chen，SHEN Hong－bin.Large－scale predictions of gram－negative bacterial protein subcellular locations [J].Journal of Proteome Research，2006，5(12):3420－3428

[11]MA J W，LIU W Q，GU H.Predicting protein subcellular locations for gram negative bacteria using neural networks ensemble [C]／／ Proceedings of CIBCB′2009.Piscataway:IEEE，2009:114－120