基于黎曼流形的蛋白质三维结构数据相似性比较

徐永红，褚泽斐，洪文学(燕山大学电气工程学院，河北秦皇岛066004)

基于黎曼流形的蛋白质三维结构数据相似性比较

徐永红*，褚泽斐，洪文学
(燕山大学电气工程学院，河北秦皇岛066004)

摘要:以NMR技术为代表的海量蛋白质空间结构数据为现代生命科学研究提供了前所未有的机遇，但后续的大数据分析却成为一大难题。充分利用已知的蛋白质三维结构信息来预测未知的蛋白质空间结构信息是研究蛋白质结构和功能关系一种重要手段。本文提出一种基于黎曼流形的蛋白质三维结构相似性比较新方法。该方法通过构建Cα坐标系和提取蛋白质结构具有旋转和平移不变性的几何特征量，将蛋白质的三维坐标序列转换为一维序列，采用黎曼距离作为三维结构相似度指标。本方法不需要对蛋白质结构做旋转和平移变换，避免了主流的RMSD方法中两蛋白质通过最小二乘拟合进行配准时产生的误差，并且完全不依赖于一级结构序列信息，对不具备序列相似性的蛋白质之间的相似性比较具有现实意义。本文分别针对不同相似度的蛋白质、Fischer提出的10个较难识别的蛋白质结构对、HOMSTRAD数据库中的700个数据这3组数据，对本文算法进行了验证。实验结果表明，与其他方法相比，本文方法的匹配精度均得到了较大提升。

关键词:蛋白质;三次样条插值; Cα坐标系;黎曼流形;结构比较

0　引言

生物大分子中的蛋白质和很多非编码RNA的功能主要取决于它们的空间结构。到目前为止，已经有超过六万个生物大分子的空间结构被测定，如何有效地比较它们之间的相似性成了生命科学中的一个重要课题［1］。蛋白质三维结构的直接获取一直是一个瓶颈问题，尽管蛋白质序列的测定已基本完成，但大量序列已知的蛋白质的三维结构尚未被实验方法测定出来，在这种情况下，充分利用一级序列信息和已知蛋白质的空间结构信息来研究预测未知蛋白质的空间结构，为结构生物学中研究蛋白质结构和功能关系的主要手段。由于用实验方法得到目的蛋白的结构困难较大，而且需要较长时间，蛋白质三维结构相似性比较可构建实验需要的结构模型，用于提出关于蛋白质功能的假设并指导进一步的实验工作［2］，因此蛋白质结构相似性比较是蛋白质功能分析、分类管理、检测等重要方法之一。

目前针对蛋白质结构比较已经有很多的研究方法与软件工具，如Dali［3］、CE［4］、VAST［5］、STRUCTAL［6］、SSM［7］、TM-align［8］等，为蛋白质空间结构的研究提供了多种分析手段。蛋白质结构比较方法主要分为3类:基于氨基酸间距离矩阵的比较(Dali、CE) ;基于蛋白质空间几何结构的比较(STRUCTAL、TM-align) ;基于蛋白质二级结构的匹配(VAST、SSM)。

传统的蛋白质相似性比较方法通常依赖于蛋白质的一级结构序列，但是两个一级结构序列不同而空间结构相同的蛋白质往往有着相同的性质。针对传统方法的不足，本文从几何角度出发，提出了一种基于黎曼流形的蛋白质结构相似性比较的新方法。通过微分流形的处理，用曲率k和挠率r这两个几何量代替蛋白质一级结构中和φ两个二面角，从而把蛋白质结构分析转变为一个纯数学问题。

1982年，Louie和Somorjai［9］最初将微分流形引入到了蛋白质结构研究中，他们将蛋白质的主链看成是连续的空间曲线，并用参数(b，ω)描述(其中α螺旋对应螺旋面，β折叠对应悬链面，b表示曲面的尺寸，ω表示曲面的扭曲程度)。通过这种描述方法，Louie和Somorjai一方面给出了识别蛋白质结构的指标;另一方面给出了蛋白质主链结构的连续描述。在此基础上，文献［10-11］对蛋白质结构三维曲线进行弹性形状分析提取弹性度量，将弹性度量用于蛋白质结构的比较中。文献［12］进一步提出弹性黎曼度量，将蛋白质结构从流形的角度进行比对，与其他结构比较方法相比效果更好，但是计算复杂度高。

