基于生物信息学、网络药理学和分子对接探讨芦丁抗新冠肺炎和乳腺癌的分子机制

叶宝倩，鲁可，罗慧燕，严倩，王雄文

（1.广州中医药大学第一临床医学院，广东 广州 510405；2.广州中医药大学第一附属医院肿瘤中心，广东 广州 510405）

新型冠状病毒肺炎（Corona Virus Disease 2019,COVID-19）是一种致命的流行病，从2019年开始在世界范围内大流行，截至2021年12月26日，全球报告的病例已经超过2.78 亿，死亡人数将近540 万[1]。目前，全球COVID-19 病例仍在激增，尽管迫切需要更有效的治疗策略，但仍然没有针对COVID-19 的特定抗病毒治疗方法，且COVID-19 的治疗指南因国家而异[2]。我们仍然急需筛选和验证针对COVID-19 的潜在抑制剂，以应对COVID-19 全球疫情的挑战。此外，有研究发现，乳腺癌（breast cancer,BC）、肺癌和结直肠癌等癌症患者感染COVID-19 的风险可能高于非癌症患者，且预后较差[3,4]。癌症已被证明是COVID-19相关死亡率的独立不良预后因素。值得关注的是，根据世界卫生组织国际癌症研究机构（International Agency for Research on Cancer,IARC）发布的2020年最新数据显示，乳腺癌新增人数达226 万，肺癌为220 万。这意味着乳腺癌正式取代肺癌，成为全球最常见癌症[5]。虽然到目前为止，可用数据无法让我们确定是否存在更容易感染COVID-19 的癌症类型，因为病例系列具有广泛的异质性[6]。但一份来自COVID-19 和癌症联盟注册数据库的队列研究报告数据显示，乳腺癌（21%）是2020年3月至2020年4月累积的1018 例COVID-19 病例中最常见的癌症类型[7]。巴黎的一项研究表明，乳腺癌合并COVID-19 患者的死亡率更多地取决于非癌症合并症，而不是先前的放射治疗或当前的抗癌治疗[8]。乳腺癌的治疗包括化疗、放疗、靶向治疗和免疫治疗等，会削弱免疫系统，并可能导致肺部问题。如果在此基础上感染病毒，发生并发症的风险将会更高[9]。中国湖北的一项多中心回顾性研究表明，在7 d内接受化疗的乳腺癌合并COVID-19 患者发生重症的风险更高，推测原因是化疗可能会引起骨髓抑制和继发感染，从而导致病情加重和COVID-19预后不良[10]。因此，有必要重视癌症特别是乳腺癌患者潜在的COVID-19 暴露风险可能，及早探讨可能的干预治疗措施，避免疾病向危重症发展。

芦丁（rutin）也称为芸香苷、维生素P，是一种植物来源的天然强抗氧化剂，广泛存在于自然界，几乎所有的芸香科和石楠科植物中均有之，也是毛冬青、连翘和槐角等常用中药材的有效成分之一，具有多种生物活性，如抗氧化、抗炎、抗菌、抗病毒和抗癌等作用[11,12]。黄曦文等[13]总结发现，芦丁在诱导乳腺癌细胞周期阻滞、调控相关蛋白表达水平以抑制乳腺癌侵袭转移、提高乳腺癌细胞对化疗的敏感性等方面具有重要作用。Elsayed Heba E 等[14]通过实验证明，芦丁可以作为抗三阴性乳腺癌的新型c-Met 受体抑制先导物。此外，多项研究已经证明，芦丁可以稳定且有效地抑制COVID-19 的主要蛋白酶，是对COVID-19 有治疗潜力的候选药物[15-17]。本研究基于生物信息学、网络药理学和分子对接，对芦丁进行深入的研究和预测，为进一步研究中药对抗COVID-19 和乳腺癌的作用机制提供有意义和有用的信息。

1 材料与方法

1.1 资料库及研究流程

具体研究流程见图1。

图1 芦丁抗新冠肺炎和乳腺癌的关键分子机制研究流程Fig.1 Flow chart of key molecular mechanisms of rutin in the treatment of COVID-19 and breast cancer

