GraphBind: 通过层次图神经网络学习蛋白质局部结构嵌入规则以用于识别核酸结合残基

《GraphBind: protein structural context embedded rules learned by hierarchical graph neural networks for recognizing nucleic-acid-binding residues》
作者：夏莹，春秋夏，潘晓勇，沈宏斌
单位：上海交通大学
发表时间：2021年2月12日
发表期刊：Nucleic Acids Research
paper
code and dataset

ABSTRACT

了解蛋白质和核酸之间的相互作用是了解各种生物活性和设计新药的基础。如何准确识别核酸结合残基仍然是一个具有挑战性的任务。在本文中，我们提出了一种基于端到端图神经网络的准确预测器GraphBind，用于识别蛋白质上的核酸结合残基。考虑到结合位点通常在局部三级结构上表现为高度保守的模式，我们首先根据目标残基结构上下文背景和空间领域构造图，利用层次图神经网络(HGNN)嵌入潜在的局部结构和生化特征模式，用于结合残基的识别。我们在DNA/RNA基准数据集上全面评估GraphBind。结果表明，GraphBind的性能优于最先进的方法。此外，将GraphBind扩展到其他配体结合残基预测，以验证其泛化能力。GraphBind的网络服务器可在http://www.csbio.sjtu.edu.cn/bioinf/GraphBind/.免费获得

INTRODUCTION

蛋白质和核酸之间的相互作用参与了多种生物学活动和过程，如基因复制和表达、信号转导、调节和代谢(1-3)。研究蛋白质和核酸之间的相互作用对于分析遗传物质、了解蛋白质功能和设计新药具有重要意义。为了研究分子间的相互作用，人们设计了许多实验方法，如X射线、核磁共振、激光拉曼光谱等。然而，它们通常既耗时又昂贵。开发可靠和准确的计算方法来大规模筛选识别核酸结合残基是非常有必要的(4)。

根据使用的数据类型，现有的识别核酸结合残基的计算方法一般可以分为两类：基于序列的方法和基于结构的方法。基于序列的方法，如ConSurf(5)，TargetDNA(6)，DRNApred(4)，Scriber(3)和d Target S(7)，使用序列衍生物理化特性特征。例如，在TargetDNA中，从蛋白质序列中提取蛋白质的进化保守信息和预测的溶剂可及性，并使用支持向量机从由滑动窗口策略确定的DNA结合残基的序列上下文中识别DNA结合残基(6)。基于序列的方法的优点是，它们可以仅根据其序列对任何蛋白质进行预测。然而，它们的性能可能是有限的，因为结合残基的潜在模式仅从它们的序列中并不明显，而在空间结构中是保守的(8，9)。因此，从蛋白质序列中捕获的特征可能不足以准确地表示残基。

与基于序列的方法不同，基于结构的方法的假设是，具有特定功能的结构基序通常在局部三级结构(8，9)上表现为高度保守的模式。基于结构的方法可分为以下两类：(I)基于模板的方法，如DR bind1(10)和TM-site(11)，通过结构比较为查询蛋白质寻找可靠的模板，并根据物理和化学原理推断蛋白质与核酸之间的相互作用；(ii)基于特征的机器学习方法，如Aarna(12)和NucleicNet(13)，其利用来自蛋白质结构的特征构造分类器。

功能位点通常由序列(14)以外的三级结构的局部模式决定。我们致力于用基于特征的机器学习方法从蛋白质结构中识别与核酸结合的残基。一个主要的挑战是如何嵌入关键的结构和生物物理化学特征来识别下游结合残基。以前的方法通常使用手工制作的特征来表示结构(12)。这些方法需要很强的领域知识，手工制作的特征可能无法捕获特定下游任务的蛋白质的关键信息。一些其他方法将蛋白质结构编码到三维(3D)欧几里德空间(15，16)。例如，DeepSite将蛋白质原子映射成3D体素来表示蛋白质结构(16)。然后使用3D卷积神经网络(3DCNNs)(17)基于3D体积表示(16)从其邻域中提取目标残基的抽象特征。蛋白质结构的三维体积表示有两个潜在的缺点：(1)残基的稀疏和不规则分布使残基的邻域信息难以表示；(2)难以保证在三维直角坐标系中旋转和平移的不变性。或者，Delia计算距离矩阵来表示残基对的距离关系。Delia将结构看作2D图像，使用固定大小的卷积核(18)从所有残基(19)的局部距离关系中学习模式，导致某些残基的邻域信息不完整，忽略了结构相邻残基之间传递的知识。

