基于黎曼流形的蛋白質(zhì)三維結構數(shù)據(jù)相似性比較

徐永紅，褚澤斐，洪文學(燕山大學電氣工程學院，河北秦皇島066004)

基于黎曼流形的蛋白質(zhì)三維結構數(shù)據(jù)相似性比較

徐永紅*，褚澤斐，洪文學
(燕山大學電氣工程學院，河北秦皇島066004)

摘要:以NMR技術為代表的海量蛋白質(zhì)空間結構數(shù)據(jù)為現(xiàn)代生命科學研究提供了前所未有的機遇，但后續(xù)的大數(shù)據(jù)分析卻成為一大難題。充分利用已知的蛋白質(zhì)三維結構信息來預測未知的蛋白質(zhì)空間結構信息是研究蛋白質(zhì)結構和功能關系一種重要手段。本文提出一種基于黎曼流形的蛋白質(zhì)三維結構相似性比較新方法。該方法通過構建Cα坐標系和提取蛋白質(zhì)結構具有旋轉(zhuǎn)和平移不變性的幾何特征量，將蛋白質(zhì)的三維坐標序列轉(zhuǎn)換為一維序列，采用黎曼距離作為三維結構相似度指標。本方法不需要對蛋白質(zhì)結構做旋轉(zhuǎn)和平移變換，避免了主流的RMSD方法中兩蛋白質(zhì)通過最小二乘擬合進行配準時產(chǎn)生的誤差，并且完全不依賴于一級結構序列信息，對不具備序列相似性的蛋白質(zhì)之間的相似性比較具有現(xiàn)實意義。本文分別針對不同相似度的蛋白質(zhì)、Fischer提出的10個較難識別的蛋白質(zhì)結構對、HOMSTRAD數(shù)據(jù)庫中的700個數(shù)據(jù)這3組數(shù)據(jù)，對本文算法進行了驗證。實驗結果表明，與其他方法相比，本文方法的匹配精度均得到了較大提升。

關鍵詞:蛋白質(zhì);三次樣條插值; Cα坐標系;黎曼流形;結構比較

0　引言

生物大分子中的蛋白質(zhì)和很多非編碼RNA的功能主要取決于它們的空間結構。到目前為止，已經(jīng)有超過六萬個生物大分子的空間結構被測定，如何有效地比較它們之間的相似性成了生命科學中的一個重要課題［1］。蛋白質(zhì)三維結構的直接獲取一直是一個瓶頸問題，盡管蛋白質(zhì)序列的測定已基本完成，但大量序列已知的蛋白質(zhì)的三維結構尚未被實驗方法測定出來，在這種情況下，充分利用一級序列信息和已知蛋白質(zhì)的空間結構信息來研究預測未知蛋白質(zhì)的空間結構，為結構生物學中研究蛋白質(zhì)結構和功能關系的主要手段。由于用實驗方法得到目的蛋白的結構困難較大，而且需要較長時間，蛋白質(zhì)三維結構相似性比較可構建實驗需要的結構模型，用于提出關于蛋白質(zhì)功能的假設并指導進一步的實驗工作［2］，因此蛋白質(zhì)結構相似性比較是蛋白質(zhì)功能分析、分類管理、檢測等重要方法之一。

目前針對蛋白質(zhì)結構比較已經(jīng)有很多的研究方法與軟件工具，如Dali［3］、CE［4］、VAST［5］、STRUCTAL［6］、SSM［7］、TM-align［8］等，為蛋白質(zhì)空間結構的研究提供了多種分析手段。蛋白質(zhì)結構比較方法主要分為3類:基于氨基酸間距離矩陣的比較(Dali、CE) ;基于蛋白質(zhì)空間幾何結構的比較(STRUCTAL、TM-align) ;基于蛋白質(zhì)二級結構的匹配(VAST、SSM)。

傳統(tǒng)的蛋白質(zhì)相似性比較方法通常依賴于蛋白質(zhì)的一級結構序列，但是兩個一級結構序列不同而空間結構相同的蛋白質(zhì)往往有著相同的性質(zhì)。針對傳統(tǒng)方法的不足，本文從幾何角度出發(fā)，提出了一種基于黎曼流形的蛋白質(zhì)結構相似性比較的新方法。通過微分流形的處理，用曲率k和撓率r這兩個幾何量代替蛋白質(zhì)一級結構中和φ兩個二面角，從而把蛋白質(zhì)結構分析轉(zhuǎn)變?yōu)橐粋€純數(shù)學問題。

