G6PD介导缺氧驱动干细胞特性：基于干性与缺氧的肝细胞癌预后分层新策略及其机制研究

本研究构建了干性与缺氧相关预后指数（SHRPI），整合多组学数据与机器学习算法，揭示G6PD通过调控HIF-1α形成正反馈环路驱动肝癌干细胞特性，为个体化治疗提供新型预测工具与靶向策略（BI2536）。

　　

1 Introduction

肝细胞癌（HCC）是一种高度侵袭性的恶性肿瘤，其高复发率和广泛治疗抵抗性与癌症干细胞（CSCs）的存在和功能密切相关。CSCs具有自我更新和分化可塑性，能够逃避免疫监视，在肿瘤起始、进展和治疗抵抗中发挥关键作用。这些细胞通过激活Wnt/β-catenin和Notch等胚胎发育信号通路上调药物外排泵，增强对常规治疗的抵抗性。作为肿瘤复发的储备库，CSCs严重依赖缺氧肿瘤微环境（TME）生存和维持干细胞特性。

在实体瘤中，低氧水平不仅促进血管生成和代谢重编程以支持肿瘤生长，还激活缺氧诱导因子（HIFs），转录上调Oct4和Nanog等干性相关基因，进一步维持CSCs的干细胞样表型，导致HCC侵袭性和治疗抵抗。此外，缺氧将TME重编程为免疫抑制生态位，放大CSC驱动的恶性程度。瘤内缺氧诱导TGF-β和IL-6等细胞因子分泌，促进调节性T细胞（Tregs）和髓源性抑制细胞（MDSCs）招募，这些免疫抑制成分协同抑制细胞毒性T细胞（CD8⁺ T细胞）活性，保护CSCs免受免疫清除。缺氧驱动的PD-L1和CTLA-4等免疫检查点上调加剧免疫逃逸。随着肿瘤进展，氧消耗增加加剧TME缺氧，形成自我强化的恶性循环，增强肿瘤侵袭性。

尽管干性-缺氧轴在HCC进展中的作用已被认可，但仍缺乏系统、严谨且有效的预后指数整合干性和缺氧特征以分层患者和识别高危亚组来指导精准治疗策略。当前HCC预后模型主要基于单一分子特征（如干性指数或缺氧评分）或临床病理参数，这些模型的有限维度限制了预测患者结局的准确性。因此，亟需开发综合预后工具整合多维分子特征（包含干性和缺氧）与临床病理参数以改进风险分层、提高预测准确性并为高危亚组的个体化治疗策略提供理论依据。

2 Materials and methods

本研究从TCGA-LIHC队列获取HCC患者发现队列，包括基因表达谱、拷贝数变异（CNV）数据、体细胞突变数据和临床信息。验证队列包含从GEO数据库检索的GSE14520数据集的基因表达谱和临床信息。经过全面筛选，本研究纳入340名TCGA-LIHC患者和242名GSE14520患者。转录组数据标准化应用log2(FPKM + 0.001)转换，并使用“SVA”R包中的Combat算法进行批次效应校正。此外，还纳入180名肝移植（LT）患者和147名GSE104580数据集患者以评估基因签名对无复发生存（RFS）和经导管动脉化疗栓塞（TACE）治疗反应的相关性。

干性特征使用单类逻辑回归（OCLR）机器学习算法确定，随后计算干性特征权重值与每个样本基因表达水平之间的相关系数，最后通过将Spearman相关系数缩放到0到1之间得到干性指数。差异表达分析使用Wilcoxon秩和检验，错误发现率（FDR）< 0.05和|log2(倍数变化)| > 1的基因被认为具有统计学意义。为增强风险模型的准确性，应用更严格的选择阈值，设定FDR < 0.01和|log2(倍数变化)| > 2。

缺氧特征评分从cBioPortal数据库获取，缺氧相关基因使用加权基因共表达网络分析（WGCNA）识别。mRNAsi相关差异表达基因（DEGs）与缺氧相关基因之间的重叠基因被共同定义为干性和缺氧相关基因（SHRGs）。使用“ConsensusClusterPlus”R包通过无监督共识聚类对SHRGs进行分类，聚类过程进行1000次以80%重采样以增强分类稳定性。

