肺腺癌EMT相关基因筛选及其ceRNA网络构建*

马楚涵,孙家琰,赵雨馨,赵幼敏,刘奕墨,赵子涵,王子剑,尹智华△

中国医科大学:1.第一临床学院;2.第四临床学院;3.公共卫生学院;4.健康管理学院,辽宁沈阳 110122

肺癌是我国发病率与死亡率最高的恶性肿瘤[1]。腺癌是肺癌最常见的病理类型[2]。中晚期腺癌患者5年生存率较低,其原因主要是肿瘤细胞的淋巴及血道转移导致术后复发。因此,对于肺腺癌侵袭、转移的研究具有重要的临床意义。上皮间质转化(EMT)是指上皮来源的细胞向间充质细胞转化的过程,上皮细胞通过EMT失去极性及细胞连接,获得富有运动能力的间充质表型,从而突破基膜进一步侵袭转移。因此EMT与肿瘤转移密切相关[3]。本研究利用基因表达图谱(GEO)及肿瘤基因组图谱(TCGA)数据库,筛选肺腺癌EMT相关基因,并对其进行了生物信息学分析及关键基因竞争性内源RNA(ceRNA)网络构建,从基因组层面系统描述了EMT过程,为相关研究提供了新思路。

1 材料与方法

1.1标本来源 在GEO数据库(http://www.ncbi.nlm.nih.gov/geo/)中下载GSE10072数据集进行分析,该数据集包含58例肺腺癌标本及49例正常组织标本,采用Affymetrix Human Genome U133A Array芯片检测全基因组基因。在TCGA数据库中下载肺腺癌项目数据,包括516例肺腺癌标本、56例正常组织标本的基因表达数据,使用R语言进行后续分析。本研究所选标本均来源于44岁以上患者,从经组织病理学证实的肺腺癌Ⅰ～Ⅳ期患者术中切除所得,未区分性别与是否吸烟。

1.2方法

1.2.1差异表达分析 本研究使用GEO2R筛选GEO数据库中的差异表达基因;构建TCGA数据库基因表达矩阵,采用R语言中的limma包进行差异分析(P<0.05且|logFC|>1),取两次分析结果的交集作为最终的差异表达基因。

1.2.2筛选肺腺癌EMT相关基因 将差异表达基因的Gene Symbol导入GenClip3(http://cismu.net/genclip3/analysis.php)。软件将从Pubmed中检索相关文献,并基于关键词对于导入的基因进行聚类分析[4]。

1.2.3基因本体论(GO)及京都基因和基因组百科全书(KEGG)富集分析 采用Metascape(http://metascape.org/)对EMT相关基因进行功能注释及GO与KEGG富集分析[5]。

1.2.4构建蛋白互作网络(PPI)及筛选关键基因 采用STRING在线工具(https://cn.string-db.org/)分析EMT相关基因编码的蛋白质之间的相互作用[6],得到PPI,将其中结合分数值>0.4的蛋白相互作用对导入Cytoscape软件中,去除离散蛋白,通过cytoNCA分析蛋白之间的相关程度。根据degree排序,选择top 10的基因作为关键基因,并针对其绘制PPI网络图。

1.2.5预测EMT关键基因潜在治疗中药 将关键基因与肺腺癌、EMT共同导入Coremine medical数据库(http://www.coremine.com/),以P<0.05为标准筛选治疗性中药。将针对10个关键基因的中药取交集,并绘制中药与关键基因相关性的气泡图。

1.2.6Kaplan-Meier分析 通过Kaplan-Meier(www.kmplot.com)对GEO、欧洲基因组表型档案(EGA)及TCGA中的患者数据进行分析,以50个月的生存率评估患者预后[7]。最后,选取Log-rankP<0.05的基因作进一步研究。

1.2.7根据ceRNA假说预测关键基因上游微小RNA(miRNA)及miRNA上游的长链非编码RNA(lncRNA) 根据ceRNA假说,lncRNA能够竞争性抑制miRNA对靶基因的调控作用,因此,靶基因的表达水平与miRNA呈负相关而与lncRNA呈正相关。使用miRTarbase(https:// mirtarbase.cuhk.edu.cn/)预测关键基因上游miRNA,为提高研究结果可信度,只选择研究报道经蛋白质印迹法、实时荧光定量PCR证实的miRNA纳入研究[8]。使用ENCORI(http://starbase.sysu.edu.cn/)计算logFC,评估miRNA及关键基因在正常组织与肿瘤组织中的差异表达情况,并绘制生存曲线,以Log-rankP<0.05判定基因表达与患者预后的关系。选择与患者预后相关并与靶基因在肿瘤组织与正常组织间差异表达趋势相反的miRNA进行下一步研究。使用miRNet database(https://www.mirnet.ca/)预测miRNA上游lncRNA[9],通过ENCORI选取与患者预后相关并与靶基因在肿瘤组织与正常组织间差异表达趋势相同的lncRNA作进一步研究。

1.2.8构建ceRNA网络 根据上一步预测的miRNA及lncRNA,使用ENCORI绘制散点图,通过Pearson相关性检验判断基因表达的相关性,并绘制关键基因的ceRNA调控网络图[10]。

1.3统计学处理 本研究采用R软件(4.1.0版)和Microsoft Excel 2019对TCGA数据进行整理。应用R软件(4.1.2版)对RNA表达数据进行差异分析(P<0.05且|logFC|>1),并将其与GEO2R所得差异表达基因取交集获得最终结果。然后,使用GenClip3软件进行聚类分析,并采用Metascape对目标基因进行GO和KEGG富集分析。同时,采用STRING在线工具预测蛋白质之间的相互作用,通过Cytoscape软件选择连接度最高的10个基因作为关键基因。在此基础上使用Coremine medical数据库筛选治疗性中药(P<0.05),并且结合生存数据进行Kaplan-Meier(Log-rankP<0.05)分析。最后,使用在线预测网站构建ceRNA,通过Pearson相关性检验判断基因表达的相关性。

2 结 果

2.1差异表达分析 使用GEO2R在线分析工具分析,得到差异表达基因662个;在TCGA数据库中使用limma包进行差异表达分析,得到差异表达基因8 295个。取交集得到差异表达基因565个,见图1。

图1 GEO和TCGA数据库差异表达基因的韦恩图

2.2EMT相关基因筛选 将565个差异表达基因导入GenClip3中,结果显示所有差异表达基因均可在文献中检索到,使用GenClip3对差异表达基因进行聚类分析,得到156个EMT相关基因。

2.3GO、KEGG富集分析 GO分析发现,这些基因主要与血管发育、激素应答、生长因子刺激的细胞应答、细胞迁移的正向调节、细胞外基质等生物过程、细胞组件、分子功能相关,见图2。KEGG通路分析发现,这些基因主要与癌通路、糖尿病并发症中的高级糖基化终末产物-受体(AGE-RAGE)通路、缺氧诱导因子-1(HIF-1)信号通路、肿瘤中的miRNA通路、细胞黏附分子通路等相关,见图3。

图2 EMT相关基因的GO富集图

图3 EMT相关基因的KEGG富集图

2.4PPI构建 根据degree排序,得到的10个关键基因,分别为钙黏蛋白1(CDH1)、白细胞介素-6(IL-6)、基质金属蛋白酶9(MMP9)、血小板内皮细胞黏附分子(PECAM1)、细胞周期蛋白依赖性激酶抑制剂2A(CDKN2A)、α1-Ⅰ型胶原基因(COL1A1)、分泌型磷酸蛋白1(SPP1)、TIMP基质金属蛋白酶抑制因子1(TIMP1)、小窝蛋白1(CAV1)、zeste 基因增强子同源物(EZH2),见表1。从PPI网络图中可以看出,关键基因编码的蛋白之间相互作用明显,见图4。

表1 关键基因及degree

图4 关键基因的PPI

2.5关键基因潜在治疗性中药 10个关键基因预测到紫草、温莪术、片姜黄等多味中药,见表2。将10个关键基因的潜在治疗性中药取交集,得到丹参、人参芦、人参、人参叶、人参花共五味中药,以P值为参考,使用R软件的ggpubr包绘制关键基因与中药相关性的气泡图,见图5。

表2 针对关键基因的部分治疗性中药

图5 中药与关键基因的相关性气泡图

2.6关键基因潜在治疗性中药的Kaplan-Meier分析 结果显示,10个关键基因均与患者预后相关(Log-rankP<0.05),其中高表达CDH1、PECAM1、CAV1的生存曲线与对照的风险比(HR)<1,表明这些基因的表达可以降低患者死亡风险;高表达IL-6、MMP9、CDKN2A、COL1A1、SPP1、TIMP1的生存曲线与对照的HR>1,表明这些基因的表达增加患者的死亡风险。见图6。

图6 关键基因的Kaplan-Meier分析结果森林图

2.7根据ceRNA假说预测关键基因上游miRNA及miRNA上游lncRNA 差异分析结果显示,IL-6(logFC=-2.64)在肿瘤组织中呈低表达,而MMP9(logFC=2.14)、COL1A1(logFC=3.05)、SPP1(logFC=4.95)、EZH2(logFC=2.70)在肿瘤组织中呈高表达,同理可知IL-6的上游miRNA hsa-miR-9-5p,MMP9的上游miRNA hsa-miR-211-5p,COL1A1的上游miRNA hsa-miR-133a-3p,SSP1的上游miRNA hsa-miR-299-5p,EZH2的上游miRNA hsa-miR-101-3p、hsa-miR-30d-5p与其靶基因的差异表达趋势相反,符合ceRNA假说;并且上述miRNA的Log-rankP<0.05,与患者预后相关,因此将其纳入下一步研究。由logFC可以看出,hsa-miR-9-5p的上游lncRNA CCDC13-AS1、LINC00996、SH3BP5-AS1、LINC00921,hsa-miR-211-5p的上游lncRNA EBLN3P、DHRS4-AS1、FAM30A、ZNF674-AS1,hsa-miR-133a-3p的上游lncRNA LINC00847、LINC02535、MIR4697HG,hsa-miR-299-5p的上游lncRNA GUSBP11,hsa-miR-101-3p的上游lncRNA LINC01579、EBLN3P、GSEC,hsa-miR-30d-5p的上游lncRNA HCG18、LINC02535与靶基因的差异表达趋势相同,而与上游miRNA差异表达趋势相反,符合ceRNA假说;并且上述lncRNA的Log-rankP<0.05,与患者预后相关,因此纳入下一步研究。见图7～9。

图7 mRNA-miRNA-lncRNA相互作用的桑基图

图8 基因差异表达分析的偏差图

图9 基因HR的棒棒糖图

2.8构建ceRNA网络 在预测的miRNA及lncRNA中,同时满足mRNA的表达水平与miRNA呈负相关而与lncRNA呈正相关,lncRNA的表达水平与miRNA呈负相关的基因被纳入ceRNA网络。结果表明,EZH2/hsa-miR-101-3p/GSEC符合ceRNA假说,见图10～16。

图10 EZH2在肿瘤组织和正常组织中表达水平差异的箱线图

图11 EZH2的Kaplan-Meier分析曲线

图12 hsa-miR-101-3p在肿瘤组织和正常组织中表达水平差异的箱线图

图13 hsa-miR-101-3p的生存曲线

图14 GSEC在肿瘤组织和正常组织中表达水平差异的箱线图

图15 GSEC的生存曲线

图16 ceRNA网络及各组分之间表达水平的相关性

3 讨 论

本研究获得了565个肺腺癌差异表达基因及156个EMT相关基因,并对其进行了功能富集分析。随后绘制了EMT相关基因的PPI网络图,并根据degree获得10个EMT关键基因。其中,CDH1编码E-钙黏蛋白(E-cadherin),E-cadherin下调被广泛报道为EMT的标志,其表达下降与肿瘤患者的不良预后密切相关,在正常细胞中,E-cadherin可以阻止β-catenin,与参与wnt信号通路的淋巴样增强因子/T细胞因子结合发挥肿瘤抑制作用[11]。IL-6可以通过Janus激酶/信号转导子和转录激活子信号通路诱导EMT[12],但也有研究显示IL-6可以动员免疫反应抑制肿瘤生长。MMP9通过参与细胞外基质分解促进肿瘤侵袭转移,有研究显示STAT3过表达时促进了EMT过程,同时MMP9的表达水平也有所升高,提示MMP9在EMT中发挥一定的作用[13]。PECAM1又名CD31,其编码的蛋白质存在于单核细胞、血小板、中性粒细胞和某些类型的T细胞表面,并构成细胞连接,CD31可以促进EMT从而促进结肠癌腹膜转移[14],也有研究将CD31作为血管生成的标志物[15],说明CD31也可能参与肿瘤血道转移。CDKN2A是一种抑癌基因,但在结肠癌中CDKN2A的表达增加促进EMT及肿瘤转移[16],并且可以与趋化素样因子超家族8、整合素连接激酶等通过转化生长因子β通路、Notch通路等形成免疫抑制微环境,从而影响患者的预后[17]。COL1A1是Ⅰ型胶原蛋白的主要成分,通过结合integrin的β1亚基激活局部黏着斑激酶(FAK)。磷酸化的FAK促进E-cadherin/catenin复合物的分解,从而激活EMT[18]。SPP1上调在前列腺癌、结直肠癌中都有激活EMT,促进肿瘤转移的作用,可能与SPP1激活细胞外信号相关激酶1/2信号通路有关[19-20]。TIMP1为MMP1的抑制剂,但其激活EMT的作用是通过与CD63结合诱导TWIST1的表达,与MMP抑制功能无关。CAV1可以通过Wnt/β-catenin通路诱导EMT过程,促进肿瘤转移。EZH2促进绝大多数肿瘤的侵袭转移,但有时也发挥抗肿瘤作用,这种矛盾现象与EMT密切相关[21]。对10个关键基因的Kaplan-Meier分析显示,IL-6、MMP9、CDKN2A、COL1A1、SPP1、TIMP1的表达与患者预后不良有关,而上述研究显示这些基因对于EMT均起促进作用,可能是其影响患者生存期的原因之一。并使用Coremine medical数据库预测了10个关键基因的治疗性中药,人参、人参叶、人参花、人参芦、丹参这五味中药对于所有关键基因均有一定效果,本研究还根据P值绘制了这五味中药与关键基因的相关性气泡图。从气泡图中分析,人参、人参叶、人参花、人参芦与EMT的关系较丹参更为密切,而人参、人参叶、人参花、人参芦来源于同一种植物,因此人参中某种成分可能有效抑制EMT从而延缓肿瘤进展。有研究表明,人参皂苷可以抑制EMT从而抑制肿瘤生长与转移[22],可能为其潜在机制之一。一项针对随机对照临床试验的Meta分析表明,人参的主要成分人参皂苷和多糖作为化疗辅助用药能够提高非小细胞肺癌患者的客观缓解率和疾病控制率,降低患者不良反应发生率[23]。根据ceRNA假说,本研究构建了EZH2/hsa-miR-101-3p/GSEC调控网络。研究表明hsa-miR-101-3p可以通过调控EZH2的表达抑制肿瘤进展,在肺鳞癌中hsa-miR-101-3p可以通过下调EZH2抑制肿瘤活性,促进肿瘤细胞凋亡[24];在肾癌中,hsa-miR-101-3p可以通过调控EZH2抑制肿瘤的侵袭转移[25];hsa-miR-101-3p还可以与Syn-cal14.1a协同作用抑制EZH诱导的乳腺癌进展[26]。GSEC/hsa-miR-101-3p与肺腺癌细胞的铁死亡相关,上调GSEC或下调hsa-miR-101-3p的表达可以抑制肺腺癌细胞的增殖与迁移能力[27]。然而EZH2与GSEC之间的相互作用及EZH2/hsa-miR-101-3p/GSEC在调控EMT中的作用仍需进一步研究。

综上所述,本研究利用GEO、TCGA数据库及GenClip3文献挖掘平台等分析工具,筛选肺腺癌EMT过程的相关基因并对其进行了功能富集分析,以及PPI的构建,从基因组层面对EMT过程进行了系统描述,加深了对于EMT过程的认识。同时,筛选出了EMT过程中的关键基因并预测了人参等针对关键基因的治疗性中药,从患者预后角度构建了EZH2/hsa-miR-101-3p/GSEC调控网络,为肺腺癌侵袭、转移的治疗提供了新思路,为肺腺癌预后标志物的研究提供了参考。