蛋白质结构域边界网络流预测方法

余众泽,彭春祥,张贵军

(浙江工业大学 信息工程学院,杭州 310023)

1 引 言

蛋白质结构域是蛋白质结构、折叠、功能、进化和设计的基本单元.80%以上的真核蛋白和67%以上的原核蛋白含有多个结构域,同时许多生物学功能也依赖于不同结构域之间的相互作用[1].因此,推断蛋白质结构域边界是推断蛋白质折叠机理、理解生物功能和注释进化机制的重要步骤[2].此外,随着AlphaFold2在蛋白质结构预测中的革命性成功[3],在生物信息领域内一致认为单域蛋白质的结构预测几乎得到解决[4].在蛋白质结构预测关键评估比赛CASP14中,AlphaFold2所预测的蛋白质结构主链准确性中位数为0.96Å[3].然而,在缺乏共进化信息的情况下,具有多域蛋白的全链建模仍然是一个挑战.当前的多域蛋白结构建模主要有两种方法[5]:第1种是直接从蛋白质序列来预测全链模型,例如AlphaFold[6]、RaptorX[7]和I-TASSER[8];第2种是将多域蛋白先划分为单个域,然后通过域组装方法将单独建模的单域结构组装成全链模型,例如DEMO[1]和AIDA[9].随着蛋白质序列长度的增加,直接预测全链模型变得非常困难且效率低下,然而组装方法几乎不受蛋白质全链长度的影响.因此,如何准确预测结构域边界是域组装方法的基础和关键.通常来说,结构域边界预测方法可以分为基于结构和基于序列的两类方法.

基于结构的方法需要实验测定或预测的蛋白质结构来识别结构域.CATHEDRAL[10]将目标蛋白质结构与来自CATH[11]数据库的模板结构进行比较来检测结构域.PDP[12]和DDOMAIN[13]根据域内接触多于域间接触的假设将蛋白质划分为多个结构域.DHcL[14]通过计算蛋白质的范德华模型来分解蛋白质结构域.SWORD[15]通过蛋白质单元的分层合并来分配结构域,其中的蛋白质单元是进化保守的子结构,它是结构域和二级结构之间的一个中间级别,用来描述蛋白质结构.此外,还有一些是基于蛋白质结构预测的方法.例如,RosettaDom[16]使用Rosetta从头结构预测方法来构建蛋白质三维模型,然后应用Taylor方法将预测的结构解析为多个结构域.这类从蛋白质结构推断结构域的方法通常具有较高的准确性,但是它高度依赖于蛋白质模型的精度.然而,由于结构预测能力有限和生物实验测定的局限性,它们只能应用于小蛋白或具有已知实验结构的蛋白质[2].与结构信息相比,蛋白质序列信息更容易获得[5].因此,从序列中预测结构域边界的方法更有意义和更具挑战性[2].

基于序列的结构域识别方法通常可以分为两类.第一类主要是基于同源性的方法,它通过将目标序列与具有已知结构域的同源序列进行序列比对来检测结构域.例如,Pfam[17]、CHOP[18]和FIEFDOM[19],它们通过隐马尔科夫(HMM)或PSI-BLAST程序在已知的蛋白质结构或家族数据库中搜索目标序列的同源信息,然后从同源模板或家族中获取结构域边界信息.ThreaDom采用基于穿线的算法来提高远同源模板检测[20].它首先使用LOMETS[21]程序在PDB中搜索目标序列的同源模板,然后基于目标序列构建多序列比对(MSA).根据这些MSA,计算结构域保守分数来衡量每个残基的保守水平,并用于进一步判断边界区域.在此基础上,ThreaDomEx[22]通过结构域片段组装来分配不连续结构域.这些基于同源性的方法在识别到质量高的模板时,可以达到很高的预测准确度,但是目标和模板之间的序列一致性较低时,准确度会急剧下降[20].另一类是从头预测方法,它可以在一定程度上克服这种限制[5],代表性的方法包括DOMpro[23]、DoBo[24]、ConDo[25]和DNN-dom[26].DOMpro[23]训练递归神经网络用于域识别,其输入特征包括序列谱、预测的二级结构和溶剂可及性.DoBo[24]将同源蛋白质中包含的进化结构域信息引入到蛋白质结构域边界预测中.ConDo[25]除了传统的局部窗口特征外,还利用远程的共进化特征训练神经网络来检测结构域.DNN-dom[26]结合卷积神经网络和双向门循环单元模型来预测蛋白质的结构域边界.尽管基于序列的从头预测方法更加实用,尤其是在没有同源序列的情况下,但它的准确性通常不高,需要进一步提升[2].