1982年，Louie和Somorjai［9］最初將微分流形引入到了蛋白質(zhì)結構研究中，他們將蛋白質(zhì)的主鏈看成是連續(xù)的空間曲線，并用參數(shù)(b，ω)描述(其中α螺旋對應螺旋面，β折疊對應懸鏈面，b表示曲面的尺寸，ω表示曲面的扭曲程度)。通過這種描述方法，Louie和Somorjai一方面給出了識別蛋白質(zhì)結構的指標;另一方面給出了蛋白質(zhì)主鏈結構的連續(xù)描述。在此基礎上，文獻［10-11］對蛋白質(zhì)結構三維曲線進行彈性形狀分析提取彈性度量，將彈性度量用于蛋白質(zhì)結構的比較中。文獻［12］進一步提出彈性黎曼度量，將蛋白質(zhì)結構從流形的角度進行比對，與其他結構比較方法相比效果更好，但是計算復雜度高。

本文通過構建Cα坐標系，對蛋白質(zhì)主鏈進行曲線插值提取幾何特征，用黎曼度量序列表示蛋白質(zhì)結構數(shù)據(jù)，從而將蛋白質(zhì)空間幾何結構的相關信息嵌入到黎曼度量序列表示的矩陣流形中。通過這種流形表示和嵌入，可以引入黎曼幾何的數(shù)學理論來提取有效的幾何結構特征。該方法能夠更直觀的體現(xiàn)出蛋白質(zhì)的結構特征，與傳統(tǒng)的方法相比幾何意義明確，為蛋白質(zhì)結構相似性比較提供了一種新的思路。

1　本文方法

本文所采用數(shù)據(jù)均來自PDB［13］(Protein Data Bank)數(shù)據(jù)庫，該數(shù)據(jù)庫是美國Brookhaven實驗室于1971年建立的大分子蛋白質(zhì)晶體結構資料數(shù)據(jù)庫，是互連網(wǎng)上唯一有生物高分子3維結構的數(shù)據(jù)庫，它的內(nèi)容主要是根據(jù)X射線結晶和核磁共振得到的實驗數(shù)據(jù)。

本文的關鍵是對蛋白質(zhì)主鏈進行曲線插值，從蛋白質(zhì)的三維幾何結構中提取具有空間平移和旋轉(zhuǎn)不變性的特征量，并將其嵌入到矩陣流形中。本文方法的具體步驟如下:

1)數(shù)據(jù)預處理

首先從數(shù)據(jù)庫中根據(jù)唯一標識的PDB ID獲取蛋白質(zhì)的PDB文件。忽略掉H原子，氨基酸的骨架是一個頂端為Cα原子的四面體，另外3個頂點分別為氨基、羧基和R基。圖1為蛋白質(zhì)1CRN中兩個氨基酸的基本原子框架，其中R基為CBSG，可以看到每個氨基酸的骨架原子為N-Cα-C的組合。從PDB文件中提取每個氨基酸的骨架原子——Cα、C、N原子的空間坐標。

2)構建Cα坐標系

弗萊納(Frenet-Serret)公式常用來描述粒子在連續(xù)可微的曲線上運動時曲線的切向、法向、副法方向之間的關系，本文在弗萊納公式的基礎上根據(jù)蛋白質(zhì)結構數(shù)據(jù)中Cα、C、N原子的空間坐標構建蛋白質(zhì)的Cα坐標系［14］:

其中，Cα、C、N分別代表Cα、C、N原子的空間坐標，X、Y、Z分別為列向量，因此得到3×3的單位正交矩陣構成的Cα坐標系序列［F1F2…Fi…Fn］，其中Fi= (XiYiZi)。每個Cα坐標系Fi均表征了對應氨基酸骨架的空間取向，每個氨基酸對應一個單位正交矩陣，則將蛋白質(zhì)氨基酸序列轉(zhuǎn)化為了單位正交矩陣序列。

圖1　氨基酸的基本結構Fig.1　Amino acid neighboring structure

3)曲線插值

