基于TCGA数据库构建肝细胞癌中RNA结合蛋白预后评估模型

杨文娟， 马 莉， 马青梅， 李 婧， 张娇弟， 马梦婷， 张 倩， 刘欣跃

兰州大学第二医院检验医学中心，甘肃 兰州 730030

肝细胞癌(hepatocellular carcinoma，HCC)是一种常见的原发性肿瘤，占原发性肝癌的85%～90%，严重影响全球人民的健康。2018年世界卫生组织的癌症统计报告显示，在十大癌症中，HCC的发病率居第六位，死亡率居第四位，且男性高于女性[1]。我国癌症流行病学数据显示，2014年和2015年HCC的新发病例数分别为36.5万、37.0万，死亡病例数分别为31.8万、32.6万，表明我国HCC的诊治并未取得较好的效果，这与HCC的病因和致病机制不明有着重要的联系[2-3]。HCC的筛查关键是早发现、早诊断、早治疗，且多部指南推荐的筛查方法为肝脏超声和血清甲胎蛋白(AFP)，但该筛查方法仅推荐用于肝硬化、病毒性肝炎、过度饮酒、长期食用含黄曲霉毒素的食物、非酒精性脂肪性肝病等高风险人群[4-8]。由于筛查手段的适用人群具有一定的局限性，使得部分患者未能及时确诊而错过最佳的手术期和化疗期[9]。此外，发现常用的索拉非尼、仑伐替尼和系统化疗等可能引起腹泻、高血压和肝肾毒性等不良反应。因此，有必要探索HCC的生物学标志物，帮助HCC患者的早期诊断、早期治疗，研发安全性更高的化学治疗药物，延长患者的生存时间。

随着RNA测序技术的发展，RNA结合蛋白(RNA-binding proteins，RBPs)被发现是一种重要的转录后调控因子[10-11]。RBPs包括mRNA、lncRNAs、miRNAs和rRNAs等，参与剪接、修饰、转运、稳定性、定位和翻译等一系列生物学过程[12]。目前，已经发现超过1 500种RBPs，然而只有小部分RBPs的功能得到了研究[13]。研究[14]表明，包括癌症在内的各种疾病与RBPs的表达异常有着密切联系。在结肠癌中，RNA结合蛋白CEBP3表达下调影响结肠癌细胞的增殖和迁移，与患者不良的生存预后有关[15]。Shen等[16]发现，RNA结合蛋白RBM47通过调控AXIN1 mRNA稳定性和Wnt/β-catentin信号通路，从而抑制非小细胞肺癌转移。此外，RNA结合蛋白IGF2BP3被发现与ELAVL1共同调控癌细胞生长、增殖[17]。上述研究均表明，RBPs在各种肿瘤的发生发展及预后起着一定的作用。然而，有关RBPs在HCC中的作用，还需进一步深入探讨。

本研究主要通过TCGA数据库进行RNA测序数据的分析，挖掘对HCC预后生存有影响的RBPs，这将为HCC发病机制的探索、治疗靶点和预后生物学标志物的研发提供新思路。

1 材料与方法

