一种抗感染性疾病药物的启发式发现方法及其在治疗新型冠状病毒肺炎药物发现中的应用初探

高 敏，徐睿峰，全 源，梁峰吉，朱月星，熊江辉2，，

〔1.哈尔滨工业大学（深圳）计算机科学与技术学院，广东 深圳 518055；2.深圳市绿航星际太空科技研究院医学健康部表观遗传学实验室，广东 深圳 518117；3.中国航天员科研训练中心航天医学基础与应用国家重点实验室，北京 100094；4.北京银谷芳香科技有限公司智慧医疗研究中心，北京 100080〕

新型冠状病毒肺炎（corona virus disease 2019，COVID-19）的暴发使抗冠状病毒药物研发再次成为世界瞩目的焦点，亟需筛选有效的药物以抗新型冠状病毒（severe acute respiratory syndrome coronavirus-2，SARS-Cov-2）。回顾性研究发现，大量的重症急性呼吸综合征（severe acute respira⁃tory syndrome，SARS）和COVID-19患者使用经验性的临床治疗策略，如核苷类药物利巴韦林、激素类药物皮质类固醇和Ⅰ型干扰素等，但这些药物不良反应非常严重［1-4］。2003年，在中国大陆SARS暴发时，约50%患者除西药外还接受了中草药作为辅助疗法，中西医结合的方法对SARS患者产生了一些积极效果［5］。

自2003年以来，研究者们探索了大量天然产物成分对冠状病毒的抑制活性，为此，本研究对相关领域进行文献调研。尽管已认识到天然产物对冠状病毒治疗的重要性，但直接通过体外实验从数目庞大的天然产物成分中筛选出有效的抗病毒成分仍然是一个巨大的挑战。运用药物重定位（老药新用）技术，通过计算工具快速筛选出最具潜力的天然产物再进行实验验证，将会极大加快药物研发速度。

药物重定位被认为是最有前途的转化医学策略之一，不仅能够使现有药物的价值得到更好的发挥，还可以避免药物的吸收、分布、代谢、排泄和毒性性质缺陷带来的药物研发失败［6］。寻找现有药物的新用途通常依赖于文本挖掘［7］、化学遗传学［8］和网络分析［9］等方法。重大感染性疾病暴发期的特点是时间紧迫且数据有限。基于病毒靶标的药物发现路径需要精细了解病毒结构与功能，而且临床疗效和副作用验证需要大量时间。基于表型的药物重定位技术是一种重要的药物发现新思路［10］，但是在感染性疾病暴发早期，很难短期获得足量的基因组学和表型组学数据。

针对这种在数据极其有限的情况下需要快速产生可验证药物假说并进行量化评价的需求，本研究借鉴人工智能领域常用的启发式搜索策略，提出一种抗感染性疾病药物的启发式发现方法（agile discovery method of drugs or natural products for epidemics），暂命名为aCODE方法。本研究的核心假设是，如果少数几个药物在细胞学或临床实践中呈现一定的临床获益，则以这组药物作为“探针”，可以归纳这组药物靶向的生物学过程（基因模块），进而对大量候选化合物的潜在功效进行预测；进一步对候选化合物（天然产物和芳香分子）进行模式匹配与统计检验排序。

本研究以艾滋病、流感和细菌感染等疾病数据建立aCODE量化评价体系，随后应用于治疗COVID-19天然产物成分筛选。本方法可随临床实践的进展将初始药物列表进行迭代更新，从而优化系统的精准度和预测能力；输出的天然产物或中药、药食同源成分可实施快速测试，形成敏捷研发测试闭环。

1 数据资源与分析方法

1.1 获取药物-疾病网络数据

药物-疾病关系来自SCG-Drug（http：//zhan⁃glab.hzau.edu.cn/scgdrug）［11］。美国食品药品监督管理局（Food and Drug Administration，FDA）批准药物来自DrugBank数据库（https：//www.drug⁃bank.ca/）［12-13］。

1.2 获取基因模块数据与药物靶基因数据

人类基因模块数据源下载自Molecular Signa⁃tures Database（MSigDB）数据库（https：//www.gsea-msigdb.org/gsea/msigdb/）［14］。 MSigDB 数据库包含常用的GO Biological Process，GO Cellular Component和GO Molecular Function等基因模块。通过STITCH数据库（ChemicalAssociationNetworks，http：//stitch.embl.de）［15］获取美国FDA批准的药物相互作用基因（STITCHscore≥200）。

1.3 药物-基因模块相关性计算