1.2 BC/COVID-19相关疾病靶点获取与差异基因筛选

于2021年12月4日从癌症基因组图谱（the cancer genome atlas,TCGA）数据库（https://portal.gdc.cancer.gov/）下载BC 的转录组数据。另外，以“COVID-19”为关键词，从Genecards 数据库（https://www.genecards.org/）、OMIM 数据库（https://www.omim.org/）、KEGG 数据库（https://www.genome.jp/kegg/）和NCBI 数据库基因组模块(https://www.ncbi.nlm.nih.gov/genome/)下载与COVID-19 相关疾病靶点，将从上述4 个数据库获得的基因靶点合并后去重。使用R 语言对BC 的转录组数据进行预处理，以基因表达数据大于0 作为过滤条件，然后与COVID-19 的潜在靶点取交集，提取COVID-19相关基因表达数据。再使用R语言Bioconductor的“limma”包筛选差异基因，将显著差异基因的筛选条件设定为错误发现率（false discovery rate,FDR）＜0.05和|log FC|＞1[18,19]，最终获得表达有显著性差异的BC/COVID-19 相关疾病靶点，并将结果以火山图形式输出。

1.3 临床相关性分析

通过TCGA 数据库下载BC 的临床数据，共获取1 097 条包括性别、年龄、生存时间、生存状态、病理分期等与之对应的临床信息用于后续分析。对数据进行预处理，删除生存时间≤30 d 的数据。将生存时间与R 语言的“limma”包筛选得到的差异基因的表达数据进行合并。使用R 语言的“survival”包进行生存分析和建模。对筛选得到的差异基因进行单因素Cox 回归生存分析（P＜0.05），并绘制森林图。对单因素Cox回归分析得到的相关基因进行多因素逐步回归分析，构建预测BC/COVID-19 患者预后风险的多因素Cox模型。将风险值的中位数作为阈值，将患者分为低风险组和高风险组，使用R 语言的“survival”“survminer”包进行生存分析，比较高低风险组生存是否具有差异，并绘制生存曲线[20]。进一步分析每个BC 患者的风险评分（risk score）、生存状态和预后基因的表达分布，并使用R 语言绘制出风险曲线、生存状态图和风险热图。对模型进行单因素与多因素独立预后分析，绘制临床相关性森林图。使用R 语言的“survivalROC”包绘制受试者工作特征曲线（receiver operating characteristic curve,ROC）曲线，利用ROC 曲线评价模型性能。

1.4 芦丁相关药物靶点搜集与筛选

以“Rutin”为关键词，通过可访问的在线平台与软件工具，包括中药系统药理学数据库和分析平台（Traditional Chinese Medicine Systems Pharmacology Database andAnalysisPlatform,TCMSP）（https://www.tcmsp-e.com/）、Swiss Target Prediction 数据库（http://www.swisstargetprediction.ch/）、TargetNet 数据库（http://targetnet.scbdd.com/）、Batman 数据库（http://bionet.ncpsb.org/batman-tcm/）和Drugbank数据库（https://www.drugbank.com/）[21,22]，搜集和筛选芦丁相关药物靶点。将搜集到的靶点合并后去重。

1.5 GO功能富集分析与KEGG通路富集分析

将R语言筛选后得到的表达有显著性差异的BC/COVID-19 相关疾病靶点与合并后去重得到的芦丁相关药物靶点取交集，获得药物 疾病共同靶点。使用R语言包，包括“clusterProfiler”“org.Hs.eg.db”“enrichplot”“ ggplot2”和“pathview”，对药物 疾病共同靶点进行GO 富集分析和KEGG 富集分析[23-25]。首先安装R 语言Bioconductor 的“org.Hs.eg.db”包，将药物 疾病共同靶点转换成entrezID。安装“clusterProfiler”包[26]和“pathview”包[27]，将已转换的entrezID，以P＜0.05，Q＜0.05 为过滤条件，对药物-疾病共同靶基因进行GO 与KEGG 富集分析，最后用“enrichplot”“ggplot2”和“pathview”包，将结果以气泡图形式输出。

1.6 蛋白质相互作用网络分析与核心靶点筛选