三次樣條插值方法解決了計算幾何中自由型曲線曲面形狀的描述問題，在整體表示的同時保持了局部性質(zhì)，因此對蛋白質(zhì)Cα主鏈結構進行了三次樣條曲線插值。

三次樣條函數(shù)定義如下:

對于［a，b］上的分劃: a = x0＜x1＜…＜xn= b，則

利用三次樣條函數(shù)進行曲線插值，即取插值函數(shù)為三次樣條函數(shù)，稱為三次樣條插值。

如圖2所示，(b)和(c)分別為蛋白質(zhì)1CRN三次樣條插值前后的結構曲線圖，數(shù)字為氨基酸的序號，‘○’處為Cα原子所在位置。可以觀察到，經(jīng)插值后蛋白質(zhì)主鏈結構轉(zhuǎn)化為光滑曲線，兩個Cα原子間的每一段曲線都對應一個三次多項式。

圖2　蛋白質(zhì)1CRN主鏈結構曲線插值Fig.2　Curve interpolation of protein 1CRN backbone structure

4)幾何特征量的提取

曲率和撓率分別反映了空間曲線的彎曲程度和扭曲程度，它們具有空間旋轉(zhuǎn)平移不變性，能夠充分地描述曲線的幾何特征，因此本文選取蛋白質(zhì)結構曲線中每個Cα原子處的曲率、撓率和相鄰Cα原子間的距離作為特征描述子，對蛋白質(zhì)的幾何結構進行描述。

若已知空間曲線的參數(shù)方程

則在參數(shù)t處曲線的曲率kt和撓率τt計算公式分別為

其中，x′，x″，x′″分別為x對t的一、二、三階導，同理可知其他。

根據(jù)式(4)計算蛋白質(zhì)結構曲線上每個Cα點處的曲率k和撓率τ，并與相鄰Cα原子間的距離d的絕對值共同構成對角陣

其中，ki、τi為第i個Cα原子處的曲率和撓率，di，(xi， yi，zi)和(xi +1，yi +1，zi +1)分別為序號為i和i +1的Cα原子的空間坐標，這樣便得到了包含蛋白質(zhì)結構曲線特征描述子的對角陣序列。

5)構造黎曼度量序列

黎曼度量［15］指的是空間上的幾何學應基于無限臨近兩點(x1，x2，…，xn)與(x1+ d1，x2+ d2，…，xn+ dn)之間的距離，用微分弧長度平方所確定的正定二次型的理解度量，亦即是由函數(shù)構成的正定對稱矩陣，這便是黎曼度量。在三維空間中可以用一個3×3的實對稱矩陣來表示空間某點的黎曼度量:

式中，det(S)＞0，λi＞0，正定對稱矩陣的分解是可逆的。

根據(jù)式(2)中正定對稱矩陣的性質(zhì)，利用步驟(2)、(3)中得到的單位正交矩陣序列和對角陣序列構造黎曼度量序列:

Si為序號為i的氨基酸所對應的黎曼度量。

6)計算兩蛋白質(zhì)對應黎曼度量之間的黎曼距離

通過將黎曼度量S分解為單位正交矩陣F和對角陣Λ，再經(jīng)過計算Fi、Fi +1和Λi、Λi +1間的距離d(Fi，F(xiàn)i +1)和d(Λi，Λi +1)，最終得到Si、Si +2之間的黎曼距離d(Si，Si +1)

其中

λ1，i、λ2，i分別為對角陣Λi、Λi +1中第i行i列的值。

系數(shù)k(Λi，Λi +1)為加權因子［16］，范圍在(0～1)之間:

λmax、λmin分別為對角陣Λ中的最大值和最小值。

7)設定閾值，比較分析

由步驟(5)得到了兩蛋白質(zhì)間的黎曼距離序列［ds1，ds2，…，dsi，…，dsn］，黎曼距離dsi的值越小，表示兩個包含蛋白質(zhì)結構幾何信息的黎曼度量在黎曼流形上離得越近，故對應的蛋白質(zhì)結構越相似。經(jīng)多次實驗驗證后取黎曼距離的閾值:

k = 0.6mean(ds) + 0.3std(ds)，(10) mean(ds)、std(ds)分別表示黎曼距離序列的均值和方差。小于閾值k則認為在兩蛋白質(zhì)在殘基i處具有相同的結構。