本文通过构建Cα坐标系，对蛋白质主链进行曲线插值提取几何特征，用黎曼度量序列表示蛋白质结构数据，从而将蛋白质空间几何结构的相关信息嵌入到黎曼度量序列表示的矩阵流形中。通过这种流形表示和嵌入，可以引入黎曼几何的数学理论来提取有效的几何结构特征。该方法能够更直观的体现出蛋白质的结构特征，与传统的方法相比几何意义明确，为蛋白质结构相似性比较提供了一种新的思路。

1　本文方法

本文所采用数据均来自PDB［13］(Protein Data Bank)数据库，该数据库是美国Brookhaven实验室于1971年建立的大分子蛋白质晶体结构资料数据库，是互连网上唯一有生物高分子3维结构的数据库，它的内容主要是根据X射线结晶和核磁共振得到的实验数据。

本文的关键是对蛋白质主链进行曲线插值，从蛋白质的三维几何结构中提取具有空间平移和旋转不变性的特征量，并将其嵌入到矩阵流形中。本文方法的具体步骤如下:

1)数据预处理

首先从数据库中根据唯一标识的PDB ID获取蛋白质的PDB文件。忽略掉H原子，氨基酸的骨架是一个顶端为Cα原子的四面体，另外3个顶点分别为氨基、羧基和R基。图1为蛋白质1CRN中两个氨基酸的基本原子框架，其中R基为CBSG，可以看到每个氨基酸的骨架原子为N-Cα-C的组合。从PDB文件中提取每个氨基酸的骨架原子——Cα、C、N原子的空间坐标。

2)构建Cα坐标系

弗莱纳(Frenet-Serret)公式常用来描述粒子在连续可微的曲线上运动时曲线的切向、法向、副法方向之间的关系，本文在弗莱纳公式的基础上根据蛋白质结构数据中Cα、C、N原子的空间坐标构建蛋白质的Cα坐标系［14］:

其中，Cα、C、N分别代表Cα、C、N原子的空间坐标，X、Y、Z分别为列向量，因此得到3×3的单位正交矩阵构成的Cα坐标系序列［F1F2…Fi…Fn］，其中Fi= (XiYiZi)。每个Cα坐标系Fi均表征了对应氨基酸骨架的空间取向，每个氨基酸对应一个单位正交矩阵，则将蛋白质氨基酸序列转化为了单位正交矩阵序列。

图1　氨基酸的基本结构Fig.1　Amino acid neighboring structure

3)曲线插值

三次样条插值方法解决了计算几何中自由型曲线曲面形状的描述问题，在整体表示的同时保持了局部性质，因此对蛋白质Cα主链结构进行了三次样条曲线插值。

三次样条函数定义如下:

对于［a，b］上的分划: a = x0＜x1＜…＜xn= b，则

利用三次样条函数进行曲线插值，即取插值函数为三次样条函数，称为三次样条插值。

如图2所示，(b)和(c)分别为蛋白质1CRN三次样条插值前后的结构曲线图，数字为氨基酸的序号，‘○’处为Cα原子所在位置。可以观察到，经插值后蛋白质主链结构转化为光滑曲线，两个Cα原子间的每一段曲线都对应一个三次多项式。

图2　蛋白质1CRN主链结构曲线插值Fig.2　Curve interpolation of protein 1CRN backbone structure

4)几何特征量的提取

曲率和挠率分别反映了空间曲线的弯曲程度和扭曲程度，它们具有空间旋转平移不变性，能够充分地描述曲线的几何特征，因此本文选取蛋白质结构曲线中每个Cα原子处的曲率、挠率和相邻Cα原子间的距离作为特征描述子，对蛋白质的几何结构进行描述。

若已知空间曲线的参数方程

则在参数t处曲线的曲率kt和挠率τt计算公式分别为

其中，x′，x″，x′″分别为x对t的一、二、三阶导，同理可知其他。

根据式(4)计算蛋白质结构曲线上每个Cα点处的曲率k和挠率τ，并与相邻Cα原子间的距离d的绝对值共同构成对角阵

