用于药物反应预测的自适应大小的子结构学习

子结构学习的引言

药物化学的理论表明，药物是一个简单的实体，其由不同的官能团（或者化学子结构）组成，这些子结构决定了药代动力学特性（生物体如何处理它）和药效学特性（它如何影响生物体），以及最终的交互反应。为此，研究人员提出了基于子结构学习的DDI预测方法。

• 首先，药物在分子图表示的基础上被表示成一个graph，为了从药物的graph中提取子结构，作者提出一个消息传递网络，其中，边具有可学习的权重，权重值约束在0到1之间。这些权重可以被认为是划分子结构的门限（gates），具有产生灵活尺寸和不规则形状子结构的效果。
• 其次，两种药物之间的 DDI 预测是基于它们学习到的子结构之间的相互作用分数，每个相互作用分数都使用共同注意机制（co-attention）由可学习的权重加权得到。对于副产物，作者提出的方法可以提示哪些子结构可能是导致 DDI 发生的原因。
• 最后，作者在两个数据集上评估了方法：1.Drugbank数据集，包含药物1704种，191808个DDI对（DDI pairs），86种反应类型；2.Twosides数据集，645种药物，4649441个DDI对，1317种反应类型。
  实验在两种设置下进行，一种是transductive，DDI对的训练和测试集共享相同的药物；另一种是inductive（也称为冷启动），测试集和训练集中的药物是不同的（测试集里的药物不在训练集中出现，训练集中的药物也不在测试集中出现）。Inductive更困难更具挑战性。
  两个数据集都被使用到transductive设置中，对于inductive设置，只使用Drugbank进行实验。

相关工作的回顾

相关工作主要是回顾以往的DDI预测任务的两个方面：1.药物是如何表示的；2.实际预测是如何进行的。

药物的表示

过去大多数DDI预测方法都是利用分子指纹表示药物，或者是其他的药物概况，比如副作用，结合靶标，转运蛋白，酶，通路或这些特征中的两个或多个的组合。分子指纹是二进制向量，其元素表示特定化学子结构的存在 (1) 或不存在 (0)。药物概况类似地表示为二进制向量，表示特定概况的存在或不存在，例如，特定的副作用或结合靶标。

一些方法通过将药物表示为相似性向量来执行进一步的预处理，该向量使用相似性度量（例如余弦相似性，Jaccard similiart）表示它与同一表示空间下的其他药物的相似程度。这是由相似/不同的药物可能具有相似/不同的生物活性的假设指导的。这种表达的缺点是它们是人工设计的，受限于专家的知识水平，缺乏可扩展性，不能编码新的药物。

近年来，GNN被用于化学分子的可学习表示，并提高了分子相关任务的性能。然而，在不过多考虑DDI主要因素（官能团或化学子结构）情况下，这种方法通常用于表示一个整体的药物。最近的一些方法在DDI中考虑了子结构信息。然而，节点在每个GNN层的隐层表达被离散地视为药物的子结构。这种方法产生的子结构形状规则，其大小由GNN的层数以及感受野决定（比如SSI-DDI）。

DDI预测

预测DDI的方法可以分为两类：一类是药物来自于知识图谱，另一类是药物被认为是相互独立的。在后者，通常我们将两个药物的表达进行聚合（求和或者拼接），然后使用分类器进行预测。对于network based方法（第一类），假设药物形成一个相互连接的系统，其中药物是节点，边表示节点之间的 DDI 或药物之间的相似性。我们通常将不同的图处理方法包括标签传播、矩阵分解、图自编码器应用于这些药物知识图以进行预测或首先学习药物的低维表示，然后执行 DDI 预测。

基于药物知识图的方法可以考虑药物互联系统中的拓扑信息，但是这种方法不能用于inductive设置。作者认为应该把药物考虑为独立的实体才能用于inductive设置下的实验。

方法

方法的整体框架如下图：

给定一个DDI的元组$(d_x,r,d_y)$，$d_x,d_y$分别用分子图$G_x=(V_x,E_x),G_y=(V_y,E_y)$表示，$r$与可学习的矩阵$M_r\in R^{b\times b}$建立联系，在子结构被提取后，$G_x$生成$|V_x|$个子结构，节点代表每个子结构的中心。$G_y$遵循同样的过程。一个交叉的子结构注意力${\gamma }_{ij}$被计算，用于表示子结构之间的相关性。每个子结构都经过线性变换，最终的DDI预测为每对药物$d_x,d_y$的子结构之间的相互作用得分之和，每对子结构的相互作用的权重为注意力${\gamma }_{ij}$。

