The value of a metabolic reprogramming-related gene signature for pancreatic adenocarcinoma prognosis prediction

Identification of DEGs from public datasets

The results of this study are summarized in a flow chart (Figure 1). The detailed information in three eligible PDAC datasets in the GEO database (GSE15471, GSE16515, and GSE32676) met our criteria. Analysis of these PDAC datasets revealed that 237 differentially expressed genes (DEGs) were shared among the 3 series of comparisons between tumor and adjacent para-cancerous tissues in the Venn and UpSet diagrams; these DEGs were regarded as credible DEGs (Figure 2A, 2B). In addition, KEGG pathway enrichment analyses were performed for these overlapping up- or down-regulated genes. As shown in Figure 2C, the DEGs were most enriched in the ECM-receptor interaction pathway, the PI3K-Akt signaling pathway and numerous energy metabolism-related pathways, such as the proteoglycans in cancer pathway and the central carbon metabolism in cancer pathway.

Prognostic signature construction

Based on these 237 DEGs and clinical features from the TCGA database, 100 genes significantly related to OS time (P < 0.05) were identified by the log-rank test, as shown in Figure 3A. Twelve glycolysis and lipid metabolism-related genes significantly associated with survival of PDAC were screened (Table 1). To avoid overfitting during the subsequent model construction, LASSO regression was performed to build a prognostic signature, which included twelve genes from the 237 previously identified DEGs. The LASSO coefficient profiles of the twelve genes are presented in Figure 3B. A formula to calculate the risk score for PDAC was derived based on individual three-gene expression levels weighted by regression coefficients, as follows:

$\begin{array}{c}\text{Risk}\text{score}=(1.0167070\times \text{expression level of MET})\\ -(0.4313794\times \text{expression level of ENO3})\\ -(0.1764747\times \text{expression level of CD36}).\end{array}$

With the median risk score as the cutoff, PDAC patients were classified into a low-risk group (n = 89) and a high-risk group (n = 88) (Figure 3C). Survival analysis confirmed that the survival times of the low risk patients were distinctly longer than those of the high-risk patients in the TCGA discovery cohort (P < 0.0001; Figure 3D). The finding was subsequently validated in ICGC validation datasets; the results obtained for these datasets agreed with those obtained for the discovery cohort (P = 0.038; Figure 3F). The same trend was also verified in the FUSCC validation cohort (P = 0.017; Figure 3H). Furthermore, to evaluate the performance of the three-gene signature in predicting the prognoses of PDAC patients, time-dependent ROC curves with respect to the 1-year, 3-year and 5-year survival rates were constructed for the TCGA, ICGC and FUSCC datasets, and the area under the curve (AUC) values for these curves were evaluated (Figure 3E, 3G, 3I). The AUC at 1, 3, and 5 years was 0.709, 0.749, and 0.836 in the TCGA set, respectively. The AUC at 1, 3, and 5 years was 0.644, 0.621, and 0.955 in the ICGC set, respectively. The AUC at 1 and 3 years was 0.601, 0.745 in the FUSCC set, respectively.

Validation of the expression of and alterations in the three genes

