GraphDTA论文阅读小白笔记（附代码注释和复现流程）_dgraphdta论文

摘要

具体代码复现以及代码注释可以查看https://github.com/zhaolongNCU/Demo-GraphDTA-

背景

1.发展前景：

• 新药设计需要花费2.6billion，17years
• 药物再利用已被用于现实的疾病中
• 为了有效地重新调整药物的用途，了解哪些蛋白质是哪些药物的靶标是有用的
• 高通量筛选方法高消耗，并且彻底地完成搜索空间的搜索几乎不可能
• 因此有强烈的动机构建新的计算方法去基于以前的药物靶标实验评估新的药物靶标对的结合作用强度

2.先前的计算方法：

• molecule docking

  • 优点在于提供更多信息，缺点在于需要晶体3D结构

• collaborative filtering

  • 协同过滤（Collaborative filtering）是一种常见的推荐系统算法，它利用用户之间的相似性和物品之间的相似性，预测用户对未评价物品的评分或偏好。在协同过滤中，用户可以为物品评分或提供反馈，该反馈信息将用于计算与其他用户或物品之间的相似性，从而生成推荐列表。

    协同过滤可以分为两种类型：基于用户的协同过滤和基于物品的协同过滤。基于用户的协同过滤通过计算用户之间的相似性来预测一个用户对一个物品的评分或偏好。基于物品的协同过滤则计算物品之间的相似性来推荐与用户评分过的物品相似的其他物品。

    协同过滤算法的主要优点是能够在不需要关于物品的元数据或用户个人信息的情况下生成推荐列表。它也可以应用于各种类型的物品和各种类型的用户，使其在电商、社交媒体和音乐等领域得到广泛应用

  • SimBoost：

    • 计算药物和药物，靶标和靶标的相似性并将其作为输入到梯度增强中，以预测未知药物靶点结合亲和力

    • 具体来说，SimBoost使用了两个核矩阵，分别用于表示化合物和蛋白质的特征。然后，它使用Boosting算法来递归地构建一个加权的核矩阵，该加权核矩阵能够更好地描述化合物和蛋白质之间的相互作用

    • 在每个Boosting迭代中，SimBoost将化合物和蛋白质的核矩阵进行组合，并计算出加权核矩阵。然后，SimBoost根据加权核矩阵的结果对训练样本进行重新加权，以便更好地处理样本之间的相关性。

    • 那么这里面的核矩阵是怎么定义的呢？

      • 对于化合物核矩阵也可以叫化合物相似度矩阵，其中每个元素表示两个化合物之间的相似度，这个相似度可以基于分子描述符，结构或者其他化学性质进行计算。一般来说化合物相似度矩阵越高，表明两个化合物越相似。

        • 举个简单的例子：

          假设我们有三个化合物A、B、C，我们可以使用分子指纹来计算化合物之间的相似度，我们假设用长度为3的哈希分子指纹来计算相似度。我们可以用一个3位的二进制来表示每个分子指纹，其中1表示该分子指纹存在，0表示不存在对于化合物A它的分子指纹为001，B为011，通过汉明距离计算出化合物之间的相似度，汉明距离值得是两个字符串之间不同字符的数量。如上AB之间只有一个位置不同，汉明距离为1，通过计算上化合物之间的分子指纹相似度，我们可以得到一个化合物-化合物相似度矩阵

                A    B    C
             A  1    1/3  0
             B 1/3   1    1/3
             C  0   1/3   1
          
          1
          2
          3
          4
          5

          这个矩阵中，每个元素表示两个化合物之间的相似度，值越大表示相似度越高。例如，矩阵中的第一个元素1表示化合物A和化合物A之间的相似度为1，即两个化合物完全相同。而矩阵中的第二个元素1/3表示化合物A和化合物B之间的相似度为1/3，即两个化合物之间只有一个分子指纹相同。

        • SimBoost模型中运用的是Tanimoto相似度来计算化合物之间的相似度

          • Tanimoto相似度（也称为Jaccard相似系数）是一种用于比较两个集合相似度的指标，通常应用于化学分子结构相似性的计算中。假设有两个集合A和B，它们的Tanimoto相似度可以通过以下公式计算：
            $$T(A,B)=\frac{|A\cap B|}{|A\cup B|}$$
            其中，|A ∩ B|表示A和B的交集元素个数，|A ∪ B|表示A和B的并集元素个数。

            在化学领域中，Tanimoto相似度通常用于计算化合物指纹之间的相似性。化合物指纹是一种将化合物结构信息编码为二进制向量的方法。Tanimoto相似度可用于计算两个化合物指纹之间的相似性。Tanimoto相似度的值介于0和1之间，数值越大表示两个化合物指纹越相似。通常，Tanimoto相似度大于0.5时可以认为化合物之间存在结构相似性。

          • Gower相似度（Gower coefficient）是一种用于比较两个向量之间的相似度的指标，它适用于不同类型的特征（如数值型、类别型和二元型特征）同时存在的情况。Gower相似度计算的思想是将两个向量中的每个特征（或属性）按其类型进行不同的权重计算，并将所有的权重加起来得到一个综合的相似度指标。
            $$Gower(A,B)=\frac{{\sum }_{i=1}^n{w}_i\cdot {\delta }_i}{{\sum }_{i=1}^n{w}_i}$$
            谷本系数是Gower相似度的一种特殊情况，即当所有特征都为二元型特征时，Gower相似度就等同于谷本系数。在化学领域中，Tanimoto相似度通常用于计算化合物指纹之间的相似性，而Gower相似度通常用于比较化学分子之间不同类型特征的相似度。

      • 蛋白质-蛋白质相似矩阵通过蛋白质序列相似性来度量

        • 本文使用的是Smith-Waterman相似性来度量，这个算法能比较两个蛋白质序列之间的相似性，并给出一个相似性得分。

          • Smith-Waterman算法是一种用于比对两个序列相似性的局部比对算法。与全局比对算法（如Needleman-Wunsch算法）不同，局部比对算法能够在两个序列的一部分区域中找到最佳的匹配。

            假设有两个序列：ACGTGA和CGTACG，现在需要找到它们之间的相似片段并计算比对得分。以下是Smith-Waterman算法的计算过程：

            1. 创建一个得分矩阵，用于记录两个序列的每个位置的比对得分。矩阵的第一行和第一列分别填充为0，表示序列1或序列2的前缀序列与空序列的比对得分为0。

               		A	C	G	T	G	A
               	0	0	0	0	0	0	0
               C	0
               G	0
               T	0
               A	0
               C	0
               G	0

            2. 对于每个矩阵中的非边界位置（即第一行和第一列以外的位置），计算该位置的得分。计算公式如下：
               $$Score(i,j)=max(0,Score(i-1,j-1)+S(x_i,y_j),Score(i-1,j)+d,Score(i,j-1)+d)$$
               其中，Score(i, j)表示矩阵中第i行、第j列的得分；S(x_i, y_j)表示序列1中第i个字符和序列2中第j个字符的匹配得分；d表示删除/插入一个字符的得分惩罚。

               对于本例中的第二行第二列的位置（即第一个字符的比对），有以下的比对得分：
               $$Score(2,2)=max(0,Score(1,1)+S(C,C),Score(1,2)+d,Score(2,1)+d)=max(0,0+2,0-1,0-1)=2$$

               		A	C	G	T	G	A
               	0	0	0	0	0	0	0
               C	0	0	1	0	0	0	0
               G	0	0	0	2	1	3	1
               T	0	0	0	1	3	2	0
               A	0	1	0	0	2	1	4
               C	0	0	2	1	1	1	2
               G	0	0	1	3	2	4	3

  • KronRLS：

    • 先来简单记录下Kronecker积：

      具体来说，设A为mxn矩阵，B为pxq矩阵，它们的Kronecker积记为：
      $$A\otimes B$$
      则其为一个mpxnq的矩阵，其满足：
      A ⊗ B = [ a 11 B ⋯ a 1 n B   ⋮ ⋱ ⋮   a m 1 B ⋯ a m n B ] A\otimes B=

      $$\begin{bmatrix} a_{11}B & \cdots & a_{1n}B \ \vdots & \ddots & \vdots \ a_{m1}B & \cdots & a_{mn}B \end{bmatrix}$$ A⊗B=[a11​B​⋯​a1n​B ⋮​⋱​⋮ am1​B​⋯​amn​B​]
      可以看出Kronecker积是用A的每个元素和B相乘

      下面举个简单的例子来说明Kronecker积：
      假设有矩阵 A = [ 1 2   3 4 ] 和矩阵 B = [ 0 5   6 7 ] 假设有矩阵A=

      $$\begin{bmatrix} 1 & 2 \ 3 & 4 \end{bmatrix}$$ 和矩阵B= $$\begin{bmatrix} 0 & 5 \ 6 & 7 \end{bmatrix}$$ 假设有矩阵A=[1​2 3​4​]和矩阵B=[0​5 6​7​]

      他们的 K r o n e c k e r 积为： A ⊗ B = [ 1 × 0 1 × 5 2 × 0 2 × 5   1 × 6 1 × 7 2 × 6 2 × 7   3 × 0 3 × 5 4 × 0 4 × 5   3 × 6 3 × 7 4 × 6 4 × 7 ] = [ 0 5 0 10   6 7 12 14   0 15 0 20   18 21 24 28 ] 他们的Kronecker积为： A\otimes B =

      $$\begin{bmatrix} 1\times 0 & 1\times 5 & 2\times 0 & 2\times 5 \ 1\times 6 & 1\times 7 & 2\times 6 & 2\times 7 \ 3\times 0 & 3\times 5 & 4\times 0 & 4\times 5 \ 3\times 6 & 3\times 7 & 4\times 6 & 4\times 7 \end{bmatrix}$$ = $$\begin{bmatrix} 0 & 5 & 0 & 10 \ 6 & 7 & 12 & 14 \ 0 & 15 & 0 & 20 \ 18 & 21 & 24 & 28 \end{bmatrix}$$ 他们的Kronecker积为：A⊗B=[1×0​1×5​2×0​2×5 1×6​1×7​2×6​2×7 3×0​3×5​4×0​4×5 3×6​3×7​4×6​4×7​]=[0​5​0​10 6​7​12​14 0​15​0​20 18​21​24​28​]

    • 具体来说，KronRLS模型首先使用药物和蛋白质的分子描述符来构建核，即两个分子之间的相似性度量。这些核可以是任何相似性度量，例如径向基函数（RBF）核。然后，使用Kronecker积来计算药物-药物和蛋白质-蛋白质核之间的关系，生成一个DTI预测的核矩阵K。KronRLS模型在这个核矩阵K上执行RLS回归，以预测药物-靶点相互作用。

    • KronRLS模型的优点是能够有效地处理大规模的DTI预测问题，同时具有良好的预测性能。该模型的缺点是需要大量的计算资源，因为Kronecker积计算的复杂度很高。