问题描述

给一个药物集合$D$，一个反应类型集合$R$，DDI预测任务被看作找一个函数$f$使得：$$D\times R\times D\Rightarrow [0,1]$$换言之，给一个包含两个药物和一个具体反应类型的三元组，我们要预测这个DDI事件发生的概率。我们需要找到函数$f$使得预测结果与数据集 M = { ( d x , r , d y ) i } i = 1 N ⊂ D × R × D M=\left\{(d_{x},r,d_{y})_{i}\right\}_{i=1}^{N}\subset D\times R\times D M={(dx​,r,dy​)i​}i=1N​⊂D×R×D近似。

药物被表示为去氢的无向图$G=(V,E)$，其中$V$为节点集合，表示原子，$E\subset V\times V$为边集合，表示原子之间的键。每个节点$v_i$都有一个对应的特征${\text{x}}_i\in R^d$，每个边$e_{ij}=(v_i,v_j)$也对应一个向量${\text{x}}_{ij}\in R^{d^{\prime }}$，由于初始阶段，边是无向的，所以$e_{ij}=e_{ji},{\text{x}}_{ij}={\text{x}}_{ji}$。

门限消息传递网络提取子结构

给一个常规的图，以节点$v_i$为中心的子结构的表达可以看作是所有可到达$v_i$的节点群的聚合，为了精细表示子结构，我们将图中的边看作门限，即，值在[0,1]的权重。通往节点$v_i$的路径上的一个节点$v_j$的加权值是连接它到$v_i$的路径上的权值的乘积。

我们假设$v_j$到$v_i$有多条路径，这种情况出现在环状的子图上，我们应该单独考虑每条路径。这里提出的用边值的乘积来衡量节点，最终会产生不同大小和形状的子结构。如果在路径上，有一条边的权重为0，则包含该边的路径权重乘积将为0，代表截断，从而得到不同形状的子结构。

图中的每个节点都可以作为一个子结构的中心，以便从图中提取尽可能多的子结构。这需要将每条边转换成双向边（每个方向都要自己的权重），因为一个节点既可以是目标节点（中心）也可以是源节点。下图是子结构提取的示例。

以图c为例，为了提取以节点v5（红色）为中心的子结构，生成从end-nodes开始到v5结束的所有路径（绿色），例子中，v5有4个end-nodes（v1，v6，v9，v10），即4条路径，在图d中，v5的向量s5表示为四条路径与节点本身的信息聚合。

现在，我们需要把无向的分子图表示$G=(V,E)$转化为有向图，边$e_{ij}$和边$e_{ji}$成为两条分开的边。为了突出这个特点，我们重命名为$e_{i\to j}$和$e_{j\to i}$，每个有向边都分配一个可学习的权重，其值在[0,1]内。

