基于生物信息学的肝细胞癌预后模型的构建

李青玲，李树德，卫小娟,2，黄映光，李思熳*

(1.昆明医科大学 基础医学院 生物化学与分子生物学系，云南 昆明 650500；2.中国科学院 昆明动物研究所，云南 昆明 650201;3.云南省第一人民医院 普外一科，云南 昆明 650034)

肝癌发病率高居世界第六，是癌症相关死亡第二大原因.肝细胞癌(Hepatocellular carcinoma，HCC)占原发性肝癌的75%～85%[1].由于HCC病情进展迅速，并发症多，大多数患者确诊时已达中晚期[2]，且大量早期HCC患者在手术治疗后还会复发.目前对HCC的治疗决策多依赖于病理及影像学证据，缺乏分子学指标.而基于HCC特异性分子构建的预后预测模型则可进一步辅助临床医师进行HCC诊治[3].肿瘤突变负荷(tumor mutation burden，TMB)作为衡量免疫治疗效果的可靠指标，在HCC预后评估中至关重要.TMB表示外显子编码区每兆碱基中替换、基因插入或突变的总数[4]，是预测免疫检查点抑制剂治疗多种肿瘤(如肺癌、子宫内膜癌、乳腺癌和结肠直肠癌)疗效的有效指标[5-9].多项研究[5-9]中已初步揭示TMB与免疫微环境之间的联系及TMB与HCC中的免疫治疗之间的联系，但由于其测序的复杂性和高成本，不太可能广泛应用于免疫治疗疗效的临床评价.Zhang等[10]基于2020年公共数据进行了HCC中TMB相关分析，然而，近2年the Cancer Genome Atlas(TCGA)数据库更新后鲜有相关研究报道.

本研究旨在基于(TCGA)数据库探索HCC中肿瘤突变概况，然后通过基因差异分析及预后分析筛选出TMB相关的预后基因，并基于筛选的预后基因构建HCC预后预测模型.

1 材料与方法

1.1 数据来源

从TCGA数据库(http：www.ncbi.nlm.nih. gov/geo)下载与HCC相关的基因组和临床病理数据，基因组数据格式为FPKM.由于数据丢失，部分样本被删除，最终选择357个肿瘤样本进行后续分析.同样，从TCGA数据库中获取HCC体细胞突变数据，并使用R软件中的“maftool”包进行分析[11].

1.2 TMB计算和分组

运用Perl脚本计算TMB.并通过X-tile结合患者总体生存率(overall survival，OS)找到TMB的最佳截断值，然后依据该值将患者分为TMB-H组和TMB-L组.

1.3 获取差异基因

使用R包“limma”处理数据，筛选出与TMB相关的差异表达基因(differentially expressed genes，DEGs).采用|log 2 FC|>1,false discovery rate(FDR)<0.05筛选出所有差异基因，并使用“pheatmap”包对差异基因进行分层聚类.

运用R包“org.Hs.eg.db”进行Gene Ontology(GO)富集分析和Kyoto Encyclopedia of Genes and Genomes(KEGG)通路分析，“clusterProfiler”和“ggplot2”包用于注释和绘图[12].采用P<0.05筛选出潜在通路.

1.4 TMB预后指数的构建

