基于对枸杞、姜黄、甘草的生物信息学分析探究GRIN3A在前列腺癌中的治疗意义

王庆哲，张恩崇，宋永胜

0 引言

前列腺癌是发生于前列腺的上皮来源的恶性肿瘤，在男性恶性肿瘤中发病率最高[1]。前列腺癌具有高度异质性，表现形式多种多样，包括低度恶性、局部侵袭性以及远端转移性。虽然该疾病在绝大多数情况下表现为低度恶性，但是对患者的长期健康具有不可忽视的危害[2-4]。对于局部前列腺癌，主要的治疗手段有主动监测、手术以及放化疗。为了减少复发，术后患者将接受雄激素去势治疗 (Androgen deprivation therapy,ADT)或者放疗[5]。对于发生转移的前列腺癌患者，早期使用化疗之后进行ADT治疗是主要的临床选择[6]。然而，当患者接受ADT治疗12个月之后病情往往会发展为去势抵抗性前列腺癌 (Castration-resistant prostate cancer，CRPC)[7]。目前，阿比特龙和恩杂鲁胺是全球公认用来治疗CRPC的药物，但是治疗效果并不理想[8]。因此，针对目前前列腺癌治疗过程中面临的挑战，探索新型的治疗靶点具有十分重要的意义[9]。

枸杞可以抑制PC-3细胞在裸鼠体内的生长，促进PC-3、DU-145细胞凋亡[10]。甘草可诱导LNCaP细胞中的细胞凋亡蛋白酶依赖和自噬相关细胞的凋亡[11]。姜黄素是姜黄的一种天然成分，可以通过抑制H3K4me3以及JNK通路从而诱导LNCaP细胞的凋亡[12]。

本研究使用BATMAN-TCM (a Bioinformatics Analysis Tool for Molecular mechANism of Traditional Chinese Medicine) 数据库分别预测了枸杞、甘草、姜黄3种药材所含化合物作用的靶蛋白[13]。从TCGA (The Cancer Genome Atlas Program) 数据库下载前列腺癌患者的RNA-seq数据，筛选出肿瘤和癌旁组织间的差异表达基因(Differentially expressed genes, DEGs)。然后对之前3类药材都作用的靶蛋白基因和DEGs取交集，作为研究的候选核心基因。根据作用化合物的多少在Cytoscape 3.7.1中筛选出核心基因[14]。随后对核心基因进行生存分析，发现GRIN3A基因高表达提示患者的无病生存时间(Disease-free survival，DFS)更短。之后通过单基因GSEA (Gene Set Enrichment Analysis) 分析发现，当GRIN3A高表达时，PI3K-Akt信号通路被激活[15-16]。该研究一方面从分子机制方面再次证实了枸杞、姜黄、甘草对前列腺癌治疗的有效性，另一方面，也为前列腺癌的精准靶向治疗提供了一个新的研究靶点。

1 材料和方法

1.1 预测枸杞、甘草、姜黄的共同作用靶向蛋白 BATMAN-TCM数据库是第1个研究传统中药分子机制的在线生物信息学分析工具。该数据库可以预测中药的潜在靶向蛋白，并对预测到的靶蛋白基因进行Gene Ontology (GO) 分析以及 KEGG 通路功能富集分析[13]。在分别查询枸杞、姜黄、甘草3类药材的靶基因后，我们选取三者共有的靶基因进行后续核心基因的研究。

