基于通路扰动相似性建模的抗类风湿关节炎药对分子协同组合发现及机制分析＊

孙继佳，王睿瑞，张磊，刘保成＊＊，吕爱平，3，4＊＊

（1.上海中医药大学中医健康协同创新中心上海201203；2.上海中医药大学数学与物理教研室上海201203；3.上海中医药大学中西医结合研究院上海201203；4.香港浸会大学中医药学院中国香港999077）

类风湿关节炎（Rheumatoid Arthritis，RA）是一种发生在滑膜关节及其他器官系统的慢性、全身性、炎症性疾病，也是一种慢性进行性自身免疫性疾病[1]。近年来随着分子生物学、免疫学、遗传学的研究发展，显示RA的发病机制与遗传因素、细菌和病毒因素、T淋巴细胞、B淋巴细胞、细胞因子以及其他免疫细胞等的改变和影响有关[2]。RA的发病率随着年龄的增长而增加，全世界每年约有0.3%-1.0%的人受其影响[3]，以35-50岁成年人多见，且女性发病率为男性的3倍左右[4]。目前，临床上治疗RA的药物仍以西药为主，但往往存在毒副作用较大[5-7]、治疗费用较高等特点，因此进一步提高抗RA药物的有效性、安全性以及降低治疗成本，是当前亟需解决的问题。

中医药治疗类风湿关节炎历史悠久，积累了丰富的临床应用经验，RA属于中医“痹证”的范畴，具有痹症共有的病因[8-9]。《素问·痹证论》曰：“风寒湿三气杂至，合而为痹也”，指出风寒湿是引起痹证的主要因素。RA中医辨证分型主要有脾肾亏虚、风寒湿痹、痰湿阻络、湿瘀痹阻等4种证型[10-11]，并认为肝肾不足、气血亏虚是RA发生的内在因素，而风寒湿是RA发生的诱发因素。目前，中药及复方在治疗RA的过程中，许多有效成分之间存在的组合机制尚不完全清楚，有效成分之间的协同性很少能够从分子调控网络角度进行更为精准地进行定量描述；因此，需要建立更为适合的研究方法和技术，发现治疗中起决定作用的主要成分，同时能够筛选出具有协同治疗作用的有效成分组合，以实现对中药及复方成分的进一步“优化”，进而为中药新组分开发提供理论和方法学上的支持。

药对是中医临床遣药组方常用的配伍形式，是历代医家逐渐积累经验而来的精妙应用形式，药对通过协同增效、相制减毒、相反相成等形式应用，介于中药和方剂之间，是一个值得研究的群体[12-14]。那么，如何从分子机制角度深入探索药对中不同成分之间的协同组合效应并进一步分析其作用机制，是亟需解决的实际问题。鉴于此，本研究从治疗RA的典型药对出发，综合运用生物信息学、化学信息学、网络药理学等理论和技术，并创新性地提出一种基于通路谱扰动相似性的计算方法，阐明不同药对中有效成分在抗RA中的协同组合特征；同时，进一步挖掘和分析药对中有效成分作用的分子机制及功能，旨在为中药药对作用分子机理阐明及新组分开发提供有力的方法学上的支持。

1 材料与方法

1.1 类风湿关节炎（RA）差异基因筛选

在GEO（Gene Expression Omnibus，https://www.ncbi.nlm.nih.gov/geo/）[15]数 据库中，输 入“Rheumatoid Arthritis”查询和收集类风湿关节炎相关的基因表达谱芯片，获得GSE55235、GSE55457和GSE77298这3个与RA有关的表达谱数据；其中，GSE55235的平台为GPL96[HG-U133A]Affymetrix Human Genome U133A Array，包含10个正常样本与10个RA样本；GSE55457的平台为GPL96[HG-U133A]Affymetrix Human Genome U133A Array，包含10个正常样本和13个RA样本；GSE77298的平台为GPL570[HG-U133_Plus_2]Affymetrix Human Genome U133 Plus 2.0 Array，包含7个健康样本和16个RA样本。

利用R软件GEOquery包下载所有3个表达谱数据芯片，并对原始数据进行标准化、校正和基因名注释等预处理；使用limma包对GSE55235、GSE55457和GSE77298分别进行差异表达基因分析（differentially expressed genes，DEGs），并根据|logFC|>1.0及校正adj.p.value<0.05筛选出每组芯片数据中的上调和下调差异表达基因。

1.2 类风湿关节炎（RA）治疗靶点数据收集