其中，ki、τi为第i个Cα原子处的曲率和挠率，di，(xi， yi，zi)和(xi +1，yi +1，zi +1)分别为序号为i和i +1的Cα原子的空间坐标，这样便得到了包含蛋白质结构曲线特征描述子的对角阵序列。

5)构造黎曼度量序列

黎曼度量［15］指的是空间上的几何学应基于无限临近两点(x1，x2，…，xn)与(x1+ d1，x2+ d2，…，xn+ dn)之间的距离，用微分弧长度平方所确定的正定二次型的理解度量，亦即是由函数构成的正定对称矩阵，这便是黎曼度量。在三维空间中可以用一个3×3的实对称矩阵来表示空间某点的黎曼度量:

式中，det(S)＞0，λi＞0，正定对称矩阵的分解是可逆的。

根据式(2)中正定对称矩阵的性质，利用步骤(2)、(3)中得到的单位正交矩阵序列和对角阵序列构造黎曼度量序列:

Si为序号为i的氨基酸所对应的黎曼度量。

6)计算两蛋白质对应黎曼度量之间的黎曼距离

通过将黎曼度量S分解为单位正交矩阵F和对角阵Λ，再经过计算Fi、Fi +1和Λi、Λi +1间的距离d(Fi，Fi +1)和d(Λi，Λi +1)，最终得到Si、Si +2之间的黎曼距离d(Si，Si +1)

其中

λ1，i、λ2，i分别为对角阵Λi、Λi +1中第i行i列的值。

系数k(Λi，Λi +1)为加权因子［16］，范围在(0～1)之间:

λmax、λmin分别为对角阵Λ中的最大值和最小值。

7)设定阈值，比较分析

由步骤(5)得到了两蛋白质间的黎曼距离序列［ds1，ds2，…，dsi，…，dsn］，黎曼距离dsi的值越小，表示两个包含蛋白质结构几何信息的黎曼度量在黎曼流形上离得越近，故对应的蛋白质结构越相似。经多次实验验证后取黎曼距离的阈值:

k = 0.6mean(ds) + 0.3std(ds)，(10) mean(ds)、std(ds)分别表示黎曼距离序列的均值和方差。小于阈值k则认为在两蛋白质在残基i处具有相同的结构。

2　实验结果与分析

为验证本文方法的可行性及有效性，进行了3组实验:第1组在SCOP数据库中分别选取家族、超家族、非超家族的蛋白质数据进行实验，验证本文方法的比对结果是否符合生物意义;第2组对Fischer［17］给出的10个较难识别结构对进行实验，从“相同残基比”方面验证本文方法的可行性;第3组基于HOMSTRAD结构比对参考库，验证本文方法在同源识别方面是否具有有效性。

2.1不同相似度的蛋白质结构比对

蛋白质结构分类数据库SCOP［18］(Structural Classification of Proteins)是提供关于已知结构蛋白质之间的结构和进化关系信息，所涉及的蛋白质包括结构数据库PDB中的所有条目。SCOP从总体上将蛋白质分为全α型，全β型，以平行折叠为主的α/β型，以反平行折叠为主的α+β型等。从SCOP数据库中选取PDB ID分别为101M、102M、1UVY、1C7Y的蛋白质数据，他们分别属于相同家族、相同超家族、不同超家族，相似度依次降低。通常，两蛋白质的关系越近，结构越相似度越高，黎曼度量序列间的黎曼距离越小。

表1为蛋白质101M与101M、102M、1UVY、1C7Y的结构比对结果，Riemann_dis为黎曼度量序列间的黎曼距离的平均值，并与传统蛋白质结构比对评价指标RMSD［19］、TM-score［20］和Z-score［21］进行对比。RMSD(Root mean square deviation)方法为蛋白质骨架直接配准的比较，该方法提出最早、使用最广泛，RMSD值越小，表示两蛋白质结构越相似。与之趋势相反的是TM-score与Z-score，值越大表示结构相似度越高。从表1中可以看出，随着蛋白质对之间关系的变远，4个参数结果均表示蛋白质结构相似度降低。图3为4对蛋白质间的黎曼距离曲线，从图中观察到随着蛋白质间的关系变远，蛋白质结构相似度降低，黎曼距离值明显减小。

