基于中心体相关基因构建肝细胞癌预后模型及治疗方法分析

屈翔宇,戴恒文,童 旭,张文静,陆 进

(蚌埠医科大学 1.临床医学院2020级;2.临床医学院2022级;3.基础医学院,安徽 蚌埠 233030)

肝癌是一种全球常见的消化道恶性肿瘤,也是癌症的主要死亡原因之一[1]。在所有肝癌病例中,肝细胞癌(hepatocellular carcinoma,HCC)占80%以上,是原发性肝癌中最常见的一种亚型[2]。HCC具有隐蔽性强,预后差的特点。早期无明显症状,而到了晚期就仅有3个月的生存时间,并且5年生存率不足20%[3]。目前,HCC的治疗方法包括肝切除,肝移植,消融和联合治疗,但这些治疗方法的复发率较高,并且无法应用于晚期患者[4]。化疗和免疫治疗可应用于晚期患者,但由于肿瘤耐药性的发展,这些治疗也可能会无效[4]。因此,探索能够预测患者预后和治疗反应的预后标志物十分重要。

中心体是一种非膜质细胞器,由一对中心粒构成。在人体内,中心体主要参与有丝分裂纺锤体的形成、染色体的分离、细胞黏附、运动和纤毛的形成。最近的研究表明,中心体异常会引起细胞分裂失败以及染色体的不稳定,进而导致肿瘤的发生和进展。目前,在许多类型的癌症中均观察到中心体异常[5]。因此,对于中心体相关基因(centrosome related genes,CRGs)的研究将有利于寻找新的HCC预后标志物,并为HCC的治疗提供新的思路。

本研究从MSigDB和已报到文献中收集CRGs,使用LASSO Cox回归构建预后模型。国际癌症组联盟(International Cancer Group Consortium,ICGC)数据库用于验证模型的准确性。此外,还进一步探究了高风险组和低风险组的生物学功能、免疫浸润以及对化疗和免疫治疗的反应。

1 材料与方法