我 们 依 次 从TTD（http://db.idrblab.net/ttd/）[16]、DrugBank（https://www.drugbank.ca/）[17]数据库收集与类风湿关节炎相关的靶点数据；在TTD中输入“Rheumatoid Arthritis”输入查找找出与RA相关的基因；在DrugBank中也同样通过输入“Rheumatoid Arthritis”搜索与治疗RA相关的已认证药物及其作用的靶点。

再将基于GEO的差异基因分析筛选后的结果和通过数据库查询得到的靶点基因进行合并、去除重复；通过UniProt（https://www.uniprot.org/）[18]数据库对所有收集到靶点信息进行确认，组成一个RA疾病相关靶点基因集合SRA。

1.3 抗类风湿关节炎（RA）典型药对数据收集

根据课题组以往的文献整理和分析，选择了治疗类风湿关节炎的两种较为常见的药对，即：肝肾不足型的典型药对，独活-桑寄生、湿热痹阻型的典型药对，黄柏-苍术，以此进行下一步研究和分析。

从TCMID 2.0（http://www.megabionet.org/tcmid/）[19]、ETCM数 据 库（http://www.nrc.ac.cn:9090/ETCM/index.php/Home/Index/）[20]分别收集到与这4种中药相关的化学成分和靶点数据，所有的中药成分化学信息以及Canonical SMILES通过PubChem（https://pubchem.ncbi.nlm.nih.gov/）获得。

1.4 药对成分潜在靶点预测与识别

药对成分的靶点数据主要由SEA（http://sea.bkslab.org/）[21]、HitPick（http://mips.helmholtz-muenchen.de/proj/hitpick）[22]以及ETCM数据库来确定。根据药对成分所对应的Canonical SMILES，通过SEA和HitPick在线系统对抗RA药对有效成分进行潜在靶点预测，其中，选取SEA数据库中评分值Max Tc>0.6，HitPick在线预测系统Precision值>50%的预测结果作为药对成分的潜在作用靶点；ETCM中根据潜在候选靶点（Candidate Target Genes）评分值大于0.9，筛选出相应成分的靶点，根据靶点基因集SRA，使用Venny 2.1（https://bioinfogp.cnb.csic.es/tools/venny/index.html）工具进行取交集，识别出药对小分子潜在的抗RA靶点。

1.5 药物小分子结构相似性计算

一般认为，如果两个药物小分子A和B的化学结构越相似，则具有协同组合治疗的可能性就越大；因此，我们先通过DRAGON 7（http://www.talete.mi.it/index.htm）软件计算出药物小分子1024维的分子指纹，再利用Tanimoto相似性系数（Tanimoto Coefficientbased Similarity）计算药物小分子之间的结构相似性即：

其中，FP(A)和FP(B)分别为药物小分子A和B的分子指纹，SFPd(A,B)为分子结构相似度。

1.6 基于通路扰动相似性的中药小分子协同组合机制发现

如果两种药物所作用的靶点属于同一条信号通路，那么这两种药物分子协同组合的可能性就越高；根据现有方法，可以将每个药物小分子所作用的靶点富集到所对应的信号通路上，构建药物作用靶点的通路向量，以此来计算不同药物分子作用通路的协同值。

虽然，这种评估方法相对简单，但没有考虑到小分子对每个通路的扰动程度；因此，本文我们先利用clusterProfiler包[23]将药物小分子作用的靶点进行KEGG通路富集分析，按p.value<0.05和q.value<0.05得到相关通路；根据注释得到的相应通路集{TPk∣k=1,2,…,n}（n为通路个数）构建一个由包含所有注释通路所构成的向量，即：通路向量TP=(tp1,tp2,…,tpn)，再根据以下流程（如图1所示）构建协同组合评价模型：

图1 基于通路扰动相似性方法的小分子协同机制发现方法示意图

如果药物小分子A作用的靶点集T(A)中的靶点属于某一条通路TPk，那么通路谱TP对应位置上对这条通路的药物扰动强度值为tpk(A)，否则为0，其中，tpk由下式计算得到，即：其中，l为通路TPk的靶点个数，m为靶点集T(A)的个数，dij为靶点之间在RA疾病相关PPI网络上的最短路径距离；由此得到药物小分子A的靶点通路扰动向量TP(A)；同理，得到药物小分子B作用的靶点通路扰动向量TP(B)。

根据所得到的向量TP(A)和TP(B)，就可以计算两种药物分子之间的通路谱扰动相似性评分STPd(A,B)，另外，为防止出现分母为0的情况，利用改进余弦距离公式进行计算，即：

1.7 GO功能注释与KEGG Pathway富集分析

