简介

肝细胞癌 (HCC) 是一种异质性疾病，预后不佳。然而，近年来，随着基因测序技术的飞速发展，我们对HCC分子发病机制的认识有了显着提高。 对大量样本进行高通量分析的累积数据表明，它可用于识别与 HCC 进展相关的关键生物标志物。然而，已知与 HCC 预后相关的生物标志物数量有限。

大约 18% 到 40% 的 HCC 患者报告了 CTNNB1 突变。 由突变的 CTNNB1 激活的致癌 Wnt/β-catenin 通路在肝脏代谢调节中起关键作用。 CTNNB1 突变的 HCC 具有独特的代谢形态类型，通常为胆汁淤积且不常出现脂肪变性。 最近的一份报告显示，参与不同代谢活动的蛋白质，如药物代谢、氨基酸代谢、糖酵解和糖异生，在 CTNNB1 突变肿瘤中富集， 11 表明它与代谢重编程密切相关。然而，潜在的机制尚不完全清楚。

代谢重编程是指细胞癌变过程中发生的代谢模式的显着变化，涉及糖酵解、氧化磷酸化、三羧酸循环、氨基酸代谢、核酸代谢和脂肪酸代谢等多个过程。肿瘤发生是一个多步骤过程，涉及对促进不受控制的增殖和消除细胞死亡的途径进行修改，这需要代谢重编程以提供用于细胞存活、生长和迁移的大分子。近年来，靶向代谢重编程已成为一种很有前景的新型治疗策略，但HCC的代谢重编程尚未被破译。

本研究通过癌症基因组图谱（TCGA）数据库对 CTNNB1 突变型 HCC 和 CTNNB1 野生型 HCC 的基因集进行富集分析，发现 CTNNB1 突变型 HCC 中显着上调的基因集均与代谢有关，进一步证实了近CTNNB1突变与代谢重编程的关系。我们进一步探讨了这些基因在 HCC 中的预后价值，并构建了一个由五个代谢基因组成的预后风险模型。多个数据集（包括 TCGA、ICGC、GSE14520和GSE116174）的结果; 共 876 名 HCC 患者）表明该风险评分在评估 HCC 预后方面非常准确。此外，我们发现该预后模型可能反映肿瘤的免疫微环境，因此具有很高的临床应用潜力。

结果

鉴定具有和不具有 CTNNB1 突变的 HCC 样本之间的差异代谢基因组