理论上，如果某药物的靶基因显著富集于某基因模块，则可认为该药物靶向该基因模块。因此，通过如下方式计算药物-基因模块相关性：某个药物i，其相互作用基因集合与基因模块j作交集，并计算超几何分布的P值，取以10为底P值的负对数（-lgP）作为药物i与基因模块j的相关性指数，记为Rij。

超几何分布检验公式如下：

其中，x是药物i靶基因与基因模块j基因集合交集的基因数量，M是全基因组基因数量，K是药物i靶基因数量，N是基因模块j基因数量。

1.4 药物在基因模块维度空间的特征向量计算

每个药物选择性地靶向某些基因模块。因此，每个药物在全局所有基因模块维度上的靶向模式，将形成一个类似指纹图谱的特征，作为该药物作用于人体潜在影响的表征。为了进行不同药物在基因模块维度空间上的相似性搜索，计算其特征向量V：药物i在药物-基因模块相关性矩阵R中第i行所有列形成的行向量，即其特征向量，记为Vi。

对每种疾病，为了实现从少数已知药物扩展产生大量候选药物的假设，取若干初始药物（“种子药物”），并计算种子药物集合的特征向量：对M个种子药物，取其特征向量的中值（median），作为该组种子药物的特征向量Vm。假设，当M=3时，种子药物为药物1，药物2和药物3，种子药物集合的特征向量Vm=中值（V1，V2，V3），其中V1，V2和V3分别为药物1、药物2和药物3在药物-基因模块矩阵中的特征向量。

1.5 5种抗感染性疾病药物分析

选择美国FDA批准药物数≥40个的5种感染性疾病——艾滋病、流行性感冒、副黏病毒感染、细菌性感染和百日咳，对于每种疾病，在其阳性药物集合中随机抽取（500次）种子药物，种子数量从2取到20，观察预测药物与该疾病阳性药物的交集。抽取种子药物后，计算种子药物集合的特征向量Vm。枚举所有候选药物，设为药物j，取其特征向量Vj，并计算Vm与Vj的皮尔森（Pearson）相关性，得到相关系数的Rho值和P值。按Rho值降序排序，按Rho值取前100个为预测药物。观察预测药物与阳性药物的交集，计算预测药物与FDA批准的该疾病药物阳性集合的交集大小，通过超几何分布计算二者交集的显著性。

通过预测结果的阳性率和其超几何检验P值2个指标来评估该方法的可靠性，其中阳性率计算公式如下：

其中，K是前100个药物的个数，即K=100；N是美国FDA批准的该疾病阳性药物的个数。

阳性率P值计算公式如下：

其中，M是纳入计算的候选药物总数；K是前100个药物的个数，即K=100；N是FDA批准的该疾病阳性药物个数；x是前100个药物与FDA批准的该疾病阳性药物的交集个数。

重复上述抽取种子、计算阳性率及其P值的计算过程，输出500次计算得到的阳性率的平均值及其P值的中值。

同时，设置对照组以验证该方法的可信度。其中，阴性对照组A用所有FDA批准的适应证非当前疾病的抗感染类药物代替“当前疾病阳性药集合”和同样计算500次的阳性率及其P值中值；阴性对照组B用所有适应证为非感染类疾病药物，代替“当前疾病阳性药集合”，同样计算500次的阳性率及其P值中值。

1.6 aCODE方法在抗COVID-19药物发现应用中的探索

选取《新型冠状病毒感染的肺炎诊疗方案（试行第七版）》［16］中的4个抗病毒药物，即洛匹那韦（lopinavir），利巴韦林（ribavirin），利托那韦（ritonavir），磷酸氯喹（chloroquine）作为治疗COVID-19的种子药物。选取天然产物（中药）化学成分数据库（Traditional Chinese Medicine Integrated Data⁃base，TCMID，http：//www.megabionet.org/tcmid/）［17］的天然产物成分作为候选药物。使用aCODE方法进行药物预测，预测流程如下：①通过STITCH数据库获得种子药物和候选药物的靶基因；②对于每个种子药物，计算其特征向量Vi（药物i在药物-基因模块相关性矩阵R中第i行所有列形成的行向量，即其特征向量）。③计算种子药物集合的特征向量：取4个种子药物特征向量的中值作为该组种子药物的特征向量Vm，Vm=中值（V1，V2，V3，V4）。④ 枚举所有候选药物，设为药物j，取其特征向量Vj，并计算Vm与Vj的皮尔森相关性，得到相关系数的Rho值和P值。按Rho值降序排序，按Rho值取前100个为预测药物。