将上述药物 疾病共同靶点导入STRING 数据库（https://cn.string-db.org/）进行分析，将物种限定为“Homo sapiens”，根据所需设置靶点间的要求和得分（置信度设为0.9），随后导出蛋白质相互作用网络（protein-protein interaction network,PPI）图及其对应的相互作用表格（TSV 格式）。将TSV 表格导入Cytoscape软件，使用Cytoscape 软件的Network-Analyzer 工具计算网络相关属性，结合R 语言，进一步筛选出PPI网络核心靶点[22]。

1.7 分子对接

分子对接能有效确定与靶受体部位空间、作用力和电性特征匹配的小分子化合物，预测与靶点结合的可能性[28]。为观察芦丁进入机体后，能否和机体的靶蛋白结合进而发挥药效，将上述获得的核心靶点及COVID-19 与芦丁进行分子对接。通过PubChem 数据库（https://pubchem.ncbi.nlm.nih.gov/）[29]获得芦丁的分子结构。通过PDB 数据库（https://www.rcsb.org/）[30]检索核心靶点及COVID-19 的蛋白质三级结构，根据配体结构下载pdb 文件。使用Discovery Studio（DS）软件处理蛋白质pdb 文件，用AutoDock Tools 1.5.6进行分子对接，PyMOL 处理分子对接结果[31]。

2 结果与分析

2.1 BC/COVID-19相关疾病靶点获取与差异基因筛选结果

首先，使用TCGA 数据库下载了BC 的转录组数据，其中正常或癌旁组有113 例样本，肿瘤组有1 109例样本。其次，从Genecards 数据库、OMIM 数据库、KEGG 数据库和NCBI 数据库基因组模块，共收集了2 612 个与COVID-19 相关的基因。对BC 的转录组数据进行预处理后，将这两个基因簇取交集，提取了2 539 个COVID-19 相关基因的表达数据。使用R 语言的“limma”包筛选差异基因，最终获得394 个表达有显著性差异的BC/COVID-19 相关疾病靶点，其中223 个基因表达上调，171 个基因表达下调（图2）。

图2 BC/COVID-19 相关差异基因表达火山图Fig.2 Volcano map of BC/COVID-19 related differential genes

2.2 BC/COVID 19 相关基因的临床相关性分析结果

为了解BC/COVID-19 相关基因与BC/COVID-19临床因素的相关性，我们对394 个差异基因进行了单因素和多因素Cox 回归分析。首先，单因素回归Cox分析确定了42 个基因与BC/COVID-19 显著相关（P＜0.05）（图3，表1）。而多因素回归Cox 分析确定了42 个基因中的20 个靶基因，包括DPP4、IL18、PDCD1、IDO1、ANO6、SDC1、PF4、ATP7B、CXCL9、NT5E、GSN、WLS、LIMCH1、EEF1A2、EZR、BAMBI、SLC35A2、KIAA0319、SLC16A6 和HBA1（表2）。进行风险预后模型的构建，即risk score=（0.040 775 916 DPP4）-（0.039 682 033 IL18）+（0.167 157 414 PDCD1）-（0.027 784 36 IDO1）+（0.046 577 695 ANO6）+（0.001 965 139 SDC1）+（1.420 090 905 PF4）-（0.047 268 456 ATP7B）-（0.006 119 714 CXCL9）+（0.075 132 227 NT5E）-（0.009 871 433 GSN）-（0.039 016 82 WLS）+（0.033 558 393 LIMCH1）+（0.001 216 977 EEF1A2）+（0.002 765 01 EZR）+（0.004 978 992 BAMBI）+（0.021 846 763 SLC35A2）+（0.290 057 506 KIAA0319）-（0.014 289 454 SLC16A6）-（0.544 899 518 HBA1）。

表1 单因素Cox 回归分析结果Table 1 Results of the univariate Cox regression analysis

表2 多因素Cox 回归分析结果Table 2 Results of the multivariate Cox regression analysis

图3 单因素Cox 回归分析森林图Fig.3 Forest plot of univariate Cox regression analysis results