1.1 数据来源本研究从癌症基因组图谱TCGA(https://portal.gdc.cancer.gov)数据库下载HCC的FPKM文件，包括374例HCC患者和50名正常肝脏组织标本的RNA测序数据；从国际癌症基因组联合会ICGC(https://icgc.org/)下载240例HCC样本的基因表达数据。

1.2 数据处理使用人类HG38表达文件，对RNA测序信息进行基因注释，用Perl软件获取基因表达矩阵文件。对已知的1 542个RBPs表达文件和HCC基因表达矩阵文件进行分析，得到HCC和正常癌旁患者RBPs的表达数据。按照|log2FC|>1且FDR<0.05的筛选标准，用R软件(3.6.1)Limma包筛选差异表达RBPs，采用“ggplot”和“pheatmap”包绘制差异RBPs的热图和火山图。

1.3 GO功能注释和KEGG富集分析GO功能注释主要包括基因执行的分子功能(molecular function)、所处的细胞组分(cellular component)及参与的生物学过程(biological process)等三方面的内容。GO功能注释和KEGG通路分析主要是明确差异RBPs在HCC发生发展中的作用。利用R软件对差异RBPs进行GO功能注释和KEGG通路富集分析。

1.4 PPI网络和功能模块分析将差异表达的RBPs输入到STRING数据库(http://string-db.org/cgi/input.pl)[18]，剔除游离的RBPs后，输出PPI相互作用的关系文件；并用Cytoscape软件的MCODE插件进行功能模块分析。

1.5 预后评估模型构建及验证从TCGA中下载患者的临床信息，使用Perl软件将患者对应的生存时间及差异RBPs的表达数据进行合并；提取ICGC数据库中患者生存时间及基因表达数据的矩阵文件；用R语言对TCGA与ICGC数据库的预后文件进行交叉合并，获得两大数据库共有的RBPs预后生存文件。对TCGA数据库中共有的RBPs进行单因素和多因素Cox分析，将多因素Cox分析获得的基因确定为关键RBPs。根据风险评估公式：风险评分=Exp1×表达系数值+Exp2×表达系数值+Expi×表达系数值(Exp表示不同关键RBPs的表达水平，i表示关键RBPs的数目)，对多因素Cox分析获得的关键RBPs的表达及生存数据建立预后模型。将TCGA数据库中的关键RBPs数据作为试验组，ICGC数据库中关键RBPs数据作为验证组。根据中位风险评分分别将试验组和验证组中的HCC患者划分为高风险和低风险组，Log-Rank检验比较两组间的总体生存时间差异。此外，对试验组获得的预后风险评估模型进行验证，分别绘制两组的生存曲线和生存状态图。

1.6 临床病理特征的预后分析通过R语言对患者年龄、性别、分期、分级、风险评分等进行单因素和多因素Cox分析，探究其是否对HCC预后生存有影响。

1.7 关键RBPs的筛选和验证采用Kaplan Meier-plotter(https://kmplot.com/analysis)和HEPIA(https://www.proteinatlas.org/)数据库对关键RBPs在HCC中的预后及表达进行验证，进一步说明所构建预后模型中RBPs的预后价值。

2 结果

2.1 HCC中差异RBPs的筛选从TCGA数据库中下载374例HCC患者和50名正常肝脏组织的RNA测序数据和临床数据，用R软件Limma包筛选差异表达RBPs。以|log2FC|>1且FDR<0.05为标准，筛选出82个差异表达RBPs，包括55个上调和27个下调的RBPs，并对差异表达RBPs作火山图和热图(见图1A～1B)。

图1 HCC中差异RBPs的筛选 A：差异RBPs聚类图；B：差异RBPs火山图Fig 1 Screening of differential RBPs in HCC A: heatmap analysis of differential RBPs; B: volcanogram of differential RBPs

2.2 差异表达RBPs的GO功能富集分析和KEGG通路分析将差异表达RBPs分为上调组(log2FC>1)和下调组(log2FC<-1)，使用R语言分析它们的GO功能富集分析和KEGG通路分析。结果表明，上调RBPs参与 mRNA的代谢、mRNA和RNA的分解代谢以及大分子甲基化等生物学过程；执行的功能为与mRNA 3′-UTR结合、RNA的催化活性、翻译调控活性和RNA解旋酶活性；也构成细胞质核糖核蛋白体、核糖核蛋白颗粒、细胞质颗粒和端粒酶等的细胞学组分。下调RBPs主要参与RNA代谢、翻译调控、细胞酰胺代谢和核酸磷酸二酯键水解等生物学过程；执行RNA催化活性、mRNA 3′-UTR、核糖核酸梅的活性和核酸酶的活性等分子生物学功能；构成mRNA帽状结合物、RNA帽状结合物、CCR4非复合物和尖端突起等的细胞学组分。此外，上调RBPs参与mRNA监控通路、DNA的复制、RNA的转录和降解等KEGG通路，下调的RBPs主要参与丙型肝炎、mRNA的监控通路和真核生物核糖体的发生等KEGG通路(见表1)。

表1 差异表达RBPs的GO功能富集分析

2.3 PPI网络和关键模块筛选将RBPs的PPI相互作用关系文件导入Cytoscape软件进行可视化，包含64个节点(RBPs)和126条边(相互作用)(见图2)。此外，用Cytoscape软件MCODE插件筛选相互作用关系的关键模块(见图3)，GO功能富集分析显示，该模块中RBPs主要参与病毒防御、Ⅰ型干扰素的反应、双链RNA的结合、核糖体的成分等，也参与丙肝、核糖体和mRNA的监视等KEGG通路。

图2 PPI相互作用的网络图Fig 2 The PPI network

图3 PPI网络中的关键模块Fig 3 Critical module of PPI network

2.4 HCC预后相关RBPs的筛选及生存模型的构建

ICGC和TCGA数据库中HCC的RBPs的表达数据和临床数据进行交叉分析，获得两大数据库共有的64个RBPs。对64个RBPs进行单因素Cox分析，得到22个RBPs(见图4)。然后，对22个RBPs进行多因素Cox分析，结果表明仅有6个RBPs是影响HCC预后总体生存时间(OS)的独立因素(见图5)。根据风险评估公式：风险评分=EZH2×表达系数值+SMG5×表达系数值+ANG×表达系数值+PPARGC1A×表达系数值+LARP1B×表达系数值+NR0B1×表达系数值，建立HCC中RBPs的预后评估模型。根据中位风险评分将生存时间和状态信息全面的370例HCC患者分为高风险组和低风险组(见图6)，Log-Rank检验结果表明高风险组患者的预后OS显著短于低风险组。为了更好地判断该模型的预后意义，绘制6个RBPs的生存曲线(ROC)，且生存曲线下面积(AUC)为0.737，表明该模型有较好的预后意义(见图7)。此外，绘制高、低风险组中6个RBPs的表达热图、风险评分和生存状态图(图8A～8C)。此外，使用ICGC数据库HCC中上述6个RBPs表达和生存数据进行验证，结果提示高风险组较低风险组预后OS短(见图9)，且AUC为0.787(见图10)。高、低风险组中6个RBPs的表达聚类图、风险评分和生存状态图如图11所示。上述结果表明，该模型对HCC具有很好的预后指导意义。

图4 TCGA队列中关键基因的单因素Cox分析;图5 TCGA队列中关键基因的多因素Cox分析Fig 4 Univariate Cox regression analysis of hub RBPs in the TCGA cohort; Fig 5 Multivariate Cox regression analysis of hub RBPs in the TCGA cohort

图6 TCGA队列高低风险组的生存曲线;图7 TCGA队列OS的AUC曲线

图8 TCGA数据库中6个基因预后模型的风险评分 A：表达热图；B：风险评分分布图；C：生存状态图Fig 8 Risk scores of 6 genes prognosis models in TCGA cohort A: expression heatmap; B: risk score distribution map; C: survival status diagram

图9 ICGC队列高、低风险组的生存曲线；图10 ICGC队列OS的AUC曲线

图11 ICGC队列中6个基因预后模型的风险评分 A：表达热图；B：风险评分分布图；C：生存状态图Fig 11 Risk scores of 6 genes prognosis models in ICGC cohort A: expression heatmap; B: risk score distribution map;C: survival status diagram

2.5 关键RBPs的列线图为了明确上述6个RBPs对HCC患者定量预后意义，建立6个基因的生存列线图，根据总得分和每个基因的预后得分可以明确HCC患者1年、2年和3年的生存几率(见图12)。此外，通过Cox分析HCC患者的临床病理特征对预后的影响，单因素分析表明，肿瘤分期和风险评分对HCC患者的OS具有影响(P<0.05)(见表2)，多因素分析表明，肿瘤分期和风险评分是HCC患者OS的独立预后因素(P<0.05)(见表2)。

图12 6个基因预后1年、2年、3年OS的列线图

表2 临床因素的预后影响

2.6 关键RBPs的预后和表达验证用Kaplan Meier-plotter数据库验证6个RBPs(包括EZH2、SMG5、ANG、PPARGC1A、LARP1B和NR0B1)与HCC患者预后OS的关系，结果表明EZH2、SMG5和NR0B1低风险组的OS较高风险组长(P<0.05,见图13A～13C)，PPARGC1A高风险组的OS较低风险组长(P<0.05,见图13D),而ANG和LARP1B的OS两组间差异无统计学意义(P>0.05,见图13E～13F)。数据分析表明，HCC患者中EZH2、SMG5、NR0B1为高表达，ANG为低表达，这也说明上述RBPs的表达变化可能影响HCC患者的生存预后。使用HPA数据库进一步验证关键RBPs的表达情况，免疫组化图中EZH2基因是高表达，这与本次研究结论一致；然而，LARP1B基因在肝癌和正常组织中的表达差异无统计学意义(见图14)。此外，HPA数据库中尚无SMG5、ANG、PPARGC1A和NR0B1的HCC免疫组化图。

图13 6种关键RBPs的生存分析 A：EZH2；B：NR0B1；C：SMG5；D：PPARGC1A；E：LARP1B；F：ANGFig 13 Survival analysis of 6 key RBPs A: EZH2; B: NR0B1; C: SMG5; D: PPARGC1A; E: LARP1B; F: ANG

图14 HPA数据库中正常组织和肝癌组织中关键基因表达验证Fig 14 Validation of hub genes expression in normal tissues and liver cancer tissues in HPA database

3 讨论

本研究通过对TCGA数据库中HCC组织及正常肝脏癌旁组织的RNA测序数据进行差异分析，成功构建HCC相关RBPs预后评估模型，这有利于HCC患者预后生存情况的预测。

通过GO功能富集分析对RBPs在HCC中所起的作用有了更深入的认识，表明RBPs参与mRNA和RNA的代谢、翻译调控、核酸磷酸二酯键的水解等过程。已有研究对RNA结合蛋白在不同肿瘤发生发展中的作用进行了阐述。如RNA结合蛋白NONO通过转录后调控SKP2与E2F8的表达进而影响乳腺癌细胞的增殖[19]。另外，研究证实在前列腺癌细胞中，RNA结合蛋白FXR1通过与FBXO4的3′UTR相结合调控FBXO4转录和mRNA代谢[20]；HUR也可与CDC6 3′-UTR相结合上调CDC6的表达[21]。RBPs也构成细胞质核糖核蛋白体、细胞质挤压颗粒和端粒酶等细胞学组分，研究证实了RNA结合蛋白Pur-alpha是构成胞质内挤压颗粒的成分[22]。KEGG通路分析表明RBP参与mRNA的监控通路、RNA的转录和降解及核糖体合成等过程，从而在HCC的发生发展中起作用。如RNA结合蛋白RBM38可以间接地影响病毒DNA的复制过程[23]。此外，RBPs相互作用中的关键模块关系图表明，在丙肝和Ⅰ型干扰素的反应中也起着作用。然而，有关RBPs在肿瘤中相关通路的研究相对较少，还需要大量的研究进行论证。

本研究构建的HCC预后评估模型的6个关键RBPs(EZH2、SMG5、ANG、PPARGC1A、LARP1B、NR0B1)，已在先前研究中明确了它们在不同肿瘤预后中所起的作用。EZH2在细胞增殖、衰老和凋亡中起着重要作用，研究表明EZH2可能在HCC的发生发展中起着作用[24]，且为多种肿瘤治疗的药物靶点[25]。LARP家族中研究较多的是LARP1，发现高表达LARP1的HCC患者，其生存时间较短[26]，然而，本次研究发现，高表达LARP1B的HCC患者较低表达者的总体生存时间更长，这可能是与不同生理条件下两者的表达水平有关。SMG5是一种介导衰老的调控因子，目前有关SMG5在肿瘤中的作用探讨，生物信息学相关的研究较多。Li等[27]发现，SMG5在HCC中较癌旁组织表达更高，这与本研究结论一致；Li等[28]研究结果表明SMG5可能对胃癌有较好的预后价值。此外，口腔癌细胞的迁移和侵袭与ANG的表达有着密切联系[29]；Hu等[30]揭示了ANG表达上调可能是胶质母细胞瘤患者预后不良的影响因素。众所周知，PPARGC1A在维持细胞内糖类和脂质等代谢中起着重要作用。PPARGC1A在miR-93-5p的调控下，影响肝癌的发展[31]；研究表明肺癌的转移可能与PPARGC1A的高表达有关[32]；另有多项研究探讨了PPARGC1A 的单核苷酸多态性与多种肿瘤如HCC、结直肠癌和乳腺癌等易感性的关系[33-35]。而Seo等[36]发现肺癌细胞系中NR0B1表达较高；NR0B1的表达受到在microRNA-106a靶向调控，进而促进乳腺癌的侵袭、迁移和增殖[37]。因此，从上述多项研究的结果中，可以确定构建HCC风险预后模型的6个基因在肿瘤的发生发展及预后中起着重要的作用。

本研究初步探索了在HCC预后中可能起作用的6个关键RBPs，旨在寻找有预后价值的生物标志物，对将来HCC生存时间的预测和靶向药物的研发具有较好的指导意义。