2.210个较难识别蛋白质结构对的比对

在文献［23］中Fischer给出了结构难以比较的10个蛋白质对，利用本文方法对这10个蛋白质对进行结构比对，并与Dali［3］、TM-align［8］、SPalign［22］这3种传统方法进行了结果比较结果如表2所示。其中，“Equ”表示进行比较的残基个数，“ResNum”表示相同的残基个数，“Riemann_dis”表示本文方法中黎曼距离的平均值。通常认为，相同残基比(Res-Num/Equ)越大，比对的效果越好。

表1　蛋白质101M与不同相似度蛋白质的结构比对结果Tab.1　Structure comparison results of different similarity of protein and protein 101M

图3　蛋白质101M与不同相似度的蛋白质之间的黎曼距离曲线Fig.3　Curve of riemann distance between protein 101M and different similarity of proteins

如表2所示，本文方法的比例值高出效果最好的Dali方法约0.5个百分点，比效果较差的SPalign方法高出了约2.8个百分点，充分验证了本文方法不仅可行，而且效果明显优于另外3种方法。

表2　10个较难识别蛋白质结构对的比对平均结果比较Tab.2　Structure comparison average results of 10 more difficult to identify protein structures

2.3基于结构比对参考库的验证

HOMSTRAD(Homologous Structure Alignment Database)数据库为一个蛋白质同源结构比对结果参考库，其结果是在MNYFIT、STAMP和COMPARER 3个程序计算结果的基础上，经过人工调整所得。当前版本有1 032个家族，每个家族有2 到41个成员。每个家族中平均蛋白质的长度为17 到855个残基，序列一致性为8%到94%。该数据库的对比结果中保存了结构对齐的序列、对齐的二级结构、变化后的原子坐标等信息。因此，采用HOMSTRAD的对齐结构作为参考，来验证自动实现蛋白质结构比对方法的可行性及准确性。

利用本文方法，对HOMSTRAD中的700个蛋白质双结构实例进行了计算，平均匹配精度达到了89%。从这700个双结构中任选出10对比对结果，与Dali、TM-align和SPalign这3种方法的结构比对精度结果进行了对比。

表3中，Match表示几种方法结构比对方法对齐残基对与HOMSTRA对齐残基对的相同部分占HOMSTRA对齐残基对的百分比。对于相似性较高的蛋白质，几种方法得到的结果几乎一致;相比之下，相似性较低的蛋白质，几种方法得到的结果有所差别。1C20: A与1IG6: A比对中，本文方法配准精度达到73%，而TM-align、SPalign分别为53%、49%，明显优于这两种方法。从平均值来看，本文方法得到的比对结果较好。

表3　4种方法对HOMSTRAD中10对实例的结果比对精度比较Tab.3　Accuracy comparison results of 4 methods for 10 examples in HOMSTRAD

3　结论

本文首先利用蛋白质结构数据构建表征氨基酸取向的Cα坐标系，并与蛋白质结构曲线中Cα原子处的曲率、挠率、相邻Cα间的距离这3个空间旋转平移不变量，共同构成黎曼度量。这样一来就将蛋白质结构的几何特征嵌入到了黎曼度量表示的矩阵流形中，因此黎曼度量完整地保留了蛋白质三维结构的所有信息。其次，使用黎曼度量作为蛋白质三维结构相似性比较的特征描述子，用黎曼距离作为相似性度量，计算各黎曼度量之间的距离，作为蛋白质结构相似度的衡量指标。最终，本文方法在3组实验数据上从不同的角度进行验证，结果均表明，本文方法得到的蛋白质三维结构相似性比较结果较其他方法效果更为显著，其中，HOMSTRAD中的700个蛋白质双结构实例计算结果的平均匹配精度达到了89%。

参考文献

［1］梁毅.结构生物学［M］.北京:科学出版社，2005.

［2］彭群生，胡敏.蛋白质三维结构相似性比较方法综述［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，Feng 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，Function，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，Function，and Bioinformatics，2012，80 (8) : 2080-2088.

［23］Elofsson A，Fischer 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-)，男，四川犍为人，博士，教授，主要研究方向为医学信息处理，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