基于代表患者风险的多因素Cox 比例风险回归分析的Coef 值（表2），我们将患者分为高风险组和低风险组。在总体生存分析中，高风险组和低风险组之间观察到显著差异（图4A）。我们进一步分析了每个BC患者的风险评分、生存状态和上述20 个基因的表达分布，并绘制出风险曲线（图4B）、生存状态图（图4C）和风险热图（图4D）。结果表明，患者的风险值越大，风险评分越高（图4B)。

图4 多因素Cox 回归分析探讨BC/COVID-19 相关基因的预后价值Fig.4 Multivariate Cox regression analysis to explore the prognostic value of BC/COVID-19 related genes

此外，我们对这20 个基因进行了单因素与多因素独立预后分析（图5）。结果表明，年龄、分期和风险评分在单因素与多因素独立预后分析中都具有显著统计学差异（P＜0.05），均可作为独立预后因子。三者的风险比率（hazard ratio,HR）均大于1，提示均为高风险因素。用ROC 曲线解析预测模型（图6），结果显示，多因素Cox 模型构建的疾病风险评分得到的曲线下面积（area under the curve,AUC）值为0.811，提示该模型具有较高的临床应用价值。基于AUC 值，与其他临床性状构建的模型进行比较，可以得出本研究构建的预后模型为最优模型。

图5 独立预后分析森林图Fig.5 Forest plot of the results of the independent prognostic analysis

图6 基于多因素Cox 比例风险回归模型的ROC 曲线Fig.6 ROC curve based on the multivariate Cox proportional hazard regression model

临床相关性分析表明（表3），IL18 的表达与原发肿瘤的范围相关，在T1-T2中的表达高于T3-T4。DPP4、IDO1、SDC1、EZR、SLC35A2、SLC16A6的表达与淋巴结转移和淋巴结区扩散的数量和范围有关。IDO1 的表达水平与BC远处转移呈明显的负相关。同时，IDO1、IL18 的表达还与分期相关，在I 和II 期中的表达高于III和IV期（图7）。此外，CXCL9、DPP4、EZR、GSN、LIMCH1、PDCD1、PF4、SDC1、SLC35A2 和WLS 的表达与年龄相关（图8）。

图7 20 个基因的临床预后分析（1）Fig.7 Results of clinical prognostic analysis of 20 genes(1)

图8 20 个基因的临床预后分析（2）Fig.8 Results of clinical prognostic analysis of 20 genes (2)

表3 临床相关性分析Table 3 Clinical correlation analysis

2.3 芦丁相关药物靶点搜集与筛选结果

使用TCMSP、Swiss Target Prediction、TargetNet、Batman 和Drugbank 数据库获取芦丁相关作用靶点。去掉重复基因后，共获得16 392 个芦丁相关靶标。将394 个表达有显著性差异的BC/COVID-19 相关疾病靶点与芦丁相关作用靶点取交集，获得332 个药物 疾病共同靶点（图9）。

图9 药物 疾病共同靶点韦恩图Fig.9 Venn diagram of drug disease common targets

2.4 GO功能富集分析与KEGG通路富集分析结果

对药物-疾病共同靶点进行GO 和KEGG 富集分析。GO功能富集分析包括3 个部分：细胞组分（cellular component,CC）、分子功能（molecular function,MF）和生物过程（biological process,BP）。GO 功能富集分析结果表明，药物 疾病共同靶点在生物过程中，主要参与白细胞迁移、炎症反应的调节、调节白细胞活化等。在细胞组分层面上，主要富集在细胞质囊泡腔、囊泡腔、细胞外基质等上。在分子功能层面上，主要与受体调节活性、受体配体活性、细胞因子活性、细胞因子受体结合等有关。按照Q 值大小（Q＜0.05），将排名前10 位的条目以气泡图形式输出（图10A）。KEGG通路富集分析结果显示，共同靶点主要富集于2019 新型冠状病毒、细胞因子 细胞因子受体相互作用和补体和凝血级联途径等相关通路上。按照Q 值大小（Q＜0.05），绘制前30 条信号通路气泡图（图10B）。气泡的大小表示该条目下差异基因的个数，点越大表示基因数越多。颜色表示Q 值，越红表示Q 值越小，富集越明显。

