一种基于网络推理的药物通路预测方法及其应用

1.本发明涉及计算机辅助药物设计及药物信息学领域，具体涉及一种基于网络推理的药物通路预测方法及其应用。


背景技术：

2.传统的新药研发方法一般遵循着“一个药物，对应一个靶标，对应一种疾病”的原则，主要包括基于靶标的药物设计、基于结构的药物设计和基于性质的药物设计等。尽管这些方法设计出了一些表现较好的候选药物，但在临床应用过程中，往往会产生脱靶效应，存在一定的局限性。近几年，越来越多的研究表明许多药物是通过调节多个靶标或疾病相关的信号通路来发挥治疗作用的，而不是仅通过单个靶标。例如：索拉菲尼通过抑制raf/mek/erk信号通路和受体酪氨酸激酶的活性来抑制肿瘤的生长和血管生成；共同抑制pi3k-akt信号通路和her激酶的活性有助于提高肿瘤患者的治疗指数。因此在传统药物设计方法的基础上需要发展新的药物设计方法，以提高药物的治疗效应和安全性、降低新药研发过程中较高的成本并提高新药研发的成功率。
3.为加快新药研发的效率，基于通路的药物发现(或基于通路的药物重定位)策略为研究人员开辟了新的视野。该策略的主要研究目的是发现新的药物与通路的关系。在疾病发生发展过程中，其表型主要体现在细胞内的分子通路水平，如基因表达水平等。药物可以通过改变疾病相关通路中的基因表达水平发挥治疗作用。这种药物通过作用于疾病相关通路来发挥治疗作用的关系，即药物-通路关系。然而，通过实验的方法来确定这种关系是非常耗时、耗力的。随着多组学技术的发展，药物基因组学可以呈现出在不同药物处理下，细胞内通路水平发生的全景式基因的变化，为从计算层次研究药物与通路的关系提供了夯实的数据基础。
4.在之前的研究中，基于矩阵分解、关联规则、机器学习等方法应用到预测药物与通路的关系。尽管这些方法已证明具有一定的可靠性，但随着药物和通路有关数据数量和类型的快速增长，远远超出了当前数据的表示和算法分析的能力。此外，目前的方法不能广泛地为一些老药、临床失败的药物和新的化学实体小分子预测潜在的通路关系。近几年来，由本课题组所研发的一系列基于网络推理的算法(nbi，sdtnbi，bsdtnbi)在药物研发领域中受到广泛关注，并应用到多重网络关系的预测，包括药物-靶标关系预测、药物-不良药物事件关系预测、药物-microrna关系预测、以及药物-atc编码关系的预测，并且该系列算法还具有不依赖阴性样本的显著优势。因此，有必要设计一种基于网络推理的药物通路预测方法来为一些老药、临床失败的药物和新的化学实体小分子预测潜在的通路关系，并使用该方法进行药物发现的重定位研究。


技术实现要素：

