Prognostic Value and Correlation With Tumor Immune Infiltration of a Novel Metabolism-Related Gene Signature in Pancreatic Cancer

Background: Energy metabolism has been considered as one of the novel features of neoplasms. This study aimed to establish the prognostic signature for pancreatic cancer (PC) based on metabolism-related genes (MRGs).

Methods: We obtained MRGs from the Molecular Signatures Database (MSigDB) and gene sequence data in the Cancer Genome Atlas (TCGA) databases. Then, differentially expressed MRGs (DE-MRGs) were identified utilizing the R software. We built the prognostic model via multivariate Cox regression. Moreover, external validation of the prognostic signature was also performed. Nomogram was created to predict the overall survival (OS). Next, this study analyzed the prognostic value, clinical relationship, and metabolism-related signaling pathways of the prognostic signature. The role in tumor infiltration was further evaluated. Eventually, the expression level of the three MRGs along with the function of NT5E was validated.

Results: Twenty-two MRGs were chosen, eight of which were identified to be most significantly correlated with the prognosis of PC. Meanwhile, a 3-MRG prognostic signature was established, and we verified this prognostic model in two separate external cohorts. What is more, the nomogram was used to predict 1-/2-/3-year OS of PC patients. In addition, the immune cell infiltration and expression of immune checkpoint were significantly influenced by the risk score. Finally, three MRGs were highly expressed in PC cell lines, and NT5E was associated with the proliferation and migration ability of PC.

Conclusion: To sum up, the study established and validated a 3-MRG prognostic signature for PC, and the signature could be utilized to predict the prognosis and assist the individualized clinical management of patients with PC.

Introduction

Pancreatic cancer (PC) remains a rather challenging and devastating disease that seriously threatens human life and health, and recently, published data reveal that it is now the sixth and fourth greatest contributor to death owing to cancer in China and the USA, respectively (1, 2). Although chemotherapy has made great strides in enhancing survival rates, it is dismaying that the 5-year survival rate for PC is only 9%, the lowest of all cancers (2). Therefore, challenges are still existing in the treatment of PC. Given the special anatomical location of the pancreas, PC is commonly symptom-free in the early stage, which only becomes apparent when the surrounding tissues are invaded or distant organs are metastasized (3). In the current case, the exploration of novel prognostic indicators to precisely predict the prognostic outcomes and guide appropriate treatment for patients suffering from cancer, including PC, is considerably urgent (4).

It is well known that tumorigenesis is strongly associated with the activation of oncogenes and the deletion of cancer suppressor genes, which is directly related to alterations in tumor metabolism (5). Otto Warburg discovered that, even in an oxygenated environment, tumor suppressor genes still chose to restrict their energy metabolism primarily to glycolysis, fermenting glucose into lactic acid to meet their energy needs (6). After that, cellular metabolism was believed to be one emerging feature of neoplasms (7). Many metabolism-related pathways are activated and crosstalked in cancer cells to promote cell survival, division, and unlimited growth (8). In addition, the impact of metabolism on the tumor microenvironment (TME) and immune regulation should not be overlooked (9). Altered tumor metabolism leads to the accumulation of specific metabolites in the TME and the increase in extracellular acidification, which facilitates the formation of a pre-metastatic ecological niche and creates a favorable environment for the metastasis of tumors, including PC (5). Therefore, numerous studies have developed novel metabolism-related prognostic signatures for predicting the prognosis of cancer (10–13). Unfortunately, no metabolism-related prognostic model is available to systematically predict the outcome of patients with PC. Thus, the construction of metabolic-based prognostic signatures that can reliably predict the outcomes of PC has remarkable clinical values.

Cancer cells obtain the ability to flourish by disrupting immunosuppressive pathways and using multiple mechanisms, such as recruiting various suppressive immune cells, secreting immunosuppressive cytokines, and producing immunosuppressive metabolites, to induce and enhance immune escape (14, 15). This saves the tumor from being cleared by the organism’s surveillance system, leading to the development of this malignancy from its origin to the abyss. Immunotherapy is an emerging cancer intervention that alters an individual’s immune system to defeat cancer through triggering a direct rejection response or interrupting inhibitory pathways (16). The therapeutic paradigm of targeting immune checkpoint by using monoclonal antibodies to inhibit cytotoxic T-lymphocyte-associated antigen 4 (CTLA-4), programmed cell death protein-1 (PD-1), and programmed cell death protein Ligand-1 (PD-L1) has achieved encouraging results in melanoma, lung cancer, and other malignancies (17). Nevertheless, most phase I and phase II clinical trials failed to demonstrate any clinical benefits on PC (18). The lacking efficacy of routine immunotherapy may be caused by the inherently aggressive biology of PC, abundant demyelinating stroma, poor immunogenicity, and poor infiltration of effector T cells, which pose serious challenges to the success of immunotherapy (19). The hypothesis has been put forward that the future of PC immunotherapy may be combined with chemotherapy, radiotherapy, or co-inhibition of different immune checkpoints (20–22). Of course, whether this will benefit a considerable number of patients with PC demands plenty of in-depth studies to confirm.

The transcriptomic data of common MRGs were utilized to build a prognostic signature, which was validated in two separate external cohorts. Subsequently, nomograms of the predictive model predicting 1-/2-/3-year OS were plotted, and the risk score was an independent prognostic factor for patients of PC. Besides, we operated the Gene Set Enrichment Analysis (GSEA) to find the metabolism-related signaling pathways primarily enriched by the prognostic signature. The study showed distinct differences in the tumor-related immune cell infiltration and immune checkpoints levels within the lower- and upper-risk groups. After that, the expression levels of three MRGs in PC cell lines were confirmed to be upregulated unanimously. Additionally, the NT5E was shown to affect the proliferative and migrative properties of PC cells. In short, this research focused on MRGs to establish one novel metabolism-related prognostic signature in PC, which could be used to guide survival risk stratification of PC to more precisely manage patients with PC.

Materials and Methods

Data Collection

The transcriptomic data, simple nucleotide variant (SNV), and matching clinical information for PC were obtained from TCGA (https://portal.gdc.cancer.gov/), GEO (https://www.ncbi.nlm.nih.gov/geo/), and the International Cancer Genome Consortium (ICGC, https://dcc.icgc.org/) database. Then, to increase the comparability and the scientific validity of the TCGA data, the Genotype-Tissue Expression (GTEx) dataset of normal pancreas samples was obtained for differentially expressed analysis on the UCSC Xena (https://xenabrowser.net/) website. Meanwhile, MRGs were extracted from the “c2.cp.kegg.v7.4.symbols” set of metabolic pathway-related genes on the MSigDB (https://www.gsea-msigdb.org/gsea/msigdb/). Subsequently, MRGs that were coexistent in the TCGA and GSE62452 datasets were extracted for subsequent analysis. To keep the gene expression levels in the training and testing sets at the same standard, a log₂(x + 1) transformation was performed in the datasets from TCGA with ICGC cohort (PACA-AU and PACA-CA) and the datasets from TCGA with GSE57495 for the construction and validation of the prognostic signature (23).

Enrichment Analysis of MRGs

DE-MRGs in PC compared to normal pancreas tissue were obtained using the R (version 4.0.5) package “limma.” Besides, the threshold for DE-MRGs screening was set to |logFC| > 1, adjusted p < 0.05. The volcano map of MRGs in TCGA was plotted using the “ggplot2” package, and the heatmap of MRGs appearing simultaneously in the TCGA and GSE62452 datasets were mapped using the “pheatmap” package. To have a clearer knowledge of the functional properties of the common MRGs and the potential signaling pathways involved, the genes were enriched with Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) utilizing the “clusterProfiler” package.

Constructing the MRGs Prognostic Signature

Samples with <30 days of follow-up and unclear clinical characteristics were removed to avoid compromising the reliability of this study. One hundred seventy-one PC samples were retained in the TCGA dataset for the construction of the prognostic signature. First, DE-MRGs significantly associated with OS of PC were identified by univariate Cox analysis, setting p < 0.01 as statistically significant. The eligible MRGs were applied to the next step of prognostic signature construction. In this study, the least absolute shrinkage and selection operator (LASSO) Cox regression method, running via the “glmnet” package, was employed to avoid the prognostic signature overfitting. Then, further multivariate Cox regression analysis of the MRGs detected by the LASSO algorithm was conducted to construct a 3-MRG prognostic signature. For per patient, the risk score is the product of the expression level and the corresponding coefficient of the prognostic signature genes, i.e., risk score = ∑(coefficient_i* expression of signature gene_i), coefficient, representing the weight of the corresponding MRG. What is more, we grouped the patients as lower- and upper-risk groups according to median risk values. We used log-rank to test the Kaplan–Meier (KM) survival analysis to compare outcome differences within low- and high-risk groups. Additionally, receiver operating characteristic (ROC) curves were plotted by the “timeROC” package to estimate the predictive accuracy of the prognostic signature.