图10 药物 疾病共同靶点的富集分析Fig.10 Gene ontology analysis and Kyoto Encyclopedia of Genes and Genomes pathway of drug-disease common targets

2.5 蛋白质相互作用网络分析与核心靶点筛选结果

将上述药物 疾病共同靶点导入STRING 数据库进行PPI 网络分析，并构建PPI 图（图11）。使用Cytoscape 软件的Network-Analyzer 工具计算网络相关属性。因介度中心性（betweenness centrality,BC）、接近中心性（closeness centrality,CC）、度中心性（degree centrality,DC）、特征向量中心性（eigenvector centrality,EC）、局部边连通性（local average connectivity,LAC）和网络中心性（network centrality,NC）等6 个参数代表整个网络中每个单独节点的拓扑重要性，其量化值越大，说明节点在该网络中越重要。因此以这6 个参数的中位值作为筛选条件，使用R 语言将数据过滤2次，每次仅选取所有参数均大于中位值的靶点进行下一次筛选，末次筛选所剩的靶点则为核心靶点[32-34]。第1 次筛选，以BC≥0、CC≥0.020 080 321、DC≥2、EC≥0.001 318 530 5、LAC≥0、NC≥0 为筛选条件，筛选后得到靶点39 个（图12A）。第2 次筛选，以BC≥17.033 333 33、CC≥0.408 602 151、DC≥4、EC≥0.080 259 524、LAC≥1.714 285 714、NC≥3 为筛选条件，最终确定了12 个核心基因靶点（图12B）。核心基因靶点包括AGT、CCL2、IL6、CSF3、STAT1、CXCL10、PIK3R1、SOCS3、CXCL2、IL12B、CCL11 和CXCR4。

图11 药物 疾病共同靶点PPI 网络图Fig.11 PPI network diagram of drug disease common targets

图12 核心基因靶点筛选Fig.12 Screening of the core target genes

2.6 分子对接结果

COVID-19 主蛋白酶在COVID-19 病毒复制过程中起关键作用，且具有不同于人体蛋白酶的特异切割位点，被认为是开发抗COVID-19 药物的理想靶标[35,36]。为了确定芦丁与12 个核心基因靶点及COVID-19 的可能结合以及确定的核心靶点，进行了分子对接分析。核心基因靶点及COVID-19 主蛋白酶的晶体结构是从PDB 数据库中收集的，PDB ID 为见表4，用于与芦丁的对接分析。对接的结果见图13，表4。分子对接得分代表分子对接吻合度，分数越小代表配体与受体蛋白结合性越好，发生作用的可能性越大。一般认为，结合能低于 20.92 kJ/mol表示对接结果比较稳定[37]。我们的结果显示，芦丁与13 个靶点的最低结合均能小于或者等于 20.92 kJ/mol，说明芦丁与以上核心基因靶点结合活性较好。其中，芦丁与CXCR4 和SOCS3的对接效果最好。

图13 芦丁与12 个核心基因靶点及COVID-19 的分子对接结果Fig.13 Results of molecular docking of rutin with 12 core gene targets and COVID-19

表4 芦丁与12 个核心基因靶点及COVID-19 的分子对接信息Table 4 Information on the molecular docking of rutin with 12 core gene targets and COVID-19

3 讨论

COVID-19 目前仍在全球范围内大流行，感染率和死亡率不断上升，病毒变异毒株的出现使其传播力更为显著，许多国家因此而经历了COVID-19 的第二次或第三次爆发[38]。疫苗接种是一项重要的预防性健康措施，可预防有症状和严重的COVID-19，但仅靠疫苗无法摆脱这场危机[39,40]。我们需要探索和开发出更安全、更特异的针对COVID-19 的疗法。此外，癌症患者被认为是COVID-19 的高度易感人群，且多项临床研究结果表明，由于肿瘤生长直接引起和抗癌治疗间接引起的免疫功能低下，患有COVID-19 的恶性肿瘤患者更易发展为重症且预后较差[7,41-43]。乳腺癌已上升为全球第一大癌，其防控形势仍然严峻。加之在COVID-19 流行期间，乳腺癌患者面临着更多的风险因素。Perrine Vuagnat[8]等人的一项临床研究显示，乳腺癌合并COVID-19 感染的患者，死亡率更多地取决于COVID-19 合并症。我们以乳腺癌患者作为切入点，研究合并COVID-19 感染的预防性干预措施，有助于防治肿瘤合并COVID-19 患者重症化，对抗击疫情有积极的意义。

