基于糖酵解相关基因模型的乳腺癌患者预后及免疫功能综合分析

Comprehensive analysis of prognosis and immune function of breast cancer patients based on glycolysis related gene model

目的探讨肿瘤有氧糖酵解对乳腺癌患者预后及其免疫微环境的影响。方法从癌症基因组图谱(TCGA)数据库中获取1988年1月至2013年12月1 097例乳腺癌组织及112例癌旁正常组织的mRNA表达谱和相应临床数据。利用基因集富集分析(GSEA)软件筛选具有显著差异性的糖酵解相关基因(GRG)。利用单/多因素Cox回归分析及最小绝对收缩和选择算法(LASSO)筛选出与OS相关的GRG后生成6-GRG模型并绘制列线图。按模型算得中位风险值划分高/低风险组,提取不同组患者的相关GRG表达量后生成热图并进行了Kaplan-Meier生存分析。受试者操作特征(ROC)曲线分析及校准图用以验证模型准确性。用实时荧光定量PCR分析三阴性乳腺癌细胞系MDA-MB-231及正常人乳腺上皮细胞系MCF-10A中相关GRG的表达量。最后将不同风险组的差异表达基因(DEG)进行基因本体论(GO)及京都基因与基因组百科全书(KEGG)分析及免疫微环境分析。结果共有6个GRG与OS相关,包括CACNA1H(HR=1.12,95%CI:1.03~1.20,P=0.009)、SDC1(HR=1.22,95%CI:1.06~1.40,P=0.008)、SDC3(HR=0.74,95%CI:0.61~0.90,P=0.003)、NUP43(HR=1.42,95%CI:1.07~1.90,P=0.015)、PGK1(HR=1.74,95%CI:1.36~2.20,P<0.001),CHST1(HR=1.14,95%CI:1.03~1.30,P=0.009),从而构建了6-GRG预后模型。根据模型测算的中位风险值,分为低风险组(n=576)和高风险组(n=521)。基因表达热图发现高风险组中SDC3表达减少,其余5个基因表达增加。生存分析发现低风险组患者的总生存率高于高风险组(χ^(2)=7.314,P<0.001),低风险组与高风险组患者的总生存期分别为(17.23±0.89)年和(8.75±1.71)年。校准图显示该模型对患者1、3、5年OS的预测准确率曲线与理想曲线(45度角灰线)贴合。ROC曲线显示,乳腺癌患者1、3、5年OS的曲线下面积分别为0.68、0.71和0.72,提示本模型具有较好的预测精准性。相较于正常乳腺细胞MCF-10A,乳腺癌MDA-MB-231细胞中CACNA1H、SDC1、NUP43、CHST1、PGK1均表达增高(F=15.36、30.73、1.08、14.92、12.93,P均<0.050),而SDC3的表达减少(F=17.50,P=0.038)。GO及KEGG结果显示不同风险组中DEG多在免疫学功能或通路显著富集。低风险组患者的正性免疫细胞比例高于高风险组,包括初始B细胞[0.037(0.011,0.095)比0.031(0.004,0.069),Z=-3.928,P=0.012]、记忆性B细胞[0.010(0.002,0.104)比0.004(0.001,0.411),Z=-5.175,P<0.001]、CD8^(+)T细胞[0.103(0.077,0.329)比0.073(0.012,0.136),Z=-4.904,P<0.001]、滤泡辅助性T细胞[0.068(0.000,0.117)比0.057(0.001,0.128),Z=-2.363,P<0.001]、γδT细胞[0.017(0.000,0.180)比0.010(0.000,0.140),Z=-1.491,P=0.001]、活化自然杀伤细胞[0.031(0.000,0.141)比0.021(0.000,0.099),Z=-1.667,P<0.001]、单核细胞[0.017(0.000,0.101)比0.015(0.000,0.085),Z=-1.093,P=0.047]、中性粒细胞[0.048(0.011,0.122)比0.021(0.008,0.069),Z=-2.776,P<0.001]。结论6-GRG预后模型具有良好的预测效能,乳腺癌患者的高糖酵解水平与不良预后及抗肿瘤免疫的功能下降密切相关。