3.现在的一些常用方法：

• DeepDTA：通过1D卷积层来提取蛋白质和化合物特征，再进行回归预测

• WideDTA：是DeepDTA的扩展，其中化合物和蛋白质被总结为更高阶的特征，例如药物最常见的子结构Ligand Maximum Common Substructures（LMCS），蛋白质最保守的子序列the Protein Domain profiles or Motifs （PDM），其他和DeepDTA基本一致
  

  • Ligand Maximum Common Substructures (LMCS)：是一种用于表示化合物分子结构的方法。该方法将一组分子中最常见的子结构组合成一个共同的子结构，以此来表示这组分子的结构特征。这种表示方法能够使得机器学习算法更好地捕捉到分子之间的相似性和差异性，有助于预测分子与蛋白质的相互作用。在药物发现领域，LMCS常用于描述分子之间的相似性和药效团特征，是一种重要的分子表示方法。

    举个例子：

    Molecule A: CC(=O)Nc1ccc(cc1)C(O)CNC(=O)C2=C(O)C=CC=C2
    Molecule B: CC(C)(C)Nc1cc(C(N)=O)ccc1C(C)(C)C
    1
    2

    在这个例子中，Molecule A和Molecule B都包含了苯环，且它们都连接到了一个氨基上，这是它们的共同特征。因此，可以将这个共同特征作为它们的LMCS。如果将所有的药物分子都找到它们的LMCS，那么这些LMCS就可以用来表示每个药物分子的结构。在表示药物分子的过程中，LMCS可以用来减少描述药物分子所需的信息量，从而简化模型并提高模型的训练速度。

  • the Protein Domain profiles or Motifs （PDM）：蛋白质结构由许多相对独立的功能和结构单元——蛋白质域组成。蛋白质域的序列和结构特征是它们的特定功能和相互作用的基础。Prosite是一种用于识别蛋白质序列中域、重复单元和其他功能单元的资源。该数据库通过将域和模体信息存储为PDM文件来描述它们。每个PDM文件包含有关单个蛋白质域的域边界、保守序列区域以及与特定蛋白质功能相关的其他信息。可以使用这些文件来比较新蛋白质序列和已知的域和模体，以确定新序列中是否存在类似的功能单元。

