Identification and enrichment analysis of DE-FRGs
The gene expression profiles and clinical information of 182 cases, including 178 tumor cases and 4 normal cases, were downloaded from TCGA database. Patients with a survival time <= 30 days were excluded from subsequent analysis. The clinical information of pancreatic cancer patients is presented in Table S2. A total of 232 FRGs were selected from the FerrDb database. After extracting the gene expression of FRGs from TCGA samples, we performed the differential analysis. According to the standards mentioned above, 18 DE-FRGs were acquired by the Wilcoxon test, among which 10 were upregulated, while 8 were downregulated (Fig. 1A-1B). Next, we performed GO and KEGG pathway analyses to further study the biological functions of ferroptosis in the development of pancreatic cancer. GO analysis indicated that these DE-FRGs were enriched in biological processes related to ferroptosis, such as cellular response to oxidative stress, chemical stress, toxic stress, response to oxidative stress and liver development. The cellular components and molecular functions of these DE-FRGs are shown in Fig. 1C. KEGG pathway analysis suggested that most of the DE-FRG pathways were significantly linked to ferroptosis, leishmaniasis, HIF-1 signaling pathway, fluid shear stress and atherosclerosis and phagosomes (Fig. 1D). Using STRING and Cytoscape software, we constructed a PPI network for these DE-FRGs (Fig. 1E). The correlation of these DE-FRGs is presented in Fig. 1F.
Construction of FRGs-based prognostic model
All FRGs were included in univariate Cox regression analysis for construction of the prognostic model. The results of univariate Cox analysis showed that 13 FRGs were associated with the prognosis of pancreatic cancer patients (Fig. 2A). The relative expression levels of the 13 genes in tumor and normal samples are shown in Fig. 2B. The 13 genes were subsequently input into LASSO regression to narrow the genes and effectively avoid overfitting (Fig. 2C-2D). Finally, the selected genes were analyzed by multivariate Cox regression analysis, and 4 genes were obtained.
Risk score = (-0.30685*ENPP2) + (-0.85263*ATG4D) + (0.119446*SLC2A1) + (0.636851*MAP3K5).
The median risk score was applied to divide the pancreatic cancer patients into a high-risk group or a low-risk group. Fig. 3A shows the risk score of each pancreatic cancer patient and Fig. 3B shows the relationship between the risk score and survival time. Fig. 3C shows the gene expression profiles of the four genes in the high-risk group and low-risk group by heatmap.
Prognostic model as an independent risk factor
We evaluated the clinical and pathological parameters and risk score for prognostic value by univariate and multivariate analyses. Fig. 4A demonstrates that age, N stage and risk score were associated with the OS, and Fig. 4B shows that the risk score can become an independent risk factor (HR 1.648, 95% CI 1.335-2.035, p < 0.001). In addition, the KM curve showed that the high-risk group had a poorer OS rate than did the low-risk group (Fig. 4C). ROC analysis was performed to explore the sensitivity and specificity of the model and the area under the curve (AUC) was 0.673, which was higher than that of other clinical and pathological parameters (Fig. 4D). In addition, the patients were stratified based on clinicopathologic parameters for further survival analysis. In almost all stratified analyses, this risk model can achieve good performance, except for stage III&IV and T1&2, which may be attributable to the low number of cases (Fig.S1).
Correlation between prognostic model and clinicopathological characteristics
Based on the association of these FRGs with OS in pancreatic cancer, we further examined the correlation between the gene expression levels of FRGs and clinicopathological characteristics, such as age (>= 65/< 65), gender (male/female), grade (G3&4/G1&2), stage (stage III&IV/stage I&II), T stage (T3&4/T1&2) and N stage (N1/N0). Only ATG4D was significantly associated with stage and N stage (Fig. S2A-S2B), and the risk score was closely related to T stage (Fig. S2C).
Construction of the nomogram for prediction of prognostic risk
To better determine the survival risk of pancreatic cancer patients, we integrated the prognostic model and clinicopathological features to construct a nomogram. A nomogram was constructed to predict the 1-year, 3-year and 5-year survival possibilities (Fig. 5A). The calibration curve showed that when compared with the actual 1-year, 2-year and 3-year survival rates, the survival rates predicted by the nomogram were able to achieve good agreement (Fig. 5B-5D).
Validation of the prognostic model
To validate the prognostic model, a total of 88 pancreatic cancer tissue samples from our center were used for further experimental verification. The KM curve showed that the OS of the low-risk group was notably better than that of the high-risk group (Fig. 6A). Fig. 6B shows that the AUC of the prognostic model was 0.647. Fig. 6C shows that the higher the risk score was, the lower the survival time was. The heatmap of the 4 genes is also shown in Fig. 6C. The nomogram was validated, and the 1-year and 1.5-year calibration curves are shown in Fig. 6D.
In addition, GSE57495 and GSE62452 datasets were downloaded for further assessment. Fig. S3 validates the prognostic model in GSE57495 and GSE62452.
Enrichment analyses in the high-risk group
To further study the enriched pathways in the high- and low-risk groups, the differential analysis was performed and the differential expressed genes were input into the online tool Metascape (Fig.7A). The enrichment analysis showed that the genes enriched in several biological functions, such as chemical synaptic transmission, regulation of ion transport, regulation of system process, neuronal system and inorganic cation transmembrane transport (Fig. 7B-7D). In addition, GSEA showed that the pathways enriched in the high-risk group were TGFβ signaling, HIF signaling pathway and adherens junction (Fig. 7E).
Immune landscape between the high- and low-risk groups
We assessed the infiltration degree of 22 types of immune cells in the high-risk and low-risk groups, and calculated the immune score and immune microenvironment score. Fig. 8A showed the immune cells, immune score, stromal score and estimate score in pancreatic cancer patients and risk score was significantly associated with immune score. Fig. 8B revealed that some immune cells differed significantly between high- and low-risk groups, indicating that immune cells played an important role. Macrophage M0 and M1 showed higher infiltration in high-risk group. To further examinate the relationship between the risk score and immune cell infiltration, the four risk genes were firstly studied by TIMER. ENPP2 expression was associated with B cell, CD8+T cell, CD4+T cell, macrophage, neutrophil and dendritic cell, and the other three genes were also related to some of the immune cells (Fig. 8C). Figure 8D-8I revealed that the risk score was negatively correlated with CD4+T cell (p=0.009, R=-0.199) and positively correlated with CD8+T cell (p=0.033, R=0.163). In addition, the relationship between SCNA of the four genes and immune cell infiltration was explored by TIMER. Fig. S4 revealed that SCNA may result in changes of immune infiltration levels of B cell and CD4+T cell.
Relationship between immune checkpoints and risk score
Immunotherapy played an important role in pancreatic cancer. In this study, we evaluated the relationship between risk score and classical immune checkpoints, such as PDCD1, CTLA4 and CD274. As shown in Fig. 9, immune checkpoints were correlated to risk score, as were the four genes included in the prognostic model. Fig. 9C revealed that CD274 may be a key immune checkpoint in pancreatic cancer treatment.