基因本体（Gene Ontology，GO）数据库是GO组织（GO Consortium）在2000年构建的一个结构化的标准生物模型，涵盖了基因的生物过程（Biological Process，BP）、分子功能（Molecular Function，MF）和细胞分组（Cellular Component，CC）。一个生物过程或通路通常是由一组基因共同参与，而不是由单个基因完成的，富集分析的主要依据是：如果一个生物学过程或通路在已知的研究中发生异常，则共同发挥功能的基因极可能被选择出来作为与这一过程或通路相关的基因集合。

本文，我们进一步利用clusterProfiler工具包对所有药对的潜在作用RA疾病的靶点集进行GO功能和KEGG Pathway富集分析，按照p.value<0.05和q.value<0.05进行筛选。

2 结果

2.1 类风湿关节炎疾病基因收集及PPI网络构建

经过分析，GSE55235中筛选得到1071个差异基因，其中，上调基因599个、下调基因472个；GSE55457中筛选得到312个差异基因，其中，上调基因189个、下调基因123个；GSE77298中筛选得到432个差异基因，其中，上调基因237个、下调基因195个。另外，将只要同时出现在任意两个芯片表达数据中的差异基因都列入RA疾病基因集，由此，共获得与RA疾病相关的267个差异基因（图2）。

从DrugBank数据库中收集到与RA治疗相关的73个已证实的药物及其对应的193个靶点基因；又从TTD数据库中收集得到139个与RA疾病相关的靶点基因；再将两个数据库所收集到的靶点集进行合并去除重复后，共获得309个与RA相关的靶点基因。

最后，再将上述两种方法（基于差异基因分析和靶点数据库收集）获得的所有基因进行合并后，总共得到555个与RA相关的靶点基因，即靶点基因集SRA；进一步再将SRA所有555个靶点基因导入STRING（https://string-db.org/）数据库得到蛋白相互作用关系，其中，参数设置：Organism设为Homo Sapiens，Combine Score阈值设为0.7；同时，采用Gephi0.9.2软件构建靶点蛋白的PPI网络（图2）。

图2 基于差异基因分析和药物靶点数据库的抗RA靶点基因收集及PPI网络构建

2.2 药对数据收集及潜在靶点预测与识别

从TCMID、ETCM数据库以及PubChem数据库分别收集和整理获得独活、桑寄生、黄柏、苍术、等4种中药成分及其相关的化学成分数据；其中，独活148个、桑寄生16个、黄柏62个、苍术95个；说明：为了便于后续分析，我们还对独活-桑寄生药对164个成分进行一一编号DS1-DS164，黄柏-苍术药对157个成分也依次编号为HC1-HC157。

利用SEA、HitPick以及ETCM，经过靶点预测和阈值筛选得到相关成分的潜在作用靶点，进一步，通过靶点基因集SRA识别出抗RA有效成分及其潜在作用靶点；其中，独活-桑寄生药对有53个成分及其作用的52个潜在RA靶点，黄柏-苍术药对有32个成分及其作用的28个潜在RA靶点，合并得到67个潜在作用靶点信息（表1）；不同药对的成分-靶点作用网络（图3）。

图3 （B）黄柏-苍术药对成分-作用靶点网络

图3 （A）独活-桑寄生药对成分-作用靶点网络

表1 两种药对所有67个潜在作用靶点信息

2.3 小分子协同组合评分结果

利用1.5所述，分别对两种药对中的小分子协同组合进行计算分析，根据协同分值Sd(A,B)值的大小来判断两种小分子是否存在可能协同作用及其强度；Sd(A,B)值越大，则对这对小分子对PPI网络的扰动程度就越大、协同组合治疗特征也越强；本文中，我们规定协同评分值Sd(A,B)≥0.15的小分子组合被认为是具有潜在的协同治疗作用，两种药对中具有协同治疗特性的小分子组合结果，如图4所示，所有具有协同治疗作用的中药小分子信息（表2）。

表2 具有协同组合治疗的2种药对小分子信息

图4 （B）黄柏-苍术药对协同小分子组合

图4 （A）独活-桑寄生药对协同小分子组合

在独活-桑寄生药对中，我们发现存在15对具有通路协同治疗特征的中药小分子组合，其中，独活中有4个小分子，分别为：DS51、DS75、DS120、DS137，桑寄生中有11个小分子，分别为：DS150、DS151、DS152、DS155、DS156、DS157、DS159、DS160、DS161、DS163、DS164。