功能富集分析使用“clusterProfiler”R包进行GO、KEGG和GSEA分析，而GSVA分析使用无监督“GSVA”R包进行。GSEA和GSVA的背景基因集从分子签名数据库（MSigDB）获取。使用“maftools”R包分析TCGA-LIHC患者的CNV和体细胞突变数据以检查不同簇之间的遗传改变。使用PROGENy算法比较不同簇之间的14条致癌通路。使用“CIBERSORT”工具量化肿瘤浸润免疫细胞的丰度以评估肿瘤免疫微环境（TIME）。

为评估TCGA-LIHC患者的免疫治疗反应，分析了先前研究中确定的68个免疫检查点相关基因的表达。随后，从TIDE门户检索肿瘤免疫功能障碍和排除（TIDE）评分、T细胞功能障碍评分、T细胞排除评分、INFG水平和MDSC水平。然后应用单样本基因集富集分析（ssGSEA）计算29个免疫相关特征的富集分数并探索指数与免疫调节之间的关联。此外，计算了来自MSigDB的三个干细胞相关基因集的ssGSEA衍生富集分数以研究指数与干性之间的关联。

为确定SHRGs与患者生存结局之间的关系，应用单变量Cox回归、LASSO回归和随机森林模型过滤SHRGs。排除绝对相关性为0.8的属性后，总共选择419个基因作为输入变量。最后，确定四个最关键基因并纳入多变量Cox回归模型以构建风险预测模型，称为SHRPI。该模型的公式如下：风险评分 = ∑(系数(i) × 表达(i))。

为评估SHRPI的稳健性，患者最初根据中位SHRPI值分为两组。使用Kaplan-Meier生存分析评估SHRPI的预后意义。通过接收者操作特征（ROC）曲线分析进一步评估SHRPI的预测准确性，使用“timeROC”R包计算曲线下面积（AUC）。为增强其临床适用性，TCGA-LIHC患者根据最佳截断值进一步分为高危（HRG）和低危（LRG）组，随后进行全面的免疫分析和药物敏感性分析。

将SHRPI评分与肿瘤分期、年龄、性别、肿瘤分级、血管侵犯状态、Child-Pugh分级、肝脏炎症状态、肝硬化状态、复发状态、BMI和AFP水平一起纳入单变量Cox回归分析。使用“survival”R包中的Cox比例风险回归模型计算每个变量的风险比（HRs）。为确定独立预后因素，进行多变量Cox回归分析，并使用“RMS”R包基于发现开发列线图。通过Schoenfeld残差和偏差残差评估模型稳定性。通过ROC分析、校准曲线和C指数评估列线图的预测性能，通过1000次自助重采样迭代计算。此外，应用决策曲线分析（DCA）评估预测模型的临床适用性。

从癌症药物敏感性基因组学（GDSC2 v8.5）、癌症治疗反应门户（CTRP v2.0）和PRISM重定位数据集（20Q2）获取基因表达数据以及相应半数最大抑制浓度（IC50）值和剂量反应曲线下面积（AUC）。使用“oncoPredict”R包预测每个样本的药物敏感性。从PubChem Compound检索化合物的分子结构，从PDB获取G6PD（PDB ID: 7UAG）的3D坐标。所有蛋白质和分子文件转换为PDBQT格式，去除水分子并添加极性氢原子以提高对接准确性。使用AutoDock Vina 1.2.2进行分子对接模拟，使用PyMol可视化所得蛋白质-配体复合物。结合能用于评估每个化合物的治疗潜力。

从GEO数据集GSE149614过滤HCC样本的单细胞RNA测序数据以去除低质量细胞。使用Seurat（v5.3.0）对基因表达进行对数标准化，随后通过PCA进行降维并通过FindNeighbors和FindClusters函数进行聚类。使用“SingleR”R包注释细胞簇。SHRPI计算为HMMR、UBE2S、G6PD和NEIL3的z得分标准化表达的加权和。根据中位SHRPI评分将细胞分为低和高SHRPI组。在t-SNE图上可视化SHRPI分布，并在单独的点图中显示每个组成基因跨细胞簇的表达。

人HCC细胞系（HuH-7、PLC/PRF/5、Hep-3B和Li-7）和HEK293从中国科学院上海细胞生物学研究所获取。所有细胞的真实性通过短串联重复（STR）分析确认。细胞在适当培养基中培养，补充10%胎牛血清和1%青霉素-链霉素，在37°C、5% CO₂培养箱中维持。Hep-3B和PLC/PRF/5细胞在MEM培养基中维持，HEK293和HuH-7细胞在高葡萄糖DMEM中培养，Li-7细胞在RPMI 1640培养基中维持。所有细胞系定期通过PCR确认无支原体。慢病毒载体构建和感染的详细信息在补充文件中提供。