中医药在抗击COVID-19 疫情中发挥了重要的作用。从有效中成药和方剂中寻找活性成分并进一步揭示其作用机制，具有重要的科学和应用价值。研究发现，中药可以通过多种成分作用于多种途径的多个靶点，在COVID-19 的治疗中发挥抗病毒、抗炎和免疫调节和靶器官保护的作用。其中，芦丁在多项研究中均被提及可以作为治疗COVID-19的有力候选先导物[44-49]。且研究证实，芦丁在抗乳腺癌中也同样发挥着重要的作用[14,50,51]。因此，我们假设芦丁可能对合并COVID-19 感染的乳腺癌患者发挥有效的药理活性。

在本研究中，我们首先筛选并确定了394 个BC/COVID-19 相关疾病靶点。差异表达基因分析显示，有223 个基因表达上调，171 个基因表达下调。对疾病相关靶点进行了临床相关性分析。根据独立的预后和生存分析，最终确定了20 个重要的差异基因，包括DPP4、IL18、PDCD1、IDO1、ANO6、SDC1、PF4、ATP7B、CXCL9、NT5E、GSN、WLS、LIMCH1、EEF1A2、EZR、BAMBI、SLC35A2、KIAA0319、SLC16A6 和HBA1，可作为预测乳腺癌合并COVID-19 感染患者预后的生物标志物。综上所述，这394 个疾病相关靶点可能为潜在的治疗靶标。此外，我们搜集和筛选芦丁相关药物靶点，并与上述疾病相关靶点取交集，获得332 个药物 疾病共同靶点。将332 个靶点进行PPI 网络分析，得出核心基因靶点CXCR4、SOCS3 和PIK3R1 等共12 个。CXCR4 又称趋化性细胞受体因子4，其在超过40% 的乳腺癌组织中高表达[52]，与乳腺癌的转移密切相关。有研究证明，CXCR4 的过表达可导致LL-37 在乳腺癌细胞中诱导的促迁移信号传导和迁移率的增强[53]。SOCS3 即为细胞因子信号转导抑制分子3，在涉及多种病毒的病毒免疫逃逸中发挥重要作用，SOCS1/3 拮抗剂有作为COVID-19 治疗剂的可能[54]。PIK3R1 是炎症基因，Shun Li 等人利用单细胞测序技术揭示，无论是COVID-19 重度或危重患者还是有中度症状患者，PIK3R1 在总T 细胞、亚群CD4+T 细胞和CD8+T细胞中均为高表达[55]。由此可见，这12 个核心基因靶点可能是芦丁治疗乳腺癌和COVID-19 的有效药理学靶点。结合分子对接的结果，也进一步说明芦丁在乳腺癌和COVID-19 的治疗中主要是通过CXCR4、SOCS3 和PIK3R1 等靶点发挥作用，且可以有效地与COVID-19 中的特定蛋白质结合。药物-疾病共同靶点的GO 和KEGG 富集分析结果表明，芦丁发挥的抗乳腺癌和抗COVID-19 作用主要与白细胞迁移与活化的调节和炎症反应的调节等相关，且相关靶点在COVID-19、细胞因子-细胞因子受体相互作用和补体和凝血级联途径等相关通路显著富集。

4 结论

本研究采以芦丁的靶蛋白、信号通路以及生物学功能的调控网络为基础，通过生物信息学、网络药理学及分子对接分析，揭示芦丁多靶点、多通路抗乳腺癌和COVID-19 的分子机制，为COVID-19 疫情下乳腺癌患者的预防性治疗提供更广泛、更深入的研究。结合上述研究结果可知，辅助补充芦丁可能会提高合并COVID-19 感染的乳腺癌患者的临床疗效。