基于m7G相关miRNA的肝癌风险预后模型构建

邱康,朱谦,金晨,孙传东

(青岛大学,山东 青岛 266071 1 医学部; 2 附属医院肝胆胰外科; 3 附属医院小儿外科)

肝癌是世界上是第七常见的癌症类型,5年生存率仅有18%,是致死率第二高的肿瘤,仅次于胰腺癌[1]。随着环境中致癌因素的增加,肝癌在我国的发病率逐年上升[2]。尽管目前的临床管理已经极大地改善了肝癌病人的生存率,但由于转移率高,肝癌病人的一般预后仍然极差[3]。因此,了解肝癌的发病机制,找到一种评估病人预后以及指导临床治疗的方法十分重要。microRNA(miRNA)作为短的单链RNA分子可以通过翻译抑制调节蛋白质编码基因[4],肝组织miRNA在肝癌中表达失调,参与了肝癌的发生和发展[5]。肝癌细胞通过多种方式影响miRNA的功能,包括染色体改变、聚腺苷基化修饰及转录因子的异常表达等[6]。N7-甲基鸟苷(m7G)具有特定修饰依赖性化学性质,普遍存在于哺乳生物中并且修饰特定转录位置。m7G存在于真核生物mRNA的5’端帽中,且在转运RNA(tRNA)、核糖体RNA(rRNA)、小核RNA(snRNA)和小核仁RNA(snoRNA)等其他种类的RNA中也发现了类似情况[7]。已有研究表明,介导内部m7G修饰的最具特异性酶甲基转移酶样1(METTL1)和其辅助因子WDR4,参与调节miRNA的m7G修饰,从而影响生物发生和细胞迁移[8]。本研究基于生物信息学方法,应用肿瘤基因组图谱(TCGA)数据库构建了m7G修饰相关的miRNA的风险预后模型,并利用Estimate数据库对免疫细胞浸润和免疫功能进行分析,旨在构建一个模型来评价肝癌病人的预后,为肝癌病人的治疗和预后评价提供参考。现将结果报告如下。

1 材料和方法

1.1 数据下载与处理