本文提出了一种基于序列的蛋白质结构域边界从头预测算法GraphDom.该算法根据预测的残基间距离分布将结构域划分问题转化为网络流分割问题,其中残基表示为网络中的一个节点,残基间的距离表示为具有一个特定容量的边,通过找到网络的一个最小切割来解决双结构域划分问题,再使用递归划分的方法来解决复杂结构域划分问题.在120个非冗余测试蛋白的实验结果表明,GraphDom在预测精度上优于ThreaDomEx、ConDo和DoBo,是一种有效的蛋白质结构域边界预测方法.

2 GraphDom算法设计

首先介绍蛋白质的网络表示,并基于Ford-Fulkerson网络流算法,提出了一种结构域划分算法GraphDom,最后详细描述GraphDom算法中的参数.

2.1 蛋白质网络的定义

一个网络是由一组节点和一组边构成的有向图.其中有两个特殊的节点,分别是源节点s和汇节点t.每条边连接两个节点,并且边的容量大于零.对于一个蛋白质,本文使用一个节点来表示蛋白质的一个残基,并使用一条边来表示所连接的两个残基在空间上是接触的,即两个残基间的距离小于或等于8Å[27].同时本文通过定义一条边的容量来衡量两个残基在空间上的紧密程度.蛋白质网络如图1所示,其中源节点s和汇节点t是两个人为定义的节点.

一个s-t切割是一组边,将这组边去除,则没有从源s到汇t的路径.例如在图1中,边{(s,1),(s,2),(s,3),(s,4),(s,5)}形成一个s-t切割.而一个最小s-t切割是指具有最小总边容量的s-t切割,例如边{(1,6),(2,6),(2,11),(7,8),(10,14),(10,15)}形成图1中的最小s-t切割.

根据最大流/最小割定理[28],可以通过找到从源s到汇t的最大流来计算最小s-t切割.在此使用f(u,v)来表示分配给边(u,v)的流值,用c(u,v)来表示分配给边(u,v)的容量值.如果分配给有向网络边的一组值满足以下3个条件,则形成一个s-t流.

1)容量限制:对于每条边(u,v),f(u,v)≤c(u,v)

2)斜对称:对于每条边(u,v),f(u,v)=-f(v,u)

3)流守恒:对于除源s和汇t以外的每个节点u,节点u的总输入流应该等于总输出流.即Σvf(u,v)=0,其中∑v表示对所有节点求和.

为了将预测的残基间距离信息转化为蛋白质网络中的边容量,本文设计了边容量公式,如公式(1)所示:

(1)

2.2 Ford-Fulkerson算法

对于一个给定的流f,边(u,v)的剩余容量定义为cf(u,v)=c(u,v)-f(u,v).根据上述流的定义,可知cf(u,v)是≥0是≥0.蛋白质网络表示图的剩余容量图如图2所示,对于一个流定义如公式(2)所示:

图2 蛋白质网络剩余容量图Fig.2 Protein network residual capacity diagram

(2)

其中x→y表示从节点x到节点y的有向边,需要注意的是剩余容量可能大于容量,因为流量可能具有负值.在图2中,具有零剩余容量的边没有绘制.