2　實驗結果與分析

為驗證本文方法的可行性及有效性，進行了3組實驗:第1組在SCOP數(shù)據(jù)庫中分別選取家族、超家族、非超家族的蛋白質(zhì)數(shù)據(jù)進行實驗，驗證本文方法的比對結果是否符合生物意義;第2組對Fischer［17］給出的10個較難識別結構對進行實驗，從“相同殘基比”方面驗證本文方法的可行性;第3組基于HOMSTRAD結構比對參考庫，驗證本文方法在同源識別方面是否具有有效性。

2.1不同相似度的蛋白質(zhì)結構比對

蛋白質(zhì)結構分類數(shù)據(jù)庫SCOP［18］(Structural Classification of Proteins)是提供關于已知結構蛋白質(zhì)之間的結構和進化關系信息，所涉及的蛋白質(zhì)包括結構數(shù)據(jù)庫PDB中的所有條目。SCOP從總體上將蛋白質(zhì)分為全α型，全β型，以平行折疊為主的α/β型，以反平行折疊為主的α+β型等。從SCOP數(shù)據(jù)庫中選取PDB ID分別為101M、102M、1UVY、1C7Y的蛋白質(zhì)數(shù)據(jù)，他們分別屬于相同家族、相同超家族、不同超家族，相似度依次降低。通常，兩蛋白質(zhì)的關系越近，結構越相似度越高，黎曼度量序列間的黎曼距離越小。

表1為蛋白質(zhì)101M與101M、102M、1UVY、1C7Y的結構比對結果，Riemann_dis為黎曼度量序列間的黎曼距離的平均值，并與傳統(tǒng)蛋白質(zhì)結構比對評價指標RMSD［19］、TM-score［20］和Z-score［21］進行對比。RMSD(Root mean square deviation)方法為蛋白質(zhì)骨架直接配準的比較，該方法提出最早、使用最廣泛，RMSD值越小，表示兩蛋白質(zhì)結構越相似。與之趨勢相反的是TM-score與Z-score，值越大表示結構相似度越高。從表1中可以看出，隨著蛋白質(zhì)對之間關系的變遠，4個參數(shù)結果均表示蛋白質(zhì)結構相似度降低。圖3為4對蛋白質(zhì)間的黎曼距離曲線，從圖中觀察到隨著蛋白質(zhì)間的關系變遠，蛋白質(zhì)結構相似度降低，黎曼距離值明顯減小。

2.210個較難識別蛋白質(zhì)結構對的比對

在文獻［23］中Fischer給出了結構難以比較的10個蛋白質(zhì)對，利用本文方法對這10個蛋白質(zhì)對進行結構比對，并與Dali［3］、TM-align［8］、SPalign［22］這3種傳統(tǒng)方法進行了結果比較結果如表2所示。其中，“Equ”表示進行比較的殘基個數(shù)，“ResNum”表示相同的殘基個數(shù)，“Riemann_dis”表示本文方法中黎曼距離的平均值。通常認為，相同殘基比(Res-Num/Equ)越大，比對的效果越好。

表1　蛋白質(zhì)101M與不同相似度蛋白質(zhì)的結構比對結果Tab.1　Structure comparison results of different similarity of protein and protein 101M

圖3　蛋白質(zhì)101M與不同相似度的蛋白質(zhì)之間的黎曼距離曲線Fig.3　Curve of riemann distance between protein 101M and different similarity of proteins

如表2所示，本文方法的比例值高出效果最好的Dali方法約0.5個百分點，比效果較差的SPalign方法高出了約2.8個百分點，充分驗證了本文方法不僅可行，而且效果明顯優(yōu)于另外3種方法。

表2　10個較難識別蛋白質(zhì)結構對的比對平均結果比較Tab.2　Structure comparison average results of 10 more difficult to identify protein structures

2.3基于結構比對參考庫的驗證