将细胞（2×10⁴）接种在聚碳酸酯膜插入物的上室中，使用200 μL无FBS培养基。下室填充800 μL含20% FBS的培养基。孵育24-48小时后，迁移通过膜的细胞洗涤，用1%多聚甲醛固定并用结晶紫染色。拍摄四个随机选择视野的照片并计数迁移细胞数。实验重复三次。

将细胞（2×10³）铺板到6孔超低附着板上，并在特殊培养基中培养，该培养基由DMEM/F12组成，补充4 μg/mL胰岛素、B27、20 ng/mL EGF和20 ng/mL碱性FGF。孵育10天后，在显微镜下拍摄直径大于75 μm的球体并计数。

Western blot分析、共免疫沉淀和定量PCR分析的详细信息在补充文件中提供。抗体的具体细节在补充表4中呈现。引物序列在补充表5中提供。

在雄性野生型C57BL/6J小鼠（6-8周）中通过水动力尾静脉注射共过表达激活AKT和c-Met生成转基因HCC小鼠模型。简而言之，将质粒m-G6PD pT3-EF1α-MYC或pT3-EF1α-MYC（MCS）（20 μg）、pT3-myr-AKT-HA（20 μg）和pT3-EF1α-c-Met（20 μg）与pCMV(CAT)T7-SB100（2.4 μg）以12.5:12.5:12.5:1.5的比例稀释在2 ml盐水（0.9% NaCl）中，通过0.22 μm过滤器过滤并在5-7秒内注射到小鼠侧尾静脉中。注射后20天，通过腹腔注射戊巴比妥钠（150 mg/kg）人道安乐死小鼠。随后收集肝肿瘤进行分析。动物实验严格按照相关指南进行并经浙江实验动物中心机构动物护理和使用委员会（IACUC, ZJCLA；批准号：ZJCLA-IACUC-20011186）批准。本研究遵循ARRIVE指南。

学生t检验或Wilcoxon秩和检验（Mann-Whitney U检验）用于评估连续分布数值数据。根据数据分布使用Pearson或Spearman相关检验进行相关分析。使用Kaplan-Meier方法生成生存曲线并使用对数秩检验比较。所有统计分析在GraphPad Prism（v9.0）和R（v4.4.1）中进行。双尾P值< 0.05被认为具有统计学意义。统计学显著性注释如下：P < 0.05，P < 0.01，P < 0.001，P < 0.0001；ns = 不显著。

3 Results

3.1 Stemness indices and hypoxia scores in HCC

使用TCGA-LIHC队列的mRNA表达和DNA甲基化数据计算了五个干性指数，从cBioPortal数据库获取缺氧评分。基于mRNA表达的干性指数（mRNAsi）和Buffa缺氧评分与HCC中的OS显著相关。具有较高mRNAsi的患者表现出显著较差的肿瘤分化、增加的血管侵犯和升高的AFP水平。类似地，较高的Buffa缺氧评分与晚期肿瘤分期、较差的肿瘤分化、增加的血管侵犯和升高的AFP水平相关。因此，选择mRNAsi量化干性特征，而Buffa缺氧评分用于评估肿瘤缺氧水平。

3.2 Identification of stemness- and hypoxia-related clusters in HCC

基于中位mRNAsi值将HCC患者分为高和低干性组， resulting in the identification of 1,341 mRNAsi相关DEGs。WGCNA揭示了11个非灰色模块，其中蓝色模块与Buffa缺氧评分显示出最强相关性。整合mRNAsi相关DEGs与缺氧相关基因确定了75个重叠基因， classified as SHRGs。对75个SHRGs进行共识聚类，确定最佳分类为k = 2。HCC患者随后分为两个簇（簇1和簇2）。簇2中的患者表现出显著较短的中位OS和较低的生存概率，以及较高的mRNAsi和缺氧评分。

进行更严格的差异表达分析使用|log2FC| > 2和FDR < 0.01作为选择标准。通过GO和KEGG对DEGs的功能富集分析揭示了细胞周期调控和DNA修复通路的显著富集。GSEA进一步证明，与簇1相比，簇2表现出E2F靶标、G2/M检查点和KRAS信号DN的显著激活，而氧化磷酸化、胆汁酸代谢、脂肪酸代谢和脂肪生成被抑制。GSVA揭示了增殖相关通路的显著上调，包括簇2中的E2F靶标、G2/M检查点、DNA修复和mTOR信号。相反，簇1主要富集于脂质和胆汁酸代谢、凝血和炎症反应。总体而言，这些结果表明簇2以高增殖表型为特征，具有增强的DNA修复和致癌信号，而簇1与代谢重编程和炎症肿瘤微环境相关。