尽管据报道 CTNNB1 突变 HCC 具有独特的代谢特征，但相关的代谢基因仍然未知。 3 GSEA 显示 CTNNB1 突变型 HCC 中显着上调的基因组均与代谢途径相关（NOM P  < 0.05）（图 1）。这些包括“细胞色素 P450 对异生物质的代谢”（NES = 1.914，NOM P  = 0.002），“药物代谢 - 细胞色素 P450”（NES = 1.813，NOM P  = 0.006），“过氧化物酶体”（NES = 1.767，NOM P  = 0.02), '酪氨酸代谢' (NES = 1.758, NOM P  = 0.01), '卟啉和叶绿素代谢' (NES = 1.745, NOM P  = 0.002), '类固醇激素生物合成' (NES = 1.725, NOM P  = 0.015) , '初级胆汁酸生物合成' (NES = 1.710, NOM P  = 0.01), '氨酰基-tRNA 生物合成' (NES = 1.707, NOM P  = 0.02), '脂肪酸代谢' (NES = 1.707, NOM P  = 0.04) , '抗坏血酸和醛酸代谢' (NES = 1.700, NOMP  = 0.01），“苯丙氨酸代谢”（NES = 1.62，NOM P  = 0.02）和“精氨酸和脯氨酸代谢”（NES = 162，NOM P  = 0.04）（表 1）。相反，在野生型 CTNNB1 HCC 中，没有与代谢途径相关的基因组显着富集。

3.2. 在 TCGA 队列中构建预后模型

尽管 CTNNB1 已被证实是 HCC中 突变频率最高的基因之一，但CTNNB1突变 型 HCC 中 显着富集的代谢基因组是否影响 HCC 的预后尚不清楚 1、2、3、4、7、18 。通过单变量 Cox 回归分析和 K-M 生存分析，我们确定了 14 个与 HCC 总生存期最相关的基因（表 2）。我们使用 Lasso 回归和多 Cox 回归进一步缩小模型基因的范围以优化模型。最后，构建了由五个代谢基因组成的风险评分：风险评分 (RS) = CYP3A5 的标准化表达水平 * -0.00248 + ALAS1 的标准化表达水平 * -0.00445 + PRDX1 的标准化表达水平 * 0.003267 + GOT2 的标准化表达水平* -0.00887 + AMD1 的标准化表达水平 * 0.078805。TCGA 队列的中位风险评分被用作划分高风险组（RS > 1）和低风险组（RS < 1）的统一分界点。与低风险组相比，高风险组的 OS 明显降低（图 2A），预后模型在0.5、1、3和5年的ROC曲线下面积分别为0.834、0.779、0.724和0.731（图 2B）。ALAS1、CYP3A5 和 GOT2 的表达水平在低风险组中较高，而 AMD1 和 PRDX1 的表达水平在高风险组中较高（图 2C）。HCC 患者的死亡风险随着风险评分的增加而增加（图 2D-E）。单变量和多变量 Cox 回归分析发现，风险评分可以独立预测预后（图 2F-G）。

 TCGA 队列中预后模型的内部验证

患者根据AFP水平、血管侵犯、组织学分级、AJCC TNM分期、初始治疗后新发肿瘤及个体肿瘤状态等临床特征分为12个亚组进行生存分析。在每个亚组中，高风险组的 OS 率似乎低于低风险组（图 3A-F）。此外，该预后模型准确评估了134例复发患者的预后（其中肝内复发65例，局部复发46例，肝外复发23例）。ROC曲线下面积在0.5、1、3和5年分别为0.824、0.692、0.746和0.780（图 3H）。

在三个不同的独立队列中对预后模型进行外部验证

为了验证预后模型是否稳健，我们使用了三个独立的队列（ICGC（n = 230）、GSE14520（n = 239）和GSE116174（n = 64））进行外部验证。使用从 TCGA-HCC 队列获得的相同公式和截止值，将患者分配到高风险组和低风险组。与 TCGA 的结果一致，高风险组的 OS 显着低于低风险组（图 4A-C）。在 ICGC 队列中，0.5、1、3 和 5 年的 ROC 曲线下面积分别为 0.705、0.759、0.718 和 0.750（图 4D）。在GSE14520队列中，0.5、1、3 和 5 年的 ROC 曲线下面积分别为 0.762、0.712、0.687 和 0.663（图 4E）。在GSE116174队列中，0.5、1、3 和 5 年的 ROC 曲线下面积分别为 0.722、0.662、0.601 和 0.647（图 4F）。如单变量和多变量 Cox 回归分析所示，风险评分可以独立预测每个队列的预后（图 4G-I）。这些结果表明我们构建的预后模型能够普遍应用。

风险评分与临床病理特征的关系

既往研究表明，AFP 水平、组织学分级、临床分期（TNM、BCLC 和 CLIP）、肿瘤直径和血管侵犯与 HCC 的预后相关。 19 , 20 , 21 , 22 本研究利用四个独立队列的临床数据分析了风险评分与预后因素之间的相关性。结果表明，较高的风险评分与较高的 AFP 水平、较高的组织学分级、较大的肿瘤直径、血管侵犯和晚期临床分期（TNM、BCLC 和 CLIP）显着相关（图 5）。我们还对不同风险组进行了卡方检验，以进行临床特征的相关性分析。结果显示，两组在临床分期、组织学分级、肿瘤数量、AFP水平和生存状态方面存在显着差异（表 3,,4,4,,5,5,,66）。

基因本体功能富集分析

我们使用 Gene Ontology (GO) 功能富集分析来注释 DEG 的功能（图 6A) 在高风险组和低风险组中。结果表明，在高风险组中显着上调的基因与多种免疫调节过程相关，包括调节淋巴细胞活化、调节免疫效应过程和体液免疫反应（图 6B）。TIMER法对免疫细胞浸润的估计显示，高危组6种免疫细胞浸润水平均高于低危组（图 6C），表明预后特征可能通过调节肿瘤的免疫微环境影响HCC患者的预后。

风险评分与免疫细胞浸润的关系

我们进一步分析了预后模型与免疫细胞浸润的相关性，发现风险评分与六种免疫细胞浸润呈正相关：CD4 T 细胞（r  = 0.175，P  = 0.001）、CD8 T细胞 ( r  = 0.295, P  < 0.001)、B 细胞 ( r  = 0.256, P  < 0.001)、巨噬细胞 ( r  = 0.472, P  < 0.001)、树突细胞 ( r  = 0.412, P  < 0.001) 和中性粒细胞 ( r  = 0.448, P  < 0.001) (图 7A）。CTNNB1 突变的 HCC 以免疫排斥为特征， 23 最近的一项临床研究表明，该亚组对免疫检查点抑制剂的治疗不敏感。 8 我们的研究发现，风险评分与六个主要免疫检查点（PDL1、LAG3、CTLA4、CD276、PDCD1 和 TIGIT）的表达水平呈正相关（图 7B)，这与 CTNNB1 突变亚组的风险评分一致，低于非突变亚组的风险评分（图 7C）。