Ford-Fulkerson算法的基本思想是重复寻找从源s到汇t的有向路径,它由具有Cf(u,v)>0的有向边(u,v)组成,然后通过p的Cf(u,v)最小值来增加路径p上的每条边的流值,这个过程将持续进行,直到找不到这样的路径为止,而一开始所有有向边的流值f(u,v)都设置为0.按照这个步骤,可以检查出图1中网络的最大s-t流值为24,如图2所示当找到最大s-t流时,标记有剩余容量的蛋白质网络.显然,在图2中没有从源s到汇t的有向路径.

然而最小s-t切割并不是唯一的,即存在多于一个切割的最小s-t切割,其边的总容量为24.例如在图2中,边{(1,6),(2,6),(2,11),(7,8),(5,10),(9,10)}和边{(1,6),(2,6),(2,11),(2,8),(3,9),(4,9),(4,5),(s,5)}也形成了最小s-t切割.因此当蛋白质网络的最小s-t切割不唯一时,GraphDom将枚举所有最小切割.

2.3 枚举所有最小切割

枚举所有最小切割,等价于枚举所有可能形成的分区.然而,对于大的蛋白质而言,直接在剩余容量图上枚举出所有可能的分区时间复杂度太高.因此,在这里本文利用强连接分量来降低时间复杂度.

首先将蛋白质剩余容量图中的强连接分量提取出来形成强连接分量图,然后再在强连接分量上枚举所有可能的最小切割.一个有向图的强连接分量是指,强连接分量中的每个节点可以通过有向边到达强连接分量的每一个其他节点,即强连接分量中的所有节点都是联通的,一个点可以到达任意其他节点.在图2中,节点{5,8,9}形成一个强连接分量,但是节点{5,9,10}并没有形成强连接分量,因为节点9不能到达节点10,节点5也不能到达节点10.关于最小s-t切割形成分区的一个简单观察是,如果s分区包含某个节点u,则s分区必须包含含有u的强连接分量,对于t分区也是如此.因此可以将强连接分量视为一个节点,合并后的蛋白质强连接分量图如图3所示.

图3 蛋白质网络强连接分量图Fig.3 Strongly connected component diagram of protein network

图3中的V1={s,1,2,3,4},V2={5,8,9},V3={10},V4={6,7,11,12,13,14,15,t},如果分别属于两个强连接分量的两个节点之间存在有向边,则在两个强连接分量之间放置有向边.

在本文中,将Spart定义为包含源s的强连接分量分区,将Tpart定义为包含汇t的强连接分量分区.通过回溯来考虑分配给Spart和Tpart的强连接分量组合,并在最小切割定义的约束下,即没有从源s到达汇t的路径,排除不合理的组合情况.以下列出图3的所有合理的s-t分区情况:

1)Spart={V1}和Tpart={V2,V3,V4}

2)Spart={V1,V2}和Tpart={V3,V4}

3)Spart={V1,V2,V3}和Tpart={V4}

需要注意的是Spart={V1,V3}和Tpart={V2,V4}形成一个不合理的s-t分区,因为Spart分区的节点V3可以到达Tpart分区的节点V2,这与最小切割的定义相驳,即没有从源s到汇t的路径.在结构域边界评估过程中,GraphDom将根据结构域的一般特性和结构域边界评估函数对不同分区情况进行评估和判断.

2.4 结构域边界评估函数

根据蛋白质结构域的一般特性设计结构域边界评估函数来评估划分的分区或区域,其遵循的原理是域内残基密度大域间残基密度小,该函数定义如公式(3)所示:

(3)

(4)

(5)

其中u、v、i分别表示不同的残基,SOu是Spart分区的残基u相对于Tpart分区的外部紧密度分数,SIu是Spart分区的残基u相对于Spart分区的内部紧密度分数,dui是残基u和残基i的距离.在本文中,结构域边界评估函数值越低表示划分的分区或区域越好.

在GraphDom中,如果划分的分区结果同时满足每个分区的残基数量大于等于40个且结构域边界分数值(DS)小于阈值DScutofff,则接收该次划分.