HOMSTRAD(Homologous Structure Alignment Database)數(shù)據(jù)庫為一個蛋白質(zhì)同源結構比對結果參考庫，其結果是在MNYFIT、STAMP和COMPARER 3個程序計算結果的基礎上，經(jīng)過人工調(diào)整所得。當前版本有1 032個家族，每個家族有2 到41個成員。每個家族中平均蛋白質(zhì)的長度為17 到855個殘基，序列一致性為8%到94%。該數(shù)據(jù)庫的對比結果中保存了結構對齊的序列、對齊的二級結構、變化后的原子坐標等信息。因此，采用HOMSTRAD的對齊結構作為參考，來驗證自動實現(xiàn)蛋白質(zhì)結構比對方法的可行性及準確性。

利用本文方法，對HOMSTRAD中的700個蛋白質(zhì)雙結構實例進行了計算，平均匹配精度達到了89%。從這700個雙結構中任選出10對比對結果，與Dali、TM-align和SPalign這3種方法的結構比對精度結果進行了對比。

表3中，Match表示幾種方法結構比對方法對齊殘基對與HOMSTRA對齊殘基對的相同部分占HOMSTRA對齊殘基對的百分比。對于相似性較高的蛋白質(zhì)，幾種方法得到的結果幾乎一致;相比之下，相似性較低的蛋白質(zhì)，幾種方法得到的結果有所差別。1C20: A與1IG6: A比對中，本文方法配準精度達到73%，而TM-align、SPalign分別為53%、49%，明顯優(yōu)于這兩種方法。從平均值來看，本文方法得到的比對結果較好。

表3　4種方法對HOMSTRAD中10對實例的結果比對精度比較Tab.3　Accuracy comparison results of 4 methods for 10 examples in HOMSTRAD

3　結論

本文首先利用蛋白質(zhì)結構數(shù)據(jù)構建表征氨基酸取向的Cα坐標系，并與蛋白質(zhì)結構曲線中Cα原子處的曲率、撓率、相鄰Cα間的距離這3個空間旋轉(zhuǎn)平移不變量，共同構成黎曼度量。這樣一來就將蛋白質(zhì)結構的幾何特征嵌入到了黎曼度量表示的矩陣流形中，因此黎曼度量完整地保留了蛋白質(zhì)三維結構的所有信息。其次，使用黎曼度量作為蛋白質(zhì)三維結構相似性比較的特征描述子，用黎曼距離作為相似性度量，計算各黎曼度量之間的距離，作為蛋白質(zhì)結構相似度的衡量指標。最終，本文方法在3組實驗數(shù)據(jù)上從不同的角度進行驗證，結果均表明，本文方法得到的蛋白質(zhì)三維結構相似性比較結果較其他方法效果更為顯著，其中，HOMSTRAD中的700個蛋白質(zhì)雙結構實例計算結果的平均匹配精度達到了89%。

參考文獻

［1］梁毅.結構生物學［M］.北京:科學出版社，2005.

［2］彭群生，胡敏.蛋白質(zhì)三維結構相似性比較方法綜述［J］.計算機輔助設計與圖形學學報，2006，18(10) : 1465-1471.

［3］Holm L，Sander C.Protein structure comparison by alignment of distance matrices［J］.Journal of Molecular Biology，1993，233(1):123-138.

［4］Shindyalov I N，Bourne P E.Protein structure alignment by incremental combinatorial extension(CE) of the optimal path［J］.Protein Engineering，1998，11(9) : 739-747.

［5］Gibrat J F，Madej T，Bryant S H.Surprising similarities in structure comparison［J］.Current Opinion in Structural Biology，1996，6(3):377-385.

［6］Levitt M.STRUCTAL.A structural alignment program［J］.1994.

［7］Krissinel E，Henrick K.Secondary-structure matching(SSM)，a new tool for fast protein structure alignment in three dimensions［J］.Acta Crystallographica Section D: Biological Crystallography，2004，60(12) : 2256-2268.

［8］Zhang Y，Skolnick J.TM-align: a protein structure alignment algorithm based on the TM-score［J］.Nucleic Acids Research，2005，33(7) : 2302-2309.

［9］Louie A H，Somorjai R L.Differential geometry of proteins: a structural and dynamical representation of patterns［J］.Journal of Theoretical Biology，1982，98(2) : 189-209.

［10］Joshi S H，Klassen E，Srivastava A，et al.A novel representation for Riemannian analysis of elastic curves in Rn［C］//IEEE Conference on Computer Vision and Pattern Recognition，Minneapolis，MN，2007: 1-7.