5.本发明针对现有技术的不足，提出一种基于网络推理的药物通路预测方法及其应用，该方法较传统基于矩阵分解、关联规则、机器学习等方法的药物通路预测而言，可以为
新的化学实体小分子预测潜在的通路关系，并且不依赖于阴性样本；本发明在使用基于网络推理算法的基础上，能够充分利用子结构-药物-通路异构网络中的信息，提高模型预测性能，简单有效，易于实施。
6.本发明的第一方面，提供了一种基于网络推理的药物通路预测方法，包括以下步骤：
7.一、构建子结构-药物-通路异构网络：利用公开的药物基因组学数据，计算药物诱导的基因标签，使用通路富集分析，构建已知的药物-通路网络；通过计算药物-通路网络中药物的各种化学子结构信息，构建子结构-药物网络；最后整合药物-通路网络和子结构-药物网络，构建子结构-药物-通路异构网络；
8.二、扩展基于网络推理算法到该异构网络：根据已知的药物-通路网络和子结构-药物网络，对于任意一个药物，其与之连接的通路节点和子结构节点各自分配一个单位的初始资源，构建基于网络推理算法的初始资源矩阵；然后在每一步资源扩散过程中，网络中拥有初始资源的子结构节点和通路节点，都会把节点本身的资源平均分配给与之相连接的邻居节点，再根据资源扩散的次数构建基于网络推理算法的转移矩阵；
9.三、预测新的药物-通路关系：对于给定的药物或新的化学实体小分子，依据所构建的初始资源矩阵和扩散后的资源转移矩阵，对于网络中的任意一个通路节点而言，其拥有资源的多少表明药物与通路之间存在关联性的强弱，即通路节点拥有的资源越多，其分数越高，药物与该通路之间存在关联性的可能性越大。
10.以下对本发明方法进行详细说明。
11.一、构建子结构-药物-通路异构网络
12.1.1利用公开的药物基因组学数据，计算药物诱导的基因标签，通过使用通路富集分析，构建已知的药物-通路网络；在该网络中，药物集合通路网络；在该网络中，药物集合表示有n种药物，通路集合表示有n种通路；这样，药物-通路网络的矩阵表示形式为：
[0013][0014]
其中，i∈(0,nd],j∈(0,n
p
]是两个正整数。
[0015]
1.2基于药物的化学结构信息，计算药物的子结构，构建子结构-药物网络；在该网络中，药物集合表示有n种药物，子结构集合表示有n种药物，子结构集合表示有n种子结构；这样，子结构-药物网络的矩阵表示形式为：
[0016][0017]
其中，i∈(0,nd],j∈(0,ns]是两个正整数。
[0018]
1.3整合上述两个网络，构建子结构-药物-通路网络；这样，子结构-药物-通路网络的矩阵表示形式为：
[0019][0020]
二、扩展基于网络推理算法到该异构网络
[0021]
2.1构建初始资源矩阵：根据基于网络推理的算法，我们在资源扩散的过程中引入了一个可调节参数α∈[0,1)，用于调节不同节点类型的初始资源分配。在初始资源分配过程中，我们对子结构-药物-通路网络中的任意一个药物节点的所有邻居节点分配一个单位的初始资源。其中，各个子结构节点赋予平均分享总量为α的初始资源，各个通路节点则赋予平均分享总量为1-α的初始资源。在资源扩散过程中，通过调节参数α，可以分析出在初始资源分配环节中占据主导地位的节点类型，有助于提升预测模型的表现能力。这样，药物-通路网络和子结构-药物网络的矩阵表示形式分别为：
[0022][0023]
其中，i∈(0,nd],j∈(0,n
p
]是两个正整数。
[0024][0025]
其中，i∈(0,nd],j∈(0,ns]是两个正整数。
[0026]
这样，子结构-药物-通路异构网络的初始资源矩阵表示形式为：
[0027][0028]
2.2构建资源转移矩阵：在构建资源转移矩阵之前，我们需要重新定义在资源扩散过程中子结构-药物网络和药物-通路网络的矩阵表示形式，仅仅是为了排除那些没有已知通路节点的全新化学实体小分子的影响，具体如下：
[0029][0030]
其中，i∈(0,nd],j∈(0,n
p
]是两个正整数。
[0031][0032]
其中，i∈(0,nd],j∈(0,ns]是两个正整数。
[0033]
根据基于网络推理的算法，我们在资源扩散的过程中还引入了另外两个可调节参数β∈[0,1)和γ∈(-∞, ∞)。参数β用于调节不同边的类型在资源扩散过程中对模型造成的影响，参数γ用于调节在资源扩散过程中枢纽节点对模型造成的影响。这样，每一步资源扩散过程的转移矩阵表示形式为：
[0034][0035][0036][0037]
其中，i,j∈(0,nd ns n
p
]是两个正整数。
[0038]
三、实现在该异构网络上的推理，预测新的药物-通路关系
[0039]
基于所构建的子结构-药物-通路异构网络的初始资源矩阵a
′
，以及每一步资源扩散过程中的转移矩阵w，迭代地在该异构网络中进行推理；假设资源扩散的次数为k，那么最终的资源扩散转移矩阵的计算公式如下：
[0040]
f＝a'
×
wk#(12)
[0041]
在该矩阵中，f(i,nd ns j)的值即是基于网络推理的算法预测药物-通路相互作用的打分，其中，i∈(0,nd],j∈(0,n
p
]是两个正整数。
[0042]
关于上述三个调节参数以及循环次数的数值确定，本发明使用不同子结构信息类型，在不同的参数α,β,γ下，在不同的扩散次数k下，进行训练模型。首先，在扩散次数k＝2和γ＝0下，寻找最优参数α和β；其次，在寻找出最优参数α和β后，进一步寻找最优参数γ和最佳扩散次数k；最终对训练后的最优模型进行评价和预测应用。通过优化，α为0.1，β为0.1，γ为-0.8，k为2。
[0043]
本发明的第二方面，提供了一种基于上述预测方法构建的药物通路预测模型，至少包括输入显示模块、药物通路关系计算模块以及存储模块。其中，输入显示模块用于输入疾病-通路网络、与疾病相关的化学实体以及待分析药物组合物，并显示分析结果；药物通路关系计算模块依据药物通路预测方法对药物与通路之间存在关联性进行计算，筛选得到潜在通路，并基于该通路对待分析药物组合物与通路之间的相互作用进行打分；存储模块实时记录药物通路关系计算模块的结果并进行存储。
[0044]
优选的，所述潜在通路为与个数最多的疾病相关的化学实体存在关系的通路。
[0045]
与现有技术相比，本发明的有益效果是：
[0046]
(1)本发明提供了一种基于网络推理的药物通路预测方法，第一次将基于网络推理的算法应用于药物通路预测领域。通过在标准数据集、外部验证集上测试表明，本发明在药物通路预测方面及基于通路的药物重定位方面具有较好的预测性能。
[0047]
(2)现有的药物通路预测方法，大部分依赖于常规的机器学习算法，需要在构建模型前，建立假定阴性样本数据集。而本发明是基于网络推理算法构建的药物通路预测框架，具有不依赖于阴性样本的显著优势。
[0048]
(3)现有的药物通路预测方法，无法为全新的化学实体小分子预测潜在的药物通路关系。而本发明通过整合子结构-药物网络和药物-通路网络，构建子结构-药物-通路异构网络，首次将子结构作为“桥梁”，为网络外的药物或者全新的化学实体小分子预测潜在
的药物通路关系，能够更加充分地帮助研究人员进行大规模的药物通路预测，为基于通路的药物发现或重定位研究提供了有效的工具。
附图说明
[0049]
为了更清楚地说明本发明现有技术中的技术方案或实施例，下面对本发明技术方案和实施例中所需要使用的附图作简单地介绍。
[0050]
图1为本发明流程图。
[0051]
图2对于预测已知通路的候选药物，基于十折交叉验证评价本发明在不同参数α和β下，不同分子指纹对应的模型表现能力。
[0052]
图3对于预测已知通路的候选药物，基于十折交叉验证评价本发明在参数γ下，不同分子指纹对应的模型表现能力。
[0053]
图4对于预测已知通路的候选药物，基于十折交叉验证评价本发明在不同的资源扩散次数k下，不同分子指纹对应的模型表现能力。
[0054]
图5为本发明实施例中肿瘤相关通路的网络图。不同颜色的矩形节点表示不同类型的肿瘤；不同大小的圆形节点表示不同通路的度值。
[0055]
图6为本发明实施例中肿瘤相关通路pi3k-akt信号通路作为基于通路的抗肿瘤药物重定位的重要特征。ln(ic
50
)是ic
50
以自然对数e为底的对数值。
具体实施方式
[0056]
以下结合具体优选的实施例用于说明本发明，但不用来限制本发明的范围。在不背离本发明精神和实质的情况下，对本发明方法、步骤和条件所作的修改或替换，均属于本发明的范围。若未特别指出，实施例中所用的技术手段为本领域技术人员所熟知的常规手段。
[0057]
基于网络推理的药物通路预测方法建立模型，其构建流程如图1所示，具体实现过程如下：
[0058]
一、构建子结构-药物-通路异构网络
[0059]
本方法所应用的数据集，包括通路集合、药物集合、子结构集合、药物-通路关系数据、子结构-药物关联数据、药物基因组学数据集。
[0060]
首先，构建药物-通路网络：使用公开的lincs药物基因组学数据集和cmap药物基因组学数据集构建药物诱导的基因标签，将lincs药物基因组学数据集按照8:2划分为pan-cancer数据集和外部验证集1；将cmap药物基因组学数据集作为外部验证集2。
[0061]
使用r包clusterprofiler对每个药物诱导的基因标签进行通路富集分析，得到药物通路关系数据，最终构建药物-通路网络(如表1所示)。对构建的药物通路网络进行筛选，仅保留符合如下标准的药物集合、通路集合及其关系条目：
①
200≤分子量≤800；
②
碳原子数》3；
③
细胞给药时间为6h，给药浓度为10μm；
④
同一药物在不同细胞上诱导的基因标签作并集处理；
⑤
除kegg pathway数据库中的疾病信号通路外，选择208条与代谢、环境信息处理、细胞过程、生物系统等相关的通路作为通路集合；
⑥
使用药物诱导的基因标签进行通路富集分析时，将校正后的p《0.01视为药物通路关系。
[0062]
表1.不同数据集中的药物-通路网络概况
[0063][0064]
注：表格中，n
drug
表示网络中药物的数量；n
pathway
表示网络中通路的数量；n
dpi
表示网络中药物通路作用关系的数量；稀疏度等于n
dpi
/(n
drug
×npathway
)。
[0065]
其次，构建子结构-药物网络：本发明使用分子指纹(fp)作为化学子结构信息的表征。本发明涉及的所有药物和小分子的化学结构转换为经典的smiles，并使用pipeline pilot软件进行标准化。使用padel-descriptor软件生成四种类型的分子指纹(fp)，包括klekota-roth fp，maccs fp，pubchem fp和substructure fp。使用python包rdkit生成morgan fp，并采用不同的原子半径(radius＝1,2)和长度(1024和2048字节)进行生成。这样，共有五种类型的分子指纹表征化学子结构信息，用于构建子结构-药物网络(如表2所示)。
[0066]
表2.不同数据集中的子结构-药物网络概况
[0067][0068]
注：表格中，n
drug
表示网络中药物的数量；ns表示网络中子结构的数量；n
dsi
表示网络中子结构药物关联的数量；稀疏度等于n
dsi
/(n
drug
×nfp
)。
[0069]
最后，构建子结构-药物-通路异构网络：本发明通过整合上述子结构-药物网络和药物-通路网络，构建用于训练和评价模型的子结构-药物-通路异构网络。
[0070]
二、扩展基于网络推理算法到该异构网络
[0071]
本发明使用基于网络推理算法构建药物通路预测方法。该方法使用不同子结构信息类型，在不同的参数α,β,γ下，在不同的扩散次数k下，进行训练模型。首先，在扩散次数k＝2和γ＝0下，寻找最优参数α和β；其次，在寻找出最优参数α和β后，进一步寻找最优参数γ和最佳扩散次数k；最终对训练后的最优模型进行评价和预测应用。
[0072]
三、实现在该异构网络中的推理，预测新的药物通路关系
[0073]
对于给定的药物drug，基于所构建的子结构-药物-通路异构网络和训练后的最优模型，进行预测潜在的通路关系。在最终的资源扩散转移矩阵中，f(i,nd ns j)的值即是基于网络推理的算法预测药物-通路相互作用的打分，其中，i∈(0,nd],j∈(0,n
p
]是两个正整数。其f打分值越高，排名越靠前，表示该药物与通路存在作用关系的可能性越大。本发明根据f为给定药物drug预测新的药物通路关系。
[0074]
构建后的模型在计算机等载体上运行，至少包括输入显示模块、药物通路关系计算模块以及存储模块。其中，输入显示模块用于输入疾病-通路网络、与疾病相关的化学实体以及待分析药物组合物，并显示分析结果；药物通路关系计算模块依据药物通路预测方法对药物与通路之间存在关联性进行计算，筛选得到潜在通路，并基于该通路对待分析药物组合物与通路之间的相互作用进行打分；存储模块实时记录药物通路关系计算模块的结果并进行存储。
[0075]
四、实验验证
[0076]
4.1评价指标
[0077]
本发明使用受试者工作特性曲线下面积(areas under a receiver operating characteristic curve，简称auc)作为评价指标。对于给定药物drug，根据预测每一个通路得到的f值，按照降序排列。给定一个阈值l，将排名在前l位的预测的药物-通路关系视为阳性，之后的视为阴性。将给定药物drug预测的阳性或阴性的药物-通路关系，与drug已知的药物-通路作用关系进行比较，可以计算得到drug真阳性数目(tp)，假阳性数目(fp)，真阴性数目(tn)，假阴性数目(fn)，进一步可以计算真阳性率(true positive rate，简称tpr)和假阳性率(false positive rate，简称fpr)。这样，通过给定一系列不同的阈值l，就可以得到一系列tpr和fpr。将tpr对fpr作图，即可得到受试者工作特性曲线(receiver operating characteristic curve，简称roc)，auc值即是该曲线下面积。除了auc值，在基于网络推理算法的预测结果中，排名越靠前的药物-通路关系在实际应用中越重要。
[0078]
本发明采用十折交叉验证的方法对模型进行训练和评价。这是一种被广泛应用于基于网络方法研究中的模型评价方式。在每一次的十折交叉验证过程中，首先会将子结构-药物-通路异构网络中的药物-通路的关系(边)随机划分为十等份；然后依次将其中的一份(10％)作为测试集，剩余的九份(90％)作为训练集。这样，即可得到十对不同的“训练集-测试集”。对每一对的“训练集-测试集”都可以计算一组评价指标，为了降低随机因素的影响，十折交叉验证重复进行了10次。使用“平均值
±
标准偏差”对评价指标进行表示。
[0079]
4.2模型评价
[0080]
本发明使用pan-cancer数据集和不同分子指纹类型所构建的子结构-药物网络进行训练，使用外部验证集1和2对训练好的最优模型进行评价和验证。如图2所示，在资源扩散次数k＝2和参数γ＝0时，对pan-cancer数据集中的药物-通路网络和不同分子指纹类型表征的子结构-药物网络进行训练，结果表明在参数α＝0.1和β＝0.1各个模型的表现能力较为突出。这表明给予模型较小的参数α值，意味着在初始资源分配环节中，每个药物节点与之连接的通路节点分配更多的初始资源，而不是子结构节点；给予模型较小的参数β值，意味着药物-通路的边设置更大的权重值，而不是子结构-药物关联的边。
[0081]
在确定模型参数α＝0.1和β＝0.1后，对模型参数γ进行优化。如图3所示，在资源
扩散次数k＝2和参数α＝0.1，β＝0.1时，对pan-cancer数据集中的药物-通路网络和不同分子指纹类型表征的子结构-药物网络进行训练，结果表明在参数γ＝-0.8大多数模型的表现能力较为突出。其中，分子指纹类型为morgan(1,1024)对应的模型表现最优(auc＝0.9358
±
0.0015)。这表明给予模型较小的参数γ值，意味着适当削弱子结构-药物-通路异构网络中刷枢纽节点的影响，对模型的表现能力起到促进作用。
[0082]
最后，在确定模型参数α＝0.1，β＝0.1和γ＝-0.8后，对模型资源扩散次数k进行评价。如图4所示，随着子结构-药物-通路异构网络中资源扩散次数的增加，各个模型的表现能力逐渐下降，并且这一结果与分子指纹表征的子结构类型无关。因此，资源扩散次数k＝2是一个合适的参数取值。根据上述结果，选择分子指纹类型为morgan(1,1024)表征子结构，资源扩散次数k＝2，参数α＝0.1，β＝0.1和γ＝-0.8时为最优模型。
[0083]
接下来，使用两个外部验证集对模型进行评价(如表3所示)，结果表明两个外部验证集在该模型上都取得较好的表现能力，auc分别为0.8519和0.7494。因此，通过在外部验证集上测试表明，本发明在药物通路预测方面具有较好的预测性能，能够更加充分地帮助研究人员进行大规模的药物通路预测，为基于通路的药物发现或重定位提供了有效的工具。
[0084]
表3.外部验证集的模型评价
[0085][0086]
五、案例分析：基于通路的药物重定位研究
[0087]
根据上述我们使用pan-cancer数据集中的药物-通路网络及其子结构-药物网络所构建的最优模型，对基于通路的药物重定位展开研究。具体地，发现具有潜在抗肿瘤效应的非肿瘤药物，主要分为以下三个步骤：
[0088]
5.1分析肿瘤相关的信号通路：我们从tcga数据库中下载了8,628个肿瘤样本的转录组学数据，包括了21种肿瘤类型。通过分别使用r包deseq2和clusterprofiler进行差异基因表达分析和通路富集分析，最终构建了肿瘤-通路网络。如图5所示，该网络包括了21种肿瘤类型、110条通路及其289条关系数据。比如：乳腺癌(brca)主要与pi3k-akt通路(hsa04151)、mapk通路(hsa04010)等通路密切相关。使用cytoscape软件中的networkanalyzer插件对该网络进行分析发现，一些通路在该网络中具有很高的度值(degree)，如：focal adhesion通路(hsa04510，degree＝14)、cell cycle通路(hsa04110，degree＝12)、pi3k-akt通路(hsa04151，degree＝10)等。因此，这些肿瘤密切相关的信号通路为基于通路的抗肿瘤药物发现或重定位研究提供了依据。
[0089]
5.2结合上述肿瘤-通路网络的分析，根据我们所构建的最优模型，确定抗肿瘤药物的通路特征。我们从文献中收集到46个全新化学实体小分子在a549上的抗肿瘤活性数据(ic
50
)，使用我们所构建的最优模型预测这些小分子的潜在通路关系(预测长度为l＝6)，通过秩和检验(wilcoxon rank sum test)的方法对预测出来的通路进行分析。根据预测的小分子与通路的关系，备择假设定义为：与通路a连接的小分子，其抗肿瘤活性数据(ic
50
)小于未与通路a连接的小分子的抗肿瘤活性数据(ic
50
)，那么通路a被认为是抗肿瘤药物的关键
通路特征。如图6所示，我们可以发现与pi3k-akt通路相连的小分子的抗肿瘤活性数据(ic
50
)要小于那些未与pi3k-akt通路相连的小分子(wilcoxon test:pvalue＝0.0013)。并且这些小分子的抗肿瘤活性数据(ic
50
)与pi3k-akt通路之间的皮尔森相关系数为r2＝-0.6489(pvalue＝1.08e-06)，而与小分子相连的其它通路并没有呈现出显著性差异。此外，通过对ccle数据库中的24个抗肿瘤药物进行预测分析，发现其中19个抗肿瘤药物与pi3k-akt通路存在潜在的关系。因此，在基于通路的抗肿瘤药物重定位策略中，我们优先考虑将pi3k-akt通路作为抗肿瘤药物的通路特征，并根据该特征发现上市药物中非肿瘤药物的潜在抗肿瘤效应。
[0090]
5.3基于pi3k-akt通路，发现非肿瘤药物中的潜在的抗肿瘤药物。使用本发明开发的药物通路预测方法，通过对drugbank数据库中上市的抗肿瘤药物进行预测分析(预测长度为l＝6)，近80％的抗肿瘤药物都预测与pi3k-akt通路具有相互作用关系；通过对drugbank数据库中上市的非肿瘤药物进行预测分析，与pi3k-akt通路具有相互作用关系的部分非肿瘤药物呈现出一定的抗肿瘤效应(如表4所示)，如：降血脂药物洛伐他汀。同时，临床前和临床研究也证明了普伐他汀、辛伐他汀可以明显抑制肿瘤细胞的生长。此外，螺内酯、西格列汀、阿卡波糖、地塞米松、己酮可可碱等非肿瘤药物也已正在相关研究中证实具有一定的抗肿瘤活性。因此，案例分析结果表明本发明的预测结果将对生物学实验具有一定的指导作用，为抗肿瘤药物的重定位或基于网络的虚拟筛选提供了一个强有力的工具。
[0091]
表4.基于pi3k-akt通路的非肿瘤药物的抗肿瘤效应
[0092]
[0093][0094]
注：抗肿瘤活性数据(ic
50
)来源于prism drug repurposing source。