1.2 通过TCGA数据库分析差异表达基因 TCGA数据库 (https://portal.gdc.cancer.gov/)是一个癌症基因组学项目，从基因分子水平上展示了33种癌症的2万多例原发性肿瘤及其配对癌旁组织样本[13]。笔者从TCGA数据库下载了前列腺癌患者的RNA-seq数据，然后在R 3.5.2 中使用DESeq2 R包处理RNA-seq数据，获得肿瘤样本和正常癌旁组织间的差异表达基因(DEGs)[17]。

1.3 对DEGs进行功能富集分析 为了更好地理解前列腺癌的生物学行为，我们对前列腺癌中分析获得的DEGs进行了GO分析和KEGG通路富集分析。作为21世纪以来具有代表性的生物信息学分析工具，GO分析和KEGG通路富集分析通过大量的基因集合锁定和该基因集相关的生物学过程[18]。GO分析包括细胞组件 (Cellular component,CC)、生物过程 (Biological process,BP)、分子功能 (Molecular function,MF)。KEGG分析主要从信号通路方面解释问题，通过R 3.5.2 中的 clusterProfiler 软件包实现[19]。我们选择GO分析以及KEGG分析中排名前10位的项目。

1.4 通过Cytoscape筛选候选核心基因 既出现在DEGs又是枸杞、姜黄、甘草共同靶点的基因具有深度探索的生物学意义，因此，笔者对DEGs集合和3种药材共同靶基因集合选取交集，将这类基因称作候选核心基因。将这些候选核心基因及通过BATMAN-TCM数据库预测作用它们的化合物输入到 Cytoscape 3.7.1 软件中[14]。通过该软件中的 CentiScaPe 2.2 插件计算出每一个基因的中心度，在网络分析研究中，具有较高中心度的节点被认为在网络中具有重要意义。因此，我们选择中心度排名前3位的基因作为候选核心基因。

1.5 利用GEPIA(Gene Expression Profiling Interactive Analysis)数据库(http://gepia.cancer-pku.cn/)筛选核心基因 GEPIA数据库是一个交互在线网页，收集了TCGA项目和GTEx项目的基因组数据并利用广泛认可的生物信息学方法分析数据。我们将候选基因输入到GEPIA数据库中，进行生存分析，检验这些基因对患者的无病生存期 (Disease-free survival，DFS)是否具有提示作用，从而确定要研究的核心基因。

1.6 单基因GSEA分析确定核心基因相关的信号通路 为进一步探索核心基因的生物学行为，我们通过单基因GSEA分析探索核心基因相关的信号通路。使用R 3.5.2 中的 clusterProfiler 软件包进行GSEA分析[19]，并用R 3.5.2 的DESeq2 软件包获取每个基因在肿瘤组织和正常癌旁组织间的差异倍数的对数 (Logarithmic fold change,LFC)，将每个基因的LFC作为GSEA分析过程中的排序依据。从Molecular Signatures Database v7.0 (http://software.broadin stitute.org/gsea/msigdb/index.jsp) 中选择C2 (curated gene sets) 作为分析用的基因集[15,20]。

1.7 COX回归检验核心基因能否成为前列腺癌患者的独立预后因素 为了探究核心基因能否成为前列腺癌患者的独立预后因素，首先对纳入的各项临床因素及核心基因的表达量进行单因素COX回归，然后选择在单基因COX回归中具有统计学意义的指标进行多因素COX回归，从而确定要研究的核心基因的生物学意义。

1.8 数据统计分析 在使用靶蛋白基因对枸杞、甘草、姜黄进行GO分析以及 KEGG 通路功能富集分析时，设定P<0.05为差异有统计学意义。在筛选DEGs过程中，设定P<0.05、|LFC|>2为筛选标准。对DEGs进行GO分析和KEGG通路富集分析时，设定P<0.05 为差异有统计学意义。生存分析、单因素COX回归、多因素COX回归以及单基因GSEA分析中，均设定P<0.05为差异有统计学意义。

2 结果

2.1 BATMAN-TCM数据库对枸杞、姜黄、甘草的分析结果 本研究分别对枸杞、甘草、姜黄在BATMAN-TCM数据库进行靶蛋白的预测，发现枸杞含有90种化合物，956个靶蛋白；姜黄含有55种化合物，677个靶蛋白；甘草含有172种化合物，685个靶蛋白。3种药材共有的靶蛋白为180个，见图1。3种药物的GO分析以及KEGG通路富集分析结果见图2。从整体来看，结果表明，这3类药物与细胞的增殖以及死亡相关。

图1 甘草、枸杞和姜黄靶向蛋白

2.2 TCGA数据库筛选DEGs 从TCGA数据库下载了539个样本的RNA-seq数据，其中肿瘤样本488个，癌旁样本51个。提取出样本中19 466个蛋白编码基因的表达数据。经过R 3.5.2 的DESeq2 R包处理后，得到649个DEGs。

2.3 DEGs的功能富集分析 对筛选得到的649个DEGs进行GO分析和KEGG通路富集分析，选取每项排名前10位的项目，见图2B。

GO分析中：BP方面提示DEGs主要与内肽酶活性的负调控、甲状腺激素代谢过程、肽酶活性的负调控、甲状腺激素生成、调节甲状腺激素的产生等方面相关；CC方面说明，DEGs主要分布在乳糜微粒、质膜顶端、高密度脂蛋白颗粒、血浆脂蛋白颗粒、血液微粒、蛋白脂质复合体、极低密度脂蛋白颗粒、富含三酰甘油血浆脂蛋白颗粒、细胞顶端；MF方面揭示，DEGs主要与肽链内切酶抑制剂活性、肽酶抑制剂活性、肽链内切酶调节物活性、肽酶调节活性、近端启动子DNA结合转录激活子活性、酶抑制剂活性、半胱氨酸型内肽酶抑制剂的活性、钾离子跨膜转运活性、转录因子活性相关。综上所述，通过GO分析，可以推测前列腺癌的发生过程主要与肽酶活性的抑制、内分泌系统的激活以及DNA转录的调控相关，与前列腺癌相关的DEGs主要分布在血液脂蛋白颗粒中。KEGG通路富集分析表明，DEGs主要富集在药物代谢通路(细胞色素P450及其他酶)、甾体类激素生物合成、花生四烯酸代谢、视黄醇代谢、戊糖、葡萄糖醛酸转换、谷胱甘肽代谢、抗坏血酸和醛酸代谢、卟啉和叶绿素代谢等细胞通路。

2.4 通过化合物及作用靶蛋白网络筛选候选核心基因 对3种药物共同作用的靶蛋白基因和DEGs取交集，获得了11个基因。这些基因既是3种药物共同作用的靶基因，同时，其在肿瘤组织和癌旁组织之间的表达又存在显著性差异，认为这11个基因有继续挖掘的价值。将这11个基因(PTGS2,GRIN3A,CACNA1D,PTGS1,SLC18A3,SCN11A,CYP19A1,CRP,CHRM2,CHRNA4,GATA3)及作用它们的化合物输入到Cytoscape 3.7.1 中构建分析网络，然后得到每个基因在网络中的中心度，结果见表1。根据结果，选择Degree排名前3位的基因(PTGS2,GRIN3A,CACNA1D)作为候选核心基因。

表1 候选核心基因在相互作用网络中的中心度

图2 对于DEGs行GO分析和KEGG通路富集分析

2.5 生存分析筛选核心基因 在GEPIA数据库中，对PTGS2、GRIN3A、CACNA1D3个候选核心基因进行生存分析。结果见图3。当GRIN3A基因高表达时，前列腺癌患者的无病生存时间显著缩短，提示GRIN3A基因对于患者不良预后具有促进作用。结果显示,其他基因与患者生存结局的关系并无统计学意义，因此，选定GRIN3A为本研究的核心基因。

2.6 对GRIN3A基因进行单基因GSEA分析 为了探索核心基因GRIN3A相关的生物学事件，使用来自TCGA数据库539个前列腺癌样本的RNA-seq数据，其中肿瘤样本488个，癌旁样本51个。根据GRIN3A的表达量水平，选取表达量前1/4的样本作为高表达组，表达量后1/4的样本作为低表达组。然后利用R 3.5.2 中的DESeq2 软件包分析了19 466个蛋白编码基因在高表达组与低表达组之间差异表达的LFC。LFC作为单基因GSEA分析中基因列表排列顺序，从Molecular Signatures Database v7.0 (http://software.broadin stitute.org/gsea/msigdb/index.jsp) 中选择C2 (curated gene sets) 作为分析用的基因集。见图4。结果显示，脂肪酸代谢、钙离子信号通路、肿瘤中的信号通路在高表达组中激活；RNA降解、凋亡信号通路被抑制。该结果进一步证实了GRIN3A在前列腺癌发生发展过程中的重要性。

2.7 COX回归检验GRIN3A基因是否可以成为患者预后的独立预测指标 为了进一步探索GRIN3A基因的生物学意义，我们获取了TCGA数据库485位前列腺癌患者的临床病理数据，见表2。在单因素COX回归中，Gleason评分(HR=2.12，95%CI：1.02～4.41，P=0.045)、PSA(HR=1.04，95%CI:1.02～1.06，P<0.001)、病理肿瘤大小分级(HR=0.17，95%CI:0.03～0.89，P=0.036)、GRIN3A(HR=4.90，95%CI：1.41～31.96，P=0.037)这4项指标对患者的生存时间具有统计学相关性。其中Gleason评分、PSA水平、GRIN3A值越高提示患者生存预后越差，病理肿瘤大小分级值越高提示患者预后越佳。为了探索这些指标能否成为前列腺癌患者预后的独立预测因素，我们进行了多因素COX回归分析，结果显示，PSA水平(HR=1.02，95%CI：1.00～1.05，P=0.041)以及GRIN3A(HR=2.13，95%CI：0.26～17.66，P=0.048)可以成为患者预后的独立预测因素。上述单因素及多因素COX回归结果见表3。综上所述，PSA水平及患者GRIN3A基因的表达量可以是患者预后的危险因素，且具有统计学意义。

3 讨论

中药在现代医学的治疗中具有重要作用，但是缺乏明确的分子机制一直是中药在治疗中存在的问题。笔者首先通过查阅文献发现，枸杞、姜黄、甘草3类药材在前列腺癌相关研究中具有明显的抑癌作用。因此，本研究通过BATMAN-TCM数据库分别预测了枸杞、甘草、姜黄3种药材所含的化合物及其作用的靶蛋白。

表2 本研究中患者的临床病理指标

本研究筛选得到PTGS2、GRIN3A、CACNA1D、PTGS1、SLC18A3、SCN11A、CYP19A1、CRP、CHRM2、CHRNA4、GATA3这11个基因。将这11个基因及作用它们的化合物输入到Cytoscape 3.7.1 中构建分析网络，根据结果，选择Degree排名前3位的基因(PTGS2,GRIN3A,CACNA1D)作为候选核心基因。PTGS2在直肠腺癌组织中的表达量是正常组织的8～9倍[21]。CACNA1D基因突变可导致肿瘤细胞的形成[22]。在GEPIA数据库中，对PTGS2、GRIN3A、CACNA1D3个候选核心基因进行生存分析，当GRIN3A基因高表达时，前列腺癌患者的无病生存时间显著下降，提示GRIN3A基因对于患者不良预后具有促进作用。结果显示，其他基因和患者的生存结局并无统计学关系。因此选定GRIN3A为本研究的核心基因。

为了探索核心基因GRIN3A相关的生物学事件，本研究进行了单基因GSEA分析，发现脂肪酸代谢、钙离子信号通路、肿瘤中的信号通路在高表达组中激活；RNA降解、凋亡信号通路被抑制。该结果进一步证实了GRIN3A在前列腺癌发生发展过程中的重要性。

图3 无病生存时间生存分析

图4 GRIN3A单基因GSEA分析

表3 临床病理指标及GRIN3A表达量的单因素、多因素COX回归结果

临床指标单因素COX回归HR95%CIP值多因素Cox回归HR95%CIP值年龄0.990.88～1.100.8201.020.89～1.160.760Gleason评分2.121.02～4.410.0451.710.77～3.780.188PSA1.041.02～1.06<0.0011.021.00～1.050.041病理淋巴结转移分级3.810.85～17.030.0813.210.48～21.370.227病理肿瘤大小分级0.170.03～0.890.0360.110.01～1.010.051肿瘤质量1.000.99～1.000.4101.000.99～1.010.527GRIN3A表达量4.902.21～31.960.0372.130.26～17.660.048

综上所述，通过对枸杞、姜黄、甘草进行一系列的生物信息学分析，三者共同作用的核心基因GRIN3A可作为预测前列腺癌患者不良预后的独立预测因素，并且在这3类药材对前列腺癌细胞的治疗作用中发挥很大作用。GRIN3A参与激活脂肪酸代谢、钙离子信号通路、肿瘤中的信号通路抑制RNA降解、凋亡信号通路。