在黄柏-苍术药对中，我们发现存在17对具有通路协同治疗特征的中药小分子组合，其中，黄柏中有14个小分子，分别为：HC9、HC16、HC23、HC24、HC33、HC39、HC60、HC9、HC16、HC23、HC24、HC33、HC39、HC60，苍术中则有8个成分，分别是：HC66、HC70、HC71、HC126、HC66、HC70、HC71、HC126。

2.4 GO功能注释与KEGG Pathway富集分析结果

利用clusterProfiler包进一步对两种药对的作用靶点集进行GO功能注释和KEGG Pathway富集分析，其中，p.value<0.05和q.value<0.05，两种药对作用靶点的功能注释与富集分析结果（图5）。

图5 （B）黄柏-苍术药对的作用靶点GO功能与KEGG通路富集分析结果

图5 （A）独活-桑寄生药对的作用靶点GO功能与KEGG通路富集分析结果

独活-桑寄生药对作用靶点的GO功能注释结果显示，所涉及的BP功能有526个，主要包括：有机阴离子运输（GO:0015711）、类固醇激素反应（GO:0048545）、调节炎症反应（GO:0050727）、营养水平反应（GO:0031667）、有机酸运输（GO:0015849）等；所涉及的MF功能有74个，主要包括：跨膜转运蛋白活性（GO:0022804）、核受体活性（GO:0004879）、转录因子活性（GO:0098531）、类固醇激素受体活性（GO:0003707）、一元羧酸结合性（GO:0033293）等；所涉及的CC功能有19个，主要包括：转录因子复合体（GO:0005667）、神经细胞体（GO:0043025）、RNA聚合酶II转录因子复合物（GO:0090575）、核转录因子复合体（GO:0044798）、膜筏（GO:0045121）等。KEGG通路富集分析显示，独活-桑寄生药对作用靶点主要富集在14个相关信号通路，分别为：胆汁分泌通路（hsa04976）、乙肝通路（hsa05161）、ABC转运通路（hsa02010）、花生四烯酸代谢通路（hsa00590）、PPAR信号通路（hsa03320）、利什曼 病 通 路（hsa05140）、PD-L1和PD-1通 路（hsa05235）、GnRH信号通路（hsa04912）、IL-17信号通路（hsa04657）、胰腺分泌通路（hsa04972）、甲状旁腺激素通路（hsa04928）、Th17细胞分化信号通路（hsa04659）、叶酸生物合成通路（hsa00790）和抗叶酸通路（hsa01523）等。

续表

黄柏-苍术药对作用靶点的GO功能注释结果显示，所涉及的BP功能有373个，主要包括：细胞对外部刺激反应（GO:0071496）、化学突触传递的调节（GO:0050804）、突触信号转导调控（GO:0099177）、有机阴离子转运（GO:0015711）、多细胞生物体内平衡（GO:0048871）等；所涉及的MF功能有36个，主要包括：嘌呤 能 受 体 活 性（GO:0035586）、激 素 结 合（GO:0042562）、ATPase偶联活性（GO:0042623）、跨膜转运蛋白活性（GO:0022804）、ATP酶活性（GO:0016887）等；所涉及的CC功能有27个，主要包括：神经细胞体（GO:0043025）、膜筏（GO:0045121）、膜微区（GO:0098857）、膜区（GO:0098589）、树突棘（GO:0043197）等。KEGG通路富集分析显示，黄柏-苍术药对作用靶点主要富集在10个相关信号通路，分别为：神经活性配体-受体相互作用通路（hsa04080）、胆汁分泌通路（hsa04976）、鞘脂信号通路（hsa04071）、cAMP信号通路（hsa04024）、cGMP-PKG信 号 通路（hsa04022）、VEGF信 号 通 路（hsa04370）、肾 细 胞 癌 通 路（hsa05211）、胰腺癌通路（hsa05212）、氮素代谢通路（hsa00910）、近端小管碳酸氢盐回收通路（hsa04964）等。

3 讨论

近年来，许多学者从不同角度，尤其是基于对文献的整理[24]、统计分析[25]、关联规则挖掘[26]的方法来总结和归纳治疗RA的核心药对组合，发现川乌-桂枝、当归-川芎、当归-黄芪、黄柏-苍术、独活-桑寄生等药对都具有较高使用频率的用药规律；此外，也有研究者利用网络药理学[27]、实验手段[28,29]来分析和阐述抗RA药对的药理作用机制，进一步揭示了药对配伍的特点。