对于一个图，上述的路径生成计算效率比较低，为此，作者提出门控消息传递网络G-MPNN模拟这个过程。MPNN是一种多层空间卷积GNN框架，每一层都包括三部分：消息传递，聚合，更新。$${\text{m}}_{j\to i}^{(l)}=M^{(l)}({\text{h}}_j^{(l-1)},{\text{h}}_i^{(l-1)},{\text{q}}_{j\to i}),\forall j:v_j\in N(v_i)$$$${\text{a}}_i^{(l)}=A^{(l)}({\left\{{\text{m}}_{j\to i}^{(l)}\right\}}_{j:v_j\in N(v_i)})$$$${\text{h}}_i^{(l)}=U^{(l)}({\text{h}}_i^{(l-1)},{\text{a}}_i^{(l)})$$其中，${\text{h}}_i^{(l)}$为节点$v_i$在第$l$层的特征，${\text{q}}_{j\to i}$是节点$v_j$到$v_i$的边特征。简而言之，在每次迭代中，节点会从相邻节点传递消息。这样做的效果是，在迭代GNN第$l$次时，一个节点将被更新的信息包括了在长度为$l$的跳数范围内可到达的所有节点特征。这与子结构提取中沿着路径中节点聚合信息的思想是符合的。图的遍历具有访问节点的冗余（一个节点可以在遍历中出现多次，但在一条路径中最多出现一次，参考上图d），为了更接近图d的路径生成过程，作者使用了directed MPNN，即D-MPNN。为了减少MPNN节点计算的冗余，消息在边之间传递而不是节点之间，D-MPNN包含三个部分： m k → j ( l ) = M ( l ) ( h k , h j , q k → j ( l − 1 ) ) , ∀ k : v k ∈ N ( v j ) \ { v i } \textbf{m}_{k\rightarrow j}^{(l)}=M^{(l)}(\textbf{h}_{k},\textbf{h}_{j},\textbf{q}_{k\rightarrow j}^{(l-1)}),\forall k:v_{k}\in N(v_{j})\backslash\left\{v_{i}\right\} mk→j(l)​=M(l)(hk​,hj​,qk→j(l−1)​),∀k:vk​∈N(vj​)\{vi​} a j → i ( l ) = A ( l ) ( { m k → j ( l ) } k : v k ∈ N ( v j ) \ { v i } ) \textbf{a}_{j\rightarrow i}^{(l)}=A^{(l)}(\left\{\textbf{m}_{k\rightarrow j}^{(l)}\right\}_{k:v_{k}\in N(v_{j})\backslash\left\{v_{i}\right\}}) aj→i(l)​=A(l)({mk→j(l)​}k:vk​∈N(vj​)\{vi​}​)$${\text{q}}_{j\to i}^{(l)}=U^{(l)}({\text{q}}_{j\to i}^{(l-1)},{\text{a}}_{j\to i}^{(l)})$$MPNN和D-MPNN的区别在于前者更新节点特征，后者更新边的特征，比如为了更新边$e_{j\to i}$，第一个公式从相邻边$e_{k\to j}$，其中 k : v k ∈ N ( v j ) \ { v i } k:v_{k}\in N(v_{j})\backslash\left\{v_{i}\right\} k:vk​∈N(vj​)\{vi​}，$v_j$被视为公共节点。在边特征更新后的最后一次迭代$L$之后，节点被表示为所有进入边特征的集合。比如节点$v_i$的最终特征${\text{s}}_i$为：$${\text{s}}_i=\sum _{j:v_j\in N(v_i)}{\text{q}}_{j\to i}^{(L)}$$现在，给出GMPNN的计算步骤，给一个图$G=(V,E)$表示药物：

• 对节点特征先进行非线性变换获得新的特征表达：$${\text{h}}_i=ML{P}_{init\_n}({\text{x}}_i),\forall v_i\in V$$其中，${\text{h}}_i\in R^f$是变换后的节点特征。
• 边的特征同样经过非线性变换处理：$${\text{h}}_{j\to i}=ML{P}_{init\_e}({\text{x}}_{ji}),\forall e_{j\to i}$$$${\text{h}}_{i\to j}=ML{P}_{init\_e}({\text{x}}_{ij}),\forall e_{i\to j}$$其中，${\text{h}}_{j\to i},{\text{h}}_{i\to j}\in R^m$分别是边$e_{j\to i},e_{i\to j}$的新特征，注意，尽管边$e_{j\to i},e_{i\to j}$是不同的边，但它们的特征目前是相同的，因为初始无向图限制了${\text{x}}_{ji}={\text{x}}_{ij}$，边$e_{j\to i}$的门限$w_{j\to i}\in [0,1]$是根据关联节点进行初始化的：$$o_{j\to i}=\frac{1}{c}({\text{h}}_{j\to i}^{T}ML{P}_w({\text{h}}_j||{\text{h}}_i)),w_{j\to i}=\sigma (o_{j\to i})$$其中，$||$表示拼接，$c$是一个常数。
• 现在，我们有节点特征${\text{h}}_i$和边的权重$w_{j\to i}$，我们应用D-MPNN机制模拟子结构的路径生成。消息在边之间传递而不是节点之间，边$e_{j\to i}$的特征表示为：$${\text{q}}_{j\to i}^{(0)}=w_{j\to i}{\text{h}}_j$$其中，${\text{q}}_{j\to i}^{(0)}\in R^f$，我们需要与前面提到的${\text{h}}_{j\to i}$进行区别，${\text{h}}_{j\to i}$只是用于计算边的权重，而${\text{q}}_{j\to i}^{(0)}$将在消息传递中使用，因为它包含了节点信息和边的权重信息。
• 应用D-MPNN机制： m k → j ( l ) = w j → i q k → j ( l − 1 ) , ∀ k : v k ∈ N ( v j ) \ { v i } \textbf{m}_{k\rightarrow j}^{(l)}=w_{j\rightarrow i}\textbf{q}_{k\rightarrow j}^{(l-1)},\forall k:v_{k}\in N(v_{j})\backslash\left\{v_{i}\right\} mk→j(l)​=wj→i​qk→j(l−1)​,∀k:vk​∈N(vj​)\{vi​} a j → i ( l ) = ∑ k : v k ∈ N ( v j ) \ { v i } m k → j ( l ) \textbf{a}_{j\rightarrow i}^{(l)}=\sum_{k:v_{k}\in N(v_{j})\backslash\left\{v_{i}\right\}}\textbf{m}_{k\rightarrow j}^{(l)} aj→i(l)​=k:vk​∈N(vj​)\{vi​}∑​mk→j(l)​$${\text{q}}_{j\to i}^{(l)}={\text{q}}_{j\to i}^{(l-1)}+{\text{a}}_{j\to i}^{(l)}$$注意到，边的特征每次都会根据权重进行缩放，这样的效果是实现了节点特征乘以连接它与子结构中心节点的边的权值的乘积。
• 最后，得到节点$v_i$为中心的子结构信息：$${\text{s}}_i=f_{sub}({\text{h}}_i+\sum _{j:v_j\in N(v_i)}{\text{q}}_{j\to i}^{(L)})$$其中，$f_{sub}$是一个非线性函数。