Validation of the Prognostic Signature

For validating the universal applicability of this MRGs prognostic signature, two separate external cohorts, the ICGC cohort (PACA-AU and PACA-CA) and the GSE57495 dataset were selected for external validation. The risk score equation from the training cohort was utilized to work out risk scores for individual patients in the validation cohorts. That is to say, individuals from the validation cohorts were classified as low- and high-risk groups according to the threshold of risk values within the training cohort, drawing KM survival curves and ROC curves immediately afterward as described previously.

Relationship of Risk Score With Clinicopathological Features

This study further verified the expression of three MRGs with the Human Protein Atlas (HPA) online database (https://www.proteinatlas.org/). Next, the univariate and multivariate Cox regression analyses were introduced to determine if this risk score might be related to the prognosis and whether it can act as an independent prognostic factor independently of clinical features (age, gender, histological grade, and pathological stage). Forest plots were drawn with the “forestplot” R package, displaying the p-value, HR, and 95% CI for every variable. On the other hand, a nomogram incorporating the risk score with clinical features was constructed using the R package “rms” to predict the 1-/2-/3-year OS of a given patient. The decision curve analysis (DCA) was plotted with “ggDCA” package. In addition, the correlation between the expression of these MRGs together with the risk score and several main clinicopathological features has been explored in subsequent work. Furthermore, time-dependent ROC curves were conducted to assess the prediction performance via comparing the risk score with clinicopathological characteristics. In the end, we generated a heatmap indicating the link between the risk score with clinicopathological features and the expression of three MRGs.

Metabolism-Related Pathways Enriched by Risk Score

Based on mRNA data from 171 PC tissues in the TCGA database, the three MRGs were ranked according to the degree of differential expression in samples from low- and high-risk groups, respectively, using GSEA software (version 4.1.0) to identify potential up- and downregulated metabolism-related signaling pathways. What is more, the sensitivity of PC patients to clinically used chemotherapeutic agents for PC was analyzed using the R software package “pRRophetic.”

Mutation Landscape Map

The “maftools” R package was utilized to calculate the mutation frequencies of all genes in the SNV data from the TCGA. The comparison of whether the top 20 genes with the highest mutation frequencies were significantly different among the separate risk groups was presented in the form of a waterfall plot. Next, patients with four classical PC driver genes KRAS, CDKN2A, SMAD4, and TP53 (24) wild or mutant type were compared to determine whether their risk scores differed. Aside from this, we compared whether there were differences in tumor mutation burden (TMB) in patients with different clinical characteristics and risks. In the end, the “survminer” package was applied to analyze the survival status of patients with different TMB.

Risk Score With Tumor Immune Infiltration and Immune Checkpoint Expression

Seven methods including TIMER, CIBERSORT, CIBERSORT-ABS, QUANTISEQ, MCPCOUNTER, XCELL, and EPIC were employed to evaluate cancer immune infiltration to reveal the association between risk score with cancer immune infiltration (25). Pearson correlation analysis was used to evaluate the associations within infiltrating immune cells. The distinctions of immune infiltration cells between the low- and high-risk groups were then compared with the help of the “vioplot” package. Previous studies have suggested that the expression pattern of immune checkpoints may influence the efficacy of immunotherapy (26). We investigated the differentially expressed 47 immune checkpoints like PD-1 (also named PDCD1), PD-L1 (or named CD274), and CTLA-4 in low- and high-risk groups to reveal the potential role of this risk score in immunotherapy response.

Western Blotting

The steps of Western blotting were carried out as mentioned previously (27). Rabbit antihuman CYP2C18 (ProteinTech Group, Rosemont, IL, USA) and INPP4B (ZEN-Bioscience, Chengdu, China) epitopes polyclonal antibodies and mouse anti-human NT5E and β-actin (ProteinTech Group, Rosemont, IL, USA) epitopes monoclonal antibodies were selected as primary antibodies. Meanwhile, a horseradish-peroxidase-labeled secondary antibody (ProteinTech Group, Rosemont, IL, USA) was used as a secondary antibody for Western blotting.

Quantitative Real-Time PCR

RNA extraction and quantitative real-time PCR (qRT-PCR) procedures were performed following the description in a prior study (28). PCR primers for CYP2C18, INPP4B, and NT5E were purchased from SangonBiotech (Sangon, Shanghai, China). Primer sequences were shown as follows: CYP2C18, 5′ -GAG TTT TCT GGA AGA GGA A-3′ (forward) and 5′ -GCA TTG GTT TTT CTC AAC TCC T-3′ (reverse); INPP4B, 5′ -GCT GGA TTG GTT TGT GGT TTT A-3′ (forward) and 5′ -TAG CAT TCC AAT TTC ATC GCT G-3′ (reverse); NT5E, 5′ -ACA ACC TGA GAC ACA CAC GGA TG-3′ (forward) and 5′ -TTC GGG AAA GAT CAT ACA CCA CAT GG-3′ (reverse); β-actin, 5′ -CCT GGC ACC CAG CAC AAT-3′ (forward) and 5′ -GGG CCG GAC TCG TCA TAC-3′ (reverse). The expression of CYP2C18, INPP4B, and NT5E was standardized by the internal control β-actin. Meanwhile, fold changes in CYP2C18, INPP4B, and NT5E were calculated in a manner of the 2^−ΔΔCT.

Paclitaxel Sensitivity Analysis

Cells were seeded onto 96-well plates, treated with paclitaxel, and incubated for 72 h. After incubation, the CCK-8 reagent (10 μl) was added to each well and incubated for another hour. The absorbance at 450 nm of the color produced in each well was measured with a microplate reader to calculate the number of viable cells. Cells receiving the different treatments were inoculated in 96-well plates. The CCK-8 reagent was added every 24 h, respectively, until 72 h after inoculation and then incubated at 37°C for 1 h before testing the absorbance at 450 nm.

Cell Transfection

The NT5E shRNA lentiviral vector (hU6-MCS-CBh-gcGFP-IRES-puromycin) was purchased from GeneChem (Genechem Co., Ltd, Shanghai, China). The transfection of Capan-2 and MIA PaCa-2 cells was conducted according to official guidelines.

Validation In Vitro

The colony formation assays were performed to assess the differences in the proliferative capacity of the differently treated cells. The cells from different groups were digested and inoculated in six-well plates (Jet Biofilter Co., Ltd., Guangzhou, China) with 1,000 cells per well. The medium was changed every 3 days and incubated continuously for 10–14 days. After visible colonies were formed, they were immobilized with 4% paraformaldehyde and then stained with crystalline violet (Solarbio Life Sciences, China). To evaluate cellular migration and invasion capacity, the wound healing assays and the Transwell assays were performed. As for wound healing experiments, cells from the different treatment groups were digested and inoculated in six-well plates. A straight line was drawn on the surface of each well with a 100-μl sterile pipette after the cells had grown to approximately 95%. The wound area was photographed using an inverted microscope (Nikon DS-RI2, Japan) at 100× magnification for 0 and 24 h, respectively. For Transwell assays, 800 μl of medium containing 10% fetal bovine serum (Corning, USA) was added to the lower chamber, and 200 μl of serum-free medium containing 20,000 cells was loaded into the upper Matrigel-coated chamber. Cells that had crossed the membrane were fixed with 4% paraformaldehyde after 24 h incubation, washed with phosphate-buffered saline (PBS) and then stained with crystal violet. Images of the cells were taken with an inverted microscope at 200× magnification.

Animal Experiments

Eight-week-old female BALB/c nude mice (HFK Bioscience Co., Ltd., Beijing, China) were fed in the animal laboratory in a specific pathogen-free level environment. A total of 5 × 10⁶ Capan-2/NT5E knockdown Capan-2 cells were suspended in 300 μl of PBS and injected into the left axilla of the nude mouse. The subcutaneous tumor volume was measured weekly with the formula: V = 0.52 × L × W², where L and W denote the long and short axes of the tumor, respectively.

Statistical Analysis

The vast majority of statistical analyses and figure outputting was processed with Perl software (version 5.34.0, https://www.perl.org/) and R version 4.0.5 (https://www.r-project.org/). Differential expression of genes was analyzed with the Wilcoxon test, and the two-tailed Student’s t-test was introduced to analyze the relative expression levels of MRGs. The log-rank test was employed for Kaplan–Meier survival curve comparisons. The p < 0.05 was accepted as considered statistical significance if not otherwise mentioned. *p < 0.05, **p < 0.01; ***p < 0.001; ns, not significant.

Results

Identification of DE-MRGs

Our workflow was summarized in the flowchart (Figure 1). To determine mRNAs differentially expressed between PC and normal samples, we integrated the datasets of TCGA (178 PC and four normal tissues) and GTEx (167 normal pancreas tissues) as the TCGA-GTEx dataset. The ICGC cohort (PACA-AU and PACA-CA) was a combination of the PACA-AU cohort (262 PC cases) and the PACA-CA cohort (84 PC cases). In addition, two GEO databases were obtained from the GEO database: GSE62452 (69 PC and 61 normal tissues) and GSE57495 (63 PC tissues). Based on the predefined screening criteria (adjust p < 0.05 and |logFC| > 1.0), 5,886 genes were significantly differentially expressed in the TCGA-GTEx dataset, consisting of 2,980 overexpressed genes and 2,906 low-expressed genes. Similarly, a total of 286 genes were found in the GSE62452 dataset comprising 188 overexpressed genes and 98 downregulated genes. Figure 2A shows that a total of 766 MRGs appeared in the TCGA-GTEx dataset, for which volcano maps were plotted (Figure 2B). For another, the metabolism-related pathway gene set contained 948 MRGs, ultimately preserving 22 common MRGs (eight upregulated MRGs and 14 downregulated MRGs) simultaneously existing in all three datasets (Figure 2C). The heatmap of the 22 common MRGs is exhibited in Figure 2D.

FIGURE 1

Figure 1 The flowchart for analyzing MRGs-prognostic signature.

FIGURE 2

Figure 2 DE-MRGs in PC. (A, B) MRGs in the TCGA database are shown in the Venn diagram and volcano plot. (C) Common MRGs in the GSE62452 dataset and TCGA database are shown in the Venn diagram. (D) Common MRGs are shown in the heatmap.

Functional Enrichment of MRGs

The potential biological functions of the 22 common MRGs were also explored through GO and KEGG enrichment analysis. Altogether, 116 GO terms along with 21 signaling pathways (adjust p < 0.05) have entered into our view. Among them, the top enriched GO terms were lipid digestion and triglyceride lipase activity (Figures 3A, B). Besides, the top 5 significant GO terms were shown in the chord charts (Figures 3C, D). On the other hand, among those enriched KEGG pathways, the common MRGs were strongly related to fat digestion and absorption, glycerolipid metabolism, pancreatic secretion, etc. (Figures 3E, F). Meanwhile, the chord diagrams presented the top 5 significantly enriched KEGG signaling pathways (Figures 3G, H).

FIGURE 3

Figure 3 Enrichment analysis. (A, B) The bar and bubble charts of GO terms are enriched in BP and MF. (C, D) The chord with cluster plots shows the top 5 GO terms and the corresponding enriched genes. (E, F) The bar and bubble plots of KEGG enriched terms. (G, H) The chord with cluster maps shows the top 3 KEGG terms and the corresponding enriched genes. The false discovery rate (FDR) <0.05 of all terms was significant.

Establishing the MRGs Prognostic Signature

To ascertain the prognostic MRGs, univariate Cox regression analysis was initially applied to analyze the expression matrix and clinical follow-up information of the TCGA cohort for 22 common MRGs. As shown in the results, nine MRGs (P4HA1, ACSL5, EPHX2, CYP2C18, HK2, MBOAT2, INPP4B, DHRS9, and NT5E) were remarkably correlated with the outcomes of PC, among which only EPHX2 was a protective factor for PC (Figure 4A; p< 0.01). MRG that may be highly intercorrelated with other MRGs was then removed via LASSO Cox regression analysis to avoid overfitting, which could confound the prediction results (Figures 4B, C). In the end, one prognostic signature was established to assess the prognostic risk for PC patients on the basis of the multivariate Cox regression analysis. As a result, the final metabolism-related prognostic signature was made up of three MRGs: cytochrome P450 family 2 subfamily C member 18 (CYP2C18), inositol polyphosphate-4-phosphatase type II B (INPP4B), and 5′-nucleotidase Ecto (NT5E, or CD73, Figure 4D). The predictive model was defined to be one combination of the expression value of the three MRGs weighted by their relative coefficients in the multivariate Cox regression in the following way:

$$risk score=(0.1426\ast expression value of CYP2C18)+(0.3876\ast expression value of INPP4B)+(0.3015\ast expression value of NT5E)$$

FIGURE 4

Figure 4 Establishing MRGs prognostic signature. (A) Recognition of the outcome-related MRGs in PC via univariate Cox regression analysis. p < 0.01 were statistically meaningful. (B) The tuning parameter (lambda) is determined at the vertical line. (C) LASSO coefficient profiles of the 5-MRGs according to the tuning parameter. (D) The 3-MRG prognostic signature was built via multivariate Cox analysis.

The risk score was calculated for each patient according to the prognostic model described above. The 171 patients were categorized by median risk score as a low-risk group (n = 86) and a high-risk group (n = 85). The scatter plots of risk score ranged from lower to higher, and survival status plots are shown in Figures 5A, B. Furthermore, the difference in expression of 3-MRG between the low- and high-risk groups was displayed in the heatmap (Figure 5C). As indicated in the survival curves, patients from the low-risk group had significantly superior OS to the high-risk group (Figure 5D). Besides, the areas under the curve (AUC) values for predicting the 1-, 2-, and 3-year OS were 0.738, 0.689, and 0.659, correspondingly (Figure 5E).

FIGURE 5

Figure 5 Risk score analysis of the 3-MRG prognostic model in the training cohort. (A–D) The risk curve, survival status of PC patients, heatmap of the 3-MRG expression, and survival curves between lower and higher-risk groups are shown. (E) Time-independent ROC curves of the risk score for predicting the OS in the training cohort.

Validation of the Prognostic Signature

To validate the metabolism-relevant prognostic signature, the ICGC cohort (PACA-AU and PACA-CA) was used as the validation cohort 1, and the GSE57495 dataset was the validation cohort 2. Patients were then classified as a low- and a high-risk group by the median-risk score of the training cohort. In the validation cohort 1, the risk score, survival status, along with heatmap of 3-MRG have been displayed in Figures 6A–C. Concordant with the outcomes generated by the training cohort, the KM curves have revealed that individuals from the higher-risk categories have exceptionally worse outcomes than those in the lower-risk ones (Figure 6D). Besides, the predicted AUC values for 1-/2-/3-year OS were 0.653, 0.659, and 0.688, respectively (Figure 6E). Similar results in validation cohort 2 are presented in Supplementary Figures S1A–C. Overall, the performance of the prognostic signature as a classifier was proved as universally applicable.

FIGURE 6

Figure 6 Risk score analysis of 3-MRG prognostic signature in the validation cohort 1. (A–D) The risk score, survival status of PC patients, heatmap of the 3-MRG expression, and survival curves between low and high-risk groups are shown. (E) Time-independent ROC analysis of risk score for predicting the OS in the ICGC cohort (PACA-AU and PACA-CA).

Expression Levels of 3-MRG

We next employed the TCGA-GTEx dataset to analyze the mRNA levels of the 3-MRG. It was shown to be greater in levels of CYP2C18, INPP4B, and NT5E from PC than from normal samples (p < 0.001, Figures 7A–C). Subsequently, 3-MRG protein levels were assessed on the HPA website. As indicated in Figure 7D, the protein level of INPP4B was remarkably greater in PC than in normal tissues. The NT5E was highly expressed in PC tissues along with exocrine glandular cells of normal pancreas tissues and moderately expressed in pancreatic endocrine cells (Figure 7E). In PC tissues, INPP4B was highly/moderately expressed in 13 samples and low/non-detected in six samples, while in pancreas tissues, INPP4B was highly/moderately expressed in five samples. In PC tissue, NT5E was highly/moderately expressed in 11 samples and low/non-detected in one sample, whereas it was highly/moderately expressed in three pancreas samples (Figures 7F, G). However, the immunohistochemical results of CYP2C18 in PC and adjacent normal tissues were not available.

FIGURE 7

Figure 7 Comparing the 3-MRG mRNA levels in PC and normal pancreas tissues. (A–C) CYP2C18, INPP4B, and NT5E. Verifying the 3-MRG expression in PC with normal tissues on the HPA database. Information about CYP2C18 was not available. (D–G) INPP4B and NT5E.

Clinical Value of Risk Score

As depicted in Figure 8A, age (p = 0.011), histological grade (p = 0.034), and risk score (p < 0.001) were remarkably associated with the outcomes of patients according to univariate Cox analysis, all of which were high-risk factors for PC. Furthermore, the results of the multivariate analysis suggested that the risk score (p < 0.001), and age (p = 0.005), can be recognized to be independent prognostic variables (Figure 8B). We then plotted the nomogram (C-index = 0.673), which was designed to score patients and predict their OS at 1, 2, and 3 years to appraise the clinical application of the risk score (Figure 8C); meanwhile, we have also corrected for the predictive performance of this nomogram, as the results suggested that the predictive performance of the nomogram was satisfactory (Figures 8D–F). In addition, the DCA curves showed that the risk score model was superior to other clinical features in predicting clinical benefit (Figure 8G). In the next step, the relationship between risk score and MRGs with clinical characteristics was investigated furthermore through assessing the potential clinical application of risk score with MRGs. There was a notable association between NT5E and risk score and clinical factors, where NT5E expression was significantly correlated with survival status (p = 0.016) and histological grade (p = 0.012, Figures 9A, B). Associations between risk score and survival status (p = 3.878e−04) along with histological grading (p = 0.048) were also significant (Figures 9C, D). On the other hand, time-dependent ROC curves were also plotted to investigate other prognostic values of the prognostic signature. The AUC values of 0.792, 0.733, 0.684, 0.675, 0.700, and 0.783 for risk score predicting 0.5-, 1-, 2-, 3-, 4-, and 5-year OS, respectively, were significantly higher than other clinical features (Figures 9E–J). The predictive value of this risk score remarkably outperformed other traditional clinicopathological characters. At last, the expression heatmap of the final 3-MRG in the lower- and higher-risk categories and the corresponding clinical information in the TCGA cohort are presented in Figure 10A.

FIGURE 8

Figure 8 Independent prognostic value of risk score. (A) Univariate analysis. (B) Multivariate analysis. (C) Nomogram based on the signature for PC patients at 1/2/3-year. (D–F) Calibration plots of the nomogram at 1, 2, and 3 years. (G) DCA of this prognostic signature.

FIGURE 9

Figure 9 Association between risk score and MRGs expression with clinical features of PC. (A, B) NT5E. (C, D) Risk score. (E–J) The ROC curve analyses of the prognostic variables in the TCGA cohort of 0.5, 1, 2, 3, 4, and 5 years, respectively.

FIGURE 10

Figure 10 (A) The heatmap showed the 3-MRG expression in the prognostic signature based on different risk, clinicopathological characteristics, and survival status in the TCGA cohort. (B) GSEA revealed the top 5 up- and downregulated metabolism-related signaling pathways enriched in low- and high-risk groups. (C–E) The IC50 showed the chemosensitivity of cisplatin, gemcitabine, and paclitaxel in the low- and high-risk groups.

Metabolism-Related Pathways Enriched by Risk Score

Figure 10B demonstrates the top 5 up- and downregulated metabolism-related signaling pathways. For one thing, the genes from the high-risk group were largely enriched in the metabolism-related pathways like glycosphingolipid biosynthesis lacto and neolacto series, glycolysis gluconeogenesis, biosynthesis of unsaturated fatty acids, fructose and mannose metabolism, and drug metabolism other enzymes. For another thing, the genes from the low-risk group were largely enriched in the metabolism-related pathways such as glycine serine and threonine metabolism, tryptophan metabolism, glycosphingolipid biosynthesis ganglio series, primary bile acid biosynthesis, and alanine aspartate and glutamate metabolism. Patients in the high-risk group were more enriched in signaling pathways such as glycolysis and lacto biosynthesis, while those in the low-risk group were more enriched in signaling pathways related to amino acid metabolism. Figures 10C–E have revealed that patients from the higher-risk group were considerably less sensitive to paclitaxel than those from the lower-risk group, suggesting a horrible outcome in higher-risk individuals may be involved in chemoresistance.

The Difference in Mutations Between Different Risk Groups

Analysis of the SNV data from TCGA showed that a total of 124 of the 158 samples harbored significant mutations. The most frequent mutations were in the driver genes of PC (TP53, KRAS, CDKN2A, and SMAD4), and the most common variant classification was the missense mutation (Figures 11A, B). In the high-risk group, 56 patients harbored KRAS mutations, 53 patients carried TP53 mutations, 17 patients had SMAD4 mutations, and 20 patients harbored CDKN2A mutations. In the low-risk group, 20 patients had KRAS mutations, 23 patients harbored TP53 mutations, 8 patients carried SMAD4 mutations, and 5 patients harbored CDKN2A mutations (Figures 11C, D). Further analysis revealed that patients with mutated KRAS, CDKN2A, and TP53 presented significantly higher risk scores than those with wild type of the relevant genes (Figures 11E–H). We also found that TMB was significantly higher in male patients or patients with high risk (Figures 11I–N). Patients carrying H-TMB had a significantly shorter survival time than those with L-TMB. Survival analysis combining the risk and TMB subgroups revealed that high-risk patients harboring H-TMB had the worst prognosis, those at low risk of having H-TMB showed the second-worst prognosis, and those at low risk of having L-TMB presented the best prognosis (Figures 11O, P).

FIGURE 11

Figure 11 (A, B) Waterfall plot for the TCGA patients. (C, D) Waterfall plot for the patients in the high- and low-risk group. (E–H) Relationship between genes (CDKN2A, KRAS, SMAD4, and TP53) status and risk score. (I–N) The association of clinical characteristics and risk scores with TMB. (O, P) Kaplan–Meier survival analysis of patients with different TMB with/or different risks.

Risk Score With Tumor Immune Infiltration and Immune Checkpoint Expression

In the PC sample of the TCGA dataset, the 22 immune cell infiltration landscapes with a marked difference in the percentage of immunocytes among samples are observed in Figure 12A. The correlation matrix showed the relationships between different infiltrating immune cells in tumor samples. From the correlation heatmap, the proportion of 22 infiltrating immune cells was correlated gradually. T cells CD8 was most strongly negatively correlated with T cells CD4 memory resting (Pearson’s correlation = −0.55), and B cells naive was moderately negatively correlated with macrophages M2 (Pearson’s correlation = −0.52), whereas T cells CD4 memory activated showed the strongest positive correlation with plasma cells (Pearson’s correlation = 0.31, Figure 12B). In addition, the results of the violin plot revealed that the low-risk group included notably higher numbers of B-cell naive, plasma cell, and T-cell CD8 infiltrates, while the high-risk group had comparably higher numbers of T-cell CD4 memory resting, dendritic cells activated, etc. (Figure 12C). Next, a heatmap of the proportions of the 22 immune cell types showed that there was an increased number of infiltrations of macrophage M0 and NK cells resting in the high-risk group, while the immune score and stroma score, along with microenvironment score, were lower, whereas the infiltrating number of B-cell, T-cell, and cytotoxicity score was higher in the low-risk group (Figure 12D). On the other hand, the immune checkpoints with low expression in the high-risk group were PDCD1, CTLA4, and CD28. Conversely, TNFSF9, CD44, CD70, TNFSF4, CD276, and HHLA2 were overexpressed in the high-risk group (Figure 12E).

FIGURE 12

Figure 12 Assessment of tumor-infiltrating cells and immunotherapy value. (A) The bar chart reveals the percentage of 22 infiltrating immune cells from the TCGA cohort. (B) The correlation matrix shows correlations between infiltrating immune cells in tumor samples. (C, D) Violin plot and heatmap present the differentially infiltrated immune cells in lower- and higher-risk groups. (E) The differential expression of immune checkpoints between low- and high-risk groups.

Experimental Validation

The CYP2C18, INPP4B, and NT5E expression levels were examined in six PC cell lines (AsPC-1, BxPC-3, Capan-2, MIA PaCa-2, PANC-1, and SW1990) and normal pancreas epithelial cells (HPC-Y5) was used as control. The levels of CYP2C18, INPP4B, and NT5E were notably higher in six PC cell lines than in the HPC-Y5 cell line (Figure 13A). On the other hand, similar results were observed at the mRNA level (Figures 13B–D). Previously, the NT5E was reported to be associated with the prognosis of PC (29). Hence, we knocked down NT5E expression to assess the effect of NT5E on the biological behavior of pancreatic cancer. Western blotting was performed to detect NT5E knockdown efficiency (Figures 13E–G). Colony formation assays revealed that the proliferative capacity of Capan-2 and MIA PaCa-2 cell lines was noticeably suppressed after NT5E knockdown (Figures 14A, B). Wound healing assays indicated that knockdown of the NT5E gene resulted in a significant reduction in migration of Capan-2 and MIA PaCa-2 cell lines (Figures 14C–E). Moreover, Transwell assays demonstrated that the invasive ability of Capan-2 and MIA PaCa-2 cell lines was also dramatically restrained upon NT5E knockdown (Figures 14F–H). The subcutaneous xenograft model was constructed to confirm the effect of NT5E on PC growth in vivo. NT5E knockdown inhibited the growth of xenograft tumor in vivo (Figures 14I–K), and Western blotting detected the relative expression of NT5E (Figure 14L). The above results suggest that NT5E may be an important promoting factor for PC. As demonstrated by drug sensitivity assays, enhanced sensitivity of Capan-2 and MIA PaCa-2 cell lines to paclitaxel could be observed after knockdown of NT5E, which indicated that the signature gene NT5E was involved in mediating resistance to paclitaxel in PC (Supplementary Figures S2A–D).

FIGURE 13

Figure 13 Validation of the 3-MRGs levels via WB and qRT-PCR. (A) The protein levels of the three MRGs in six PC cell lines and pancreatic epithelial cells. (B–D) The mRNA expression of the three MRGs in six PC cell lines and pancreatic epithelial cells. (E–G) Confirmation of NT5E knockdown efficiency by WB. *p < 0.05, **p < 0.01, ***p < 0.001, ns, not significant.

FIGURE 14

Figure 14 Validation of the function of NT5E via colony formation assay, wound-healing assay, Transwell assay, and animal experiments. (A, B) Colony formation assays. (C–E) Wound-healing assays. (F–H) Transwell assays. (I–K) Animal experiments. (L) Confirmation of NT5E knockdown efficiency in vivo by WB. N, negatively control; K, NT5E knockdown. *p < 0.05, **p < 0.01, ***p < 0.001, ns, not significant.

Discussion

As of today, pancreatic cancer remains a very disruptive malignancy at the global level with a rather poor prognosis. An accurate prediction of OS in PC would facilitate the personalized treatment regimen. However, nearly no reliable and valid biomarkers to predict the prognosis of PC are currently available, for which the need to identify credible biomarkers and predictor patterns to foresee the outcome of PC is urgent. In this study, we built and verified a novel prognostic signature with MRGs based on outcome forecasting in PC patients.

To explore the value of metabolism-related genes for prognosis prediction and treatment guidance in PC. Two research teams have now established a metabolism-related prognostic signature in PC. Huo et al. (30) developed a prognostic signature based on 15 MRGs to predict the prognosis of PC patients. The proportion of CD8+ T cells infiltrated was higher in low-risk patients from this model, which predicted longer survival. Furthermore, the expression of PD-L1 was positively correlated with the risk score, whereas high PD-L1 expression induced immune escape and a worse prognosis. However, the authors did not conduct further studies on the expression of numerous immune checkpoints to explore whether risk stratification is a guide to immunosuppressive therapy. Tan’s work (31) was more concerned with metabolic reprogramming, which is one of the hallmarks of cancer, including enhanced lipid metabolism and glycolysis (32, 33), and therefore established a three-gene prognostic signature based on genes related to glycolysis and lipid metabolism. The clinical therapeutic value of targeting the signature gene c-Met for oncogenesis and metastasis and the association of CD36 with PC resistance to gemcitabine were proposed. Unfortunately, the study did not further stratify the patients with the help of the prognostic signature to more precisely indicate the patients benefiting from c-MET-targeted therapy or the patients resistant to gemcitabine. Finally, the signatures established by both teams were based on bioinformatic level analysis with no further exploration of the impact of signature genes on the biological behavior of PC.

In the current study, a 3-MRG prognostic signature was established using LASSO Cox and multifactorial Cox regression, which could differentiate between PC patients with different risks. The risk score was an independent prognostic factor for PC, and a nomogram containing clinical features was constructed on this basis. Given that paclitaxel occupies a very critical role in the clinical management of PC, the sensitivity to paclitaxel was lower in high-risk patients, who had higher expression of NT5E. NT5E gene silencing can dramatically enhance paclitaxel sensitivity in Capan-2 and MIA PaCa-2 cells. It was hypothesized that NT5E may be involved in mediating chemoresistance in high-risk patients. According to our prognostic model, patients in the low-risk group may benefit more from treatment with paclitaxel than those in the high-risk group. One application of our model was to provide a theoretical basis to guide clinical decisions on chemotherapy regimens for patients with PC. This is a complement of this research to those two relevant publications in terms of conceptual and practical aspects. In the high-risk group, the frequency of mutations in PC driver genes was dramatically higher than those in the low-risk group. High risk was also associated with high TMB, where patients at high risk with high TMB exhibited the worst outcomes. There were considerable differences in the abundance of tumor-infiltrating immune cells and in the expression of the immune checkpoint profile among different risk groups. Ultimately, the biological role of the signature gene NT5E in PC was rigorously validated by experiments in vitro and in vivo, which took this study a step forward, which is what this study exceeds these two relevant publications in practical advance.

In recent years, a substantial number of studies have been devoted to the identification of molecular subtypes of tumors, and NT5E was found to be differentially expressed in different subtypes of PC. In Collison’s work (34), the NT5E gene was found to be upregulated in a highly aggressive quasi-mesenchymal subtype. Similarly, Bailey et al. (35) revealed that NT5E gene expression was upregulated in gene programs 2 (GP2) of the mutant TP53 knockout squamous cell subtype. The two subtypes refer to the same PC subtype (quasi-mesenchymal and squamous cell subtype) as defined by the TCGA pan-cancer analysis (35). Other metabolic subtyping/stratification studies have also yielded associations between a glycolytic metabolism in PDAC with mesenchymal features and poor prognosis, where a lipogenic phenotype is associated with epithelial features and better prognosis (36). More recently, the work by Gao et al. has identified that INPP4B was enriched in a metabolic gene enriched subtype (C1) and was more highly expressed in patients at high risk of this subtype, which had a worse prognosis (37). There is an association between human cytochrome P450 (CYP) enzyme activity with the risk of multiple cancers. However, the correlation between CYP2C18 and cancer risk has not been reported (38).

Overall, KEGG enrichment suggested that all PC patients in TCGA were predominantly enriched in pathways associated with lipogenic metabolism, which was relevant to epithelial (classical) subtype characteristics with better prognostic (36, 39). Nevertheless, there is pronounced metabolic heterogeneity in PC patients, and enrichment analysis in general for all patients may not elicit accurate conclusions (40). Therefore, patients may benefit more from further patients subtyping or stratification to locate patients with similar metabolic patterns and target their metabolic characteristics. With further risk stratification of patients, a shift in metabolic pattern towards glycolysis was noted in high-risk patients, accompanied by increased lactate production. This conversion was accompanied by an enrichment of the glycolysis-related genes LDHA and ENO₂, which are responsible for catalyzing the final step of aerobic glycolysis, facilitating the glycolytic process by transforming pyruvate into lactate (41). Acidification of the extracellular microenvironment through increasing lactate production promotes tumor invasion and metastasis (42). Such patients were mainly tightly related to the mesenchymal (QM-PDA) subtype, showing a more aggressive capability and markedly worse prognosis (36, 39). Various preclinical studies have shown that blocking LDHA to reduce glycolysis could inhibit tumor growth and metastasis, presumably indicating LDHA as a potential therapeutic target (43–46). In Rajeshkumar’s work, treatment of pancreatic cancer cells carrying mutated TP53 with FX11, a small molecule inhibitor of LDHA, was revealed to have decreased metabolic activity, increased apoptosis, and attenuated tumor growth (42). In the high-risk group, the tumors were primarily glycolytic, carrying a higher proportion of TP53 mutation frequencies, which may imply a potential metabolic vulnerability to glycolytic inhibitors (e.g., FX11) in this group of patients. This is a further confirmation of the value of our model for risk stratifying patients for treatment guidance.

Cancer is caused by the accumulation of somatic mutations in oncogenes and tumor suppressor genes (47). KRAS mediates downstream signaling of growth factor receptors, and its mutational activation drives over 90% of pancreatic cancers (48). Oncogenic mutations in KRAS are involved in regulating alterations of cell metabolism, mainly in the form of increased uptake of glucose, a shift from oxidative phosphorylation to aerobic glycolysis, and enhanced levels of lactate (49, 50). All these alterations would result in enhanced proliferation, migration, and invasion of PC cells (51, 52). This study indicated that more KRAS mutations were carried by patients from the high-risk group who were more likely to exhibit aggressive metastasis. Liang and colleagues demonstrated that SMAD4 mutations induce upregulation of phosphoglycerate kinase 1 (PgK1) expression in PC, the first ATP-generating enzyme in the glycolytic pathway (53). Nuclear PgK1 preferentially drives cell metastasis via mitochondrial oxidative phosphorylation (54). An earlier study by Oshima suggested that mutations in SMAD4 were significantly associated with tumor size, lymphatic infiltration, and lymph node metastasis, whereas CDKN2A mutations were associated with lymphatic infiltration and extensive postoperative metastasis. High-risk patients screened by our prognostic model had a higher frequency of SMAD4 and CDKN2A mutations (although the former were not significant), suggesting that this group of patients may be more susceptible to lymphatic infiltration or metastasis and extensive postoperative recurrence. Studies have indicated that TP53 inactivation was associated with the malignant progression of PC (35, 55). In another study, Morton and colleagues revealed that only TP53 mutant PC cells exhibited invasive activity compared to TP53 gene deletion (56). In conclusion, mutations in genes may allow cancer cells to acquire more aggressive properties by altering their metabolic pattern.

A correlation between infiltrating immune cells in the TME with the prognosis of many cancers (e.g., ovarian cancer, kidney cancer, colorectal cancer, breast cancer, and pancreatic cancer) has been increasingly observed (30, 57). The CD8+ and CD4+ T cells, dendritic cells (DCs), regulatory T cells (Treg), and tumor-associated macrophages (TAMs) are mainly present in TME (58). It is generally accepted that the presence of T-cell infiltration predicts a better prognosis (59, 60). This study also identified significantly lower infiltration of CD8+ T cells in high-risk patients, consistently predicting a worse prognosis for these patients. Studies have suggested that T cells entering a resting state would reduce the proliferation of CD8+ T cells in lymph nodes, preventing them from recognizing cancer cells, which effectively mitigated the body’s immune response and promoted further proliferation of cancer cells (61, 62). In the high-risk group, the proportion of CD4+ T-cells memory resting was dramatically higher, possibly reducing the proliferation of CD8+ T cells, which prevented the organism’s immune system from recognizing the cancer cells and producing an effective immune response. In this case, the cancer cells were permitted to continuously proliferate. In addition, the TAMs were suggested to be associated with poor prognosis and rapid disease progression (63, 64). Xu et al. (65) proposed that M0 macrophages were an independent predictor of poor prognosis in patients with PC after finding an accumulation of M0 macrophages in the tumor tissue. The prognostic value of M0 macrophage infiltration was further confirmed by the current study, where the level of M0 macrophage infiltration was markedly higher in high-risk patients with a worse prognosis. Besides, dendritic cells activated in this study exhibited a greater level of infiltration in the high-risk group, which may facilitate the metastasis and immune escape of cancer cells by upregulating the immunosuppressive WNT pathway (65, 66). Overall, the infiltration levels of several immune cells related to the prognosis of PC patients differed significantly in the current study, and the TME may act on tumor metastasis and immune escape through these immune cells. The identification of genes aberrantly expressed in tumors is important for the development of individualized therapies, which could improve treatment outcomes (67). The expression levels of PD-1 and CTLA4 are highly correlated with the prognosis of patients receiving immunotherapy. Liu et al. revealed that high serum PD-1 and CTLA4 levels predicted better outcomes for LIHC patients treated with cytokine-induced killer cell (CIK) immunotherapy (68). PD-1 was highly expressed in older cancer patients. Similarly, PD-1 and CTLA4 levels were higher in black patients than in Caucasian or Asian patients, who may have had better outcomes with immunotherapy (68). The combination of anti-PD-1 and CTLA4 immune checkpoint inhibitors in several cancers demonstrated better efficacy than monotherapy (69–72). Analysis of immune checkpoint expression profiles revealed that PD-1 and -CTLA4 had higher expression levels in the low-risk group of patients, inferring therapeutic opportunities for novel immune modulators in low-risk patients. This result also highlighted the importance of our model in guiding the clinical personalized treatment of patients.

Admittedly, a few limitations exist during the current work. Initially, our prognostic model was based on mRNA expression levels for predicting patient prognosis, but the protein levels of genes may be more closely aligned with the real conditions of patients in actual clinical situations. Validation of the association between protein expression levels of signature genes with patient prognosis is currently not achievable due to the unavailability of clinical specimens and corresponding survival information; therefore, this very critical component will need to be refined in later clinical work. Second, the reliability and stability of our prognostic signature have to be proven on a large PC research center with prospective studies. Third, the bio-functions of the other two MRGs in PC remain to be further verified through a series of experiments.

Conclusion

A novel prognostic signature based on 3-MRG and a prognostic nomogram to predict 1-/2-/3-year OS for PC. Molecular features based on this signature offered fresh perspectives on the malignant development of PC. The difference between immune cell infiltration and the level of immune checkpoint expression in PC could be an important guide to the prognosis and treatment of PC patients. NT5E is of crucial importance for the malignant behavior of pancreatic cancer. This prognostic signature may be of great value in the prognostic assessment and treatment guidance for PC patients.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

Ethics Statement

The animal study was reviewed and approved by the ethics committee of the Affiliated Shengjing Hospital of China Medical University.

Author Contributions

HC designed the project under the direction of XT. HC, FZ, and TZ collected data and wrote the first draft. HC, TZ, ZC, JW, PL, and ZL did the statistical analysis and output the figures. HC and FZ performed the validation experiments. LZ, HW, and HT did the interpretation of data for the work and helped with part of English writing. All authors contributed to the article and approved the submitted version.

Funding

This research was supported by the evaluation of the perioperative curative effect of main pancreatic duct jejunum bridge catheter drainage in robotic pancreaticoduodenectomy (grant number 2020JH 2/10300130) and clinical application and promotion of a novel early diagnosis method for pancreatic cancer (grant number 2020JH6/10500055).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

Thanks to Suting Huang for the suggestions and English editing on the manuscript. All authors sincerely appreciate the sharing of data from TCGA, GEO, and GTEx databases.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2021.757791/full#supplementary-material

Abbreviations

TCGA, The Cancer Genome Atlas; GEO, Gene Expression Omnibus; GTEx, Genotype-Tissue Expression; MsigDB, Molecular Signatures Database; ICGC, International Cancer Genome Consortium; PC, pancreatic cancer; MRGs, metabolism-related genes; DE-MRGs, differentially expressed metabolism-related genes; FDR, false discovery rate; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; KM, Kaplan–Meier; OS, overall survival; ROC, receiver operating characteristic; AUC, area under the curve; HR, hazard ratio; CI, confidence interval; SDS-PAGE, sodium dodecyl sulfate–polyacrylamide gel electrophoresis; PVDF, polyvinylidene fluoride; TME, tumor microenvironment.

References

1. Chen WQ, Zheng RS, Baade PD, Zhang SW, Zeng HM, Bray F, et al. Cancer Statistics in China, 2015. CA Cancer J Clin (2016) 66:115–32. doi: 10.3322/caac.21338

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Siegel RL, Miller KD, Jemal A. Cancer Statistics, 2020. CA Cancer J Clin (2020) 70:7–30. doi: 10.3322/caac.21590

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Mizrahi JD, Surana R, Valle JW, Shroff RT. Pancreatic Cancer. Lancet (2020) 395:2008–20. doi: 10.1016/S0140-6736(20)30974-0

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Feng ZY, Shi MM, Li KX, Ma Y, Jiang LX, Chen H, et al. Development and Validation of a Cancer Stem Cell-Related Signature for Prognostic Prediction in Pancreatic Ductal Adenocarcinoma. J Transl Med (2020) 18:360. doi: 10.1186/s12967-020-02527-1

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Gottfried E, Kreutz M, Mackensen A. Tumor Metabolism as Modulator of Immune Response and Tumor Progression. Semin Cancer Biol (2012) 22:335–41. doi: 10.1016/j.semcancer.2012.02.009

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Warburg O. On the Origin of Cancer Cells. Science (1956) 123:309–14. doi: 10.1126/science.123.3191.309

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Hanahan D, Weinberg RA. Hallmarks of Cancer: The Next Generation. Cell (2011) 144:646–74. doi: 10.1016/j.cell.2011.02.013

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Qin C, Yang G, Yang JS, Ren B, Wang HY, Chen GY, et al. Metabolism of Pancreatic Cancer: Paving the Way to Better Anticancer Strategies. Mol Cancer (2020) 19:50. doi: 10.1186/s12943-020-01169-7

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Ruocco MR, Avagliano A, Granato G, Vigliar E, Masone S, Montagnani S, et al. Metabolic Flexibility in Melanoma: A Potential Therapeutic Target. Semin Cancer Biol (2019) 59:187–07. doi: 10.1016/j.semcancer.2019.07.016

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Hu JJ, Yu HS, Sun LD, Yan YL, Zhang LL, Jiang GN, et al. Identification of an Individualized Metabolism Prognostic Signature and Related Therapy Regimens in Early Stage Lung Adenocarcinoma. Front Oncol (2021) 11:650853. doi: 10.3389/fonc.2021.650853

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Yuan CY, Yuan MQ, Chen MQ, Ouyang JH, Tan W, Dai FF, et al. Prognostic Implication of a Novel Metabolism-Related Gene Signature in Hepatocellular Carcinoma. Front Oncol (2021) 11:666199. doi: 10.3389/fonc.2021.666199

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Li Y, Weng YL, Pan YH, Huang ZW, Chen XC, Hong WQ, et al. A Novel Prognostic Signature Based on Metabolism-Related Genes to Predict Survival and Guide Personalized Treatment for Head and Neck Squamous Carcinoma. Front Oncol (2021) 11:685026. doi: 10.3389/fonc.2021.685026

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Ye ZS, Zheng M, Zeng Y, Wei SH, Huang H, Wang Y, et al. A 13-Gene Metabolic Prognostic Signature Is Associated With Clinical and Immune Features in Stomach Adenocarcinoma. Front Oncol (2021) 11:612952. doi: 10.3389/fonc.2021.612952

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Khalil DN, Smith EL, Brentjens RJ, Wolchok JD. The Future of Cancer Treatment: Immunomodulation, CARs and Combination Immunotherapy. Nat Rev Clin Oncol (2016) 13:394. doi: 10.1038/nrclinonc.2016.65

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Ghalamfarsa G, Kazemi MH, Mohseni SR, Masjedi A, Hojjat-Farsangi M, Azizi G, et al. CD73 as a Potential Opportunity for Cancer Immunotherapy. Expert Opin Ther Targets (2019) 23:127–42. doi: 10.1080/14728222.2019.1559829

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Allard B, Longhi MS, Robson SC, Stagg J. The Ectonucleotidases CD39 and CD73: Novel Checkpoint Inhibitor Targets. Immunol Rev (2017) 276:121–44. doi: 10.1111/imr.12528

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Brower V. Checkpoint Blockade Immunotherapy for Cancer Comes of Age. J Natl Cancer Inst (2015) 107:djv069. doi: 10.1093/jnci/djv069

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Royal RE, Levy C, Turner K, Mathur A, Hughes M, Kammula US, et al. Phase 2 Trial of Single Agent Ipilimumab (Anti-CTLA-4) for Locally Advanced or Metastatic Pancreatic Adenocarcinoma. J Immunother (2010) 33:828–33. doi: 10.1097/CJI.0b013e3181eec14c

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Torphy RJ, Schulick RD, Zhu YW. Understanding the Immune Landscape and Tumor Microenvironment of Pancreatic Cancer to Improve Immunotherapy. Mol Carcinog (2020) 59:775–82. doi: 10.1002/mc.23179

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Henriksen A, Dyhl-Polk A, Chen I, Nielsen D. Checkpoint Inhibitors in Pancreatic Cancer. Cancer Treat Rev (2019) 78:17–30. doi: 10.1016/j.ctrv.2019.06.005

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Schizas D, Charalampakis N, Kole C, Economopoulou P, Koustas E, Gkotsis E, et al. Immunotherapy for Pancreatic Cancer: A 2020 Update. Cancer Treat Rev (2020) 86:102016. doi: 10.1016/j.ctrv.2020.102016

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Ma R, Qu XJ, Che XF, Yang BW, Li C, Hou KZ, et al. Comparative Analysis and In Vitro Experiments of Signatures and Prognostic Value of Immune Checkpoint Genes in Colorectal Cancer. Onco Targets Ther (2021) 14:3517–34. doi: 10.2147/OTT.S304297

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Wu MW, Li XB, Zhang TP, Liu ZW, Zhao YP. Identification of a Nine-Gene Signature and Establishment of a Prognostic Nomogram Predicting Overall Survival of Pancreatic Cancer. Front Oncol (2019) 9:996. doi: 10.3389/fonc.2019.00996

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Kamisawa T, Wood LD, Itoi T, Takaori K. Pancreatic Cancer. Lancet (2016) 388:73–85. doi: 10.1016/S0140-6736(16)00141-0

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Xu QH, Xu H, Deng RS, Wang ZJ, Li NJ, Qi ZX, et al. Multi-Omics Analysis Reveals Prognostic Value of Tumor Mutation Burden in Hepatocellular Carcinoma. Cancer Cell Int (2021) 21:342. doi: 10.1186/s12935-021-02049-w

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Goodman A, Patel SP, Kurzrock R. PD-1-PD-L1 Immune-Checkpoint Blockade in B-Cell Lymphomas. Nat Rev Clin Oncol (2017) 14:203–20. doi: 10.1038/nrclinonc.2016.168

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Zhou L, Tan XD, Kamohara H, Wang W, Wang BS, Liu JG, et al. MEK1 and MEK2 Isoforms Regulate Distinct Functions in Pancreatic Cancer Cells. Oncol Rep (2010) 24:251–55. doi: 10.3892/or_00000853

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Tan XD, Zhou L, Wang HT, Yang YF, Sun Y, Wang ZP, et al. Differential Expression Profiles of microRNAs in Highly and Weakly Invasive/Metastatic Pancreatic Cancer Cells. Oncol Lett (2018) 16:6026–38. doi: 10.3892/ol.2018.9352

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Kong LM, Liu P, Fei X, Wu TY, Wang ZP, Zhang BH, et al. A Prognostic Prediction Model Developed Based on Four CpG Sites and Weighted Correlation Network Analysis Identified DNAJB1 as a Novel Biomarker for Pancreatic Cancer. Front Oncol (2020) 10:1716. doi: 10.3389/fonc.2020.01716

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Ino Y, Yamazaki-Itoh R, Shimada K, Iwasaki M, Kosuge T, Kanai Y, et al. Immune Cell Infiltration as an Indicator of the Immune Microenvironment of Pancreatic Cancer. Br J Cancer (2013) 108:914–23. doi: 10.1038/bjc.2013.32

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Tan Z, Lei YB, Xu J, Shi S, Hua J, Zhang B, et al. The Value of a Metabolic Reprogramming-Related Gene Signature for Pancreatic Adenocarcinoma Prognosis Prediction. Aging (Albany NY) (2020) 12:24228–41. doi: 10.18632/aging.104134

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Sunami Y, Rebelo A, Kleeff J. Lipid Metabolism and Lipid Droplets in Pancreatic Cancer and Stellate Cells. Cancers (Basel) (2017) 10:3. doi: 10.3390/cancers10010003

CrossRef Full Text | Google Scholar

33. Yang JS, Ren B, Yang G, Wang HY, Chen GY, You L, et al. The Enhancement of Glycolysis Regulates Pancreatic Cancer Metastasis. Cell Mol Life Sci (2020) 77:305–21. doi: 10.1007/s00018-019-03278-z

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Collisson EA, Sadanandam A, Olson P, Gibb WJ, Truitt M, Gu SD, et al. Subtypes of Pancreatic Ductal Adenocarcinoma and Their Differing Responses to Therapy. Nat Med (2011) 17:500–03. doi: 10.1038/nm.2344

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Bailey P, Chang DK, Nones K, Johns AL, Patch AM, Gingras MC, et al. Genomic Analyses Identify Molecular Subtypes of Pancreatic Cancer. Nature (2016) 531:47–52. doi: 10.1038/nature16965

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Daemen A, Peterson D, Sahu N, McCord R, Du XN, Liu BN, et al. Metabolite Profiling Stratifies Pancreatic Ductal Adenocarcinomas Into Subtypes With Distinct Sensitivities to Metabolic Inhibitors. Proc Natl Acad Sci USA (2015) 112:E4410–17. doi: 10.1073/pnas.1501605112

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Gao Y, Zhang EC, Fei X, Kong LM, Liu P, Tan XD. Identification of Novel Metabolism-Associated Subtypes for Pancreatic Cancer to Establish an Eighteen-Gene Risk Prediction Model. Front Cell Dev Biol (2021) 9:691161. doi: 10.3389/fcell.2021.691161

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Agundez JAG. Cytochrome P450 Gene Polymorphism and Cancer. Curr Drug Metab (2004) 5:211–24. doi: 10.2174/1389200043335621

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Espiau-Romera P, Courtois S, Parejo-Alonso B, Sancho P. Molecular and Metabolic Subtypes Correspondence for Pancreatic Ductal Adenocarcinoma Classification. J Clin Med (2020) 9:4128. doi: 10.3390/jcm9124128

CrossRef Full Text | Google Scholar

40. Liang C, Qin Y, Zhang B, Ji SR, Shi S, Xu WY, et al. Energy Sources Identify Metabolic Phenotypes in Pancreatic Cancer. Acta Biochim Biophys Sin (Shanghai) (2016) 48:969–79. doi: 10.1093/abbs/gmw097

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Feng YB, Xiong YL, Qiao TY, Li XF, Jia LT, Han Y. Lactate Dehydrogenase A: A Key Player in Carcinogenesis and Potential Target in Cancer Therapy. Cancer Med (2018) 7:6124–36. doi: 10.1002/cam4.1820

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Rajeshkumar NV, Dutta P, Yabuuchi S, de Wilde RF, Martinez GV, Le A, et al. Therapeutic Targeting of the Warburg Effect in Pancreatic Cancer Relies on an Absence of P53 Function. Cancer Res (2015) 75:3355–64. doi: 10.1158/0008-5472.CAN-15-0108

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Sheng SL, Liu JJ, Dai YH, Sun XG, Xiong XP, Huang G. Knockdown of Lactate Dehydrogenase A Suppresses Tumor Growth and Metastasis of Human Hepatocellular Carcinoma. FEBS J (2012) 279:3898–10. doi: 10.1111/j.1742-4658.2012.08748.x

PubMed Abstract | CrossRef Full Text | Google Scholar

44. McCleland ML, Adler AS, Shang YL, Hunsaker T, Truong T, Peterson D, et al. An Integrated Genomic Screen Identifies LDHB as an Essential Gene for Triple-Negative Breast Cancer. Cancer Res (2012) 72:5812–23. doi: 10.1158/0008-5472.CAN-12-1098

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Xie H, Hanai JI, Ren JG, Kats L, Burgess K, Bhargava P, et al. Targeting Lactate Dehydrogenase–a Inhibits Tumorigenesis and Tumor Progression in Mouse Models of Lung Cancer and Impacts Tumor-Initiating Cells. Cell Metab (2014) 19:795–809. doi: 10.1016/j.cmet.2014.03.003

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Rong YF, Wu WC, Ni XL, Kuang TT, Jin DY, Wang DS, et al. Lactate Dehydrogenase A Is Overexpressed in Pancreatic Cancer and Promotes the Growth of Pancreatic Cancer Cells. Tumour Biol (2013) 34:1523–30. doi: 10.1007/s13277-013-0679-1

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Vogelstein B, Papadopoulos N, Velculescu VE, Zhou SB, Diaz LA Jr, Kinzler KW. Cancer Genome Landscapes. Science (2013) 339:1546–58. doi: 10.1126/science.1235122

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Collins MA, Bednar F, Zhang YQ, Brisset JC, Galban S, Galban CJ, et al. Oncogenic Kras Is Required for Both the Initiation and Maintenance of Pancreatic Cancer in Mice. J Clin Invest (2012) 122:639–53. doi: 10.1172/JCI59227

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Jonckheere N, Vasseur R, Van Seuningen I. The Cornerstone K-RAS Mutation in Pancreatic Adenocarcinoma: From Cell Signaling Network, Target Genes, Biological Processes to Therapeutic Targeting. Crit Rev Oncol Hematol (2017) 111:7–19. doi: 10.1016/j.critrevonc.2017.01.002

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Perera RM, Bardeesy N. Pancreatic Cancer Metabolism: Breaking It Down to Build It Back Up. Cancer Discov (2015) 5:1247–61. doi: 10.1158/2159-8290.CD-15-0671

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Buscail L, Bournet B, Cordelier P. Role of Oncogenic KRAS in the Diagnosis, Prognosis and Treatment of Pancreatic Cancer. Nat Rev Gastroenterol Hepatol (2020) 17:153–68. doi: 10.1038/s41575-019-0245-4

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Nishikawa M. Reactive Oxygen Species in Tumor Metastasis. Cancer Lett (2008) 266:53–9. doi: 10.1016/j.canlet.2008.02.031

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Bernstein BE, Hol WG. Crystal Structures of Substrates and Products Bound to the Phosphoglycerate Kinase Active Site Reveal the Catalytic Mechanism. Biochemistry (1998) 37:4429–36. doi: 10.1021/bi9724117

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Liang C, Shi S, Qin Y, Meng QC, Hua J, Hu QS, et al. Localisation of PGK1 Determines Metabolic Phenotype to Balance Metastasis and Proliferation in Patients With SMAD4-Negative Pancreatic Cancer. Gut (2020) 69:888–900. doi: 10.1136/gutjnl-2018-317163

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Hosoda W, Chianchiano P, Griffin JF, Pittman ME, Brosens LA, Noe M, et al. Genetic Analyses of Isolated High-Grade Pancreatic Intraepithelial Neoplasia (HG-PanIN) Reveal Paucity of Alterations in TP53 and SMAD4. J Pathol (2017) 242:16–23. doi: 10.1002/path.4884

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Morton JP, Timpson P, Karim SA, Ridgway RA, Athineos D, Doyle B, et al. Mutant P53 Drives Metastasis and Overcomes Growth Arrest/Senescence in Pancreatic Cancer. Proc Natl Acad Sci USA (2010) 107:246–51. doi: 10.1073/pnas.0908428107

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Zhang L, Conejo-Garcia JR, Katsaros D, Gimotty PA, Massobrio M, Regnani G, et al. Intratumoral T Cells, Recurrence, and Survival in Epithelial Ovarian Cancer. N Engl J Med (2003) 348:203–13. doi: 10.1056/NEJMoa020177

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Shiao SL, Ganesan AP, Rugo HS, Coussens LM. Immune Microenvironments in Solid Tumors: New Targets for Therapy. Genes Dev (2011) 25:2559–72. doi: 10.1101/gad.169029.111

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Pages F, Berger A, Camus M, Sanchez-Cabo F, Costes A, Molidor R, et al. Effector Memory T Cells, Early Metastasis, and Survival in Colorectal Cancer. N Engl J Med (2005) 353:2654–66. doi: 10.1056/NEJMoa051424

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Sato E, Olson SH, Ahn J, Bundy B, Nishikawa H, Qian F, et al. Intraepithelial CD8+ Tumor-Infiltrating Lymphocytes and a High CD8+/regulatory T Cell Ratio Are Associated With Favorable Prognosis in Ovarian Cancer. Proc Natl Acad Sci USA (2005) 102:18538–43. doi: 10.1073/pnas.0509182102

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Chen LP, Han X. Anti-PD-1/PD-L1 Therapy of Human Cancer: Past, Present, and Future. J Clin Invest (2015) 125:3384–91. doi: 10.1172/JCI80011

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Naidoo J, Page DB, Li BT, Connell LC, Schindler K, Lacouture ME, et al. Toxicities of the Anti-PD-1 and Anti-PD-L1 Immune Checkpoint Antibodies. Ann Oncol (2015) 26:2375–91. doi: 10.1093/annonc/mdv383

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Canioni D, Salles G, Mounier N, Brousse N, Keuppens M, Morchhauser F, et al. High Numbers of Tumor-Associated Macrophages Have an Adverse Prognostic Value That can be Circumvented by Rituximab in Patients With Follicular Lymphoma Enrolled Onto the GELA-GOELAMS FL-2000 Trial. J Clin Oncol (2008) 26:440–46. doi: 10.1200/JCO.2007.12.8298

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Leek RD, Lewis CE, Whitehouse R, Greenall M, Clarke J, Harris AL. Association of Macrophage Infiltration With Angiogenesis and Prognosis in Invasive Breast Carcinoma. Cancer Res (1996) 56:4625–29.

PubMed Abstract | Google Scholar

65. Xu CM, Sui SL, Shang YR, Yu ZY, Han J, Zhang GX, et al. The Landscape of Immune Cell Infiltration and Its Clinical Implications of Pancreatic Ductal Adenocarcinoma. J Adv Res (2020) 24:139–48. doi: 10.1016/j.jare.2020.03.009

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Argentiero A, De Summa S, Di Fonte R, Iacobazzi RM, Porcelli L, Da Via M, et al. Gene Expression Comparison Between the Lymph Node-Positive and -Negative Reveals a Peculiar Immune Microenvironment Signature and a Theranostic Role for WNT Targeting in Pancreatic Ductal Adenocarcinoma: A Pilot Study. Cancers (Basel) (2019) 11:942. doi: 10.3390/cancers11070942

CrossRef Full Text | Google Scholar

67. Cao Z, Zhang SH. An Integrative and Comparative Study of Pan-Cancer Transcriptomes Reveals Distinct Cancer Common and Specific Signatures. Sci Rep (2016) 6:33398. doi: 10.1038/srep33398

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Liu JN, Kong XS, Huang T, Wang R, Li W, Chen QF. Clinical Implications of Aberrant PD-1 and CTLA4 Expression for Cancer Immunity and Prognosis: A Pan-Cancer Study. Front Immunol (2020) 11:2048. doi: 10.3389/fimmu.2020.02048

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Hammers HJ, Plimack ER, Infante JR, Rini BI, McDermott DF, Lewis LD, et al. Safety and Efficacy of Nivolumab in Combination With Ipilimumab in Metastatic Renal Cell Carcinoma: The CheckMate 016 Study. J Clin Oncol (2017) 35:3851–58. doi: 10.1200/JCO.2016.72.1985

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Hellmann MD, Rizvi NA, Goldman JW, Gettinger SN, Borghaei H, Brahmer JR, et al. Nivolumab Plus Ipilimumab as First-Line Treatment for Advanced Non-Small-Cell Lung Cancer (CheckMate 012): Results of an Open-Label, Phase 1, Multicohort Study. Lancet Oncol (2017) 18:31–41. doi: 10.1016/S1470-2045(16)30624-6

PubMed Abstract | CrossRef Full Text | Google Scholar

71. Larkin J, Chiarion-Sileni V, Gonzalez R, Grob JJ, Cowey CL, Lao CD, et al. Combined Nivolumab and Ipilimumab or Monotherapy in Untreated Melanoma. N Engl J Med (2015) 373:23–34. doi: 10.1056/NEJMoa1504030

PubMed Abstract | CrossRef Full Text | Google Scholar

72. Wolchok JD, Kluger H, Callahan MK, Postow MA, Rizvi NA, Lesokhin AM, et al. Nivolumab Plus Ipilimumab in Advanced Melanoma. N Engl J Med (2013) 369:122–33. doi: 10.1056/NEJMoa1302369

PubMed Abstract | CrossRef Full Text | Google Scholar