Combining Single-Cell and Transcriptomic Data Revealed the Prognostic Significance of Glycolysis in Pancreatic Cancer

Abstract

Background: Pancreatic cancer (PC), the most common fatal solid malignancy, has a very dismal prognosis. Clinical computerized tomography (CT) and pathological TNM staging are no longer sufficient for determining a patient's prognosis. Although numerous studies have suggested that glycolysis is important in the onset and progression of cancer, there are few publications on its impact on PC. Methods: To begin, the single-sample gene set enrichment analysis (ssGSEA) approach was used to quantify the glycolysis pathway enrichment fraction in PC patients and establish its prognostic significance. The genes most related to the glycolytic pathway were then identified using weighted gene co-expression network analysis (WGCNA). The glycolysis-associated prognostic signature in PC patients was then constructed using univariate Cox regression and lasso regression methods, which were validated in numerous external validation cohorts. Furthermore, we investigated the activation of the glycolysis pathway in PC cell subtypes at the single-cell level, performed a quasi-time series analysis on the activated cell subtypes and then detected gene changes in the signature during cell development. Finally, we constructed a decision tree and a nomogram that could divide the patients into different risk subtypes, according to the signature score and their different clinical characteristics and assessed the prognosis of PC patients. Results: Glycolysis plays a risky role in PC patients. Our glycolysis-related signature could effectively discriminate the high-risk and low-risk patients in both the trained cohort and the independent externally validated cohort. The survival analysis and multivariate Cox analysis indicated this gene signature to be an independent prognostic factor in PC. The prognostic ROC curve analysis suggested a high accuracy of this gene signature in predicting the patient prognosis in PC. The single-cell analysis suggested that the glycolytic pathway may be more activated in epithelial cells and that the genes in the signature were also mainly expressed in epithelial cells. The decision tree analysis could effectively identify patients in different risk subgroups, and the nomograms clearly show the prognostic assessment of PC patients. Conclusion: Our study developed a glycolysis-related signature, which contributes to the risk subtype assessment of patients with PC and to the individualized management of patients in the clinical setting.

Keywords: glycolysis; immune infiltration; pancreatic cancer; prognosis; single-cell.

FIGURE 1

Flow chart of our study.

FIGURE 2

ssGSEA analysis and weighted gene correlation network analysis (WGCNA). (A) ssGSEA analysis showed that the glycolysis score was obvious elevated in the dead PC patients. (B) Survival analysis revealed that the high-glycolysis group has a worse prognosis. (C) Best soft threshold of WGCNA was 7. (D) WGCNA analysis found 18 no-gray gene modules. (E) Correlation between the modules and glycolysis. The black and red modules had the strongest correlation with the glycolytic phenotype (Cor = 0.5 and p < 0.001). (F,G) Relation between module membership and gene significance in red and black modules. (H) Correlation between red and black modules and glycolysis.

FIGURE 3

Gene signature was constructed in TCGA cohort. (A,B) LASSO Cox regression was used to identify the most important genes, and the optimal lambda was 0.111. (C) GLCS was obviously elevated in the dead PC patients (p = 2.3E-10). (D) Survival analysis reveals that GLCS-high has a worse prognosis (p < 0.0001). (E) Multivariate Cox analysis reveals that GLCS was an independent prognostic factor (p < 0.001). (F) AUC of GLCS and clinical features. The AUC value of GLCS was higher than that of other clinical features.

FIGURE 4

Assessment of the gene signature in extra validation cohorts. (A) Survival analysis in the PACA-AU cohort suggested that the prognosis of the GLCS-high group was worse (p = 0.0051). (B) Survival analysis in the PACA-CA cohort suggested that the prognosis of the GLCS-high group was worse (p < 0.001). (C,D) Multivariate Cox analysis in the PACA-AU cohort and PACA-CA cohort revealed that GLCS was an independent prognostic factor. (E,F) AUC of GLCS and clinical features in the PACA-AU cohort and PACA-CA cohort.

FIGURE 5

Gene signature is suitable for different clinical patients. (A) Among the TI and TII stage PC patients, the prognosis of the high-GLCS group was worse (p = 0.002). (B) Among the TIII&TIV stage PC patients, the prognosis of the high-GLCS group was worse (p < 0.001). (C) Among the PC patients with N0, the prognosis of the high-GLCS group was worse (p < 0.001). (D) Among the PC patients with N1, the prognosis of the high-GLCS group was worse (p < 0.001). (E) Among PC patients with stage I, the prognosis of the high-GLCS group was worse (p = 0.003). (F) Among PC patients with stage II, the prognosis of the high-GLCS group was worse (p < 0.001). (G) Among PC patients with age<=65, the prognosis of the high-GLCS group was worse (p < 0.001). (H) Among PC patients with age>65, the prognosis of the high-GLCS group was worse (p < 0.001). (I) Among male PC patients, the prognosis of the high-GLCS group was worse (p < 0.001) (J) Among female PC patients, the prognosis of the high-GLCS group was worse (p < 0.001).

FIGURE 6

Exploration of the relation between GLCS and immune infiltration. (A) Immune landscape of PC patients. (B) Difference in immune infiltration levels between the two groups in the form of a box plot. (C) Differences in the immune checkpoint gene expression between high-risk and low-risk groups. (D) Differences in the expression of immunogenic cell death genes between high-risk and low-risk groups (*p < 0.05, **p < 0.01, and ***p < 0.001).

FIGURE 7

Exploration of the relation between GLCS and tumor mutation. (A) Landscape of genetic mutations in PC patients. (B) Mutation symbiosis among the top mutation genes. (C) TMB were higher in the GLCS-high group.

FIGURE 8

Single-cell analysis. (A) Fifteen samples were integrated with the CCA method. (B) Dimension reduction and cluster analysis. The cell types were shown with the umap plot. (C) Glycolysis pathway was activated in different cell types (D–K). Genes in the prognostic signature expressed differently in different cell types. (L–N) Analysis of epithelial cell locus differentiation and (O) the genes in the signature expressed differently during the epithelial cell locus differentiation.

FIGURE 9

Clinical implications of the signature. (A) Decision tree analysis could divide PC patients into five risk subtypes. (B) Five risk subtypes have different prognosis, and cluster 1 has a best prognosis. (C) Nomogram analysis showed the 1-, 3-, and 5-year mortality rates of patient TCGA-S4-A8RP. (D,E) 2- and 3-year calibration curves of the nomogram. (F) AUC of the nomogram and other clinical features to evaluate the prognosis of PC.

FIGURE 10

qRT-PCR to verify the expressions of seven model genes in PC. The seven model genes were all upregulated in PC compared with the normal adjacent tissues (*p < 0.05, **p < 0.01, and ***p < 0.001). (A) MET; (B) FAM25A; (C) LY6D; (D) FAM111B; (E) ITGB6; (F) CENPE; (G) KCTD14.