Objective To investigate the effect of tumor aerobic glycolysis on the prognosis and immune microenvironment of breast cancer patients.Methods The mRNA expression profiles of 1 097 cases of breast cancer tissues and 112 cases of normal adjacent tissues and corresponding clinical data were obtained from The Cancer Genome Atlas (TCGA) database from January 1988 to December 2013. The gene set engineering analysis (GSEA) software was used to screen out glycolysis-related genes (GRGs) with a significant difference. The single/multiple factor Cox regression analysis and the least absolute contraction and selection algorithm (LASSO) regression method were used to find GRGs related to OS, and then the 6-GRG model was generated and the nomogram was drawn. The patients were divided into high risk and low risk groups according to the median risk value of the model. The expression of related GRGs in different groups was extracted to generate heat maps and the Kaplan-Meier method was conducted for survival analysis. The receiver operating characteristic (ROC) curve analysis and calibration chart were used to verify the accuracy of the model, and single/multiple factor regression analysis is used to judge the independence of the model. The validity of the model was also verified by RT qPCR cell experiment. Finally, GO, KEGG analysis and immune microenvironment analysis were conducted for differentially expressed genes(DEGs) of different risk groups.Results Six GRGs were significantly related with OS, including CACNA1H (HR=1.12, 95%CI:1.03-1.20,P=0.009), SDC1 (HR=1.22, 95%CI:1.06-1.40,P=0.008), SDC3 (HR=0.74, 95%CI:0.61-0.90,P=0.003), NUP43 (HR=1.42, 95%CI:1.07-1.90,P=0.015), PGK1 (HR=1.74, 95%CI: 1.36-2.20,P<0.001), CHST1(HR=1.14, 95%CI: 1.03-1.30,P=0.009). Thus, we established a prognostic model of 6-GRG. According to the median risk value measured by the model, the patients were divided into low-risk group (n=576) and high-risk groups (n=521). The heat map showed that the expression of SDC3 in high-risk group decreased, and the expression of other 5 genes increased. The survival analysis showed that the overall survival rate of patients in low-risk group was significantly higher than that in high-risk group (χ^(2)=7.314,P<0.001), and the overall survival of patients was (17.23±0.89) years in low-risk group and (8.75±1.71) years and high-risk group. The calibration plot shows that the predictive curve of the model for 1-, 3- and 5-year OS fitted the ideal curve (gray line at 45-degree angle). The ROC curve showed that the area under the curve (AUC)of 1-, 3-and 5-year OS was 0.68, 0.71 and 0.72, respectively, indicating a good prediction accuracy of this model. Compared with normal breast cell line MCF-10A, the expression of CACNA1H, SDC1, NUP43, CHST1 and PGK1 in breast cancer cell line MDA-MB-231 significantly increased (F=15.36, 30.73, 1.08, 14.92, 12.93, allP<0.050), while the expression of SDC3 significantly decreased (F=17.50,P=0.038). The GO and KEGG analysis indicated that most of DEGs in different risk groups were significantly enriched in immunological functions or pathways. The patients in low-risk group had significantly higher proportion of positive immune cells compared with high-risk group, including initial B cells [0.037 (0.011, 0.095)vs0.031 (0.004, 0.069),Z=-3.928,P=0.012], memory B cells [0.010 (0.002, 0.104)vs0.004 (0.001, 0.411),Z=-5.175,P<0.001], CD8^(+)T cells [0.103 (0.077, 0.329)vs0.073 (0.012, 0.136),Z=-4.904,P<0.001], follicular helper T cells [0.068 (0.000, 0.117)vs0.057(0.001, 0.128),Z=-2.363,P<0.001], γδ T cells [0.017 (0.000, 0.180)vs0.010 (0.000, 0.140)Z=-1.491,P=0.001], activated natural killer cells [0.031 (0.000, 0.141)vs0.021 (0.000, 0.099),Z=-1.667,P<0.001], monocytes [0.017 (0.000, 0.101)vs0.015 (0.000, 0.085),Z=-1.093,P=0.047], and neutrophils [0.048 (0.011, 0.122)vs0.021 (0.008, 0.069),Z=-2.776,P<0.001].Conclusions The 6-GRG prognostic model has good predictive performance. The high glycolysis level in breast cancer patients is closely related to poor prognosis and decreased anti-tumor immunity.