为了更好地捕捉蛋白质结构信息和残基之间的空间关系，用图表示蛋白质结构，节点表示残基，边根据残基之间的空间关系定义。该图表示法不仅具有旋转和平移的不变性，而且可以处理数量不断变化的残基无序邻域。最近，图神经网络(GNNs)已经成为计算生物学中处理图数据的有力工具(20)。例如，Fout等人。提出了一种基于GNN的方法，用于从蛋白质结构中对成对残基相互作用进行分类(21)。十角形使用图卷积网络(GCNS)预测不同药物组合的副作用(22)。DimiG使用半监督GCNS在交互图上推断与microRNA相关的疾病(23)。Torng和Altman提出了一个两步图-卷积(Graph-CNN)框架，用于预测药物与靶标的相互作用(24)。以上研究表明，神经网络在处理生物和化学图数据方面是有效的。

在这项研究中，我们提出了一个准确的核酸-酸结合残基预测器，GraphBind，基于结构上下文和层次图神经网络(HGNNs)构建的图。为了从蛋白质结构中提取关键的结构和生物物理化学特征的局部模式，对于每个目标残基，我们首先基于目标残基的局部环境构建一个图。初始节点特征向量包括进化保守性、二级结构信息、其他生物物理化学特征和位置嵌入。位置嵌入是根据定义结构上下文中残基的空间关系的几何知识来计算的。初始边缘特征向量也是由几何知识导出的。然后，我们构造了一个层次图神经网络来学习潜在的局部模式，用于结合残基预测。设计了边更新模块、节点更新模块和图更新模块，学习目标残基的高层几何和生物物理化学特征，并对目标残基进行固定大小的嵌入。此外，门控递归单元(25)用于堆叠多个GNN块，利用了所有块的信息，避免了梯度消失问题。实验结果表明，GraphBind在核酸结合残基预测方面具有较好的性能。此外，我们还证明了GraphBind可以扩展到其他配体结合残基的预测，具有良好的性能。

MATERIALS AND METHODS

在本节中，构建了两个基准数据集来评估GraphBind的性能。然后介绍了HGNN的图结构和体系结构。最后，对评估方案和详细的实验环境进行了简要的总结

Benchmark datasets

为了评估GraphBind的性能，并将其与其他方法进行公平的比较，我们从BioLiP数据库(26)中构建了两个核酸结合蛋白基准数据集，并根据发布日期将它们分为训练集和测试集。基准数据集可在http://www.csBio.sjtu.edu.cn/Bioinf/GraphBind/Datets.html上获得。

DNA/RNA结合蛋白收集自2018年12月5日发布的BioLiP数据库。这个数据库是生物上相关的配体蛋白相互作用的集合，这些相互作用在复合物中结构上被解决。如果靶残基与核酸分子之间的最小原子距离小于0.5˚A加上最近的两个原子范德华半径之和，则定义为结合残基。

BioLiP包含2018年12月5日6,342个核酸-蛋白质复合体中的48133个核酸结合位点。这些复合体分为4344个DNA-蛋白质复合体(9574个DNA结合蛋白链)、1558个RNA蛋白质复合体(7693个RNA结合蛋白链)和440个DNA-RNA-蛋白质复合体。我们排除了DNARNA-蛋白质复合物以避免混淆，因为在BioLip数据库中没有做任何注释来区分DNA或RNA结合残基。根据发布日期，2016年1月6日之前发布的蛋白质链被分配到原始训练集(6731条DNA结合蛋白链和6426条RNA结合蛋白链)，其余的蛋白质链被分配到原始测试集中(2843条DNA结合蛋白链和1267条RN绑定蛋白链)。