4.思考

• 深度学习确实是DTA预测中最好的模型之一，但是这些模型将药物表示为字符串，这不是表示分子的较为自然的方法，当使用smiles时分子的结构信息会丢失，这可能会影响模型预测能力。
• 已经有图卷积神经网络被用于药物发现、相互作用预测、合成预测、de Novo分子设计和定量结构预测，但是图神经网络还未用于药物靶点相互作用预测。

5.本文所做

• 提出新的模型GraphDTA能够将化合物表示为分子图，并且预测性能比最新的深度学习模型基准都要好。
• 另外，为了更好的理解我们的基于图形的模型是如何工作的，我们对模型的潜在空间进行了多元统计分析（用来了解他们的模型如何使用图形表示来自动分配药物特征的重要性，帮助我们更好地理解数据中变量之间的关系，为模型训练和优化提供指导），发现了隐藏节点的激活与药物的一些特定化学或结构特征的描述，例如药物中具有的特定官能团的数量存在相关性，这表明图神经网络可以自动将重要性分配给具有特定化学特征的药物，而不需要任何先验的知识。
• 同时本文发现一些药物对总体预测误差的贡献不成比例，即它们的预测误差显著高于其他药物，然而，对这些药物做箱线图进行异常点分析，发现他们并不是异常点，即它们的化学特征与其他药物类似，并没有被判定为异常点，这表明该模型对于这些药物化合物的预测能力存在一定局限，需要进一步的探究和改进。
• 综上，该模型具有高度的准确性和一定的实际意义，但也有预测失败的问题。最后作者讨论这些见解可以帮助研究者更好地利用深度学习模型进行药物发现，并有望为药物研究提供新的思路方法。

数据和方法

GraphDTA概述

本文模型将药物表示为分子图，利用图卷积神经网络直接学习原子之间的键特征，蛋白质特征依旧用CNN去提取。

药物表征

1.SMILES字符串编码基本规则：