3.3 Genomic and TIME characteristics of the two stemness- and hypoxia-related clusters

检查了两个簇中的体细胞突变和CNV以研究其不同预后的潜在机制。在几个基因中检测到 recurrent mutations，包括TP53、TTN和CTNNB1，并且两个簇表现出不同的单核苷酸变异（SNV）替代模式。TP53、LRP1B、RB1、ABCB5和ZNF469在簇之间显示出显著不同的突变频率。此外，簇2表现出较高的非整倍性评分、肿瘤突变负荷（TMB）和同源重组缺陷（HRD） compared with簇1。这些发现表明具有 elevated DNA损伤的肿瘤可能 possess enhanced免疫逃避能力 and reduced对免疫治疗的反应性。

为进一步验证分类，将患者簇与先前建立的分子分类进行比较，其中“炎症”亚型（C3）与最佳预后相关。大多数C3患者被分配到预后较好的簇1，与先前发现一致。

鉴于干性、缺氧和免疫相关通路之间复杂的相互作用，进一步探索了两个簇之间TIME的差异。CIBERSORT分析显示，簇2中的患者表现出显著较高水平的活化记忆CD4⁺ T细胞、滤泡辅助T细胞和Tregs，而M0巨噬细胞、静息记忆CD4⁺ T细胞和M2巨噬细胞显著减少。接下来，分析了两个簇中先前报告的免疫检查点相关基因的表达。在簇2中，已知抑制T细胞免疫活性的基因，包括CTLA4及其配体以及PD-1及其配体，显著上调。随着免疫治疗在HCC治疗中的日益突出，使用TIDE模型评估患者的潜在反应。鉴于较高的TIDE评分表明增加的免疫逃避和减少的免疫治疗功效，观察到簇2具有显著高于簇1的TIDE评分，表明簇2中的患者可能 benefit less from免疫治疗。

使用ssGSEA评估治疗签名和29个免疫相关基因签名，涵盖免疫、基质和其他细胞过程。首先，簇2表现出与细胞周期、DNA复制和错配修复相关的基因签名的显著上调。其次，簇2显示免疫抑制细胞浸润增加，包括肿瘤相关巨噬细胞（TAMs）、MDSCs和Tregs，以及增强的肿瘤细胞增殖，导致高度免疫抑制和升高的促肿瘤免疫评分。最后，进行PROGENy分析以评估癌症相关信号通路的活性。结果表明，簇2在雌激素、缺氧、MAPK、NF-κB、p53、TNFα和WNT信号通路中表现出显著较高的活性，而簇1在VEGF信号通路中显示出相对较高的活性。

3.4 Construction of the 4-gene SHRPI and evaluation of its prognostic significance

对两个簇之间的DEGs进行Cox回归分析以识别预后基因。为减轻共线性效应，排除Pearson相关系数大于0.80的基因， resulting in 221 DEGs。然后应用LASSO回归分析， yielding五个稳定预后基因。为进一步降低假阳性率并增强模型准确性，使用随机森林模型选择平均减少基尼大于1的基因。随后，四个重叠关键基因纳入Cox回归模型， with corresponding coefficients used to construct the SHRPI：SHRPI = 0.10211669 × HMMR + 0.15981839 × UBE2S + 0.20537184 × G6PD + 0.06132584 × NEIL3。

基于中位SHRPI评分将患者分为低和高SHRPI组。差异表达分析显示所有四个关键基因在高SHRPI组中上调。在TCGA-LIHC队列中，高SHRPI组中的患者表现出显著较短的OS compared with低SHRPI组。时间依赖性ROC曲线分析证明SHRPI对生存表现出稳健且稳定的预测性能，1年、3年和5年OS的AUC值分别为0.82、0.70和0.67。单变量和多变量Cox回归分析确认分期、复发状态和SHRPI是HCC中OS的独立风险因素。在NC-LT队列中，高SHRPI组中的患者表现出显著较短的RFS compared with低SHRPI组。单变量和多变量Cox回归分析确认肿瘤直径、AFP水平和SHRPI是LT后HCC中RFS的独立风险因素。此外，SHRPI的适用性