为了验证预测结果有效性，本研究通过文献调研的方式获取已知的具有抗冠状病毒活性的天然产物为验证集进行随机取样统计。文献调研使用谷歌学术（https：//f.glgoo.top/scholar）（2020年2月8日）进行：初始检索设定时间限制为“2003年及以后”，关键词以“coronavirus”和“natural product”进行检索，检索结果按“相关性”排序，选择排序靠前的文献［1-3，5，18-26］精读记录，并以这些参考文献及引用文献中所涉及的天然产物为线索，整理一个具有抗冠状病毒活性的天然产物的初步成分。随后依次以“coronavirus”和相应成分的英文名称为关键词进行检索，从而对初步成分进行核减和扩充，扩充的新成分再次遍历，直至无新扩充天然产物成分为止。在文献检索过程中，凡综述中所述成分，均溯源其原始文献核实，核实不到的均予以核减。全部天然产物成分按其参考文献引用次数（Cita⁃tions）排序，有多篇文献的采用多篇引用次数的总和（Total）进行排序。

2 结果

2.1 基于5种抗感染性疾病药物分析的aCODE方法量化评价

抗感染性疾病药物除了靶向病毒或细菌，同时也可能靶向病原微生物与宿主的相互作用通路。为了确认是否可以通过宿主靶基因网络这个维度空间来表征抗感染性疾病药物的特征，本研究首先分析了所有美国FDA批准的抗感染类药物的特征向量。对某一疾病，将美国FDA批准的所有药物分为2类，第一类为有该疾病的适应证（阳性药），第二类为无该疾病适应证（阴性药）。将这2类药物的基因模块特征向量进行t检验，发现所有抗感染性疾病的阳性药和阴性药均存在显著差异的基因模块（表S1，http://202.38.153.236:81/Jweb_cjpt/CN/item/down⁃loadFile.jsp?filedisplay=20200811174457.xlsx）。以HIV为例，其阳性药物富集的基因模块（表S2，http://202.38.153.236:81/Jweb_cjpt/CN/item/download⁃File.jsp?filedisplay=20200811174457.xlsx）中 包 含“艾滋病病毒生命周期”（HIV life cycle）、“病毒生命周期”（viral life cycle）等重要的基因模块。差异最显著的基因模块是“辅因子代谢过程”（cofactor metabolic process），是抗逆转录病毒疗法中被认为非常有吸引力的新靶标［27］。以上结果表明，对于感染性疾病而言，其有疗效的药物对于宿主基因网络的扰动存在共性特征，且与其药物作用机制密切相关。确认上述理论基础后，本研究量化分析了通过输入种子药物、计算候选药物与种子药物特征向量相关性，进而筛选药物的可行性。选择了美国FDA批准药物≥40个的5种感染性疾病——艾滋病、流行性感冒、副黏病毒感染、细菌性感染和百日咳进行分析。对于每种疾病的实验组、阴性对照组A和阴性对照组B，随机抽取（500次）M个种子药物，M从2取到20，观察预测药物与阳性药物的交集，计算预测药物与FDA批准的该疾病药物阳性集合的交集大小，通过超几何分布计算二者交集的显著性。重复进行上述计算500次后，计算预测结果中阳性药物比例（阳性率）的平均值如图1、表1、表2和表S3-S5（http://202.38.153.236:81/Jweb_cjpt/CN/item/downloadFile.jsp?filedisplay=20200811174457.xlsx）所示，反映交集显著性的P值的中值如表1、表2和表 S3-S5（http://202.38.153.236:81/Jweb_cjpt/CN/item/downloadFile.jsp?filedisplay=20200811174457.xlsx）所示。

实验结果显示，当种子药物从2取到20时，实验组的预测结果阳性率始终高于阴性对照组A和阴性对照组B。以“HIV infections”为例，如表1所示，当输入2个阳性药物作为种子药物时，实验组阳性率为9.6%且P值显著（P=2.40E-04）；而阴性对照组A和阴性对照组B的阳性率分别为4.29%（P=0.32）和3.9%（P=0.32），P值均不显著。说明用已知的阳性药物作为种子药物进行预测是可行的。

同时，实验组与阴性对照组的显著区别随着种子药物数量的增加而增大。如图1所示，当“种子药物”数量从2变化到20时，5种感染性疾病实验组的预测结果中，阳性药物的比例随着种子药物数量增加而增加，提示种子药物数越多，由这组种子药物预测的与该疾病相关的基因模块特征越准确，预测结果中阳性药物的比例越高，最终阳性率达到一个相对稳定的水平。以上结果提示，种子药物达到一定数量后，新增种子药物所获得的疾病特征信息增益趋于平稳，具有“饱和效应”。所有5种感染性疾病的阴性对照组A和阴性对照组B，其阳性率并不随着种子药物数量的增大而提高，显现出基本持平或略有下降趋势（图1）。