风险评分与肿瘤免疫微环境的关系

我们使用 ssGSEA 进一步探索预后模型与肿瘤免疫微环境 (TIME) 之间的内部关系。在免疫细胞浸润方面，在四个独立队列中，高危组的 Tregs、aDCs、Th2 细胞和巨噬细胞的丰度显着高于低危组（图S1A-D）。在免疫功能方面，结果显示低风险组的 II 型 IFN 反应明显强于低风险组（图S1A-D）。此外，免疫检查点相关基因组的表达在高危组中上调，这与之前的结果一致（图S1A-B）。

材料和方法

2.1。从 TCGA 收集数据

我们从 TCGA 网站 ( https://portal.gdc.cancer.gov/repository ) 获得了 374 个 HCC 样本的序列数据。相应的临床资料包括总生存时间、生存状态、性别、年龄、种族、甲胎蛋白（AFP）水平、体重指数（BMI）、血管侵犯、组织学分级、AJCC TNM分期、癌症家族史、肿瘤史从 UCSC Xena 网站 ( https://xenabrowser.net/ )获得 TCGA 队列的初始治疗后新肿瘤的存在和个体肿瘤状态。CTNNB1 突变样本列表来自 cBioPortal 网站 ( https://www.cbioportal.org/ )。我们保留了平均表达值大于 1 的基因，同时去除了低丰度的 RNA 测序数据。 15 , 16 本研究符合 TCGA 的发表要求 ( http://cancergenome.nih.gov/publications/publicationguidelines )。

2.2. 从 ICGC 和 GEO 收集数据

三个独立的队列用于外部验证（ICGC-LIRI-JP、GSE14520和GSE116174）。我们获得了基因表达文件（来自Illumina HiSeq RNA-seq平台的ICGC-LIRI-JP基因表达文件，来自GPL571平台的GSE14520基因表达文件和来自GPL13158平台的GSE116174基因表达文件）和来自基因表达的相应临床数据综合 (GEO) 数据库 (https://www.ncbi.nlm.nih.gov/geo/ )和国际癌症基因组学联盟 (ICGC , https://icgc.org/ ）。通过 SVA R 包（R Core Team，R Foundation for Statistical Computing，Vienna，Austria）中包含的战斗功能消除了 RNA-seq 数据集的批次效应。下载的配置文件均符合 GEO 和 ICGC 数据访问规则。

2.3. 基因集富集分析（GSEA）

本研究利用基因集富集分析来揭示 TCGA 队列中 CTNNB1 突变组 (n = 98) 和非突变组 (n = 276) 之间代谢相关基因的差异。选择带注释的基因集文件（c2.cp.kegg.v7.0.symbols.gmt）作为参考。阈值被确认为 NOM P值 < 0.05。

2.4. 代谢相关预后模型的构建和验证

我们首先使用单变量 Cox 回归和 Kaplan-Meier (K-M) 对每个基因进行生存分析，并选择两种算法中P  < 0.05 的基因作为候选基因来构建模型。 17 接下来，Lasso 回归分析有助于进一步缩小候选基因数量。最后，我们构建了一个风险评分系统，将每个代谢基因表现出的归一化表达水平与从多变量 Cox 回归分析中获得的回归系数相乘。TCGA 队列（n = 343）的中位风险评分用于对高风险组和低风险组进行分类。Kaplan-Meier (K-M) 生存分析（对数秩检验）与时间相关的受试者工作特征曲线 (ROC) 分析一起用于评估上述四个独立队列中预后模型所表现出的预测能力（ TCGA、ICGC、GSE14520和GSE116174）。本研究进行了单变量和多变量 Cox 回归分析，以确认风险评分是否可以独立预测预后。P  < 0.05 被认为具有统计学意义。R包'survminer'中的K-M方法用于生成生存曲线，R包'survivalROC'用于生成ROC曲线。需要指出的是，为了避免其他因素对患者预后的影响，我们排除了生存时间小于1个月的患者。

2.5. 临床病理参数与风险评分的相关性分析

我们使用 Wilcoxon 符号秩检验（2 组）或 Kruskal-Wallis 检验（>=2 组）分析临床病理学（包括 AJCC TNM 分期、Barcelona 分期、CLIP 分期、主要肿瘤大小、肿瘤分化、AFP和血管肿瘤细胞类型）和风险评分。P  < 0.05 被认为具有统计学意义，箱线图是通过 R 软件中实现的 beeswarm 包生成的。

2.6. 不同风险组之间差异表达基因 (DEG) 的鉴定

使用 R 包“limma”中的 Wilcoxon 测试方法检查了两个 HCC 组的组织中的 DEG。阈值被确认为 |log2 倍变化 (FC)| > 1.0 和 FDR < 0.05。基因本体（GO）富集分析进一步用于使用 R 包“clusterProfiler”对 DEG 进行生物学功能注释。

2.7. 评估不同风险组的免疫细胞浸润

从肿瘤免疫评估资源（TIMER）网站下载HCC患者免疫细胞浸润水平的相关数据，比较TCGA数据集中不同风险组之间的免疫细胞浸润水平。采用单样本基因集富集分析（ssGSEA，由 29 个免疫相关基因集代表的免疫细胞类型、功能和通路）来量化免疫细胞、功能或通路在高、低基因组中的活性或富集水平。来自四个独立队列的风险样本。从 ssGSEA 计算的归一化富集分数 (NES) 用于“GSVA”和“GSEABase”R 包中。采用独立样本 t 检验比较高危组和低危组之间免疫浸润水平和免疫功能的差异，以及建议P值 < 0.05 具有统计学意义。

2.8. 统计分析

R 软件 v3.6.1（R Foundation for Statistical Computing，Vienna，Austria）与 GraphPad Prism v7.00（GraphPad Software Inc，USA）一起协助进行统计分析。Fisher 精确检验或 Pearson 卡方检验有助于分析定性变量。定量变量的分析依赖于非参数 Wilcoxon 秩和检验（对于未配对样本）。应用 Kruskal-Wallis 检验对多个组进行归一化。如果以上未指定，则认为P  < 0.05 具有统计学意义。

讨论

在过去的几十年中，在了解 HCC 的危险因素和分子特征方面取得了相当大的进展。 24 , 25 , 26 然而，许多国家的发病率和癌症特异性死亡率继续增加。 26 , 27 现有的预后分期系统在指导准确治疗和预测更准确的临床结果方面仍有许多局限性。 28 , 29 我们需要利用基因测序技术进一步完善现有的预后分期系统，让患者受益于精准治疗。 8 越来越多的证据表明，代谢重编程可以显着影响 HCC 的发展，并可能与晚期疾病和不良临床结果有关 30 , 31 ；因此，靶向肿瘤代谢似乎有助于有效治疗 HCC。

最近的一份报告发现，参与药物代谢、糖异生、糖酵解和氨基酸代谢等不同代谢活动的蛋白质在 CTNNB1 突变型 HCC 中富集。 11 NRF2突变是HCC发生的早期事件，被认为是促进癌前肝细胞进展为HCC的重要因素。 32 现有研究发现，NRF2 信号通路的异常转导与辛伐他汀有关。 33 , 34 相比之下，CTNNB1突变通常发生在HCC进展的后期，对其具体机制的相关研究尚缺乏。考虑到这一点，我们在 TCGA 数据库中对 mRNA 表达谱进行了 GSEA。根据CTNNB1是否突变，将374个肿瘤组织分为两组，其中突变98个，未突变276个。有趣的是，突变组中显着富集的基因与 12 条代谢途径相关。我们采用单变量Cox回归、K-M生存分析、Lasso回归和多Cox回归逐步分析代谢基因，最终构建了五种代谢基因（包括AMD1、PRDX1、ALAS1、CYP3A5和GOT2）风险评分模型. 考虑到统一的风险评分公式和风险分类的阈值，将所有纳入的患者分为高风险组和低风险组。我们首先根据 TCGA 队列的临床特征（例如初始治疗后是否存在新肿瘤、个体肿瘤状态、AJCC 分期、组织学分级、血管侵犯和 AFP 水平）对预后模型进行了内部验证。然后将患者分为12个亚组。根据 K-M 生存曲线，在每个亚组中，高风险组的 OS 明显低于低风险组。随后，我们通过使用三个独立队列ICGC（n = 230）对预后模型进行了外部验证，我们首先根据 TCGA 队列的临床特征（例如初始治疗后是否存在新肿瘤、个体肿瘤状态、AJCC 分期、组织学分级、血管侵犯和 AFP 水平）对预后模型进行了内部验证。然后将患者分为12个亚组。根据 K-M 生存曲线，在每个亚组中，高风险组的 OS 明显低于低风险组。随后，我们通过使用三个独立队列ICGC（n = 230）对预后模型进行了外部验证，我们首先根据 TCGA 队列的临床特征（例如初始治疗后是否存在新肿瘤、个体肿瘤状态、AJCC 分期、组织学分级、血管侵犯和 AFP 水平）对预后模型进行了内部验证。然后将患者分为12个亚组。根据 K-M 生存曲线，在每个亚组中，高风险组的 OS 明显低于低风险组。随后，我们通过使用三个独立队列ICGC（n = 230）对预后模型进行了外部验证，与低风险组相比，高风险组的 OS 明显降低。随后，我们通过使用三个独立队列ICGC（n = 230）对预后模型进行了外部验证，与低风险组相比，高风险组的 OS 明显降低。随后，我们通过使用三个独立队列ICGC（n = 230）对预后模型进行了外部验证，GSE14520 (n = 239) 和GSE116174（n = 64）。与 TCGA 结果一致，高风险组的预后明显低于低风险组。正如单变量和多变量 Cox 回归分析所证明的那样，该预后模型可以独立预测 HCC 的预后。这些结果支持我们的预后模型具有很强的普遍适用性。GO富集分析表明，高危组中表达上调的基因与免疫调节密切相关。具有浸润性巨噬细胞的肿瘤与风险评分之间的相关性最为突出。众所周知，巨噬细胞在肿瘤组织中含量最多，对肿瘤有明显的调节作用，能够促进肿瘤细胞的增殖， 35 , 36 , 37 正如预期的那样，风险评分与六个常用免疫检查点（PDL1、LAG3、CTLA4、CD276、PDCD1 和 TIGIT）的表达水平呈正相关。我们得出结论，免疫抑制性肿瘤微环境可能是导致高风险组预后不良的主要因素。值得注意的是，风险评分与新兴的免疫检查点 CD276 (B7-H3) 密切相关，最近有报道称其在多种人类实体癌中高度过表达，并且通常与患者的不良临床结果相关，使其成为肿瘤免疫治疗的潜在靶点。 38 , 39 此外，高危评分可反映与生存结局相关的不良临床特征，如 AFP 水平（>300 ng/mL）、肿瘤血管侵犯、低肿瘤分化、晚期 AJCC 分期、晚期 BCLC 分期、晚期 CLIP 分期、新发肿瘤初始治疗后，主要肿瘤大小> 5 cm和多发肿瘤，这也可能有助于解释高危患者预后不良的原因。

腺苷甲硫氨酸脱羧酶 1 (AMD1) 是影响多胺（包括精胺）生物合成的必需酶。Zabala-Letona 发现 AMD1 的上调可以激活 PTEN-PI3K-mTORC1 通路来维持前列腺癌细胞的生长和增殖。 40 此外，Xu 报道 AMD1 对人胃癌的预后具有致瘤作用，但 AMD1 在 HCC 中的潜在功能仍不清楚。 41 Peroxiredoxin 1 (PRDX1) 作为抗氧化酶的过氧化物酶家族成员。Fang Y 和 Sun 等发现 PRDX1 在 HCC 组织中的高表达对应于不良的临床结果，其机制可能与促进肿瘤血管生成和调节细胞迁移和侵袭有关。 42 , 43 5'-氨基乙酰丙酸合成酶 1 (ALAS1) 在血红素生物合成过程中起限速酶的作用。最近进行的研究表明，ALAS1 可以影响许多细胞功能，并对非小细胞肺癌、结直肠癌和口腔癌有重要影响。 44、45 然而， ALAS1 在 HCC 中的作用仍不为人所知。CYP3A5 是一种细胞色素 P450 蛋白，在肝脏中许多致癌物质和抗癌药物的代谢中发挥作用。江等人发现CYP3A5通过调节mTORC2/Akt信号通路在HCC中发挥抗肿瘤作用。 46 近年来，谷氨酸草酰乙酸转氨酶 2 ( GOT2 ) 被多次报道与胰腺癌的进展有关，47、48 但 其在 HCC 中的作用仍不清楚。

尽管上述研究支持我们的预后模型具有很高的临床应用潜力，但该研究面临一些局限性。尽管从大量样本中积累的高通量分析数据得到了优化应用，但还需要通过前瞻性研究进一步验证。此外，HCC 中五个基因的特定生物学功能需要通过实验来探索。