［11］Klassen E，Srivastava A，Mio W，et al.Analysis of planar shapes using geodesic paths on shape spaces［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，2004，26(3) : 372-383.

［12］Liu W，Srivastava A，Zhang J.A mathematical framework for protein structure comparison［J］.PLoS Computational Biology，2011，7(2) : e1001075.

［13］Berman H M，Westbrook J，F(xiàn)eng Z，et al.The protein data bank ［J］.Nucleic Acids Research，2000，28(1) : 235-242.

［14］Hanson A J，Thakur S.Quaternion maps of global protein structure ［J］.Journal of Molecular Graphics and Modelling，2012，38: 256-278.

［15］紀永強，許志才.微分流形與黎曼幾何［M］.西安:陜西師范大學出版社，1994.

［16］Collard A，Bonnabel S，Phillips C，et al.An anisotropy preserving metric for DTI processing［J］.arXiv: 1210.2826，2012.

［17］Fischer D，Elofsson A，Rice D W，et al.Assessing the performance of inverted protein folding methods by means of an extensive benchmark［C］//Proceeding of the First Pacific Symposium on Biocomputing，1996: 300-318.

［18］Murzin A G，Brenner S E，Hubbard T，et al.SCOP: a structural classification of proteins database for the investigation of sequences and structures［J］.Journal of Molecular Biology，1995，247(4) : 536-540.

［19］Maiorov V N，Crippen G M.Significance of root-mean-square deviation in comparing three-dimensional structures of globular proteins ［J］.Journal of Mmolecular Biology，1994，235(2):625-634.

［20］Zhang Y，Skolnick J.Scoring function for automated assessment of protein structure template quality［J］.Proteins: Structure，F(xiàn)unction，and Bioinformatics，2004，57(4) : 702-710.

［21］Shindyalov I N，Bourne P E.Protein structure alignment by incremental combinatorial extension(CE) of the optimal path［J］.Protein Engineering，1998，11(9) : 739-747.

［22］Yang Y，Zhan J，Zhao H，et al.A new size-independent score for pairwise protein structure alignment and its application to structure classification and nucleic-acid binding prediction［J］.Proteins: Structure，F(xiàn)unction，and Bioinformatics，2012，80 (8) : 2080-2088.

［23］Elofsson A，F(xiàn)ischer D，Rice D W，et al.A study of combined structure/sequence profiles［J］.Folding and Design，1996，1(6): 451-461.

Similarity comparison of 3D protein structure based on Riemannian manifold

XU Yong-hong，CHU Ze-fei，HONG Wen-xue
(College of Electrical Engineering，Yanshan University，Qinhuangdao，Hebei 066004，China)

Abstract:As the representative technology of protein spatial structure exploration，NMR technology provides an unprecedented opportunity for modern life science research.But subsequent large data analysis has become a major problem.It is an important means to study protein structure and functional relationship by known information of proteins' three－dimensional structures to predict the unknown spatial structure of proteins.A method for similarity comparison of 3D protein structures based on Riemannian manifold theory is proposed in this paper.By constructing Cα frames and extracting geometric feature of protein，3D coordinates of proteins are converted into one dimension sequences with rotation and translation invariance.The Riemann distance is used as the three－dimensional structure similarity degree index.Spatial transformation on protein structure is not needed in this method，which avoiding errors when matching two proteins in the traditional method for registration by the least squares fitting.This method is independent of sequence information completely.It has realistic significance for proteins which do not have a similarity between sequences.Three experiments are designed according to 3 sets of data: proteins of different similarity，ten pairs whose protein structures are more difficult to identify proposed by Fischer，700 proteins in the HOMSTRAD database.Compared with the traditional method，the experiment results show that the matching accuracy of this method has been greatly enhanced.

Key words:protein; cubic spline interpolation; caframe; Riemannian manifold; structural comparison

作者簡介:*徐永紅(1975-)，男，四川犍為人，博士，教授，主要研究方向為醫(yī)學信息處理，Email: xyh@ysu.edu.cn。

基金項目:國家自然科學基金資助項目(60873121)

收稿日期:2014-06-26

文章編號:1007-791X(2015) 01-0035-07

DOI:10.3969/j.issn.1007-791X.2015.01.006

文獻標識碼:A

中圖分類號:R318