由于DNA/RNA结合残基预测存在数据不平衡问题，即DNA/RNA结合残基的数量远小于非结合残基的数量，因此我们将数据增强方法应用于原始训练集。在先前的研究(3，4，27-29)之后，我们从相似的蛋白链中转移结合注释，以增加训练集中的结合残基数量，原因如下：(I)具有相似序列和结构的蛋白质，虽然可以来自不同的生物体，但可能具有相同的生物学功能；(Ii)不同的分辨率可能导致同一蛋白质的结构略有不同。为此，我们首先应用bl2seq(30)(E值<0.001)和TM-Align(31)来评估蛋白质链对之间的序列同一性和结构相似性。其次，我们对序列同一性>0.8和TM得分>0.5的链进行聚类。第三，将同一簇中蛋白质链的注释转移到残基数量最多的链中。在转移结合注释后，我们进一步用CD-HIT(32)去除多余的蛋白质链，将训练集中的序列同一性降低到30%以下。最后，我们得到了573条DNA结合蛋白链和495条RNA结合蛋白链作为训练集。数据增强使DNA结合残基和RNA结合残基分别增加了30.7%和24.3%。去除从原始DNA/RNA结合测试集合到DNA/RNA结合训练集合中的任意链的蛋白质链，所述原始DNA/RNA结合测试集合具有CD-HIT(32)测量的超过30%的序列同一性。最后，我们分别获得129个DNA结合蛋白和117个RNA结合蛋白作为DNA和RNA结合测试集。表1汇总了数据集的详细信息(未增加数据的训练集见补充表S1)。

Graph construction based on structural contexts

提取了基于序列和基于结构的多种特征，包括伪位置特征、残基原子特征、二级结构特征和进化会话特征。然后，使用3D空间中定义的滑动球来提取以残基为中心的目标残基的结构上下文。基于残基在结构上下文中的伪位置计算的邻接矩阵被用来构造图。此外，在节点和边缘特征向量中嵌入了几何知识和生物物理化学特征。图1显示了图构建的管道。

图1.GraphBind中使用的图构造管道。它包括三个模块：特征提取、结构上下文提取和图构造。(A)特征提取。从蛋白质结构中提取残基的伪位置和原子特征。利用DSSP、PSI-BLAST和HHblits从蛋白质结构和序列中提取二级结构特征和进化会话特征。(B)结构上下文提取。目标残基的结构上下文由以残基为中心的预定半径的滑动球体确定。(c )图构造。结构上下文进一步由图$G=(V,E,u,A)$表示。V、E、u和A分别表示节点的特征向量集合、边的特征向量集合、图特征向量和邻接矩阵。图中的节点表示残基。节点I的原始特征向量$V_i^{raw}\in {\mathbb{R}}^{72}$是节点I的位置嵌入和残基特征的拼接，距离矩阵是根据残基的伪位计算的。在距离矩阵上应用二进制阈值$r_v$，得到记录节点连接的邻接矩阵A。边ij的原始特征向量$e_{ij}^{raw}\in {\mathbb{R}}^2$分别由两个相邻结点之间的欧几里得距离和两个向量之间从球心到两个相邻结点的夹角θij的余弦编码

特征提取。推导出四种类型的残基水平特征如下：

第一个是伪位置。包括残基的主链和侧链原子的残基的质心被表示为该残基的伪位置，因为蛋白质和核酸之间的相互作用可以发生在主链和侧链原子(33)上。

二是残基的原子特征。对于残基，我们提取了属于残基(不包括氢原子)的每个原子的以下七种特征：原子质量、B因子、是否为残基侧链原子、电荷、与其键合的氢原子数、是否在环中以及原子的范德华半径。残基的原始原子特征表示为 { f s , t } s = 1 , … , 7 , t = 1 , … , N a \left\{f_{s, t}\right\}_{s}=1, \ldots, 7, t=1, \ldots, N_{a} { fs,t​}s​=1,…,7,t=1,…,Na​， { f s , t } \{f_{s, t}\} { fs,t​}代表第t-th原子的某些特征，$N_a$代表属于残基的原子数。由于不同的残基可能具有不同的原子数，我们将所有原子的某物特征平均为残基的加工的某物原子特征xs，这导致对于每个残基 { X s } s \{X_s\}_s { Xs​}s​=1，…，7：

$$x_s=\frac{1}{N_a}\left(\sum _{t=1}^{t=N_a}{f}_{s,t}\right)$$

最后，我们为含有$L$个残基的查询蛋白质生成一个$L\times 7$的原子特征矩阵。三是二级结构剖面。DSSP(34，35)生成$L\times 14$矩阵形式的二级结构轮廓，包括1列渣油暴露表面，5列键角和扭角，8列8种状态的单热编码二级结构。二级结构的8种状态包括B(residue in isolated β-bridge),E(extended strand, participates in β-ladder), G(310-helix),H(α-helix), I(π-helix), S(bend), T(H-bonded turn) 和其他.