从TCGA数据库(https://portal.gdc.cancer.gov)中检索并下载肝癌的miRNA表达数据和临床特征数据,其中正常组织50例,肿瘤组织375例,临床数据371例。利用软件Strawberry Perl 5.32.1.1将基因表达数据进行转换整理为表达矩阵,用于后续分析。GSEA数据库(http://www.gsea-msigdb.org/gsea/index.jsp)下载得到m7G相关基因并进行单因素Cox回归筛选,得到差异基因用于后续分析。在Estimate数据库中下载肝癌相关的373个样本的免疫评分,用于后续免疫细胞浸润和免疫功能的分析。

1.2 m7G差异基因分析

使用R语言中的“edgeR”和“Limma”包,分析肿瘤组织和正常组织间差异基因的表达情况。差异基因的筛选标准:|log2FC|≥1,P<0.05。

1.3 构建m7G风险评分模型

利用“survival”包和“glmnet”包对差异表达miRNA进行单因素和多因素回归分析,得到预后相关miRNA;同时结合临床数据计算模型的风险评分,构建m7G风险评分模型。公式为:风险评分(Riskscore)=∑回归系数×基因表达量。然后利用“Rtsne”包执行PCA分析。依据风险评分,绘制Kaplan-Meier生存曲线以及ROC曲线,进一步对该模型准确性进行评估。

1.4 风险评分与临床因素分析

利用“survival”包构建具有风险评估和临床因素的Cox模型。在临床肝癌病人中,通过单因素和多因素Cox回归分析,评估影响其预后的独立危险因素。

1.5 功能富集分析

使用R语言中的“clusterprofile”“enrichplot”“ggplot2”包进行基因本体(GO)和KEGG通路分析,以确定与差异表达的m7G相关miRNA的作用位点和信号通路。

1.6 风险评分模型与免疫浸润的相关性分析

利用“Limma”包和“edgeR”包对肝癌病人高风险组和低风险组进行基因表达差异分析(筛选标准:|log2FC|>1,FDR<0.05)。采用“GSVA”包中的基因集富集分析算法(ssGSEA)计算肝癌病人免疫细胞和免疫功能表达的基因丰度并绘制热图,进行免疫细胞和功能的相关性分析,比较高风险组和低风险组的免疫细胞表达及免疫功能。

1.7 统计学处理

利用R语言软件(4.1.3)进行统计学分析。应用单因素和多因素Cox回归分析风险模型的独立预后价值,Kaplan-Meier法进行生存分析,Mann-Whitney检验比较两组免疫浸润和免疫通路激活情况。所有生存结局均以P值和风险比(HR)表示,置信区间(CI)为95%。P<0.05为差异有显著性。

2 结 果

2.1 肝癌病人与m7G相关的miRNA差异基因

从GSEA数据库下载得到34个可能与m7G相关的基因,通过单因素Cox分析肝癌病人预后数据,结合相关文献结果[1]筛选出12个与m7G相关的基因。利用TargetScan数据库预测m7G调控的相关miRNA,并得到相应表达矩阵。结合TCGA数据库中肝癌病人的miRNA数据与得到的m7G相关miRNA进行差异分析,根据相关筛选标准(|log2FC|>1、FDR<0.05)获得363个miRNA差异基因并绘制火山图(图1A)。其中上调基因320个,下调基因43个,差异表达最显著的前20个上游miRNA热图见图1B。

A:火山图;B:热图。

2.2 Lasso回归筛选m7G相关miRNA

对363个差异表达miRNA与生存时间取交集进行单因素 Cox回归分析,结果显示共有32个miRNA与肝癌病人的预后显著相关。将候选的32个miRNA纳入Lasso回归构建模型,选取误差最低的数值,进一步筛选出11个有意义的特征基因miRNA,分别为miR-4661-5p、miR-301a-3p、miR-760、miR-9-3p、miR-7156-5p、miR-561-5p、miR-3911、miR-513c-5p、miR-4652-3p、miR-346、miR-548aq-5p(图2A、B)。根据临床生存时间进行多因素Cox回归分析显示,结果显示has-miR-5p、has-miR-9-3p、has-miR-3911、hsa-miR-513c-5p、has-miR-4652-3p、has-miR-538aq-5p差异有统计学意义(HR=1.004～1.553,P<0.05)。见图2C。

A:Lasso回归系数分析;B:Lasso回归交叉验证图;C:m7G相关miRNA多因素独立预后分析。

2.3 m7G相关miRNA的风险模型构建

应用上述11个与m7G相关的miRNA构建风险模型,并且基于风险评分将肝癌病人分为高、低风险组。PCA聚类分析结果表明,二者区分明显(图3A)。根据风险评分构建Kaplan-Meier曲线,经Cox回归分析显示,低风险组的整体生存水平明显高于高风险组(HR=2.23,P<0.001) (图3B)。ROC曲线分析结果显示,肝癌病人的模型风险评分对1、3、5年生存评估的AUC分别为0.674、0.725和0.697(图3C)。

2.4 风险评分与临床因素的模型构建

为了评估风险水平模型在临床中的预测效果,将风险水平与临床因素中的年龄、性别、肿瘤分期(G0=T1+T2,G1=T3+T4)和风险评分纳入分析。根据风险评分以及临床因素绘制列线图,对肝癌病人的预后进行预测,结果显示,该预测模型的一致性指数(C-index)为0.611(0.563～0.658),肝癌病人第1、3、5年的生存率分别为78.1%、51.8%和37.1%(图4A)。为了进一步评估预测模型是否符合实际情况,采用Bootstrap自抽样法进行检测,并绘制了预测模型的校准曲线。结果表明,3年和5年拟合线重合度较高,表明该预测模型具有较好的生存预测功效,其3年和5年生存期预测值与实际值相似(图4B)。为了鉴别预测模型的区分度,绘制了预测模型的ROC曲线,结果显示1、3、5年的AUC分别为0.613、0.701和0.673(图4C)。

A: 风险评分和临床因素预后相关列线图;B:预测模型标定曲线;C:风险评分和临床因素模型对预后区分度的ROC曲线。

2.5 风险评分与免疫评分相关基因的功能分析

对风险模型与免疫功能的相关性分析显示,在风险模型中的高风险和低风险组中发现了539个差异表达的mRNA。将来自Estimate数据库的经免疫评分分组的960个基因与风险模型的差异表达基因相交,产生183个与风险模型和免疫浸润相关的mRNA(图5A)。GO分析显示,这些基因主要参与了消化、内分泌过程、摄食行为等生物学过程(BP);细胞成分(CC)主要富集在含胶原蛋白的细胞外基质、分泌颗粒腔、细胞质囊泡腔中;分子功能(MF)富集的主要是受体配体活性、信号受体激活剂活性、生长因子活性等(图5B)。KEGG分析显示,神经活性配体-受体相互作用、胰腺分泌、蛋白质的消化吸收是主要参与的富集通路(图5B)。在众多基因中筛选出了与免疫功能和免疫细胞相关的mRNA,其中S100A8、KRT4、FCGBP与免疫浸润呈正相关,而SOHLH1、CRISP2、CTSV则与免疫浸润呈负相关(图5C)。

A:风险评分与免疫评分差异基因韦恩图;B:GO和KEGG分析;C:差异基因与免疫浸润相关性热图。

2.6 风险模型与免疫微环境之间的关系

应用ssGSEA分析肝癌病人风险评分与免疫细胞表达丰度以及免疫功能相关性结果显示,高风险和低风险组中辅助性T细胞、人类白细胞抗原(HLA)和主要组织相容性复合物Ⅰ类(MHCⅠ)在肿瘤免疫微环境中显著表达(图6A)。肿瘤浸润淋巴细胞(TIL)和CD8+T细胞在肝癌病人肿瘤免疫微环境中相关性最高(r=0.870),滤泡辅助性T细胞(Tfh)和肥大细胞的相关性最低(r=-0.004),在免疫功能中T细胞共抑制途径和免疫检查点显示最高的正相关性(r=0.930)(图6B)。高风险组活化的树突状细胞(aDCs)、巨噬细胞相较于低风险组富集更明显(P<0.05),而低风险组中B细胞、肥大细胞、中性粒细胞、NK细胞、辅助性T细胞、细胞溶解活性、亚炎症以及Ⅰ、Ⅱ型干扰素反应均显著升高(P<0.05)(图6C)。因此,肝癌高风险组病人可能更倾向于应用抗巨噬细胞或抗aDCs细胞进行免疫治疗。

A:免疫细胞与免疫功能表达丰度热图;B:免疫细胞和免疫功能之间的相关性热图;C:不同风险组中免疫浸润相关表达丰度的箱式图。

3 讨 论

m7G是最普遍的RNA修饰之一,有研究表明m7G修饰显著参与肿瘤的发生发展[9]。此外,m7G修饰不仅发生在mRNA、rRNA、tRNA上,还通过修饰miRNA介导RNA的代谢和功能[8]。但是针对m7G修饰miRNA是否同样影响了肝癌的发生发展,尚未见报道。既往研究表明,miRNA在转录后水平上调节基因的表达,在人的肝癌和正常肝脏中的表达存在显著差异,miRNA的失调提示肝脏可能发生早期癌变[10]。已有研究结果显示,m7G修饰相关基因METTL1和WRD4参与了对miRNA的调控,推测可能有更多的m7G修饰基因参与了miRNA的调控。本研究提取肝癌中具有预后价值的m7G相关基因,筛选下游相关miRNA构建风险模型,评估模型的预测价值,为临床治疗提供参考。

本研究通过肝癌的预后数据筛选得到了12个m7G相关基因,并通过Cox和Lasso回归分析构建模型,筛选出11个受m7G修饰影响的miRNA。其中血清外泌体中的miR-4661-5p可以作为早期肝癌的潜在诊断标志物[11],miR-301a-3p可以通过靶向VGLL4抑制肝癌的恶性进展,miR-760可以调控下游的NACC-1促进肝癌的发生发展[12],miR-561-5p通过影响CX3CR1/NK细胞的浸润和功能促进转移性肝癌的发生发展[13],海绵化的miR-346在circMDK促进肝癌进展的过程中发挥作用[14],而miR-9-3p则通过下调TAZ在肝癌细胞中发挥肿瘤抑制作用[15]。目前肝癌研究中尚未见miR-7156-5p、miR-3911、miR-513c-5p、miR-4652-3p、miR-548aq-5p报道,其功能和机制有待进一步研究。

本文通过12个相关miRNA与生存时间构建风险模型,计算风险评分并进行PCA聚类分析,结果显示风险评分能够较好地将病人区分为高低风险组;Kaplan-Meier分析显示,高风险和低风险组病人生存率存在差异,低风险组存活率显著高于高风险组;ROC曲线分析显示,该模型的1、3和5年生存率的AUC均大于0.6,表明预后模型具有预测价值,尤其是对肝癌病人3年生存率的预测。本文根据临床因素建立预后模型列线图分析显示,风险评分对于病人预后的评价具有显著意义,根据风险评分构建的预后模型可以更好地预测临床病人的预后情况,为临床治疗提供指导意见。

RNA甲基化修饰除了直接影响肿瘤发展外,还通过调节肿瘤免疫力来影响肿瘤治疗的有效性[16]。

作为m7G的关键基因,METTL1在肝癌放疗中可能通过m7G修饰参与了免疫微环境的调节[17]。本文利用已经构建的风险模型与免疫相关基因分析得到m7G修饰相关miRNA的下游基因,并分别进行了功能富集分析和相关免疫功能和免疫细胞的分析,结果显示,受体配体活性、信号受体激活剂活性与免疫浸润的抗肿瘤免疫力相关[18];ssGSEA分析结果显示,在m7G相关的免疫环境中,辅助性T细胞、HLA和MHCⅠ的表达水平更高。不同风险组免疫浸润相关表达丰度分析显示,高风险组肝癌病人可能更倾向于抗巨噬细胞或抗aDCs细胞的免疫治疗,提示m7G相关的miRNA可能与免疫浸润有关,这可为免疫治疗提供可能的指导方案。

综上所述,本文研究构建了与m7G修饰相关miRNA的风险预后模型,将肝癌病人分为高低风险组,基于风险评分构建预测模型对于肝癌病人的预后评价良好。风险评分可以作为预测病人预后的独立预后因素,对于病人的预后预测有着显著意义。免疫浸润相关分析提示高风险组病人的预后与巨噬细胞和aDCs细胞显著相关,可能作为免疫治疗的潜在方向。该模型的建立有助于为肝癌病人的生存和免疫治疗提供预测和指导,但是其实际价值还有待于更大样本量和相关实验的验证。