• 原子符号：使用元素符号来代表原子，如C表示碳，O表示氧，N表示氮等。
• 隐式氢原子：如果一个原子是一个简单的原子（如碳、氮、氧等），则SMILES字符串中默认情况下它会有一定数量的隐式氢原子。隐式氢原子可以用小写字母h来表示，如C表示一个碳原子，C1表示一个碳原子带有一个氢原子，C=表示一个碳原子带有两个隐式氢原子。
• 连接符号：使用连接符来表示化学键，如单键使用“-”表示，双键使用“=”，三键使用“#”。
• 环：用数字表示一个环，数字代表着该原子在环中的顺序编号。例如，C1CCCCC1代表着一个含有6个碳原子的环。
• 分支：使用小括号来表示分支结构，括号内的分子结构是分支结构，括号后面跟着数字表示分支结构连接的原子编号。例如，CC©C表示一个分子由一个中心碳原子连接着三个碳原子组成的分支结构。

​ SMLIES实现了计算机可读分子，快捷检索，还可以将其视为字符串用NLP处理表征，或者之间在卷积神经网络CNN中使用。

2.但是，本文将药物视为原子之间相互作用形成的图表，每个节点为一个原子，每个节点是表达五条信息的多位二进制特征向量：原子符号、相邻原子原子数量、相邻原子氢原子数量、原子的隐含价、原子是否在芳香结构中。具体的特征表示可以看下边代码理解：

• 原子隐含价：在化学中，每个原子都有一定的化学价（valency），即能够和多少个其他原子形成化学键。有时，原子上可能有其他的键，而这些键又不会参与化学反应，这些被称为隐含键（implicit bond）。因此，隐含价是指一个原子的实际价与它已经形成的化学键的总价之差。例如，一个氧原子通常的化学价为2，如果它已经形成了一个双键，那么它的隐含价就是0。在SMILES表示法中，这个隐含价可以表示为一个整数值。

• 芳香结构：通常由含有苯环或其它类似结构的分子组成。具有芳香性的化合物具有独特的化学性质和生物活性，因此在药物研究和设计中经常被使用。在化学上，芳香结构具有一定的共振稳定性，使得其中的电子能够比较自由地运动。这种稳定性也使得芳香结构对化学反应有一定的影响。

  • 常见的芳香结构：
    • 苯环：由六个碳原子组成，每个碳原子都与一个氢原子相连，分子中有共轭的π电子体系，呈平面六边形。
    • 噻吩环：由一个硫原子和一个氮原子组成的五元环。
    • 咔唑环：由一个氮原子和一个碳原子组成的五元环。
    • 噻唑环：由一个硫原子和一个氮原子组成的五元环。
    • 吡啶环：由一个氮原子和五个碳原子组成的六元环。
    • 哒唑环：由一个氮原子和两个碳原子组成的五元环。
    • 呋喃环：由一个氧原子和四个碳原子组成的五元环。

• 节点特征代码：

  def atom_features(atom):
      return np.array(one_of_k_encoding_unk(atom.GetSymbol(),['C', 'N', 'O', 'S', 'F', 'Si', 'P', 'Cl', 'Br', 'Mg', 'Na','Ca', 'Fe', 'As', 'Al', 'I', 'B', 'V', 'K', 'Tl', 'Yb','Sb', 'Sn', 'Ag', 'Pd', 'Co', 'Se', 'Ti', 'Zn', 'H','Li', 'Ge', 'Cu', 'Au', 'Ni', 'Cd', 'In', 'Mn', 'Zr','Cr', 'Pt', 'Hg', 'Pb', 'Unknown']) +
                      one_of_k_encoding(atom.GetDegree(), [0, 1, 2, 3, 4, 5, 6,7,8,9,10]) +
                      one_of_k_encoding_unk(atom.GetTotalNumHs(), [0, 1, 2, 3, 4, 5, 6,7,8,9,10]) +
                      one_of_k_encoding_unk(atom.GetImplicitValence(), [0, 1, 2, 3, 4, 5, 6,7,8,9,10]) +
                      [atom.GetIsAromatic()])
  """
  
  这段代码定义了一个函数 atom_features，它接受一个分子中的原子作为输入，并返回一个表示该原子特征的 NumPy 数组。
  
  具体来说，该函数使用 one_of_k_encoding_unk 函数对原子的符号、节点数、氢原子数、隐式化合价、原子是否在芳香环内等属性进行编码
  生成相应的 one-hot 向量，并将这些向量拼接在一起形成特征向量
  其中，one_of_k_encoding_unk 函数将离散属性转化为 one-hot 向量，未知属性使用 all_unk 来编码。
  
  atom_features()函数的作用是将节点（原子）的多个属性编码成一个特征向量，为后续的分子学习任务提供输入X
  
  """
  1
  2
  3
  4
  5
  6
  7
  8
  9
  10
  11
  12
  13
  14
  15
  16
  17

蛋白表征

1.虽然蛋白质也可以用图来表征，但由于其三级结构有些不是和可靠，所以本文选择用ont-hot编码来对其蛋白进行表征。

2.从uniprot数据库中获得蛋白质序列，用整数编码例如，丙氨酸(A)是1，胱氨酸©是3，天冬氨酸(D)是4，依此类推，用如下代代码进行处理为one-hot编码形式：

#蛋白氨基酸序列转化为one-hot编码
def seq_cat(prot):
    x=np.zeros(max_seq_len)                        #生成全为0的numpy数组用来存储编码结果
    for i,ch in enumerate(prot[:max_seq_len]):     #遍历蛋白质序列并生成元素和元素索引
        x[i]=seq_dict[ch]                          #并将该编码所对应的整数值赋值给x数组中对应位置的元素
    return x