G-MPNN是可以用于inductive设置的，因为分子总是有原子和键的特征表示，即我们总是可以得到节点和边的特征，从而可以使用G-MPNN提取自适应大小的子结构。

DDI预测的过程

给一个DDI元组$(d_x,r,d_y)$，预测该事件的概率：$$P(d_x,r,d_y)=\sigma (\sum _i^{|V_x|}\sum _j^{|V_y|}{\gamma }_{ij}({\hat{\text{s}}}_i^{(x)}{)}^T{M}_r{\hat{\text{s}}}_j^{(y)})$$其中，${\hat{\text{s}}}_i^{(x)},{\hat{\text{s}}}_j^{(y)}$分别是药物$d_x$的子结构$i$，和$d_y$的子结构$j$经过线性变换后的表达：$${\hat{\text{s}}}_i^{(x)}=W{\text{s}}_i^{(x)},{\hat{\text{s}}}_j^{(y)}=W{\text{s}}_j^{(y)}$$其中，$W\in R^{b\times b}$是可学习的参数。$M_r\in R^{b\times b}$是反应类型$r$的可学习的参数，并且是对角矩阵：$M_r=diag({\text{m}}_r)$，${\text{m}}_r\in R^b$是反应类型$r$的可学习向量。

交叉的子结构注意力${\gamma }_{ij}$为：$${\gamma }_{ij}=softmax(ML{P}_{\gamma }({\text{s}}_i^{(x)}||{\text{s}}_j^{(y)}))$$现在，DDI预测被分解为DDI元组的二元预测。DDI数据集$M$中的样本都是正样本，我们还需要构造负样本，通过破坏元组$(d_x,r,d_y)$即可，比如将$d_x$或$d_y$替换成别的药物。

通过最小化BCE Loss优化目标，损失函数为：$$L=-\frac{1}{|M|}\sum _{i:(d_x,r,d_y{)}_i\in M}(log(p_i)+log(1-p_i^{\prime }))$$其中，$p_i,p_i^{\prime }$分别为正负样本的预测概率。

结果

直推式学习Transductive设定下的模型比较。下表加粗了每项指标的最优值。论文所提方法对每个指标的改进情况显示在最后一行。

下图可视化展示亚硝酸戊酯和西地那非这两种药物之间的DDI prediction。$p=0.82$是预测输出值。灰色填充的两个原子是所关心的子结构的中心点。（红色虚线框中的硝酸基团在这两个药物关联之间起了重要作用，根据文献，这两个药不能同时服用）

归纳式学习Inductive设定下的模型表现比较。下表加粗了每项指标最优值。

作者从药物集合$D$中划分药物，随机选取五分之一药物作为$D_{new}$，$D_{new}$不参与训练，即对模型来说是陌生的，$M_{s1}$表示数据集中的$d_x,d_y$都来自$D_{new}$，$M_{s2}$表示数据集中的$d_x$来自$D_{new}$，$d_y$来自$D_{old}=D-D_{new}$。