2.2 aCODE方法在治疗COVID-19药物发现中应用的初步探索

分别计算4个治疗COVID-19的种子药物和候选药物的药物-基因模块相关性、及其在基因模块维度空间的特征向量Vm和Vj，计算Vm与Vj的皮尔森相关性，并以皮尔森相关系数排序。前100个候选药物为aCODE方法预测所得具有潜在治疗COVID-19活性的药物，如表S6（http://202.38.153.236:81/Jweb_cjpt/CN/item/downloadFile.jsp?filedisplay=202008 11174457.xlsx）所示。

Fig.1 Effect of different quantities of seed drugs on proportion of positive drugs in prediction results.For the above five infectious diseases，HIV infections（A），influenza（B），paramyxoviridae infections（C），bacterial infections（D）and whooping cough（E），respectively，M positive drugs approved by the United States Food and Drug Administration（FDA）were randomly select⁃ed as seed drugs，M was taken from 2 to 20，and the feature vector Vm of the seed drug set was calculated.Enumerate all candidate drugs，set it as drug j，its feature vector Vj was taken，and the Pearson correlation between Vm and Vj was calculated to obtain the Rho value of the correlation coefficient.Rho value was sort in descending order，and the Top 100 as predicted drugs were taken.The reliability of the method was evaluated by the positive rate of the prediction result.In the results，as the number of seed drugs changed from 2 to 20，the proportion of positive drugs in the experimental group gradually increased.However，the proportion of positive drugs in the prediction results of negative control group A and B did not change significantly with the increase in the number of seed drugs.

由于目前尚无研究报道确切的治疗COVID-19的天然产物，本研究采用已知的具有抗冠状病毒活性的天然产物作为验证集，从侧面验证预测得到的这100个潜在治疗COVID-19药物的阳性概率。

验证集通过文献调研获得，共得到55种具有抗冠状病毒活性的天然产物成分，调研结果如表S7（http://202.38.153.236:81/Jweb_cjpt/CN/item/down⁃loadFile.jsp?filedisplay=20200811174457.xlsx）所示，该55种天然产物成分全部为经体外实验验证具有抗冠状病毒活性的成分。

在aCODE预测所得的前100个天然产物中，与验证集重合的天然产物成分达到4种，分别是厚朴酚（magnolol）［25］、和厚朴酚（honokiol）［25］、杨梅素（myricetin）［3］和人参皂苷 Rb1（ginsenoside Rb1）［19］。其中，YU等［3］通过荧光共振能量转移的双链DNA解旋研究了杨梅素对SARS病毒解旋酶的影响。实验结果显示，杨梅素可通过影响ATP酶活性而在体外有效抑制SARS病毒解旋酶蛋白，且杨梅素对正常乳腺上皮MCF10A细胞无细胞毒性。WU等［19］和WEN等［25］均通过Vero E6细胞实验证明，和厚朴酚和人参皂苷Rb1对SARS病毒具有体外抑制作用。这预示以上预测得到的4种抗冠状病毒的天然产物成分，可能具有抗SARS-Cov-2活性，值得通过体外实验进一步验证。

为了考察该结果的统计显著性，本研究在候选药物中每次随机选取100个化合物，与验证集中的55个化合物做交集，重复5000次。其中23次与验证集交集的化合物数量≥4，计算得P=23/5000=0.0046，结果差异显著。

3 讨论

通过对抗细菌和病毒感染性疾病药物的基因模块特征向量进行分析发现，在每类疾病中，阳性药物在基因模块相关性方面存在显著的共性特征，且共性特征与疾病的发生发展机制、药物作用机制密切相关。说明抗感染性疾病药物除直接靶向病原微生物外，靶向宿主基因网络也是其重要特征，该特征正是可通过宿主基因模块实现药物相似性匹配、进而实现药物重定位与天然产物活性预测的重要基础。