seq_voc = "ABCDEFGHIKLMNOPQRSTUVWXYZ"              #常见的氨基酸序列
seq_dict = {v:(i+1) for i,v in enumerate(seq_voc)} #生成氨基酸序列对应字典，比如A对应1，B对应2，方便电脑运算
seq_dict_len = len(seq_dict)                       #计算字典中键的数量并赋值为seq_dict_len
max_seq_len = 1000                                 #氨基酸序列的最大长度设置为1000

"""
这段代码定义了一个名为 seq_cat 的函数，用于将蛋白质序列编码成一个固定长度的向量。
该函数需要在外部定义一个名为 max_seq_len 的变量，指定输出向量的长度。

在函数实现中，代码首先创建了一个长度为 max_seq_len 的全0 NumPy 数组 x，用于存储编码结果。
然后，代码使用 enumerate() 函数遍历蛋白质序列的前 max_seq_len 个氨基酸，
对于每个氨基酸，代码在 seq_dict 字典中查找其整数编码，并将该编码赋值给 x 数组中对应位置的元素。最后，函数返回编码后的数组 x。

另外，这段代码还定义了一个名为 seq_voc 的字符串，其中包含了26个字母的缩写，代表蛋白质中常见的氨基酸。
接下来，代码使用一个字典推导式将这些字母映射到整数。具体地，代码使用 enumerate() 函数将 seq_voc 中的每个字母和一个整数索引对应起来。
由于字典的键必须是唯一的，因此在这里我们使用每个字母的整数索引加1作为其在字典中的值，以避免出现重复的键。
最后，代码使用 len() 函数计算了 seq_dict 字典中键的数量，并将其赋值给变量 seq_dict_len。

例如氨基酸序列为"ABCDEFGHIKLMNOPQRSTUVWXYZ" 则生成
array([ 1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10., 11., 12., 13.,
       14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25.,  0.,....])
1行1000列的numpy数组

"""
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

分子图的深度学习

1.将CNN推广到图的过程中有两个主要的挑战：

• 数据点不是排列成欧几里得的网格，感受野怎么设置？

  • 理解下接受域：

    在传统的卷积神经网络中，输入的数据通常是排列成欧几里得网格的，例如图像数据就是二维的矩阵，文本数据可以表示为一维的序列。在这种情况下，卷积操作可以通过一个固定大小的窗口在数据上滑动，提取局部特征，这个窗口称为感受野(receptive field)。

    但是，在一些应用中，数据是以图形的形式表示，例如化学分子、社交网络等。在这些数据中，节点和边通常没有固定的排列方式，因此不能像传统的欧几里得网格数据那样使用卷积操作进行特征提取。为了解决这个问题，需要设计一种新的操作来构建感受野。

    意思就是说在CNN中可以用卷积块框数据，但是图数据无法用卷积块框数据，因为数据是节点数据，所以开发了邻居节点定义为新的感受野。

• 下采样的图形池化操作：如何将数据特征缩小为更小的表示（常见的图池化方法有max-pooling和average-pooling）,在池化操作时需要考虑保留重要的节点和边信息，避免信息丢失。

2.下图显示了本文的模型示意图，其中蛋白质先被转化为one-hot编码再通过嵌入层转化为128维向量，将每个符号表示为一个连续的向量，将蛋白质数据转化为（1000,128）维向量表示存储在矩阵中。

GCN

1.模型介绍：

• 形式化地，给定一个药物的图形表示为G = (V, E)，其中V是由N个节点组成的集合，每个节点由一个C维向量表示，E是由邻接矩阵A表示的边的集合。多层图卷积网络接受节点特征矩阵X∈RN×C（其中N = |V|，C表示每个节点的特征数量）和邻接矩阵A∈RN×N作为输入，然后生成一个节点级输出Z∈RN×F（其中F表示每个节点的输出特征数量）。并且为了保持稳定性，通常将数据进行归一化处理。传播规则为：

​ H(l+1) = σ(~D⁻¹/₂ ~A ~D⁻¹/₂ H(l) W(l))

​ 其中 H(l) 表示第 l 层节点的特征矩阵，W(l) 表示第 l 层的可训练参数，~A 是增广邻接矩阵（即加上了自连接的邻接矩阵），~D 是对角线矩阵，其第 i 个对角线元素是节点 i 的度数。

​ 在传播规则中，首先通过对 H(l) 乘以 W(l) 进行线性变换，得到了一个新的 特征矩阵。然后通过将邻接矩阵 ~A 进行归一化处理，即 ~D⁻¹/₂ ~A ~D⁻¹/ ₂，来计算节点与邻居节点的加权和。最后再通过一个激活函数 σ 进行非线 性变换，得到了下一层节点的特征矩阵 H(l+1)。

​ 通过不断地重复这个传播规则，我们就可以在 GCN 网络中学习到节点的特 征表示，这些特征表示可以用于节点分类、图分类等任务。