2.5 源汇极点的选择

所选择的源汇点代表着需要识别的结构域,不同的源汇点组合可能会导致不同的划分结果.因此,本文中选择许多不同的源汇点组合,为了获得好的结构域划分结果.在GraphDom中,采用以下两条规则来获得源汇点的集合:

1)根据半径12Å以内的邻域残基密度,选择目标序列密度前10%的残基作为候选的源汇点;

2)为了避免密度靠前的残基聚集在某一个域内,每隔20个残基,选择一个残基作为候选的源汇点.

2.6 算法全流程及参数确定

全流程的详细细节如图4所示,首先根据输入的目标序列预测其残基之间的距离分布,即预测距离图.在本文中,残基间距离分布图来自于trRosetta[29].然后,将预测距离图转换成初始蛋白质容量图,对源汇点集合中的所有极点组合使用Ford-Fulkerson算法获得蛋白质剩余容量图,再使用深度优先算法获得强连接分量和使用回溯算法枚举所有可能的最小切割,最后根据结构域的一般特性和域边界评估函数,对获得的所有s-t分区结果进行评估,并根据最优的分区结果判断是否需要继续往下递归划分.如果该次划分有效,则将输入的蛋白质容量图根据分区结果分成两个蛋白质容量图,并继续往下进行划分,划分过程的停止条件就是划分后的s-t分区不满足每个分区至少包含40个残基或结构域边界分数DS大于等于阈值DScutoff.

图4 GraphDom的管道图Fig.4 Pipeline diagram for GraphDom

在GraphDom程序中有2类参数需要确定:1)与边容量有关的参数ω1和ω2;2)与递归划分过程的停止条件有关的参数DScutoff,通过优化GraphDom在训练集上的预测性能来对这些参数进行调整.

为了训练和测试GraphDom,本文从SCOPe2.07-stable数据库[30]中收集一组具有已知结构域结构的非冗余蛋白质,在收集过程中以序列相似度<30%和序列长度>40个残基进行筛选过滤.该数据集包含240个多域蛋白质,这些蛋白质被随机分成120个训练蛋白和120个测试蛋白.在GraphDom中,ω1=5、ω2=10和DScutoff=0.45都是通过在训练蛋白上进行训练得到的.

3 实验结果与分析

为了测试GraphDom的域边界预测性能,本文在120个测试蛋白上比较GraphDom和其他3种方法.这120个测试蛋白来自SCOPe2.07-stable数据库[30],它们包括104个两域蛋白和16个三域蛋白,且序列相似度<30%.在这120个测试蛋白中还包含12个不连续域蛋白.该测试集由多种类别的蛋白质组成,可以公平地测试GraphDom和其他3种主流方法的性能.

NDO[31]和DBD[32]分数用于在CASP[33]比赛中评估结构域划分和结构域边界预测的质量,其中NDO分数是通过计算预测的结构域区域和真实的结构域区域之间的重叠区域的归一化获得,而DBD分数是通过计算预测的结构域边界和真实的结构域边界之间的距离获得,其中如果真实的结构域边界位于Loop区域时,则整个Loop区域都被视为结构域边界.

3.1 多域蛋白质的识别

单域蛋白和多域蛋白是蛋白质的两种类别,这里本文对GraphDom和其他3种先进方法的结构域识别能力进行比较,比较结果如表1所示.

表1 在120个测试蛋白上的多域识别结果Table 1 Multi-domain classification results on 120 test proteins

在120个蛋白的测试集上,GraphDom的多域蛋白召回率为0.908,分别比ThreaDomEx和ConDo高11.1%和31.2%,然而比DoBo低2.7%.测试集的平均结构域数量为2.117,GraphDom、ThreaDomEx、ConDo和DoBo的预测结果的平均结构域数量分别为2.083、2.292、1.867和3.367,与测试集平均结构域数量的差值分别为-0.034、0.175、-0.25和1.25.由此可见,DoBo对结构域有过度预测的倾向,平均每个蛋白质多预测出一个结构域,即DoBo倾向于高估多域蛋白的数量,所以在多域蛋白召回率上高于GraphDom.