有了上述理论上的可能性，继而需要探讨技术上的可行性，即通过少数几个种子药物，归纳该疾病的基因网络特征，进而对候选化合物潜在功效进行预测。本研究设计了一个模拟实验，即每次随机抽取不同数量的种子药物进行预测，评估结果的准确度；随机抽取500次，评估该方法的性能。对病毒或细菌感染性疾病的量化评价结果显示，实验组在种子药物数目区间为2到20之间，种子药物越多，所得预测准确度越高，统计学意义越显著；而阴性对照组不仅初始准确率低，且随着种子药物的增加准确率呈下降的趋势，且准确率不显著。实验组与阴性对照组在预测结果上的显著差异提示，在感染性疾病的疾病-基因相关性证据极少的情况下，通过有一定疗效的种子药物归纳该疾病的基因网络特征，继而实现药物筛选是可行的。种子药物与该疾病的匹配度越高，且种子药物数越多，预测结果准确率越高，这使得迭代式的开发成为可能。例如，当获取少数几个有一定疗效的药物时，通过种子药物在更大化合物范围内搜索潜在的活性药物进行验证，是一种可行的研发路径。

本方法在COVID-19药物发现和应用中的初步探索也进一步证实，抗病毒类药物的作用机制与其靶向宿主基因网络有密切的关系。《新型冠状病毒感染的肺炎诊疗方案（试行第七版）》指南中的药物利巴韦林、利托那韦和磷酸氯喹均可对宿主蛋白及相关通路产生影响。例如，利巴韦林抑制宿主肌苷单磷酸脱氢酶，随后耗尽GTP库是该药物作用机制之一［28］；利巴韦林还可通过消除辅助型T细胞2型炎症反应，对病毒发挥宿主的免疫调节作用，增强宿主对病毒的免疫力［29］。利托那韦可通过P-糖蛋白和多药耐药相关蛋白流出通道来限制其他蛋白酶抑制剂的细胞转运和流出［30］。磷酸氯喹可抑制宿主血管紧张素转换酶2的末端糖基化，血管紧张素转换酶2则是SARS-CoV和SARS-CoV-2靶向进入宿主细胞的受体［31-32］。

本研究结果提示，基于感染性疾病早期有限的用药经验反馈，推测并计算疾病特征进而对更大范围的候选药物进行量化分析是可行的。本研究提出的aCODE方法具有如下特点：①敏捷迭代性。本方法可以随临床实践的进展将初始药物列表进行迭代更新，从而优化系统的精准度与预测能力；②快速闭环性。输出的天然产物、中药或药食同源成分，可实施快速测试，形成敏捷研发测试闭环；③药理机制方面。可对候选药物与初始药物的药物作用的细胞分子生物学机制进行进一步的聚类和比对。

本研究旨在提出通过迭代计算逼近更理想药物的方法，重点不是对相关药物进行验证。本方法的局限性在于预测结果对种子药物的依赖性较大，对于有一定疗效但由于副作用而临床获益不显著的候选药物，本方法是否能够归纳多个药物疗效的共性特征而实现有效预测，还有待观察。应用本方法尝试COVID-19治疗药物发现时，由于目前COVID-19治疗药物研发处于早期探索阶段，我国《新型冠状病毒感染的肺炎诊疗方案（试行第七版）》中推荐的利巴韦林等4个抗病毒药物的临床获益究竟有多大，还需根据对这些药物后期研究的结果确定。有相关综述通过对9篇关于COVID-19的临床试验、体外实验和专家报告调研发现，氯喹和羟氯喹具有治疗COVID-19的活性［33］。然而，氯喹的药动学和安全性研究表明，氯喹虽然可用于急性病毒感染的治疗，但其剂量安全范围很窄，尤其是儿童［34］。因此，在尚无临床批准的治疗COVID-19药物的前提下，将我国《新型冠状病毒感染的肺炎诊疗方案（试行第七版）》中推荐的药物作为种子药物，预测结果虽然具有一定的可参考性，但对其可靠度应保持质疑。

aCODE方法在以下几方面可持续更新：①药物靶标数据库的更新。本研究目前使用的药物靶标数据来自公开的STITCH数据库，后续可以接入更精细的机构专有药物靶标数据库。②疾病相关基因的更新。如果有病毒感染的宿主细胞组学等数据，可对本研究用到的疾病相关基因模块数据进行更新。③基因模块特异性的更新。当前所用的MsigDB基因模块，大部分来自基因本体（Gene Ontology），该知识表示体系是一个层级体系，为了兼顾特异性和一般性，选择了基因数量在50～500个之间的基因模块；后续可以对该数据库进行扩充或微调，并且优先选择基因数量比较小、特异性更高的基因模块。④其他药物证据的整合能力。本方法是一个纯数据驱动的计算过程，如果该结果与其他来源的临床或文献证据相符，按贝叶斯原理，则提示该候选化合物的成药性增强，后续可以结合其他因素进一步测试其成药潜力。