1.1 数据来源与处理从癌症基因组图谱(The Cancer Genome Atlas,TCGA)数据库(https://portal.gdc.cancer.gov/)中下载TCGA-LIHC队列的转录组数据和相应的临床数据,其中包括50例正常组织样本和374例HCC样本;从ICGC数据库(https://dcc.icgc.org/)中下载LIRI-JP队列的243例HCC样本的转录组数据和临床数据。CRGs从MSigDB(http://software.broadinstitute.org/gsea/index.jsp)数据库和已报到文献[6]中获得。具体分析思路与流程图见图1。

图1 研究流程图

1.2 CRGs差异分析和预后分析使用R软件“limma”包比较CRGs在TCGA-LIHC队列正常样本和肿瘤样本之间的表达差异,以|log2FC|>1且FDR<0.05为标准筛选出差异基因。对差异基因进行单因素Cox分析,以P<0.001为标准筛选出预后基因。

1.3 聚类分析基于预后基因的表达量,使用R软件“ConsensusClusterPlus”包对TCGA-LIHC队列的HCC样本进行聚类分析,根据累积分布函数(cumulative distribution function,CDF)和CDF曲线下面积变化确定最佳聚类结果。Kaplan-Meier曲线用于探究不同聚类之间的生存差异。R软件“GSVA”包用于探究不同聚类之间生物学功能的差异。

1.4 构建预后模型为防止模型过度拟合,以TCGA-LIHC队列作为训练组,使用最小绝对值选择与收缩算子(least absolute shrinkage and selection operator,LASSO) Cox回归对预后基因进行进一步筛选并构建预后模型。根据基因表达量和回归系数(coef)计算风险评分,风险评分=Σcoef×基因表达量。依据风险评分的中位数将HCC患者分成高风险组和低风险组。Kaplan-Meier曲线用于探究两组之间的生存差异,接受者操作特性(receiver operating characteristic,ROC)曲线用于评估模型的准确性。随后,使用ICGC数据库的LIRI-JP队列进行外部验证,采用与训练组相同的公式计算风险评分,并根据训练组风险评分的中位数将患者分成高风险组和低风险组。Kaplan-Meier曲线和ROC曲线用于验证预后模型。

1.5 基因富集分析使用R软件“limma”包,以|log2FC|>1且FDR<0.05为标准筛选出高风险组和低风险组之间的差异基因。京都基因与基因组百科全书(Kyoto Encyclopedia of Genes and Genomes,KEGG)和基因本体论(gene ontology,GO)用于探究差异基因的生物学功能。

1.6 免疫浸润分析使用XCELL、TIMER、QUANTISEQ、MCPCOUNTER、EPIC、CIBERSORT-ABS、CIBERSORT七种方法分析免疫浸润与风险评分的相关性。单样本基因集富集分析(single sample gene set enrichment analysis,ssGSEA)用于计算高风险组与低风险组中免疫细胞和免疫相关功能的评分,以比较两组之间的免疫活性。随后本文探究了免疫检查点基因在高风险组和低风险组中的表达。

1.7 化疗和免疫治疗分析肿瘤免疫功能障碍和排除(tumor immune dysfunction and exclusion,TIDE)评分用于评估免疫治疗的效果,本研究计算了高风险组和低风险组的TIDE评分。IMvigor210队列用于探究高风险组和低风险组对抗PD-L1治疗的反应。IMvigor210队列的表达数据和生存数据来源于http://research pub.gene.com/IMvigor210CoreBiologies。R软件“pRRophetic”包用于计算多柔比星、埃博霉素B、丝裂霉素C在高风险组和低风险组中的半最大抑制浓度(half maximal inhibitory concentration,IC50),以评估两组对三种化疗药物的敏感性。将高风险组和低风险组之间的差异基因上传到CMAP数据库,以预测能够治疗HCC的小分子化合物。

1.8 统计分析R软件(4.2.1)用于本文的统计学分析,Wilcoxon检验用于评估两组之间的差异,P<0.05具有统计学差异。

2 结果

2.1 CRGs的差异表达和预后分析从MSigDB中的11个基因集和先前文献中共收集到132个CRGs。使用TCGA-LIHC队列正常样本和肿瘤样本的表达数据鉴定出77个差异基因(图2A～B)。随后,对差异基因进行单因素Cox分析,发现有46个基因与预后相关(图2C)。

图2 中心体相关基因的差异分析和预后分析

2.2 聚类分析基于预后基因的表达,对TCGA-LIHC队列的HCC样本进行聚类分析,根据CDF曲线及曲线下面积变化确定K=3为最佳聚类结果(图3A～C)。热图结果显示,中心体相关预后基因在A聚类中高表达,在B聚类中低表达(图3D)。生存分析结果显示,B聚类的总生存率最高,A聚类的总生存率最低(图3E)。使用GSVA富集分析探究高表达聚类(A聚类)和低表达聚类(B聚类)生物功能的差异,结果发现细胞周期、DNA复制以及同源重组等通路在两组之间具有显著差异(图3F)。

图3 肝细胞癌样本的聚类分析

2.3 预后模型的构建基于TCGA-LIHC队列的样本数据,对46个中心体相关预后基因使用LASSO Cox回归分析,共筛选出8个基因(CEP85、KIF3A、MYCN、NDRG1、NPM1、PARD3、SAC3D1、TUBG1)(全称见表1)用于构建预后模型(图4A～B)。根据LASSO Cox回归分析的结果计算风险评分:风险评分=0.051 317×CEP85+0.006 093×KIF3A+0.029 806×MYCN+0.0016 26×NDRG1+0.002 449×NPM1+0.001 187×PARD3+0.000 952×SAC3D1+0.004 556×TUBG1。根据风险评分的中位数将患者划分为高风险组和低风险组。TCGA-LIHC队列和LIRI-JP队列的Kaplan-Meier分析结果均显示高风险组的总生存率(overall survival,OS)降低(图4C,E)。ROC分析结果显示,TCGA-LIHC队列1年、2年、3年OS的曲线下面积(area under curve,AUC)为0.789,0.696,0.692(图4D),LIRI-JP队列1年、2年、3年OS的AUC值为0.689,0.750,0.735(图4F),表明模型具有良好的预测能力。TCGA-LIHC队列的生存状态分布图显示,随着风险评分的增加,患者的死亡率升高,生存率降低(图5A～B)。热图显示8个风险基因在高风险组中高表达(图5C)。而对LIRI-JP队列的分析也得到了同样的结果(图5D～F)。对临床性状和风险评分进行单因素Cox和多因素Cox分析。单因素Cox分析结果显示,TCGA-LIHC队列的临床分期、T分期、M分期、风险评分以及LIRI-JP队列的性别、临床分期和风险评分是影响预后的危险因素(图6A,C)。多因素Cox分析结果显示,TCGA-LIHC队列的风险评分、LIRI-JP队列的性别、临床分期以及风险评分是影响预后的独立因素(图6B,D)。

表1 8个风险基因

图4 预后模型的构建和验证

图5 预后模型的风险评分和生存状态分布

图6 风险评分的独立预后分析

2.4 基因富集分析为进一步探究高风险组和低风险组潜在的生物学机制,对两组之间的差异基因进行KEGG和GO分析。KEGG结果显示,差异基因在肝癌、膀胱癌以及p53信号通路等癌症相关通路中显著富集(图7A)。GO结果显示,差异基因在中心体重复、Toll样受体的结合及中心粒等通路富集(图7B)。

图7 高风险组和低风险组的生物学功能和免疫状态

2.5 免疫浸润分析免疫浸润在肿瘤进展过程中具有重要作用,且能够反映免疫治疗的效果。本研究使用ssGSEA计算高风险组和低风险组免疫细胞和免疫功能的评分。结果显示,活化的树突状细胞(activated Dendritic Cells,aDCs)、肥大细胞、自然杀伤细胞(Natural Killer cells,NK cells)等免疫细胞的评分;检查点、主要组织相容性复合体I(Major Histocompatibility Complex I,MHCI)、2型干扰素等免疫功能的评分在两组之间具有显著差异(图7C)。之后分析了风险评分和免疫浸润之间的相关性。结果显示,大部分免疫浸润与风险评分呈正相关(图7D)。人的肿瘤共有6种免疫亚型,包括C1(伤口愈合)、C2(IFN-G主导)、C3(炎症)、C4(淋巴细胞耗竭),C5(免疫沉默)、C6(TGF-B主导)。本文探究了每种免疫亚型与风险评分的关系,发现C1与高风险评分相关,而C3与低风险评分相关(图7E)。免疫检查点分析结果显示,免疫检查点基因在高风险组中高表达(图7F)。

2.6 化疗和免疫治疗分析由于高风险组的免疫检查点基因表达量更高,表明高风险组可能对免疫治疗更敏感,TIDE评分的计算结果也表明高风险组接受免疫治疗的效果更好(图8A)。为进一步探究高风险组接受免疫治疗的效果,本研究使用IMvigor210队列探究高风险组和低风险组对抗PD-L1治疗的反应。结果显示,免疫治疗反应组的风险评分更高。高风险组接受抗PD-L1治疗后拥有更高的缓解率和更长的生存时间(图8B～D)。本研究还计算了多柔比星,埃博霉素B,丝裂霉素C三种化疗药物在高风险组和低风险组中的IC50,结果发现三种化疗药物在高风险组的IC50更低,表明高风险组对于三组化疗药物更敏感(图8E～G)。此外,将高风险组和低风险组之间的差异基因上传到CMAP数据库,以评分小于-80为标准共筛选出38个具有潜在治疗效果的小分子化合物(图8H)。

图8 高风险组和低风险组对化疗和免疫治疗的反应

3 讨论

肝细胞癌是一种发病率和死亡率都很高的恶性肿瘤。由于肿瘤的异质性,临床上很难预测患者的预后,给临床治疗带来了很大的困难。目前越来越多的研究表明,中心体的异常与肿瘤的发生和发展密切相关[7]。例如,中心体蛋白70(CEP70)过表达会导致微管紊乱和多极染色体的形成,进而促进胰腺癌的进展[8]。中心体蛋白72(CEP72)通过介导SERPINE1促进膀胱癌的侵袭和迁移[9]。但CRGs在HCC中的预后尚不清楚。本研究构建了CRGs在HCC中的预后模型,能够预测HCC患者的预后及对化疗和免疫治疗的敏感性,为临床治疗提供帮助。

本研究从MSigDB和先前的文献中收集到132个CRGs,其中有77个基因具有表达差异,46个基因与预后相关。最终筛选出8个CRGs(CEP85、KIF3A、MYCN、NDRG1、NPM1、PARD3、SAC3D1、TUBG1)用于构建预后模型。训练组和验证组ROC曲线的AUC值表明模型对HCC患者预后的预测较为准确。单因素Cox和多因素Cox分析结果表明风险评分可以作为肝细胞癌的独立预后因子。

CEP85在机体中负责调节中心体的分离。研究表明,CEP85与STIL之间的相互作用能够介导PLK4从而促进癌细胞的迁移[10]。KIF3A是驱动蛋白家族成员之一,在机体中参与纤毛生成、细胞迁移和分裂[11]。研究表明,KIF3A能够通过Wnt信号通路促进前列腺癌的增殖和侵袭[12]。此外,抑制KIF3A的表达可以通过抑制Rb-E2F信号传导和上皮间充质转化来抑制三阴性乳腺癌的生长和转移[13]。MYCN是MYC家族成员之一。实验表明,MYCN在HCC中过表达,并与HCC的不良预后相关[14]。最近的研究证实MYCN可作为HCC复发的生物标志物和抗HCC治疗的靶标[15]。NDRG1是NDRG家族成员之一,在机体中能够促进脂肪的生成并维持脂肪的功能[16]。研究表明,NDRG1在HCC中过表达,并与较差的预后相关[17]。此外,NDRG1可以通过调节整合素β3以抑制HCC的增殖和侵袭[18]。NPM1是一种核仁伴侣蛋白,在机体中参与基因组稳定,中心体复制,DNA修复等多种生物学功能[19]。NPM1能够调节激活转录因子5(ATF5)进而促进HCC的增殖和存活[19]。并且,NPM1也能够介导HCC对索拉菲尼的耐药[20]。PARD3是一种支架蛋白,与多数癌症的增殖,迁移和侵袭密切相关[21]。研究表明,PARD3在HCC中高表达,并与HCC的不良预后相关[22]。敲低PARD3能够抑制HCC的增殖,自噬和血管形成。SAC3D1是一种中心体相关蛋白,在机体中参与细胞周期,中心体复制和纺锤体形成。研究表明,SAC3D1在HCC中过表达,并与HCC的不良预后相关,可以作为HCC的预后标志物[23]。TUBG1是微管蛋白超家族的成员,参与细胞周期和微管的形成。TUBG1在HCC中高表达,并与HCC的不良预后相关。TUBG1在体外能够促进癌细胞的增殖,迁移和侵袭,并能抑制细胞凋亡[24]。

为了进一步探究高风险组和低风险组的生物学功能,本研究进行了KEGG和GO分析。结果表明差异基因主要与癌症和中心体相关。ssGSEA分析结果显示,aDCs、肥大细胞、NK细胞等免疫细胞评分,检查点、MHC I、2型干扰素等免疫功能评分在两组之间具有显著差异。本研究分析了免疫浸润和风险评分的相关性,结果显示,多数的免疫浸润与风险评分呈正相关。免疫亚型与风险评分关系的分析结果显示,C1和C2与高风险评分相关,C3和C4与低风险评分相关。免疫治疗已成为肿瘤治疗中的重要部分,而免疫治疗抑制剂是肿瘤免疫治疗中最主要的部分。目前CTLA-4和PD-L1等免疫检查点的抑制剂已应用于癌症的治疗,为癌症患者带来了新的希望。本研究探究了免疫检查点基因在高风险组和低风险组中的表达,结果表明高风险组接受免疫检查点抑制剂治疗的效果更好。IMvigor 210队列是一项为评估PD-L1抗体阿特珠单抗在晚期尿路上皮癌中活性的Ⅱ期研究[25]。使用IMvigor 210队列探究抗PD-L1治疗反应组和非反应组与风险评分的关系,结果表明高风险组接受抗PD-L1治疗的效果更好。除免疫治疗,本文还探究了高风险组和低风险组对多柔比星、埃博霉素B和丝裂霉素C的敏感性,结果显示,高风险组对三种化疗药物的敏感性更高。此外,使用CMAP数据库分析高风险组和低风险组之间的差异基因,最终鉴定出38个具有潜在治疗效果的小分子化合物。

综上所述,本研究从中心体相关基因出发,构建了HCC预后模型,为预测HCC患者的预后以及对化疗和免疫治疗的反应提供帮助。但是,本研究仍然存在一定局限性,使用公共数据库构建预后模型,还需要大量的临床数据验证模型的准确性。此外,8个风险基因在肝细胞癌中的作用还需要实验的进一步验证。