最后是两种进化的对话模式。
(1)PSI-BLAST配置文件。比对工具PSI-BLAST使用启发式算法和动态规划来搜索NCBI的非冗余数据库(NR)，以查找具有三次迭代且$E_{value}<10^{-3}$(36)的同源序列。生成的特定位置评分矩阵(PSSM)的大小为$L\times 20$。PSSM中的每个元素$x$通过Sigmoid函数归一化到范围[0，1]：
$$\bar{x}=\frac{1}{1+e^{-x}}$$

(2)HHblits剖面。基于隐马尔可夫模型(HMM)的HHblits被用于针对具有默认参数的uniclust30数据库进行搜索，以生成查询序列的HMM矩阵(37)。HMM矩阵的大小为$L\times 30$。HMM矩阵由同源序列中20种氨基酸的20列观测频率、7列转换频率和3列局部差异组成。每个分数都转换为范围[0，1]：

$$\bar{x}=\frac{1}{10000}$$
PSI-BLAST和HHblits配置文件是互补的，因为它们的后端算法和搜索的数据库是不同的，这在我们随后的实验中得到了证实。

综上所述，对于一个查询蛋白质，我们得到了大小为$L\times 3$的伪位矩阵和大小为$L\times 71$的特征矩阵。对于特征矩阵中的每一列，执行最小-最大归一化以将该值线性归一化为[0，1]：

$$\bar{x}=\frac{x-x_{min}}{x_{max}-x_{min}}$$

其中，$x_{min}和x_{max}$值分别是训练集中该特征的最小值和最大值

结构化上下文提取。根据三级结构中残基的伪位，一个球体沿着多肽链滑动，以获得每个残基的结构上下文。对于目标残基，结构上下文被定义为以该残基为中心半径为$r_g$的球体。球体中的所有残基及其几何知识形成了目标残基的局部结构上下文。与蛋白质的整体结构相比，结合位点通常更多地与其局部结构环境的几何和生物物理化学性质有关(8-9，15)。

图构造。在此步骤中，残基的结构上下文进一步表示为图。图$G=(V，E，u，A)$，其中 V = { v i } i = 1 , … , N v V=\left\{\boldsymbol{v}_{i}\right\}_{i=1, \ldots, N_{v}} V={ vi​}i=1,…,Nv​​和${\mathbf{v}}_i\in {\mathbb{R}}^{D_v}$分别表示$N_v$节点的特征向量集合和节点i的特征向量。添加$N_v\times N_v$形状的邻接矩阵。 E = { e i j ∣ A i j = 1 } E=\left\{\boldsymbol{e}_{i j} \mid \boldsymbol{A}_{i j}=1\right\} E={ eij​∣Aij​=1}表示$N_e$的特征向量集合。${\mathbf{e}}_{\mathbf{i}\mathbf{j}}\in {\mathbb{R}}^{D_e}$表示节点i和j之间的边ij的特征向量。${\mathbf{e}}_{ij}\in E$ if ${\mathbf{A}}_{ij}=1,{\mathbf{e}}_{ij}\notin E$ if ${\mathbf{A}}_{ij}=0$。U代表图特征向量。在图中，残基表示为节点。由相应残基的伪位置定义的第i个节点$P_i$的位置。靶残基周围的残基可能形成特定的局部几何图案，这对结合残基的识别是有用的。基于这一观察，我们使用位置嵌入来表示目标残基与其每个上下文残基之间的位置关系，因为它包含目标残基周围的局部几何知识。将节点i的位置嵌入定义为节点$ii$与球心之间的归一化欧几里德距离，
$$P{E}_i=\frac{1}{r_g}\left|\overrightarrow{{\mathbf{p}}_o{\vec{p}}_i}\right|$$

其中$P_o和P_i$分别表示球心和节点i的位置，$\overrightarrow{{\mathbf{p}}_o{\mathbf{p}}_i}$ 是从$P_o$到$P_i$的向量。节点i的原始特征向量 v i raw  ∈ R 72 \boldsymbol{v}_{i}^{\text {raw }} \in \mathbb{R}^{72} viraw </