A four oxidative stress gene prognostic model and integrated immunity-analysis in pancreatic adenocarcinoma

Background and aims: Pancreatic adenocarcinoma (PAAD) is highly aggressive and characterized by a poor prognosis. Oxidative stress has great impacts on the occurrence and development of tumors. However, the predictive role of oxidative stress related genes on PAAD patients’ prognosis remains unclear. In this study, we aimed to construct a prognostic model for PAAD based on oxidative stress genes and to evaluate its predictive value.

Methods: The Cancer Genome Atlas (TCGA) and three Gene Expression Omnibus (GEO) datasets were used to identify differentially expressed oxidative stress genes. Univariate Cox regression, Kaplan-Meier and multivariate Cox regression analysis were used to select genes and to construct a prognosis model. According to the median value of the model’s risk score, patients were divided into high and low risk groups, and gene set enrichment analysis (GSEA), immune infiltration and immunotherapy effect, drug resistance and the expression of immune checkpoint related genes and synthetic driver genes of T cell proliferation were analyzed. Finally, the mRNA and protein levels of four genes in PAAD were verified by the clinical proteomic tumor analysis consortium (CPTAC) database and the immunostaining of patients’ tissue.

Results: 55 differentially expressed oxidative stress genes were identified, and four genes including MET, FYN, CTTN and CDK1 were selected to construct a prognosis model. GESA indicated that immune related pathways, metabolic pathways and DNA repair pathways were significantly enriched in the high risk group as compared to the low risk group. The frequency of genetic mutations was also significantly higher in high risk groups than that in low risk groups. Moreover, the infiltration level of 23 immune cells as well as the expression of immune checkpoint related and synthetic driver genes of T cell proliferation were significantly altered, with the better immunotherapy effect occurring in low risk group. In patient PAAD tissues, the mRNA and protein levels of these four genes were up-regulated.

Conclusion: We have successfully constructed a four oxidative stress gene prognostic model that has important predictive value for PAAD patients, and this model might be a promising guidance for prognostic prediction and efficacy monitoring in clinical individualized therapy.

Introduction

Pancreatic adenocarcinoma (PAAD), a common digestive system tumor, is the seventh leading cause of cancer-related death worldwide (1, 2). It is estimated that 60,430 people were diagnosed, and 48,220 deaths were ascribed to PAAD in the United States in 2021 (1). PAAD has a high degree of malignancy and usually diagnosed at advanced stages, and it has the poorest prognosis with a 5-year overall survival (OS) rate of about 10% among various solid malignancies (3). Therefore, it is an urgent need to find a new sensitive and reliable prognostic model to accurately predict patient survival and guide reasonable treatment options (4).

In recent years, the tumor microenvironment (TME) has attracted much more attention for its important role in tumor development and progression, especially immunotherapy has brought great hope to cancer patients. As a highly complex system that consists of tumor cells, immune infiltrating cells, cancer-associated stromal cells and others, the stability and balance of TME components is crucial. In contrast, the upset of the whole balance by various conditions would accelerate or slow down tumor development or progression, and thus affect patients’ prognosis (5). Among various TME disturbed conditions, oxidative stress, which is involved in each stage of tumorigenesis as well, is considered a hallmark of cancer (6).

Oxidative stress is an important factor in tumor occurrence and development. In 1985, Helmut Sies firstly introduced the concept of oxidative stress in a review entitled “biochemistry of oxidative stress”, and defined oxidative stress as a disturbance in the prooxidant-antioxidant balance, which mainly emphasizes the prooxidant (7). Upon exposed to continuous environmental stress, for example, UV, metabolic stress, and anti-cancer drugs, the intracellular biological redox steady state is broken and excessive reactive oxygen species (ROS) is produced, which then alters cell growth, cell death, immune dysfunction and signal transduction (8). More importantly, ROS could lead to the changes in base modification or rearrangement of DNA sequence, DNA damage-derived miscoding lesions and activation of the oncogene, which further jointly contributing to the development and progression of tumors (9). In fact, many cancers including brain cancer (10), breast cancer (11), pancreatic adenocarcinoma (12, 13), lung cancer and so on are all connected to oxidative stress (14, 15), and so does other aging-related diseases (16). Additionally, there has been advanced cognition about the influences of ROS on the development, function and control of dendritic cells (DC), macrophages, natural killer (NK) cells, T cells and B cells (17, 18); more importantly, oxidative stress mediators suppress effector T cell function in TME (19). Being a double-edged sword, oxidative stress was reported to boost or prevent tumor development or progression either by affecting tumor cells or by remodeling TME (20). However, the global influences of oxidative stress on PAAD patients’ prognosis and the relationship between ROS-related biomarkers and TME in PAAD progression remains to be uncovered.

In this study, we have constructed a prognostic model with four oxidative stress related genes that screened from the TCGA database and three GEO datasets, and have developed a nomogram that may precisely predict PAAD patients’ prognoses. Moreover, we have explored the level of immune infiltration and the expression of immune related genes in patients with different prognoses based on the four oxidative stress gene prognostic model. Our findings could help to elucidate the role of oxidative stress in the progression and prognosis of PAAD, and this model might become promising guidance for prognostic prediction and efficacy monitoring in clinical individualized therapy.

Materials and methods

Data collection and preprocessing