• 再来简单介绍下图中常见的几个矩阵：
  • 特征矩阵：在图数据中，特征矩阵通常是一个NXF的矩阵，其中 N表示图中节点的数量，F 表示每个节点所包含的特征的维数。每个节点的特征可以是任意类型的数据，例如节点的度、节点的标签、节点的位置、节点的向量表示等等。这些特征可以被组织成一个特征向量，并在特征矩阵的每一行中进行存储。在图卷积神经网络中，特征矩阵是输入的基本数据结构，用于表示图中节点的特征信息。
  • 邻接矩阵：在图数据中，邻接矩阵（adjacency matrix）是用于表示节点之间连接关系的矩阵。对于一个包含N个节点的图，邻接矩阵是一个NXN的矩阵A，其中第i行第j列的元素A[i,j]表示节点i和节点j之间是否存在一条边。如果存在一条边，则A[i,j]通常为1或者一个权重值；否则为0或者没有权重值。对于无向图，邻接矩阵是对称的；对于有向图，则不一定对称。在某些情况下，邻接矩阵也可以被扩展为包含自环（self-loop）的情况，即图中每个节点都与自身相连，也就是自连接。
  • 度矩阵：对于无向图来说度矩阵就是个对角矩阵对角元素为与这个节点相连边的数量。

2.本文使用三个连续的GCN层，每个用ReLU函数激活，然后添加一个全局最大池化层来获取图形的特征向量

GAT

1.图注意力网络是一种用于处理图数据的神经网络模型，它可以学习出每个节点的重要性，并利用这些重要性为每个节点赋予不同的权重，从而在图中进行信息传播和节点分类等任务。它的核心思想是使用自注意力机制，使每个节点能够动态地选择和聚合其邻居节点的特征信息，而不是简单地对所有邻居节点进行加权平均。这使得模型更加灵活和精确，能够捕捉到节点之间的复杂依赖关系，并在处理具有不同规模和结构的图时表现出更好的性能。

GAT的基本单元是图注意力层（Graph Attention Layer），该层对每个节点进行线性变换，然后通过注意力系数进行加权，得到加权和作为该节点的嵌入表示。注意力系数计算使用了多头自注意力机制，可以通过学习得到不同的权重分配，使得不同节点之间的关系得到更准确的建模。

具体来说，GAT的基本原理是将节点之间的关系抽象为图，通过在图上进行自注意力聚合，将节点的特征向量进行更新。在GAT中，每个节点都有一个自注意力向量，表示该节点与其邻居节点之间的关系权重。这个自注意力向量可以通过学习得到，其中注意力权重是由节点的特征向量和相邻节点的特征向量之间的相似度计算得出。通过自注意力聚合，每个节点的特征向量都会被更新，同时保留了相邻节点的信息。最终，通过多层GAT的叠加，可以获得更高级别的节点表示。
$$h_i^{\prime }=\text{ReLU}(\sum _{j\in \mathcal{N}i}\alpha ijW{h}_j)$$

$$其中，h_i表示节点i的输入特征，h_i^{\prime }表示节点i的输出特征，W是一个权重矩阵，{\alpha }_{ij}表示节点i与节点j之间的注意力系数，{\mathcal{N}}_i是节点i的一阶邻居节点集合。$$

图注意力网络的关键在于注意力系数的计算。常见的计算方式是通过对节点特征的内积来计算相似度，再通过 softmax 函数将其归一化，得到注意力系数。具体地，节点 i 与节点j之间的注意力系数可以表示为：
$${\alpha }_{ij}=\frac{\exp (\text{LeakyReLU}({\vec{a}}^T[W{h}_i|W{h}_j]))}{\sum _{k\in {\mathcal{N}}_i}\exp (\text{LeakyReLU}({\vec{a}}^T[W{h}_i|W{h}_k]))}$$

$$其中，\vec{a}是一个可学习的向量，|表示向量拼接操作，LeakyReLU是一个带有负半段的修正线性单元激活函数，用于处理注意力系数的非线性关系$$

2.本文使用两个GAT层，由一个Relu函数激活，然后跟随一个全局最大池化层来或得图提取的特征向量，对于第一个GAT层应用多头注意，其中头部被设置为10，并且输出特征数量被设置为输出特征数量，第二个GAT的输出特征设置为128。

GIN

1.GIN（Graph Isomorphism Network）是一种基于图同构性的图神经网络模型，由Xu等人于2019年提出。与传统的基于节点特征的图神经网络不同，GIN模型直接对节点的邻居特征进行聚合，从而更好地利用了图的局部结构信息。GIN模型的主要思路是通过对图中节点对子图同构性的判断，来学习节点的特征表示，进而实现节点分类、图分类等任务。

2.本文采用5个GIN层，每个层后边跟着一个批量归一化层，最后加入全局最大池化层得到图提取特征。

GAT-GCN

1.一个GAT加一个GCN再加一个全局最大池化

代码中是这样处理分子为图数据的：