药对小分子协同组合评价结果显示，独活-桑寄生药对中的DS51（4-Methoxyacetophenone）与DS163（Quercitrin-7-olate）具有最高的协同性评分值，且DS51与桑寄生中的其他7个成分（DS160、DS157、DS159、DS152、DS161、DS156、DS164）等都存在一定的通路协同作用机制。王友庆[30]等人研究证实，槲皮素（DS150，Quercetin）可通过抑制NF-κB激活，减弱炎症环境下软骨细胞内MMP-13产生、基质降解和细胞凋亡，起到保护软骨的作用，而通过计算结果显示，槲皮素可能与DS120这个小分子化合物具有协同治疗作用。邓莉[31]等人研究发现，TNF-α能促进成纤维样滑膜细胞IL-6及IL-1β的表达和分泌，并促进p38MAPK的磷酸化及NF-κB核转位，而齐墩果酸（DS152，Oleanolic acid）呈浓度依赖性抑制TNF-α诱导的HFLS炎症因子的表达（IL-6及IL-1β），并抑制TNF-α诱导p38MAPK的磷酸化及NF-κB核转位，而计算结果也提示4-Methoxyacetophenone可能在药对中与齐墩果酸具有组合治疗效应。黄柏-苍术药对中的H24（28-Isofucostanol）与HC126化合物具有最高的协同评分，苍术中的HC71（Icariside F2）与黄柏中的其他5个成分（HC33、HC16、HC9、HC39、HC60）均存在潜在的协同组合治疗作用。研究显示[32]，盐酸药根碱（HC9，Jatrorrhizine）能够抑制体外滑膜细胞的增殖，迁移和分泌，并改善类风湿关节炎的大鼠模型；还有研究显示[33]，Jatrorrhizine能抑制RANKL诱导的破骨细胞生成并防止磨损颗粒诱导的溶骨，已经具有发展成为治疗RA的新型治疗剂的巨大潜力。绿原酸（HC23，Chlorogenic acid）能够通过下调核因子κB配体诱导的活化T细胞c1表达的受体激活剂来抑制破骨细胞分化和骨吸收，可能是破骨细胞相关疾病伴发炎性骨破坏的潜在治疗选择[34,35]；另有研究也证明[36]，绿原酸和木犀草素（Luteolin）通过调节NF-κB和JAK/STAT信号通路的活化，协同抑制白介素1β诱导的成纤维样滑膜细胞的增殖。桑寄生中的儿茶素（HC66，Epicatechin）[37]具有抗氧化，抗炎，抗微生物，抗过敏，抗癌，抗血栓形成和保肝等治疗作用。由此，可以发现，通过我们所提出的计算和分析方法，不仅可以找药对中一些抗RA的主要活性成分，同时，也能够发现药对中的哪些其他成分与这些主要活性成分产生了“1+1>2”作用，从分子和定量角度揭示了中药药对配伍的特征。

从药对所作用靶点富集的通路来看，独活-桑寄生作用靶点所富集通路中，如已有研究报道IL-17能够调节类风湿关节炎中SHP-2和IL-17RA/STAT-3所依赖的Cyr61、IL-23和GM-CSF的表达以及RANKL介导的成纤维细胞样滑膜破骨细胞的生成，并可能揭示RA的发病机理的新机制[38]。已有报道还显示[39]，可以通过调节成骨细胞中花生四烯酸的代谢以对类风湿关节炎的起到治疗作用。黄柏-苍术药对中，例如，VEGF信号通路，目前已经明确报道与RA疾病多种发病分子机制有关[40,41]；还有研究显示[42]，可以通过调节MMP14和VEGF介导的间充质细胞表达miRNA-150-5p外来体来治疗类风湿关节炎。

此外，我们还发现独活-桑寄生药对主要涉及14个相关信号通路，黄柏-苍术药对则涉及10个相关信号通路，而两种药对中只有一个信号通路Bile secretion（hsa04976）是相同的，提示中医在RA治疗中针对不同证型采用不同的中药配伍，因此，涉及到对不同信号通路的作用调控，从一定程度上反映了中医“同病异治”的特征。

综上所述，本文我们综合了生物信息、化学信息学、网络药理学及计算建模等多种手段，通过GEO差异基因分析筛选RA关键基因，再结合已知数据库中的RA的相关靶点信息，构建了一个较为可靠的RA疾病PPI网络；选取了两种不同证型的经典药对：独活-桑寄生和黄柏-苍术，整理得到了两种药对的主要成分数据，并预测其潜在作用靶点；通过我们所提出的基于通络扰动的相似性计算建模方法，挖掘出了不同药对中具有协同组合治疗的小分子组合；本文所提出的研究方法，从一个全新角度揭示了中药药对内在成分的作用机制和规律，为后续深入开展实验研究提供技术上有力的支持。