使用Gene Expression Profiling Interactive Analysis(GEPIA)数据库(http://gepia.cancer-pku.cn/)[13]分析差异基因的OS率.然后根据logrankP<0.05筛选出关键基因.然后根据风险评分(risk score,RS)将所有患者分为高危组和低危组.采用K-M曲线表示RS水平与OS之间的关系.最后，绘制受试者生存工作特征曲线(receiver operating characteristic curve,ROC)，计算曲线下面积(caculate the area under the curve,AUC)，以评估预后模型的价值.

1.5 统计分析

使用R Studio(https://www.rstudio.com)进行统计分析[14].OS采用K-M和log-rank检验方法计算.2组连续变量比较采用t检验或单因素方差分析.P值<0.05，表示差异有统计学意义.

2 结果

2.1 TMB在HCC中的研究概况

首先，通过评估每个HCC样本的突变情况发现：单核苷酸多态性(Single nucleotide polymorphism，SNP)、错义突变和C>T突变是HCC中最常见的变异类型.其中，C>T变异最高，达 15 420.HCC中突变前10位的基因为：TP53、TTN、CTNNB1、MUC16、PCLO和ALB、RYR2、ABCA13、MUC4频率均高于10(P<0.001)(图1).

2.2 突变基因的瀑布效应和相互关系

图2(a)描绘了HCC中前30个突变基因的瀑布图.我们进一步分析前20个突变基因表达相关度，结果发现:CTNNB1与TP53、AX1N1表达负相关(P<0.05)，CTNNB1与SPTA1、OBSCN；MUC16与SPTA1、ABCA13、RYR2；PCLO与CACNA1E、HMCN1；RYR2与XIRP2；APOB与LRP1B；FLG与CSMD3、OBSCN；LRP1B与CACNA1E;HMCN1与SPTA1；TTN与LRP1B表达正相关(P<0.05)(图2(b)).

2.3 TMB与临床相关性

利用X-tile确定TMB=4.61为最佳截止值，并据此将357个样本分为TMB-H组和TMB-L组.TMB-H组和TMB-L组的平均值分别为8.19和2.79(图3(a)).生存分析结果显示TMB与HCC预后显著负相关(P<0.05).TMB高的患者预后较差[15](P<0.05)(图3(b)).图3(c)～(e)反映了TMB与临床病理特征的关系：在HCC患者中TMB与N分期(P<0.05)、年龄(P<0.05)、性别(P<0.05)相关.总的来说，老年(>65岁)无淋巴转移的男性TMB较高，但TMB水平与T/M分期、临床分期和分级无明显相关性(图3(f)～(i)).

(a)变异的分类 (b)变异类型 (c)单核苷酸多态性分类 (d)突变前10位的基因

(a)HCC中前30个突变基因的瀑布图

(b)HCC中前20个突变基因的相互关系

(a) TMB在各组间的分布 (b)TMB与HCC预后相关性 (c)TMB与年龄的关系

(d)TMB与性别的关系 (e)TMB与N期的关系 (f)TMB与T期的关系

(g)TMB与M期的关系 (h) TMB与临床分期的关系 (i) TMB与临床分级的关系

2.4 差异基因的富集分析

首先通过基因差异分析，筛选出39个与TMB相关的DEGs.图4(a)为DEGs的聚类热图.我们进一步通过GO功能富集分析研究DEGs相关的主要生物过程(biological process，BP)、分子功能(molecular function，MF)和细胞成分(cellular component，CC) (图4(b)).KEGG通路分析显示，DEGs主要富集于细胞外基质组织、细胞外结构组织、硫化合物生物合成过程和硫化合物代谢过程(图4(c)).Gene set enrichment analysis(GSEA)分析结果显示，TMB-H组富集在RNA代谢过程，TMB-L组富集在其他途径包括调控超分子纤维组织、调控细胞器组织、正调控细胞器组织、肌动蛋白丝组织、调控细胞骨架组织、肌动蛋白聚合或解聚、调控肌动蛋白丝基过程、调控肌动蛋白丝组织(图4(d)).

2.5 COX回归分析

基于生存分析，可进一步发现：SFRP4、IL7R、FBLN2、COLEC10和CHGA的表达可能是影响HCC患者预后的独立因素(图5).由此，可利用COX回归分析建立预后预测模型，并利用ROC曲线验证模型的准确性.通过R包“merge”下载了357例HCC病例的转录组数据.根据多元COX回归模型，计算TMBPI为:

PI=(0.088 945 38×SFRP4-0.009 912 19×COLEC10-0.000 690 51×CHGA-0.024 612 34×FBLN2-0.165 479 22×IL7R).

每个个体的RS通过TMBPI得到，然后可根据风险评分的中位数将样本分为高风险组和低风险组(图6(a)).该模型1,3,5年生存率的AUC分别为0.64、0.67、0.59(图6(b)).结果显示高危组患者预后较差(P<0.001)(图6(c)).

(a)DEGs的聚类热图

(b)DEGs的GO功能富集分析 (c)DEGs的KEGG通路分析

(d)TMB组间的GSEA分析

(a)SFRP4的表达与HCC患者预后相关 (b) IL7R的表达与HCC患者预后相关 (c) FBLN2的表达与HCC患者预后相关

(d)CHGA的表达与HCC患者预后相关 (e) COLEC10的表达与HCC患者预后相关

3 讨论与结论

HCC侵袭性强，且具有高度异质性.多项研究[16-18]表明肝癌细胞内信号通路与正常肝细胞迥异.例如，PI3K/AKT和IKK/NF-κB通路在HCC中被激活，从而促进细胞增殖并诱导上皮-间质转化[19,20].此外，免疫系统在肿瘤的发展过程中起着至关重要的作用.Marina等[21]观察到MYC/CTNNB1通路的激活促进了免疫逃避和肝癌药物耐药性.TMB作为免疫系统识别和攻击肝癌细胞的关键因素，在很大程度上影响患者的预后.此外，在各种癌症中可观察到TMB与客观缓解率显著相关[22].在本研究中，我们发现错义突变是HCC中最常见的类型.研究表明TP53的错义突变可将HCC分化为不同的亚型，并可能促进疾病的进展[23].在HCC中SNP变异占绝大多数，在我们的研究中，C>T变异数量达到 15 420，这与HCC的发生密切相关[24,25].

功能分析结果表明，DEGs主要与细胞外基质组织和硫化合物代谢过程有关.细胞外基质是HCC免疫微环境的组成部分，靶向给药已取得一定成效[26].硫化合物代谢可能与HCC中的自噬密切相关[27].此外，我们发现RNA代谢过程的活性与TMB之间存在明显的负相关.越来越多的证据表明，不同的非编码RNA(non-coding RNAs，ncRNAs)参与多种肝脏疾病(包括乙肝、丙肝和非酒精性脂肪肝病)代谢状态调控，并最终导致HCC[28-30].这可能是TMB-H组总体预后不佳的最主要原因之一.

SFRP4、IL7R、FBLN2、COLEC10和CHGA在HCC发生发展中发挥重要作用.SFRP4位于染色体7p14.1上，含有一个富含半胱氨酸的结构域.该结构域可通过与Wnt直接结合来调节Wnt信号通路，从而形成沉默复合体并抑制HCC[31].IL-7R在T细胞分化和淋巴细胞发育中能发挥作用[32].IL-7R上调可激活细胞内通路，诱导相关分子表达，促进肝癌细胞增殖和迁移[33].COLEC10编码肝脏胶原凝集素1，是C凝集素家族的一员[34].研究表明[35]，COLEC10受miR-452-5p的调控，促进肝癌细胞的增殖、侵袭和迁移.CHGA(或CGA)在HCC样本中高度表达，既往研究[36]表明它可以作为HCC的辅助诊断分子.本研究显示FBLN2高表达组与低表达组的预后差异有统计学意义.然而，FBLN2在肝癌中的作用及其应用价值尚未见报道，这可能是一个值得探索的靶点.

基于此，本文构建了多元COX回归模型，实现了对患者个体生存概率的预测和分析.在本研究中，该模型1、3、5年生存率的AUC分别为0.64、0.67、0.59，具有一定的准确度.该模型将为探索HCC新的预后因素提供了新的思路.此外，该模型无须识别患者的体细胞突变，而是通过直接检测TMB相关的5个基因，这使得基于关键基因的靶向测序技术更加常规.基于5个基因的预后预测模型可有效评估患者预后，这可能有助于指导肝癌患者的临床治疗.在常规临床实践中，临床医师多通过病理分期决定患者治疗手段.但HCC病理活检损伤大，多数患者依赖术前影像学(如CT、MRI)或常规生化指标(如AFP,CEA,CA199)等评估病情.本研究创新性地将TMB、DEGs与HCC临床数据相结合，通过COX回归模型预测患者预后，其可操作性强，侵入性更小.在精准医疗时代，利用分子分析结合病理分期的预后预测将比传统方法更准确.

该模型具有对HCC患者的实质性预测能力和价值，然而其局限性也应当被意识到.首先，目前研究缺乏患者的免疫治疗信息，这限制了我们对TMB介导的HCC免疫微环境变化的研究.其次，由于TCGA项目中只有一部分样本具有临床病理和TMB数据可用于分析，目前的样本量存在限制.最后，由于缺乏肝癌患者的全外显子测序信息，我们没有通过本中心临床数据来验证结论.因此，预测模型的有效性值得在未来使用独立队列进行更多的外部验证.

综上，本研究建立了一个新的基于5个基因的HCC患者预后预测模型.该模型可能是HCC预后风险分层的有效工具，并为肝癌的临床诊断结合病理和分子分析提供了新的思路.