RNA sequencing (RNA-Seq) data of TCGA were collected from the University of California Santa Cruz Xena (http://xena.ucsc.edu/), which included 170 PAAD tissue samples and 4 normal tissue samples. Along these data, somatic mutation data in the “Masked Somatic Mutation” type were processed by VarScan2 (21); clinical data including histological type, age, gender and survival data were obtained from the TCGA database (“https://portal.gdc.cancer.gov/”). Besides, datasets GSE62165, GSE62452 and GSE28735 were downloaded from the GEO database. Next, all the gene names were transformed into corresponding gene symbols; and 1399 oxidative stress related genes with the relevance score ≥ 7 were downloaded from the GeneCards database(https://www.genecards.org) (Supplementary Table 1) for further analysis.

Identification of differentially expressed mRNA in PAAD

170 RNA-Seq data of TCGA with complete clinical information were randomly divided into a training set (106 cases) and a validation set (64 cases). The original expression data of TCGA were transformed with log2, and differentially expressed genes (DEGs) were identified by the “DESeq2” package of R language (22). Further, DEGs were selected using the following criteria: (a) |log₂FC| > log₂(1.2), (b) P < 0.05. DEGs in GSE62165, GSE62452 and GSE28735 were screened under the same conditions with the “limma” package (23).

Identification of oxidative stress related genes prognostic signatures for PAAD

The DEGs between the training set of TCGA and the expression profiling analysis of GSE62165, GSE62452, and GSE28735 were intersected, the oxidative stress related genes among the above DEGs were further intersected with 1399 oxidative stress related genes, and then the differentially expressed oxidative stress related genes (DEOSGs) were obtained. Subsequently, DEOSGs were further screened with univariate Cox regression analysis by exploring the relationship between genes and the overall survival of patients using the “survival” R software package, and the genes with P < 0.05 were identified as prognosis related DEOSGs. Finally, the Kaplan-Meier method and log-rank test were applied to these genes to screen survival related DEOSGs using the R Bioconductor “survival” package.

Clinical correlation analysis

To explore the correlation between clinical information and patients’ prognosis, RNA-Seq data were divided into 2 groups according to gender (male, female), age (> 65, ≤ 65), clinical stage (stage I, stage II/III), then Kaplan-Meier method and log-rank test were applied using R Bioconductor “survival” package. P < 0.05 was considered statistically significant.

Prognostic model construction and efficacy evaluation

Based on the survival related DEOSGs, information on patients’ survival status and survival time, multivariate Cox regression analysis was applied to construct a prognostic model, and each patient was assigned a risk score. The risk score of each patient was calculated following the equation below:

$$RiskScore={\displaystyle \sum _{n=1}^n(CiEi)}$$

For short, the risk score was generated as the sum (∑) of the coefficient value (C_i) multiplied by the expression value (E_i) of each selected DEOSG. With the median value of risk score as the cut-off value, patients were classified into high and low risk groups. Additionally, the “timeROC” packages in R were used to verify the predictive accuracy of this signature. Finally, the validation set was used to test the accuracy of this model. To facilitate clinical use, patients’ information was used to construct a nomogram with the “RMS” package (24). All oxidative stress related genes identified by multivariate Cox regression analysis and survival related clinical information were used to build a nomogram to investigate the probability of 1 and 2-year OS of PAAD.

Gene set enrichment and immune infiltration analyses

Gene set enrichment analysis (GSEA) was performed to clarify the differences in prognostic correlated signaling pathways in high risk and low risk groups (25). In detail, the enriched KEGG pathways, biological processes, cellular components and molecular functions were identified by using the R package “GSEA” based on the Molecular Signatures Database v. 6.2. C2 (curated gene sets) and C5 (GO gene sets). |NES| > 1 and FDR < 0.25 or P < 0.05 were considered statistically significant. Patients were divided into high and low risk groups with the median value of risk score as the cut-off value, and the frequency of gene mutation in high and low risk groups was counted respectively. For an extension of the GSEA algorithm (26, 27), the single sample gene set enrichment analysis (ssGSEA) algorithm was used to explore the immune infiltration difference between high and low risk patients, and the infiltration of 28 immune cells was calculated by R package “ssGSEA” with reported cell makers. Herein, P < 0.05 was considered statistically significant. Given that the expression level of immune checkpoint related genes was related to the treatment responses of immune checkpoint inhibitors (ICIs) (Supplementary Table 2). The relationship between risk score and treatment responses of ICIs was explored by detecting differences in expression levels of immune checkpoint related genes between high risk and low risk groups. The similar analysis was also applied for synthetic driver genes of T cell proliferation (28) (Supplementary Table 3). Finally, Tumor Immune Dysfunction and Exclusion (TIDE) algorithms (http://tide.dfci.harvard.edu/) was further used to assess the drug resistance of ICI therapy (29).

Protein database validation, quantitative real‐time pcr and immunohistochemistry (IHC) analysis

The protein expression of CDK, MET, FYN and CTTN in PAAD were inquired in the CPTAC database (https://cptac-data-portal.georgetown.edu/). And, the mRNA levels of the above four genes were determined with quantitative real-time PCR (qRT-PCR) in 15 paired PAAD/adjacent non-tumor samples collected from patients who had undertaken surgery in Xijing Hospital (Fourth Military Medical University). The study was approved by the ethics committee of Fourth Military Medical University (XJYYLL-2015564). Briefly, total RNA was extracted from these samples using Total RNA Kit II (Omega, USA) following the manufacturer’s instructions, and then was reversely transcribed to cDNA using PrimeScript RT Reagent Kit (TaKaRa, Biotechnology, Japan). qPCR was performed using SYBR Premix Ex Taq II Kit (TaKaRa, Biotechnology, Japan), with GAPDH was using as the reference gene. The qPCR data were analyzed using the ΔΔCt method. Primer sequences were shown as follows: GAPDH forward: GTCTCCTCTGACTTCAACAGCG, reverse: ACCACCCTGTTGCTGTAGCCAA; and CDK1 forward: GGAAACCAGGAAGCCTAGCATC, reverse: GGATGATTCAGTGCCATTTTGCC; MET forward: TGAGAGCTGCACCTTGACTT, reverse: AATTTCCAGTTAAAGTAAG; CTTN forward: AGGTGTCCTCTGCCTACCAGAA, reverse: CCTGCTCTTTCTCCTTAGCGAG; FYN forward: CTGGTCACCAAAGGAAGAGTGC, reverse: GGTCCTTTTTCCAGCAGTGGATC.

In addition, 5 paired PAAD/adjacent non-tumor frozen tissues were cut into 5 μm sections, transferred to an adhesive-coated slide, and fixed in 4% (wt/vol) paraformaldehyde for 30 min at 37°C. Endogenous peroxidase was blocked using 3% H₂O₂ dissolved in methanol, and cell membrane were permeabilized using 0.25% Triton X-100. Then, slides were blocked with 5% goat serum for 2 hours at room temperature. Afterward, slides were incubated with primary anti-CDK1 antibodies (Cat No. 19532-1-AP, Proteintech, China), anti-CTTN antibodies (Cat No. 11381-1-AP, Proteintech, China), anti-FYN antibodies (Cat No. 66606-1-Ig, Proteintech, China) and anti-MET antibodies (Cat No. 25869-1-AP, Proteintech, China) at 4°C overnight. The next day, the sections were incubated with biotinylated secondary antibody for 1 hour at room temperature, followed by being visualized in DAB and observed under a light microscope, staining intensity and percentage score for IHC were assessed by two pathologists in a blinded manner. The staining intensity was scored as: 0 for negative, 1 for weak, 2 for intermediate, and 3 for strong staining. The percentage score ranges from 0 to 4, that is: 0, no immunostaining; 1, 1-35% of cells are stained; 2, 36-70% are positive; and 3, ≥70% of cells are positive. For statistical analysis, IHC score was calculated by plus staining intensity and percentage score.

Statistical analysis

R software (version 4.0.3) was used for statistical analysis in this study. Overall survival analysis was performed using the Kaplan-Meier method and the log-rank test. Univariate and multivariate Cox regression analyses were performed to calculate the prognostic significance of DEOSGs in PAAD patients. The statistical significance of differences between independent groups was calculated using student’s t-test. The Venn diagrams were drawn using Evenn (http://www.ehbio.com/test/venn/#/). If not specified above, P < 0.05 was considered statistically significant.

Results

DEOSGs identification

To screen genes related to oxidative stress, we performed our study as described in the flow chart (Figure 1). In the TCGA database (training set), a total of 2994 genes were identified as DEGs at mRNA level (Figure 2A). Three datasets in GEO were also screened and 4603, 7291 and 4428 genes were identified as DEGs in GSE28735, GSE62165 and GSE62452, respectively (Figures 2B–D, Table 1). Among the identified DEGs above, 356 DEGs are shared by the four data sets (Figure 2E), of which 55 genes were related to oxidative stress and were identified as DEOSGs (Figure 2F). Considering the gene expression level (mean value of fragments per kilobase of transcript per million mapped reads (FPKM)> 1), 52 out of 55 DEOSGs were selected for further analysis.

FIGURE 1

Figure 1 Flowchart of the study. OS, Oxidative stress; TCGA, The Cancer Genome Atlas; DEG, differentially expressed genes; DEOSGs, differentially expressed oxidative stress related genes; GSEA, Gene set enrichment analysis; ssGSEA, single sample gene set enrichment analysis; CPTAC, clinical proteomic tumor analysis consortium.

FIGURE 2

Figure 2 Identification of differently expressed OS genes. (A–D) Volcano plot of TCGA (A), GSE28735 (B), GSE62165(C), GSE62452 (D), dots in green represent down-regulated genes, dots in red represent up-regulated genes, and dots in black represent unchanged genes. (E) DEGs are genes shared by the training set of TCGA and GSE28735, GSE62165, GSE62452. (F) DEOSGs are genes shared by DEGs from the above four datasets and oxidative stress (OS) related genes.

TABLE 1

Table 1 The numbers of screened DEGs in four databases.

Construction and evaluation of prognostic model

To construct the prognostic model, 52 DEOSGs were further applied for analyzing their relationship with the prognosis of PAAD patients. Firstly, 22 genes were identified as PAAD prognostic-associated genes with P < 0.05 by univariate cox regression analysis, and the top 10 prognostic-associated DEOSGs were listed in the forest plot (Figure 3A). Further, 12 genes significantly related to overall survival were selected out by Kaplan-Meier analysis (Figure 3B); among which four genes (MET, CDK1, CTTN, FYN) with P < 0.1 were chosen by multivariate Cox regression analysis (Figure 3C). Finally, the above four genes were applied to construct the prognostic model, and the risk scores were calculated as follows: risk score = (0.8481 × MET expression) + (0.3674 × CDK1 expression) - (0.2532 × CTTN expression) - (0.7350 × FYN expression).

FIGURE 3

Figure 3 Construction of prognostic model and the nomogram for PAAD patients. (A) The forest plot of the top 10 DEOSGs in PAAD by univariate Cox regression analysis. (B) 12 genes significantly related to overall survival were selected by Kaplan-Meier analysis. (C) The forest plot of 10 DEOSGs in PAAD screened by multivariate Cox regression analysis. (D) Overall survival analysis of PAAD patients grouped by age, gender and clinical stage respectively. (E) A nomogram was constructed to predict the probability of 1 and 2-year overall survival of PAAD.

To improve clinical applicability, we jointly applied risk scores and patients’ information to build a prediction model. Among three parameters that have completely available clinical information, which including age, gender and clinical stage, the clinical stage was significantly correlated with the patients’ prognosis, while no significant correlation between patients’ prognosis and gender or age was found (Figure 3D). Together, by combing the expression level of the above four genes and the clinical stage information of PAAD patients, a nomogram which could predict the 1 and 2- year OS probabilities of PAAD patients was constructed (Figure 3E).

Furthermore, we analyzed the relationship between risk scores and patients’ survival time, survival status or expression level of the above four genes. As shown in Figure 4A, the survival probability of patients as a whole decreased significantly as the risk scores increased. As expected, the expression of risk factors, that were MET and CDK1, CTTN, showed an enhancing trend as the risk scores increased; on the contrary, the expression of protective factor, that is FYN, was declining as the risk scores increased (Figure 4A). Finally, we used log-rank test to evaluate the survival probability in high and low risk groups, and a significant prognostic difference was observed (P = 1.5 x 10^-5) (Figure 4B). In addition, time ROC analysis indicated that the prediction model was credible, as the areas under the curve (AUCs) of 1 and 2-year OS for the training set were 0.87 and 0.91, and for the validation set were 0.74 and 0.81, respectively (Figures 4C, D). These data suggested that the prognosis model constructed by these four oxidative stress genes has a good prediction ability on the prognosis of PAAD patients.

FIGURE 4

Figure 4 Evaluation and validation of the prognostic model for overall survival in PAAD patients of the training set. (A) Risk scores (up), overall survival (OS, middle) and expression profiles of four genes in high and low risk groups (down). (B) Kaplan-Meier analysis in PADD patients with high and low risk. (C) Time ROC in the training set. (D) Time ROC in the validation set.

Gene set enrichment and immune related genes expression analysis

To explore the potential signaling pathways that involved in the different prognoses of PAAD patients with high or low risk, we carried out the GSEA. Results displayed that there were several pathways significantly enriched in low risk patients as compared to that in high risk patients, which including B cell activation (NES=-1.62, P < 0.001), B cell proliferation (NES=-1.73, P < 0.001), leukocyte mediated immunity (NES=-1.52, P < 0.001), leukocyte proliferation (NES=-1.44, P < 0.001), T cell activation (NES=-1.36, P < 0.001), T cell receptor complex (NES=-1.82, P < 0.001), chemokine signaling pathway (NES=-1.63, P < 0.001), immune receptor activity (NES=-1.62, P < 0.001) and chemokine signaling pathway (NES=-1.63, P < 0.001)(Figure 5A). Besides, 3 pathways related to cell metabolism, that is glutamate receptor signaling pathway (NES=-1.90, P < 0.001), insulin secretion (NES=-1.77, p < 0.001) as well as glycine, serine and threonine metabolism (NES=-1.67, P < 0.001), were also significantly enriched in low risk patients (Figure 5B). On the contrary, 3 significant pathways including DNA repair (NES=1.52, P < 0.001), recombinational repair (NES=1.41, P < 0.0019) and cell cycle (NES=-2.01, P < 0.001) were significantly enriched in high risk patients (Figure 5C). As for gene mutation in PADD, the total mutation count in the high and low risk groups was 10,717 and 864, respectively (Figure 5D). And, the top 10 genes with the highest mutation frequency in the high risk group were TP53, KRAS, CDKN2A, SMAD4, RNF43, TTN, MUC1, ADAMTS12, FLG, and KMT2D, while in the low risk group were TP53, KRAS, CDKN2A, SMAD4, RNF43, TTN, TGFBR2, SYNE1, PKHD1, FLRT2(Figure 5E). As for six genes (TP53, KRAS, CDKN2A, SMAD4, RNF43, TTN) shared by both groups, the mutation frequency in high risk groups are also higher than that in low risk groups (Figure 5F). Together, these results indicated that immune related pathways and metabolic pathways were significantly enriched in the low risk group, while DNA repair pathways were significantly enriched in the high risk group.

FIGURE 5

Figure 5 Functional GSEA and gene mutation analysis in PAAD patients with high and low risk. (A) Immune-related signatures were significantly enriched in low risk patients. (B) Metabolic-related signatures were significantly enriched in low risk patients. (C) Recombinant repair, DNA repair and cell cycle signatures were significantly enriched in high risk patients. (D) Total mutation counts in PAAD patients with high and low risk. (E) Top 10 genes with the highest mutation frequency and their overlaps in PAAD patients between high and low risk. (F) The mutation frequency analysis of 6 common genes in PAAD patients with high and low risk.

As the majority of enriched pathways were related to the immune system, we then turn to perform the immune infiltration analysis in high and low risk groups. As shown in Figure 6A, 23 out of 28 immune cells were significantly altered between high and low risk groups. Among 23 immune cells, 19 cells (activated B cell, activated CD8⁺ T cell, activated dendritic cell, central memory CD4⁺ T cell, effector memory CD4⁺ T cell, effector memory CD8⁺ T cell, eosinophil, immature B cell, immature dendritic cell, macrophage, mast cell, MDSC, monocyte, natural killer T cell, natural killer cell, plasmacytoid dendritic cell, regulatory T cell, T follicular helper cell, and type 1 T helper cell) have significantly negative correlations with a risk score. On the contrary, 4 cells that include CD56dim natural killer cell, neutrophil, type 17 T helper cell and type 2 T helper cell have significantly positive correlations with risk scores (P < 0.05). Due to the facts that 10 out of 23 significantly altered immune cells were related to T cells, we speculated that the four oxidative stress related gene prognostic model has a close relationship with T cell immunity. Considering the recent breakthrough in checkpoint-based immunotherapy (30), we conducted a differential expression analysis of 25 immune checkpoint genes between high and low risk groups. The results showed that CD47, TNFSF9 and PVR were positively correlated with risk scores; while HLA-DQB1, CD96, SIRPA, CD48, HLA-DRB1, BTN2A2, HLA-DOA, HLA-DPA1, HLA-DPB1, HLA-DMB, HLA-DMA, BTNL9, HLA-DRA, CD27 and HLA-DOB showed an opposite trend (Figure 6B).

FIGURE 6

Figure 6 Analysis of immune infiltration and immune related genes. (A) Differential infiltration analysis of 28 immune cells based on ssGSEA algorithm between high and low risk groups. (B) Differential expression analysis of 25 immune checkpoint genes between high and low risk groups. (C) Differential expression analysis of 29 synthetic driver genes of T cell proliferation between high and low risk groups. (D) TIDE algorithm to predict the efficacy of immunotherapy in high and low risk groups. *P < 0.05, **P < 0.01, ***P < 0.001. ns P<0.05.

From GSEA and immune infiltration analysis, we noticed that there were significant differences in T cell activation between high and low risk groups. Therefore, we carried out the differential expression of 29 T cell proliferation genes between high and low risk groups. Analysis of synthetic driver genes of T cell proliferation showed that LTBR, IL1RN, BATF, AHNAK, SLC10A7, CLIC1, RAN, CDK2, CDK1 and AHCY were positively correlated with risk score; while CXCL12, NFYB, FOSB, CYP27A1, GPD1 and ITM2A have an opposite trend (Figure 6C). Taken together, our data revealed that oxidative stress may reshape tumor immune environment by activating CD4⁺ T and CD8⁺ T cell, and thus affect the prognosis of PAAD patients.

Immunotherapy using ICIs has brought great hope to PAAD patients, we then analyzed the drug resistance of ICI therapy by TIDE algorithms. As shown in Figure 6D, the TIDE score, Exclusion score and Dysfunction score in high risk group are significantly increased, indicating that the immune escape potential of high-risk group patients is increased, and thus the efficacy of ICI could be poor.

Expression validation of four prognostic DEOSGs

To validate the expression level of four prognostic DEOSGs in PAAD, we firstly detected their protein levels by digging into the CPTAC database and by performing immunohistochemistry (IHC) staining. We found that the protein levels of four prognostic DEOSGs were all significantly enhanced in PAAD as compared to normal tissues (Figure 7A). Meanwhile, the results of IHC staining further verified that these four genes were highly expressed in PAAD patients (Figures 7B, C). Besides, consistent with the protein levels, the mRNA levels of MET, CDK1, CTTN and FYN were also significantly upregulated in PAAD as compared to that of adjacent para-cancer tissues (Figure 7D). Taken together, our data showed that the expression of four prognostic DEOSGs was upregulated in PAAD; which suggesting that the jointly high expression of four prognostic DEOSGs in PAAD may associate with poor patients’ prognosis.

FIGURE 7

Figure 7 Verification of four DEOSG genes at protein and mRNA levels. (A) The protein expression of four DEOSG genes was analyzed in the CPTAC database. (B) The mRNA expression of four DEOSG genes was verified by qPCR in fifteen pairs of PAAD tissues and adjacent para-cancer tissues (n=15). (C) The protein expression and subcellular localization of four DEOSG genes were verified by immunohistochemistry in five pairs of PAAD and adjacent para-cancer tissues. Scale bars, 200 µm. (D) IHC scoring statistics of immunohistochemistry **P < 0.01, ***P < 0.001.

Discussion

In this paper, we have successfully screened out DEOSGs (MET, CDK1, CTTN, FYN) and constructed a four oxidative stress gene related prognostic model, which has good prediction ability, specifically with the AUC of 1-year survival more than 0.8 and the AUC of 2 years at least 0.74. Based on our four oxidative stress gene prognostic model, patients bearing low risk scores have significantly enriched immune related pathways, metabolic pathways, and DNA repair pathways. In addition, patients bearing low risk scores have high level infiltrated CD4⁺ T and CD8⁺ T cells, expressed synthetic driver genes of T cell proliferation, and achieved good immunotherapy effect. More importantly, significant differences in the protein and mRNA expression of the above four oxidative stress genes between cancer and para-cancerous were validated by CPTAC database and patients’ tissues sample. Collectively, our prognostic model indicated that oxidative stress has a great influence on PAAD patients’ prognosis, and has a significant correlation with immune infiltration, especially T cell infiltration during PAAD progression. Therefore, our prognostic model could predict the efficacy of ICI therapy.

Redox homeostasis is pivotal to mediating multiple physiological signaling pathways that are required to maintain cell metabolism, differentiation and proliferation. Once the oxidative stress balance is broken, the overall impact is often huge and harmful even if the changes of the intracellular environment are slight. When exceeding the ability of cellular self-repair, high levels of oxidative stress will lead to an increase in gene mutations and thus facilitate cancer initiation and progression by activating oncogenes and/or inhibiting tumor suppressors as well as activating multiple pro-tumorigenic signaling pathways. In PAAD, the influences of oxidative stress are multifaceted, including oxidative damage to DNA and protein, and dysregulation of cell cycle and apoptosis (13, 31). In this paper, we have found that total mutation frequency differed significantly (10717 vs 864) between patients with high or low risk scores, with the mutation frequency positively correlating with the risk scores. Among the top 10 genes with the highest mutation frequency, KRAS, TP53 and SMAD4, the highest frequency of mutations in high risk group, are also reported to be critical events for the initiation of PAAD (32). Therefore, oxidative stress may influence PAAD patients’ prognosis by affecting the frequency of gene mutation and repair. Besides, we also have investigated the influences of oxidative stress on PAAD from the perspective of the immune environment, and analyzed the differences in the infiltration immune cells and the expression of immune checkpoint related genes or newly reported synthetic driver genes of T cell proliferation between different risk groups, all of which have not been reported yet.

For the potential biological significance of our four oxidative stress gene related prognostic model to patients’ prognosis, we found that the risk score was negatively correlated to the immune pathway which including immune receptor activity, T cell receptor complex, T cell activation, and the immune infiltration of activated CD8⁺ T cell, central memory CD4⁺ T cell, effector memory CD4⁺ T cell, effector memory CD8⁺ T cell; as well as to the expression of T cell proliferation genes which including CXCL12, NFYB, FOSB, CYP27A1, GPD1 and ITM2A. Therefore, our prognostic model indicated that PAAD patients with low risk might benefit from the infiltration of CD4⁺ T and CD8⁺ T cells and the activation of corresponding immune genes or pathways. These findings are consistent with advanced cognitions about the influences of ROS on immune infiltration (28). Moreover, our results also showed that several immune checkpoint genes, for example, CD47, TNFSF9 and PVR, were highly enriched in patients with high risk scores, which indicating that our prognostic model based on the four oxidative stress genes might be a guide for individualized immunotherapy for PAAD patients.

In our four genes oxidative stress gene prognostic model, MET, FYN and CDK1 are all kinase families, being important proteins for controlling cell growth, proliferation, and differentiation. And, MET and FYN are also proto-oncogenes. MET mutation drives oncogene amplification and overexpression, which have been implicated in a variety of human cancers such as renal cell carcinoma (33, 34). Besides, MET also has a certain mutation probability in pancreatic adenocarcinoma. Mechanismly, MET protein transduces extracellular matrix signals into the cytoplasm by binding to hepatocyte growth factor/HGF ligand, and regulates many physiological processes including proliferation, dissemination, morphogenesis and survival, then contributing to the tumor process (35). FYN gene also plays an important role in the process of cancer development, including regulation of cell growth and survival, cell adhesion, cell signaling, cell motility, and immune response (36–38). It has been reported that FYN is associated with ROS through NADPH oxidases (NOX) and CO-releasing molecule-2 (CORM2) (39, 40). CDK1 is a member of the Ser/Thr protein kinase family, which plays a key role in controlling the eukaryotic cell cycle by regulating the centrosome cycle and mitotic initiation (41, 42). What’s more, CDK1 is also an important protein for autophagy regulation, the relative excessive accumulation of ROS could break cellular homeostasis, and induce autophagy (43, 44). CTTN can regulate intermolecular adhesions as well as cytoskeletal and cell adhesion structures that organize epithelial and cancer cells, which is relatively less researched (45–47). Among these four genes, only the MET gene has been reported to be related to the tumor microenvironment. In detail, MET affects the polarization of macrophages through regulating the HGF/c-MET pathway (48, 49). Herein, we found that FYN, MET, CTTN and CDK1 genes may play a role in the prognosis of PAAD patients by associating with the immune microenvironment, especially with the activation and proliferation of T cells, which are our new findings and deserve more attention.

As a highly lethal malignancy, few PAAD patients lived over 3 years (50), thus our model only predicts overall survival for 1 or 2 years. Due to the fact that there were less than 200 PAAD patient samples with survival data and only four samples in the normal group in the TCGA database, the accuracy of the prognostic model might be compromised. Therefore, with the purpose to enlarge data volume, the selection criteria for DEGs were relatively relaxed (|log₂FC| > log₂(1.2), P < 0.05), and thus more samples in three GEO datasets (sample for GSE28735, GSE62165 and GSE62452 datasets were 90, 131, 130, respectively) were analyzed, so does normal samples (total 123 cases). In the future, we will further optimize our model with more data and more samples

Conclusions

In summary, we have constructed a four oxidative stress related prognostic model for patients with PAAD. In our model, the risk score was significantly related to immune infiltration levels, especially the immune response, activation of T cells and the efficiency of immunotherapy. Our model might be promising a guidance for prognostic prediction and efficacy monitoring of clinical individualized therapy.

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 authors.

Ethics statement

The studies involving human participants were reviewed and approved by the First Affiliated Hospital of the air force Medical University(KY20213199-1). The patients/participants provided their written informed consent to participate in this study.

Author contributions

HW, R-FT and XL were involved in the data analyses, carried out the experiments, wrote, reviewed, and edited the manuscript. JF and Z-CD, J-JZ, X-YF contributed to prepare figures, data analyses and reviewed the manuscript. Z-NC, D-SY and LL contributed to the discussion, and reviewed the manuscript, HW, Z-NC, D-SY and LL conceived the study, designed and oversaw the study, evaluated data and revised the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the National Natural Science Foundation of China (81872349, 31571469, 81872482 and 82002940).

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.

Supplementary material

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

References

1. Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer statistics, 2021. CA Cancer J Clin (2021) 71(1):7–33. doi: 10.3322/caac.21654

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Rahib L, Smith BD, Aizenberg R, Rosenzweig AB, Fleshman JM, Matrisian LM. Projecting cancer incidence and deaths to 2030: The unexpected burden of thyroid, liver, and pancreas cancers in the united states. Cancer Res (2014) 74(11):2913–21. doi: 10.1158/0008-5472.CAN-14-0155

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Grossberg AJ, Chu LC, Deig CR, Fishman EK, Hwang WL, Maitra A, et al. Multidisciplinary standards of care and recent progress in pancreatic ductal adenocarcinoma. CA Cancer J Clin (2020) 70(5):375–403. doi: 10.3322/caac.21626

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Gorrini C, Harris IS, Mak TW. Modulation of oxidative stress as an anticancer strategy. Nat Rev Drug Discovery (2013) 12(12):931–47. doi: 10.1038/nrd4002

CrossRef Full Text | Google Scholar

5. van der Pol A, van Gilst WH, Voors AA, van der Meer P. Treating oxidative stress in heart failure: Past, present and future. Eur J Heart Fail (2019) 21(4):425–35. doi: 10.1002/ejhf.1320

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Zuo J, Zhang Z, Li M, Yang Y, Zheng B, Wang P, et al. The crosstalk between reactive oxygen species and noncoding rnas: From cancer code to drug role. Mol Cancer (2022) 21(1):30. doi: 10.1186/s12943-021-01488-3

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Sies H. Oxidative stress: A concept in redox biology and medicine. Redox Biol (2015) 4:180–3. doi: 10.1016/j.redox.2015.01.002

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Reuter S, Gupta SC, Chaturvedi MM, Aggarwal BB. Oxidative stress, inflammation, and cancer: How are they linked? Free Radic Biol Med (2010) 49(11):1603–16. doi: 10.1016/j.freeradbiomed.2010.09.006

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Czarny P, Wigner P, Galecki P, Sliwinski T. The interplay between inflammation, oxidative stress, DNA damage, DNA repair and mitochondrial dysfunction in depression. Prog Neuropsychopharmacol Biol Psychiatry (2018) 80(Pt C):309–21. doi: 10.1016/j.pnpbp.2017.06.036

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Salganik RI, Albright CD, Rodgers J, Kim J, Zeisel SH, Sivashinskiy MS, et al. Dietary antioxidant depletion: Enhancement of tumor apoptosis and inhibition of brain tumor growth in transgenic mice. Carcinogenesis (2000) 21(5):909–14. doi: 10.1093/carcin/21.5.909

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Brown NS, Bicknell R. Hypoxia and oxidative stress in breast cancer. oxidative stress: Its effects on the growth, metastatic potential and response to therapy of breast cancer. Breast Cancer Res (2001) 3(5):323–7. doi: 10.1186/bcr315

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Edderkaoui M, Hong P, Vaquero EC, Lee JK, Fischer L, Friess H, et al. Extracellular matrix stimulates reactive oxygen species production and increases pancreatic cancer cell survival through 5-lipoxygenase and nadph oxidase. Am J Physiol Gastrointest Liver Physiol (2005) 289(6):G1137–47. doi: 10.1152/ajpgi.00197.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Martinez-Useros J, Li W, Cabeza-Morales M, Garcia-Foncillas J. Oxidative stress: A new target for pancreatic cancer prognosis and treatment. J Clin Med (2017) 6(3):29. doi: 10.3390/jcm6030029

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Azad N, Rojanasakul Y, Vallyathan V. Inflammation and lung cancer: Roles of reactive Oxygen/Nitrogen species. J Toxicol Environ Health B Crit Rev (2008) 11(1):1–15. doi: 10.1080/10937400701436460

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Sosa V, Moline T, Somoza R, Paciucci R, Kondoh H, ME LL. Oxidative stress and cancer: An overview. Ageing Res Rev (2013) 12(1):376–90. doi: 10.1016/j.arr.2012.10.004

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Visconti R, Grieco D. New insights on oxidative stress in cancer. Curr Opin Drug Discovery Devel (2009) 12(2):240–5.

Google Scholar

17. Zhang M, Shao W, Yang T, Liu H, Guo S, Zhao D, et al. Conscription of immune cells by light-activatable silencing nk-derived exosome (Lasneo) for synergetic tumor eradication. Adv Sci (Weinh) (2022) 9(22):e2201135. doi: 10.1002/advs.202201135

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Lee GY, Han SN. The role of vitamin e in immunity. Nutrients (2018) 10(11):1614. doi: 10.3390/nu10111614

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Bhattacharyya S, Saha J. Tumour, oxidative stress and host T cell response: Cementing the dominance. Scand J Immunol (2015) 82(6):477–88. doi: 10.1111/sji.12350

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Zuo J, Zhang Z, Luo M, Zhou L, Nice EC, Zhang W, et al. Redox signaling at the crossroads of human health and disease. MedComm (2020) (2022) 3(2):e127. doi: 10.1002/mco2.127

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Koboldt DC, Zhang Q, Larson DE, Shen D, McLellan MD, Lin L, et al. Varscan 2: Somatic mutation and copy number alteration discovery in cancer by exome sequencing. Genome Res (2012) 22(3):568–76. doi: 10.1101/gr.129684.111

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for rna-seq data with Deseq2. Genome Biol (2014) 15(12):550. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma powers differential expression analyses for rna-sequencing and microarray studies. Nucleic Acids Res (2015) 43(7):e47. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Iasonos A, Schrag D, Raj GV, Panageas KS. How to build and interpret a nomogram for cancer prognosis. J Clin Oncol (2008) 26(8):1364–70. doi: 10.1200/JCO.2007.12.9791

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U.S.A. (2005) 102(43):15545–50. doi: 10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Hackl H, Charoentong P, Finotello F, Trajanoski Z. Computational genomics tools for dissecting tumour-immune cell interactions. Nat Rev Genet (2016) 17(8):441–58. doi: 10.1038/nrg.2016.67

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Barbie DA, Tamayo P, Boehm JS, Kim SY, Moody SE, Dunn IF, et al. Systematic rna interference reveals that oncogenic kras-driven cancers require Tbk1. Nature (2009) 462(7269):108–U22. doi: 10.1038/nature08460

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Legut M, Gajic Z, Guarino M, Daniloski Z, Rahman JA, Xue X, et al. A genome-scale screen for synthetic drivers of T cell proliferation. Nature (2022) 603(7902):728–35. doi: 10.1038/s41586-022-04494-7

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med (2018) 24(10):1550–8. doi: 10.1038/s41591-018-0136-1

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Dong ZY, Wu SP, Liao RQ, Huang SM, Wu YL. Potential biomarker for checkpoint blockade immunotherapy and treatment strategy. Tumour Biol (2016) 37(4):4251–61. doi: 10.1007/s13277-016-4812-9

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Yu JH, Kim H. Oxidative stress and cytokines in the pathogenesis of pancreatic cancer. J Cancer Prev (2014) 19(2):97–102. doi: 10.15430/jcp.2014.19.2.97

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Deramaudt T, Rustgi AK. Mutant kras in the initiation of pancreatic cancer. Biochim Biophys Acta (2005) 1756(2):97–101. doi: 10.1016/j.bbcan.2005.08.003

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Bruni D, Angell HK, Galon J. The immune contexture and immunoscore in cancer prognosis and therapeutic efficacy. Nat Rev Cancer (2020) 20(11):662–80. doi: 10.1038/s41568-020-0285-7

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Guo R, Luo J, Chang J, Rekhtman N, Arcila M, Drilon A. Met-dependent solid tumours - molecular diagnosis and targeted therapy. Nat Rev Clin Oncol (2020) 17(9):569–87. doi: 10.1038/s41571-020-0377-z

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Recondo G, Che J, Janne PA, Awad MM. Targeting met dysregulation in cancer. Cancer Discovery (2020) 10(7):922–34. doi: 10.1158/2159-8290.CD-19-1446

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Elias D, Ditzel HJ. Fyn is an important molecule in cancer pathogenesis and drug resistance. Pharmacol Res (2015) 100:250–4. doi: 10.1016/j.phrs.2015.08.010

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Schenone S, Brullo C, Musumeci F, Biava M, Falchi F, Botta M. Fyn kinase in brain diseases and cancer: The search for inhibitors. Curr Med Chem (2011) 18(19):2921–42. doi: 10.2174/092986711796150531

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Saito YD, Jensen AR, Salgia R, Posadas EM. Fyn: A novel molecular target in cancer. Cancer (2010) 116(7):1629–37. doi: 10.1002/cncr.24879

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Matsushima S, Kuroda J, Zhai P, Liu T, Ikeda S, Nagarajan N, et al. Tyrosine kinase fyn negatively regulates Nox4 in cardiac remodeling. J Clin Invest (2016) 126(9):3403–16. doi: 10.1172/JCI85624

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Uddin MJ, Jeong J, Pak ES, Ha H. Co-Releasing molecule-2 prevents acute kidney injury through suppression of ros-Fyn-Er stress signaling in mouse model. Oxid Med Cell Longev (2021) 2021:9947772. doi: 10.1155/2021/9947772

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Asghar U, Witkiewicz AK, Turner NC, Knudsen ES. The history and future of targeting cyclin-dependent kinases in cancer therapy. Nat Rev Drug Discovery (2015) 14(2):130–46. doi: 10.1038/nrd4504

CrossRef Full Text | Google Scholar

42. Malumbres M, Barbacid M. Cell cycle, cdks and cancer: A changing paradigm. Nat Rev Cancer (2009) 9(3):153–66. doi: 10.1038/nrc2602

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Odle RI, Florey O, Ktistakis NT, Cook SJ. Cdk1, the other 'Master regulator' of autophagy. Trends Cell Biol (2021) 31(2):95–107. doi: 10.1016/j.tcb.2020.11.001

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Li L, Tan J, Miao Y, Lei P, Zhang Q. Ros and autophagy: Interactions and molecular regulatory mechanisms. Cell Mol Neurobiol (2015) 35(5):615–21. doi: 10.1007/s10571-015-0166-x

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Ramos-Garcia P, Gonzalez-Moles MA, Gonzalez-Ruiz L, Ayen A, Ruiz-Avila I, Navarro-Trivino FJ, et al. An update of knowledge on cortactin as a metastatic driver and potential therapeutic target in oral squamous cell carcinoma. Oral Dis (2019) 25(4):949–71. doi: 10.1111/odi.12913

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Luo ML, Shen XM, Zhang Y, Wei F, Xu X, Cai Y, et al. Amplification and overexpression of cttn (Ems1) contribute to the metastasis of esophageal squamous cell carcinoma by promoting cell migration and anoikis resistance. Cancer Res (2006) 66(24):11690–9. doi: 10.1158/0008-5472.CAN-06-1484

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Du XL, Yang H, Liu SG, Luo ML, Hao JJ, Zhang Y, et al. Calreticulin promotes cell motility and enhances resistance to anoikis through Stat3-Cttn-Akt pathway in esophageal squamous cell carcinoma. Oncogene (2009) 28(42):3714–22. doi: 10.1038/onc.2009.237

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Giannoni P, Pietra G, Travaini G, Quarto R, Shyti G, Benelli R, et al. Chronic lymphocytic leukemia nurse-like cells express hepatocyte growth factor receptor (C-met) and indoleamine 2,3-dioxygenase and display features of immunosuppressive type 2 skewed macrophages. Haematologica (2014) 99(6):1078–87. doi: 10.3324/haematol.2013.091405

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Comoglio PM, Trusolino L, Boccaccio C. Known and novel roles of the met oncogene in cancer: A coherent approach to targeted therapy. Nat Rev Cancer (2018) 18(6):341–58. doi: 10.1038/s41568-018-0002-y

PubMed Abstract | CrossRef Full Text | Google Scholar

50. O'Reilly EM, Oh DY, Dhani N, Renouf DJ, Lee MA, Sun W, et al. Durvalumab with or without tremelimumab for patients with metastatic pancreatic ductal adenocarcinoma: A phase 2 randomized clinical trial. JAMA Oncol (2019) 5(10):1431–8. doi: 10.1001/jamaoncol.2019.1588

PubMed Abstract | CrossRef Full Text | Google Scholar