乳腺癌是女性最常见的肿瘤之一,也是仅次于肺癌的癌症相关死亡的第二大原因
[1]。尤其在女性中,2020年约有685 000例女性患者死于乳腺癌,即每6个死于癌症的女性患者中就有1个是乳腺癌,预计到2040年,乳腺癌的新发病例会增加至300万例,死亡病例会增加至100万
[2]。随着早期筛查的进步以及抗癌策略的不断发展,乳腺癌的治疗有了很大的改善。然而,乳腺癌的复发及转移率仍然很高
[3]。目前乳腺癌的治疗方法主要包括手术治疗、化疗、放疗、靶向治疗、内分泌治疗等。肿瘤的复发转移被认为是治疗失败的重要因素,早期乳腺癌已经具有相对成熟的治疗方案,而对于晚期乳腺癌则尚未有得到广泛认可的标准方案
[4]。肿瘤转移过程中的基因再生使癌细胞能够对不同的化疗和放疗产生耐药性,这是晚期乳腺癌高病死率的原因,也是乳腺癌治疗过程中的主要障碍
[5-6]。因此,寻找更有效的治疗方式成为乳腺癌目前研究的重点。
一项纳入了729例晚期乳腺癌患者的回顾性队列研究
[7]结果表明,中医辅助治疗可提高晚期乳腺癌患者的生存率。山柰酚是一种广泛存在于中草药中的黄酮类化合物,具有抗癌、抗炎及抗菌等作用
[8]。流行病学研究
[9-10]表明,山柰酚的摄入与各种癌症的发生率降低有关,包括结直肠癌、肝癌、胰腺癌、膀胱癌、卵巢癌和胃癌。既往研究
[11-16]表明许多癌症相关分子,如基质金属肽酶(matrix metallopeptidase,MMP)、促凋亡和抗凋亡蛋白、各种生长因子、细胞周期蛋白和细胞周期蛋白依赖性激酶可受到山柰酚的调节。Kim等
[17]的体外实验结果表明山柰酚可通过上调组织蛋白酶、细胞周期蛋白D1、细胞周期蛋白E的表达,同时下调Bax和P21的表达拮抗三氯生对人乳腺癌细胞MCF-7的促进增殖作用。Zhu等
[18]体内实验结果表明,山柰酚可通过抑制COX-2/PGE2通路抑制乳腺癌细胞生长并治疗乳腺癌相关抑郁症。尽管山柰酚作为抗乳腺癌药物的潜力已得到了部分验证,但由于山柰酚作用靶点及信号通路的多样性,对于山柰酚治疗乳腺癌的具体作用机制仍未完全阐明。
网络药理学以系统生物学和多药理学为基础,提供了一种“多靶点、多效应、复杂疾病”的新型网络模式
[19]。可用作了解药物-靶标相互作用的非常有用的工具。本研究基于网络药理学方法,综合基因、蛋白质等生物分子网络,探究山柰酚治疗乳腺癌可能的作用靶点和通路,并筛选预后基因,建立预后模型,为临床上使用山柰酚治疗乳腺癌提供理论依据,同时为下一步的研究提供方向。
1 材料与方法
1.1 数据来源
HERB本草组鉴数据库(
http://herb.ac.cn)是一个整合了多个中药的数据库,包含了迄今为止最全面的草药及成分清单的数据库
[20]。以“Kaempferol”为关键词在HERB数据库中收集并下载与山柰酚相关靶点。然后使用UniProt数据库(
https://www.uniprot.org)将靶点名称转变为uniprot ID及靶基因名称。利用GeneCards数据库(
https://www.genecards.org)以“breast cancer”为关键词进行检索,筛选relevance score>5的靶点作为乳腺癌靶基因。利用venny 2.1.0将山柰酚靶基因与乳腺癌靶基因进行分析,获得交集靶基因。
1.2 蛋白互作网络(protein-protein interaction,PPI)建立
将山柰酚与乳腺癌的交集靶基因导入STRING数据库(
https://string-db.org)进行PPI的建立,将置信度设置为0.4,下载分析得到的数据并将其导入Cytoscape 3.10.2软件,绘制“药物-靶点”相互作用图及PPI网络图。
1.3 基因本体(gene ontology,GO)
分析及京都基因组百科全书信号通路富集分析(Kyoto encyclopedia of genes and genomes,KEGG)将山柰酚与乳腺癌的55个交集靶基因导入R 4.4.1,使用R语言中“cluster Profiler”包进行GO分析及KEGG分析,并设定
P<0.05,FDR<0.25
[21]。将其结果进行可视化。
1.4 分子对接
通过Pubchem数据库(
https://pubchem.ncbi.nlm.nih.gov)获取山柰酚的结构,在RSCD PDB数据库(
https://www.rcsb.org)收集前5位关键靶点蛋白结构,使用PyMol软件对靶点进行清除水分子和小分子配体等优化处理并用AutoDock Tools进行加氢与电荷处理。以关键靶点为受体,以其对应的有效成分作为配体,利用PyRx软件内部的Vina进行分子对接并进行结合能的计算。最终利用PyMol软件进行结果可视化处理。
1.5 LASSO-Cox回归预测模型的构建
从TCGA数据库(
https://portal.gdc.cancer.gov)获取1 098例乳腺癌肿瘤样本和113例正常样本的mRNA表达数据及962例乳腺癌患者的临床信息(
表1),利用R语言整理每例患者的生存时间、生存状态及靶基因在每例患者样本中的表达。将整理后的数据使用LASSO-Cox方法对靶基因进行进一步筛选,导入R语言中鉴定出乳腺癌的预后基因
[22]。建立风险评分公式,取各患者风险评分平均值,根据该平均值,将患者分为高危组和低危组,并绘制Kaplan-Meier生存曲线,以评估相关风险基因的预测性能。
1.6 单基因预后分析
在LASSO-Cox回归预测模型的基础上对纳入预后基因进行在线分析,获得预后基因的表达分别对乳腺癌患者总生存(overall survival,OS)期的影响并将其结果可视化。
1.7 预后模型与肿瘤免疫的相关性
使用TCGA数据库中获取的RNAseq和相应的临床信息,使用Spearman相关分析来评估风险模型评分与免疫评分之间的相关性(Spearman相关分析对于数据中异常值不敏感,结果更具稳健性
[23])。
P<0.05为差异有统计学意义。
2 结 果
2.1 山柰酚及乳腺癌基因靶点的获取
在HERB本草组鉴数据库中以“Kaempferol”为关键词进行检索,获取65个非重复靶点。使用GeneCards数据库以“breast cancer”为关键词进行检索,筛选后获得4 084个非重复靶点。将山柰酚与乳腺癌的靶基因导入venny 2.1.0后获得55个交集基因及韦恩图(
图1)。
2.2 药物-靶点网络及PPI网络建立
利用Cytoscape 3.10.2软件对山柰酚及其预测靶点进行连接,得到药物-靶点相互作用网络图(
图2)。将山柰酚与乳腺癌的55个交集靶蛋白导入String数据库,分析PPI关系,构建山柰酚与乳腺癌共同靶点的PPI网络(
图3),该网络包含55个节点和477条边,节点的平均度值为17.3。每个节点代表1个靶点,蛋白之间线条表示靶点之间互作关系,线条越多表示关联度越大,即度值越大
[24]。筛选出度值排名前5的靶点即Akt1、Bcl2、CASP3、ESR1、AR(
表2)。
2.3 GO分析和KEGG通路富集分析
将山柰酚与乳腺癌的55个交集靶基因导入R中进行GO分析及KEGG通路富集分析。获得GO条目共1 753条,其中生物学过程(biological process,BP)1 604条,涉及对外来刺激的反应、对肽的反应、对氧化应激的反应、细胞对化学应激的反应等。细胞组成(cellular component,CC)条目20个,涉及膜微区、膜筏、转录调节复合物、转移酶复合物,转移含磷基团等。分子功能(molecular function,MF)条目121个,涉及DNA结合转录因子结合、RNA聚合酶Ⅱ特异性DNA结合转录因子结合等。分别选取BP、CC、MF前10个条目,绘制GO分析气泡图(
图4)。KEGG通路富集共获得141条信号通路,根据
P值升序排列,筛选前10条信号通路(
图5),其中山柰酚治疗乳腺癌的可能相关通路包括化学致癌—活性氧信号通路、糖尿病并发症中的AGE-RAGE信号通路、TNF信号通路、IL-17信号通路等。
2.4 分子对接
将山柰酚与关键靶点Akt1、Bcl2、CASP3、ESR1、AR分别进行分子对接验证。一般认为,分子对接结合能<0时,说明两分子具有自发结合能力,分子对接结合能<-1.2 kcal/mol(-5.0 kJ/mol)时,表明两分子结合良好
[25]。分子对接结果(
图6)(
表3)表明山柰酚与关键靶点均能自发结合,其中与ESR1、AR、Bcl2具有强烈的结合活性,其结合能分别为-9.1 kcal/mol、-8.5 kcal/mol、-7.6 kcal/mol,与ESR1结合效果最好。其结果说明山柰酚能够通过作用于乳腺癌靶点发挥乳腺癌防治作用。
2.5 LASSO-Cox回归预测模型的构建
整合TCGA数据库中每例乳腺癌患者的生存时间、生存状态及靶基因在每例患者样本中的表达,将其导入R,在R语言中进行LASSO-Cox预后回归分析对55个靶基因进行进一步筛选。结果显示,随着λ值的变化,系数越早被压缩为0,靶基因对乳腺癌的影响越小(
图7A)。交叉验证图显示LASSO回归筛选最佳变量数的过程(
图7B)。当λ=0.014 8时误差最小,以λ=0.014 8进行变量筛选获得10个预后基因:Bcl2、CYP1B1、DPP4、GSTM1、GSTM2、MMP1、NCOA2、NOS2、NR1I3、PTGS2。将预后基因中各基因回归系数带入风险评分公式:risk score=(-0.058 9)Bcl2+(-0.035 9)PTGS2+(-0.006 8)GSTM1+(0.011 5)MMP1+(-0.049 9)CYP1B1+(0.222 1)NOS2+(0.164 2)NCOA2+(0.022 9)DPP4+(-0.039 7)GSTM2+(0.060 8)NR1I3。根据公式计算出每例患者的风险评分后,取其平均值,在此基础上将患者分为高危组和低危组,并绘制生存曲线(
图8),采用Log-rank检验可得出高危组与低危组的预后差异有统计学意义(
P<0.001),其中高危组和低危组的中位生存时间分别是9.2年和17.2年。为了验证预后回归模型的准确性绘制了ROC曲线(
图9),AUC值>0.5。
为了分析不同风险评分与患者随访时间、生存状态及基因表达之间的关系,绘制了风险曲线、风险和生存散点图(
图10)和预后基因的表达热图(
图11)。随着风险评分的升高,患者OS率明显降低。预后基因的表达热图提示,随着Bcl2、CYP1B1、GSTM1、GSTM2、PTGS2表达增加,患者的风险评分逐渐降低,随着DPP4、MMP1、NCOA2、NOS2、NR1I3的表达增加,患者的风险评分逐渐升高。
2.6 单基因预后分析
通过LASSO-Cox分析获得了10个预后基因,为了评估单基因对乳腺癌患者生存状态及生存时间的影响,利用TCGA数据库中乳腺癌患者生存数据,分别绘制了各基因处于高表达和低表达状态下的生存曲线(
图12)。单基因生存曲线提示Bcl2(
HR=0.61,95%
CI=0.43~0.86,
P=0.005)、CYP1B1(
HR=0.68,95%
CI=0.49~0.94,
P=0.022)、GSTM1(
HR=0.68,95%
CI=0.47~0.98,
P=0.037)、GSTM2(
HR=0.64,95%
CI=0.46~0.90,
P=0.010)、PTGS2(
HR=0.62,95%
CI=0.44~0.86,
P=0.005)高表达组的OS期均高于低表达组。而MMP1(
HR=1.72,95%
CI=1.22~2.41,
P=0.002)、NCOA2(
HR=1.71,95%
CI=1.12~2.60,
P=0.013)、NOS2(
HR=1.67,95%
CI=1.20~2.32,
P=0.002)、NR1I3(
HR=1.69,95%
CI=1.21~2.37,
P=0.002)高表达组的OS均低于低表达组。DPP4(
P=0.087)暂不认为是影响乳腺癌预后的因素。该结果提示Bcl2、CYP1B1、GSTM1、GSTM2、PTGS2可能是乳腺癌预后的保护因素,而MMP1、NCOA2、NOS2、NR1I3可能是乳腺癌预后的危险因素。
2.7 预后模型与肿瘤免疫的相关性
为了进一步探讨模型中的基因是否影响免疫反应,利用从TCGA数据库中下载的乳腺癌患者样本RNAseq数据,进行风险模型评分与免疫评分之间的Spearman相关性分析(
图13)。由图可知,风险评分与T细胞、CD8
+ T细胞、髓样树突细胞、NK细胞、B细胞的浸润呈负相关,与单核/巨噬细胞的浸润呈正相关(均
P<0.05)。
3 讨 论
山柰酚作为一种天然黄酮类化合物,可以调节各种与癌症相关的过程和活动,如细胞周期、氧化应激、细胞凋亡、增殖、转移和血管生成
[26]。通过抑制乳腺癌细胞的生长、诱导乳腺癌细胞的凋亡及抑制乳腺癌细胞的迁移和侵袭等多种途径
[27]发挥作用,为乳腺癌的治疗提供了新思路。
网络药理学是基于多向系统药理学和生物学提出的一种新的药物研究方法,近年来已成功用于疾病-基因靶向药物网络的构建和可视化
[28]。GO分析描述了基因的BP、CC及MF,KEGG富集分析描述了药物治疗机制研究中显著改变的信号通路
[29]。本研究利用网络药理学方法,获得山柰酚与乳腺癌的55个交集靶点,并进行PPI分析,得到度值较大的靶点基因:Akt1、Bcl2、CASP3、ESR1、AR。Qin等
[30]发现山柰酚可通过抑制PI3K/Akt通路抑制胆管癌的生长和转移。El-Kott等
[31]发现山柰酚通过抑制PI3K/Akt信号通路增加卵巢癌细胞对化疗的敏感性。Da等
[32]发现山柰酚以剂量依赖性方式抑制二氢睾酮诱导的AR表达和核积累从而抑制激素敏感性前列腺癌细胞生长。Yi等
[33]发现山柰酚可能通过下调抗凋亡蛋白Bcl2表达并诱导促凋亡蛋白Bax表达抑制MCF-7乳腺癌细胞的增殖。为了更深入地了解山柰酚治疗乳腺癌的生物学功能和信号通路,本研究进行了GO和KEGG富集分析。GO分析表明,与该治疗相关的生物学功能主要包括对外来刺激的反应、对肽的反应、对氧化应激的反应、细胞对化学应激的反应等。涉及的CC包括膜微区、膜筏、转录调节复合物、转移酶复合物,转移含磷基团等。涉及的MF主要包括DNA结合转录因子结合、RNA聚合酶Ⅱ特异性DNA结合转录因子结合等。KEGG通路富集分析表明,山柰酚治疗乳腺癌可靶向多条通路,主要包括化学致癌-活性氧信号通路、糖尿病并发症中的AGE-RAGE信号通路、TNF信号通路、IL-17信号通路等。Calderón-Montaño等
[34]研究发现山柰酚在高浓度下增加抗氧化酶的表达,在低浓度下清除羟基(OH)自由基及过氧亚硝酸自由基。Zeng等
[35]研究发现山柰酚可以抑制小鼠骨髓来源的中性粒细胞中活性氧(reactive oxygen species,ROS)的产生,并通过ROS-PAD4通路减少乳腺癌的转移。基于这些结果,将山柰酚与关键靶点Akt1、Bcl2、CASP3、ESR1、AR分别进行分子对接验证。分子对接结果显示山柰酚与关键靶点均能自发结合,其中与ESR1结合效果最好,提示山柰酚与关键靶点之间具有强烈的结合活性,进一步验证了山柰酚能通过作用于关键靶点发挥治疗乳腺癌作用。
此外,为探索山柰酚治疗对乳腺癌预后的影响,应用LASSO-Cox回归分析构建了具有风险评分的LASSO-Cox预后回归模型。通过该模型获得患者的风险评分和预后生存率,发现高风险组的OS率明显低于低风险组,该模型可用于进行临床治疗策略的调整。同时分别对10个预后基因进行预后分析,结果提示Bcl2、CYP1B1、GSTM1、GSTM2、PTGS2可能是乳腺癌预后的保护因素,而MMP1、NCOA2、NOS2、NR1I3可能是乳腺癌预后的危险因素,暂不认为DPP4是乳腺癌预后的影响因素。通过Spearman相关性分析发现风险评分与T细胞、CD8+ T细胞、髓样树突细胞、NK细胞、B细胞的表达呈负相关关系,与单核/巨噬细胞的表达呈正相关关系,差异均具有统计学意义。该模型可用于指导临床的免疫治疗。
综上所述,山柰酚可通过多靶点、多通路发挥治疗乳腺癌作用。根据LASSO-Cox预后回归模型获得了10个预后相关基因,并探索了预后模型和免疫的相关性,这些研究有助于指导临床治疗并为下一步实验提供理论依据。
但本研究存在一定的局限性:本研究仅在数据库对山柰酚治疗乳腺癌的机制及预后进行了研究,样本量的限制、数据库选择及潜在模型误差等均可对研究结果产生影响,且缺乏进一步的体内外实验验证。其次,未经验证或记录的靶标可能未包含在本研究中。此外,通过网络药理学方法鉴定的靶基因在复杂的人体系统中的影响尚未得到充分的探索及解释。在后期研究中,笔者将及时更新数据资料、使用多种分析方法或多个数据库进行验证,减小研究结果误差,并设计进一步体内外试验对山柰酚治疗乳腺癌的机制进行探究。
湖南省自然科学基金资助项目(2022JJ30335)