3.2 结构域边界的预测

为了测试不同方法预测结构域边界的能力,本文在测试集上比较GraphDom和其他3种先进方法,并在表2中总结各个方法的平均NDO和DBD分数.

表2 在测试集上的多域蛋白质预测结果总结Table 2 Summary of prediction results of multi-domain proteins in the test set

如表2所示,在120个蛋白质的测试集上,GraphDom的平均NDO和DBD分数分别为0.856和0.603,高于其他3种方法的NDO和DBD分数,通过单尾成对双样本t-校验计算GraphDom预测结果和其他3种方法预测结果之间的P-值,从中可以看到这些P-值都是<0.05,这表明GraphDom在域边界预测能力上与其他3种比较方法具有显著差异.在测试集蛋白上,GraphDom的NDO分数分别比ThreaDomEx、ConDO和DoBo的NDO分数高18.9%、13.8%和36.3%,GraphDom的DBD分数分别比ThreaDomEx、ConDO和DoBo的DBD分数高45.3%、66.6%和95.8%.

除此之外,从表2中可知,ThreaDomEx和ConDo在NDO分数方面是接近相等的,它们NDO分数之间的P-值为0.09,然而在DBD分数方面,ThreaDomEx比ConDo高14.7%,说明二者的预测结构域区域与真实结构域区域重叠比例相近,但是ThreaDomEx预测结构域边界能力强于ConDo.

3.3 不连续域蛋白质的预测结果

由序列上不连续的片段组成的不连续结构域相对于连续结构域而言建模更加困难,这里本文在测试集中的12个不连续多域蛋白质上,对GraphDom和其他3种主流方法的预测性能进行了比较,具体比较结果如表3所示.

表3 在12个不连续域测试蛋白上的预测结果总结Table 3 Summary of prediction results on 12 discontinuous domain test proteins

如表3所示,在测试数据集的12个不连续多域蛋白质中,GraphDom识别出其中9个蛋白包含不连续域,召回率为75%.ThreaDomEx是基于ThreaDom开发的域边界预测方法,并针对不连续域进行优化,然而它只检测到33.4%的蛋白质包含不连续域.其他两种基于机器学习的方法ConDo和DoBo,没有检测出任何含有不连续域的蛋白质.

值得注意的是,对于连续域蛋白质和不连续域蛋白质,GraphDom的结构域边界预测性能非常接近.更具体地说,GraphDom在连续域和不连续域蛋白质的NDO和DBD分数分别为0.860/0.813和0.602/0.616.这些结果表明,GraphDom的性能并不明显依赖于多域蛋白质的结构域类型,从而突出了GraphDom在迭代划分过程中识别复杂结构域的有效性.

4 结 论

本文提出一种基于网络流的蛋白质结构域边界预测算法GraphDom.首先,根据边容量公式将预测的残基间距离分布图转换为蛋白质容量图,然后根据最大流/最小割定理通过Ford-Fulkerson算法获得目标蛋白的剩余容量图,对于具有非唯一最小切割的网络,先使用深度优先算法获得强连接分量图,再使用回溯算法枚举每一个可能的最小切割,最后基于蛋白质结构域的一般特性设计的结构域边界评估函数来对所有可能的划分进行评估,判断是否可以继续划分,如果可以划分则将输入的蛋白质容量图分成两个,每一个蛋白质容量图继续执行上述步骤,直到不满足每个分区至少包含40个残基或结构域边界分数小于阈值为止.在120个非冗余测试蛋白上的实验结果表明,GraphDom具有较强的结构域边界预测能力和复杂结构域识别能力,是一种有效的结构域边界预测算法.在下一步研究中,我们将利用已知的结构域序列信息来辅助识别目标序列的结构域数量和预测结构域边界可能存在的区域,进一步提高蛋白质结构域边界预测精度.