The mRNA expression levels of the three genes were validated using the Gene Expression Profiling Interactive Analysis (GEPIA) database and GEO. In the GEPIA dataset, only MET mRNA expression was significantly elevated in tumor samples (Figure 4A–4C). In the GEO dataset, we found that MET expression levels were significantly higher in PDAC tumor tissue than in paracancerous tissue in the GEO32676 and GEO15471 datasets (Figure 4D). In contrast, the mRNA expression of ENO3 and CD36was markedly down-regulated in tumor tissues (Figure 4E, 4F). Protein expression of the DEGs in pancreatic tumor and non-tumor tissues was evaluated by the Human Protein Atlas (https://www.proteinatlas.org/) [12]. Immunohistochemical (IHC) staining information for the three genes in PDAC from the Human Protein Atlas database is shown in Figure 4G–4I. Moreover, the OS times were also obtained from the GEPIA database (Figure 4J–4L). A total of 848 patients in four cohorts (the ICGC, Queensland Centre for Medical Genomics [QCMG], TCGA and UT Southwestern Medical Center [UTSW] cohorts) were included in this study for analysis of the mutation information of the three glycolysis- and lipid metabolism-related genes. Overall, amplification and deep deletion were the most common types of mutations in MET, CD36 and ENO3 (Figure 4M).

Validation of the risk score formula in the GEO and TCGA cohorts

ROC curve analysis was performed with AUCs to evaluate the usefulness of the risk score of the gene signature in distinguishing PDAC tissues from control samples in four GEO cohorts (Figure 5A). Tumor tissues could be reliably identified based on this risk score. Subgroup analyses were performed to assess the three-gene signature and clinical characteristics of pancreatic cancer (including American Joint Committee on Cancer [AJCC] stage and grade) from appropriate datasets. In terms of AJCC stage and grade, stage T3-4 and grade 3-4 patients had higher risk scores than stage T1-2 and grade 1-2 patients. However, the same trend was not noted for N stages (Figure 5B–5D). The Kaplan-Meier curves of these clinical features are shown in Figure 5E–5H. The survival time of the high-risk score group was significantly shorter than that of the low-risk score group for all pathological T stages and for grades 1-2. In terms of mutations, risk scores were observed to be highly correlated with the mutation states of key genes in PDAC. The risk scores for the CDKN2A, KRAS and TP53 mutant groups were significantly higher than that of the wild-type group (Figure 5I–5L).

Correlations between the three-gene signature and clinical characteristics

We examined the association of the three-gene signature (risk score) with clinical features in PDAC patients using univariate and multivariate Cox proportional hazard regression analyses. Univariate Cox proportional hazard regression showed that N stage, age, and risk score could predict poor survival of PDAC patients, as shown in Table 2 (P < 0.05). In addition, multivariate Cox proportional hazard regression showed that the risk score (P < 0.001) was an independent prognostic indicator of PDAC (Figure 6A). Based on the results acquired from multivariate Cox regression of OS in the TCGA dataset, we developed a nomogram to predict 1-, 3- and 5-year survival probability in PDAC (Figure 6B). The C-index of our nomogram in the TCGA cohort was 0.689 (95% confidence interval, 0.628-0.749), and the calibration curves for this nomogram presented good agreement between the possibility estimated by the nomogram and the actual proportion (Figure 6C–6E).

Analysis of gene set enrichment relevance of the three-gene signature

To elucidate the molecular mechanism of the three-gene signature, patients from the TCGA PDAC dataset were divided into high- and low-risk groups according to the cutoff value. GSEA compared the high- and low-risk groups. Gene ontology (GO) gathers information about molecular function, biological processes, and cellular components in a number of different organisms. These enriched gene ontology (GO) data revealed that molecular alteration in the high-risk group was closely related to the nuclear periphery, apical junction and nuclear matrix. Oncological signatures were significantly enriched, including HOXC6, SMAD2, SMAD3 and MTORC (Figure 7).

Acquisition of PDAC gene expression data and clinical data

Microarray gene expression profiles for PDAC were downloaded from the NCBI Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/), TCGA (http://cancergenome.nih.gov/) and the International Cancer Genome Consortium (ICGC) database (https://dcc.icgc.org/). The gene expression microarray datasets GSE15471, GSE16515 and GSE32676 were chosen and downloaded for differentially expressed gene (DEG) analysis. The datasets met the following criterion: they contained data on human pancreatic tissue samples from tumors and adjacent paracancerous areas. Datasets that were unanalyzable, that included few DEGs (fewer than 100), or that included defectively annotated genes were excluded. The probe names were transformed into gene symbols using the annotation files supplied by the manufacturer. The median ranking value was used as the expression value if several probes matched a single gene. Robust multiarray average (RMA)-normalized data were log2-transformed for further analysis. Normalized RNA-sequencing data in transcripts per million (TPM) format and the associated clinical information of the PDAC samples were downloaded from the TCGA dataset. For validation, we used transcriptomic data from the ICGC and recruited 82 patients who underwent surgical resection at Fudan University Shanghai Cancer Center (FUSCC) from 2010 to 2012. These patients were diagnosed according to strict pathological criteria, and total RNA was extracted from their tumor tissues. Strictly controlled postoperative follow-ups were conducted for all of the patients. All of the procedures were implemented after obtaining approval from the Clinical Research Ethics Committee of FUSCC, and informed consent was acquired from each patient.

DEG identification and bioinformatic analysis

The “limma” package of R software (version 3.6.2) was used to establish a prognostic gene signature by identifying the DEGs from the mRNA data. A |log2(fold change [FC]) | > 1, a P-value < 0.05, and a false discovery rate (FDR) < 0.05 were considered the cutoffs for DEGs between tumor and nontumor samples in the GEO database. Next, the DEGs that overlapped among 3 series, as determined with Venn and UpSet diagrams, were regarded as credible DEGs. Gene ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were applied to explore the potential biological roles of the DEGs. Significantly relevant signaling pathways were identified with the Database for Annotation, Visualization, and Integrated Discovery (DAVID) (https://david.ncifcrf.gov/) [26]. P < 0.05 was considered the threshold for statistical significance.

Establishment of the LASSO regression model and calculation of risk score

We used log-rank tests to examine the associations between gene expression and OS in the TCGA dataset. DEGs with P-values < 0.05 were considered statistically significant and were included in the following analyses. Genes obtained from the Molecular Signature Database (MSigDB) gene sets “REACTOME_GLYCOLYSIS” and “REACTOME_METABOLISM_OF_LIPIDS_AND_LIPOPROTEINS” were used as glycolytic and lipid metabolism-related genes, respectively [7, 8, 27]. Twelve glycolysis and lipid metabolism-related genes were selected to construct LASSO Cox regression analysis models with the R package “glmnet”. LASSO analysis is the most popular method for evaluating survival data and is especially appropriate for analyzing gene expression profiles. The “glmnet” package returned a series of models. For each model, the tuning parameter λ was inversely correlated with the complexity of the model and the value of deviance. Ten-fold cross-validation was performed to screen the best model from the series (the model with the minimum mean cross-validation error). The patients were divided into high-risk and low-risk groups according to a cutoff median risk score. Kaplan-Meier analysis and time-dependent receiver-operating characteristic (ROC) curve analysis were used to test predictive differences between the high- and low-risk groups in the discovery and two validation datasets. Cox regression analysis was employed to test whether the risk score was an independent prognostic factor. A nomogram and related calibration curves were created based on the TCGA cohort for further clinical application.

Gene set enrichment analysis (GSEA) and mutation analyses

GSEA was applied to elucidate the molecular mechanisms of the prognostic gene sets [28]. The TCGA dataset was divided into high- and low-risk groups according to the cutoffs. GSEA was performed in the Java program GSEA, version 3.0. The MSigDB v. 6.2. C2 collection (curated gene sets), C5 collection (GO gene sets), and H collection (hallmark gene sets) were searched to identify enriched KEGG pathways and biological information and dysfunctional oncogenic signatures associated with the high-risk group. Results with a |normalized enrichment score (NES)| > 1 and an FDR < 0.05 were considered statistically significant. Somatic mutations (SNPs and small indels) in TCGA were obtained with MuTect2 in the University of California, Santa Cruz (UCSC), with the Xena browser (https://xenabrowser.net/datapages/). Genetic alterations in the three genes associated related to pancreatic cancer were obtained from cBioPortal (https://www.cbioportal.org/) [29].

RNA extraction, reverse transcription, and qRT-PCR analysis

RNA was obtained from patients’ samples and clinically diagnosed with PDAC at the Fudan University Shanghai Cancer Center (FUSCC) from 2010 to 2012. Informed consent was obtained from each patient, and all of the experiments were performed with the approval of the Clinical Research Ethics Committee of FUSCC (ethical code: 050432-4-1212B). In the FUSCC validation set, total RNA was extracted from 82 patient samples using a MiniBEST Universal RNA Extraction Kit (9767, TaKaRa). A PrimeScript RT Reagent Kit (K1622, Thermo Scientific) was used to synthesize first-strand cDNA isolated from total RNA. Then, SYBR Green qRT-PCR was conducted on an ABI 7900HT platform (Applied Biosystems, USA). We used GAPDH mRNA as an internal reference. The primers used for the mRNA molecules tested in this study were synthesized by Sangon (Shanghai, China), and the sequences are listed in Table 3.

Statistical analysis

Statistical analysis was performed in R software (version 3.6.2, https://www.r-project.org/) and GraphPad Prism software, version 7.01 (GraphPad Software, La Jolla, CA, USA). Categorical variables were analyzed by the χ² test or Fisher’s exact test. Continuous variables were analyzed using Student’s t-test for paired samples. Kaplan-Meier curves and log-rank tests were used. Calibration plots were produced to assess whether actual outcomes matched predicted outcomes for the nomogram. All of the statistical tests were 2 sided, and a P-value < 0.05 was considered to indicate statistical significance.