'''

def smile_to_graph(smile):
    mol = Chem.MolFromSmiles(smile)                     #先转换化为可以处理的mol格式
    
    c_size = mol.GetNumAtoms()                          #计算分子原子个数
    
    features = []                                       #生成这个分子原子特征列表
    for atom in mol.GetAtoms():                         #遍历整个分子的原子
        feature = atom_features(atom)                   #计算出这个原子的特征
        features.append( feature / sum(feature) )       #归一化并加入特征里对应一个原子由78个特征

    edges = []                                          #生成这个分子化学键列表
    for bond in mol.GetBonds():                         #遍历整个分子的键
        edges.append([bond.GetBeginAtomIdx(), bond.GetEndAtomIdx()]) #获取这个键两端原子索引，然后将这两个索引保存在edges列表中
    g = nx.Graph(edges).to_directed()                   #利用networkx生成无向图，to/diected()将其变为有向图双向的
    edge_index = []                                     #生成边索引列表
    for e1, e2 in g.edges:                              #遍历整个图的边索引
        edge_index.append([e1, e2])                     #保存到边索引列表中
        
    return c_size, features, edge_index                 #返回分子原子个数，所有原子特征，边索引
"""
这段代码定义了一个名为 smile_to_graph 的函数，该函数接受一个 SMILES 格式的化合物表示字符串作为输入，并返回该分子的图形表示。
具体来说，该函数首先将输入的 SMILES 字符串转换为 RDKit 中的分子对象 mol
然后分别计算分子中每个原子的特征向量和分子中所有化学键的连接关系
最后返回一个包含分子大小、每个原子的特征向量以及连接关系的元组

具体来说，该函数的实现步骤如下：
将输入的 SMILES 字符串转换为 RDKit 中的分子对象 mol。
对于分子中的每个原子，调用 atom_features 函数计算其特征向量，然后将该向量除以其元素之和，得到归一化的特征向量，将所有归一化的特征向量保存在一个列表 features 中。

对于分子中的每个化学键，获取它的两端原子的索引，然后将这两个索引保存在一个边列表 edges 中。
使用 edges 列表创建一个有向图 g，并将图中所有的边转换为边索引列表 edge_index。

返回一个包含分子大小 c_size、每个原子的归一化特征向量列表 features 和连接关系的边索引列表 edge_index 的元组。

需要注意的是，该函数中的 atom_features 函数和 one_of_k_encoding_unk 函数是在之前的代码段中定义的，用于计算原子的特征向量和将输入值转换为 one-hot 向量。
另外，该函数使用了 NetworkX 库中的 Graph 类和 DiGraph 类来表示分子的连接关系，并使用 to_directed 方法将无向图转换为有向图。

例如
s='NC(C)C(=O)O' 
b=smile_to_graph('NC(C)C(=O)O')
生成的有向图索引为：
[[0, 1],
 [1, 0],
 [1, 2],
 [1, 3],
 [2, 1],
 [3, 1],
 [3, 4],
 [3, 5],
 [4, 3],
 [5, 3]]

"""

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56

基准

1.模型：DeepDTA、KronRLS、SimBoost、WideDTA

2.数据集：Davis（[72，442]，[5.0,10.8]），Kiba([2116，229]，[0,17.2])

3.超参数设置

模型解释

1.深度神经网络中的潜在变量（latent variables）和如何使用外部数据来分析这些变量的关系。通过图神经网络层获取了128个潜在变量，并通过冗余分析（redundancy analysis）直接分析了这些变量，以了解模型性能与领域知识之间的关系。这种多变量统计方法可以衡量外部数据源能够解释潜在变量总方差的百分比。在这种分析中，外部数据源是38个分子JoeLib特征/描述符的矩阵，每个药物都有一个这样的矩阵。同时，文章还将这些潜在变量的主成分值与每种药物的测试集误差进行比较。在这里，每种药物（或每种蛋白质）的误差是指所有包含该药物（或该蛋白质）的测试集对之间的预测DTA和真实DTA之间的绝对误差的中位数。在这些分析中，作者专注于使用GIN模型（因为它的性能更优）和Kiba数据集（因为它的药物目录更大）。

结果讨论

图模型的表现超过了其它模型

1.Davis：

MSE：所有4个图模型的MSE都是最低的，最好的是0.229比最好的基准0.261低了14%。

CI：四个图模型中只有两个2个具有最高的CI，基准模型最高CI是0.886，相比下GAT和GIN的CI分布为0.892和0.893。

2.KIBA:

MSE：四个图模型中有三个最低，这里最低的为0.139，比最好的基准模型0.179低了28.8%

3.在所有的图模型中GIN是唯一一个在两个数据集和两个性能指标上都具有最佳性能的模型，出于这个原因之后的统计分析都将把重点放在这个模型上。

4.自己复现的情况

Davis：批量化次数均设置为128，迭代次数设置为50，本文设置的是512、1000

由于KIBA数据集数据太多复现太慢算力不够暂不复现

图模型发现已知药物特性

1.图神经网络将每个药物分子抽象为一个新的特征向量，即128个潜在变量组成的特征向量。然后利用这些潜在向量与蛋白提取到的特征进行DTA预测，因此我们假设这些潜在变量代表对DTA有实质性贡献的图形特征。

2.但是比较难确定每个潜在变量对应的分子亚结构。因此，本文将学习到的潜在空间与已知的分子描述符矩阵（外部数据源是38个分子JoeLib特征/描述符的矩阵）进行回归分析，用来寻找它们之间的重叠，即所谓的冗余分析，它可以衡量外部数据源对神经网络潜在向量的解释能力。在这个例子中，外部数据源是每种药物的38个分子描述符（来源于ChemMine Tools）。分析发现，20.19%的潜在空间可以由已知描述符解释，其中“羟基化的脂肪链个数”对解释方差做出了最大贡献。此外，两个潜在变量V58和V14与该描述符之间的相关性最强，即当羟基化的脂肪链个数较多时，它们的激活程度较高。

• 冗余分析：冗余分析是一种多元统计方法，用于研究多个解释变量（在这个例子中是图神经网络中的潜在变量）对一个响应变量（在这个例子中是分子描述符）的影响程度。在本文中，冗余分析被用来测量外部数据源（即38个分子描述符）对于神经网络中的潜在变量的解释能力。

  具体操作流程如下

  1. 对于神经网络中的每个潜在变量，得到其在所有分子中的值，并将其构成一个向量。
  2. 对于每个分子，将其对应的38个分子描述符组成一个向量。
  3. 将所有分子的潜在变量和分子描述符向量分别进行归一化处理。
  4. 运行冗余分析，得到潜在变量和分子描述符向量之间的相关系数和贡献率。
  5. 根据贡献率排序，找到对潜在变量解释能力最强的分子描述符。

少数几种药物对总误差贡献不成比例

1.尽管GraphDTA模型由于其他模型，但是我们想更多的了解它的预测为什么失败，为此本文对Davis和Kiba测试集的每种药物和每种蛋白的预测误差进行了平均。下图展示了亲和力预测的误差从小到大排序，发现少数药物和少数蛋白质对总体误差的贡献率不成比例，并不是均匀的分布在药物或者蛋白质上，预测某些药物或者蛋白质的亲和力值比预测其他的蛋白质的亲和力要难。值得注意的是CHEMBL1779202(一种ALK抑制剂)、CHEMBL1765740(一种PDK1抑制剂)和蛋白质CSNK1E的MAE都大于2。

2.对药物误差分析

下面这幅图展示了每个药物的预测误差与潜在向量的前6个主成分的关系，我们可以看到难以预测的药物通常出现在原点附近。本文解释这个结果是因为具有独特分子结构的药物更容易被预测。

3.对蛋白质误差分析

分析同源蛋白对测试集性能的影响。首先，使用CLANS序列聚类算法对测试集中的目标蛋白进行聚类，得到5个主要聚类见下图。其次，计算每个聚类的平均绝对误差（见下）。第三，估计每个聚类在训练集中的代表情况。在这里，我们可以看到不同聚类之间误差有所不同。例如，第2个聚类的误差较高，而第4个聚类的误差较低。然而，我们发现训练集测试集平等地代表了这5个聚类，比例为5:1。这表明测试集的性能不能简单地解释为训练集中蛋白质家族的不对称表示。

模型解释和意义

1模型的可解释性对于研究具有重要意义，图神经网络可以在没有任何先验知识的情况下学习已知分子描述符的重要性，但是大多数学习到的潜在变量仍然无法通过已知的描述符解释。但是，模型的性能又表明这些学习的潜在变量在亲和力预测中是有用的，这也表明机器和人类专家看到的化学物质的方式既有相似又有不同。了解这种区别可能会进一步提高模型性能。

2.在测试集中存在一些“问题药物”或者“问题蛋白”它们的预测非常困难，这可以通过采集更多此类药物进行补充数据，或者利用化学知识来设计能够补充分子图的特征。另外，PCA的离群值是最容易预测的药物于，这表明需要一些额外的特征输入来区分缺乏明显分子亚结构的药物。

• “离群值”指的是在PCA降维后空间中距离其他数据点较远的数据点。原文指出这些离群点在预测中往往是最容易被预测的，这意味着当模型无法准确预测数据点时，提供更多的特征（比如离群点的附加特征）可能会有所帮助，以便更好地区分不同的数据点，进而提高预测准确率。

3.作者认为2D图像比1D字符串包含更多信息，但是模型仍然忽略了分子的立体化学信息，未来研究可以测试3D（或2D）中进行药物（或蛋白质）表征，是否会进一步提高模型性能。

4.此外作者还发现训练集中蛋白质的欠表示是”问题蛋白质“的主要原因，作者在下表中对同源蛋白对测试集性能的影响进行了分析，结果显示虽然测试集误差在聚类蛋白组中存在差异，但训练集能够平等地表示所有的蛋白质簇群，因此测试集性能的变化不能简单地解释为训练集中蛋白质组数据的不对称。

总结和未来工作

1.本文提出GraphDTA模型，性能表现优于其他模型

2.进行了潜在变量与已知分子描述符的冗余分析，侧面反映出图神经网络可以在没有先验知识的情况下学习到比较重要的分子描述符

3.分析了模型的性能发现少数药物、蛋白质对预测误差的贡献不成比例，反映出数据中未知结构对模型的影响

4.未来将蛋白质也表示为图形，例如3D结构的图形，可能会进一步提高性能，但也蛋白质3D结构难以获得、同时也会带来过拟合的风险。不过加入结构信息、结合位点、结合构象、药物分子在溶液中的环境信息（在化学反应中，溶液中的各种条件，如温度、压力、pH值等）来增强模型可能是有意义的。