Clinical and gene expression data of patients who had a pathological diagnosis of PDAC were systematically collected from TCGA, ICGC, GEO, and ArrayExpress databases. Patients who survived less than 60 days were removed from the cohorts, as these deaths might be a sequela of post-surgical complication, cardiovascular disease, and mortality other than cancer. As a result, a total of 11 PDAC cohorts contained 1337 PDAC patients were obtained. Detailed clinical information of the 11 cohorts are presented in Supplementary Table S1.

The IRRGs were collected from the “Hallmark inflammatory response” and “GOBP inflammatory response” gene sets in the Molecular Signatures database (MSigDB) (http://www.gsea-msigdb.org/) (Subramanian et al., 2005; Liberzon et al., 2015). To identify IRRGs associated with patients’ overall survival (OS) in E-TMAB-6134 and PACA-AU cohorts, the univariate Cox regression analysis was carried out with a p-value threshold of 0.001 by using the R package “survival.” The R package “VennDiagram” was then used to screen intersecting genes, which were selected as OS-related IRRGs. In order to determine which OS-related IRRGs may play a role in PDAC development, we then compared their expression levels between tumor and paired adjacent tissue in two GEO cohorts (GSE15471 and GSE16515). Taken together, genes with HR <1 and decreased expression in tumor or genes with HR >1 and increased expression in tumor were used as candidates for the signature construction.

The LASSO regression analysis (Tibshirani, 1997), which was designed for variable selection and shrinkage, was conducted by using the R package “glmnet” to construct an IRRGs prognostic signature. In LASSO, variables with a regression coefficient equal to zero after the shrinkage process are excluded from the model. A tuning parameter lambda (λ) controls the amount of shrinkage, with increased shrinkage for higher λ values. An optimal λ was chosen when the partial likelihood deviance reached its lowest, which was based on the 10-fold cross-validation. The risk score of the signature for each patient was calculated as follows: R i s k   s c o r e = ∑ i = 1 n ( E x p i ∗ β i ) where n represents the number of prognostic genes, Exp _i represents the expression value of gene i, and β _i represents the regression coefficient of gene i.

We selected E-MTAB-6134 cohort as a training set because it contains the largest number of patients among all cohorts. Other 10 cohorts were used as validation sets to validate reliability of the signature. In each cohort, patients were divided into the high- or low-risk group according to the median value of the risk score. To compare the survival difference between high- and low-risk groups, Kaplan–Meier survival and receiver operating characteristic (ROC) analyses were conducted by using the R package “survival” and “timeROC.” We also calculated the time-varying concordance index (C-index) of our signature and compared it with six previous published signatures which were based on glycolysis (Song et al., 2021), ferroptosis (Feng et al., 2021a), hypoxia (Ding et al., 2021), m6A (Meng et al., 2020), autophagy (Yu J. et al., 2021), and immune (Mao et al., 2021) related genes. The information of these six signatures are listed in Supplementary Table S2. Furthermore, to verify the IRRGs signature was independent of other clinical variables (tumor grade, tumor stage, and resection margin) for the OS prediction, univariate and multivariate Cox regression analyses were performed in cohorts with clinical variables documented.

Gemcitabine response of patients in TCGA-PAAD cohort was obtained by using R package TCGAbiolinks. Patients’ risk scores were compared among four types of gemcitabine response, including complete response (CR), partial response (PR), stable disease (SD), and progressive disease (PD). Furthermore, risk scores of 808 cancer cell lines in Genomics of Drug Sensitivity in Cancer (GDSC) (Yang et al., 2013) were calculated by the IRRGs signature. The correlation between IC50 value of gemcitabine and risk score of cancer cell lines was explored.

Expression data of human cancer cell lines (CCLs) were downloaded from the Broad Institute Cancer Cell Line Encyclopedia (CCLE) project (https://portals.broadinstitute.org/ccle/) (Ghandi et al., 2019). Drug sensitivity data of CCLs were extracted from the Cancer Therapeutics Response Portal (CTRP v.2.0, https://portals.broadinstitute.org/ctrp) (Rees et al., 2016). The area under the dose-response curve (AUC) was used as a measure of drug sensitivity in CTRP, and higher AUC values indicate decreased sensitivity to drugs. Based on drug sensitivity and gene expression data from CTRP, drug response of patients was predicted according to their gene expression levels by using the R package pRRophetic (Geeleher et al., 2014). The AUC value of each drug in each patient was estimated. Additionally, the connectivity map (CMap) analysis (Subramanian et al., 2017) was also used to predicted potential small molecule drugs for the treatment of high-risk PDAC patients.

Somatic mutations of the KRAS, TP53, and CDKN2A are the common events in PDAC patients, and KRAS and TP53 mutations have been tightly linked to tumor-promoting inflammation (Yachida et al., 2012; Uehara and Tanaka, 2018; Hamarsheh et al., 2020). Thus, we compared the frequencies of KRAS, TP53, and CDKN2A mutations between high- and low-risk PDAC patients. The Mutation Annotation Format (MAF) files of TCGA-PAAD cohort were downloaded from the GDC database (https://portal.gdc.cancer.gov/). The simple somatic mutation format files of PACA-AU and PACA-CA cohorts were obtained from ICGC database (https://dcc.icgc.org/) and converted to MAF files by using R package “matfools” (Mayakonda et al., 2018). Likewise, all the MAF files were analyzed by using “matfools” package. Additionally, we extracted the KRAS, TP53, and CDKN2A mutation status of patients in E-MTAB-6134 cohort, which were documented in a clinical information file.

Next, to gain insight into the difference in tumor immune microenvironment between high- and low-risk groups, CIBERSORT algorithm (https://cibersortx.stanford.edu/) (Chen et al., 2018) was used to estimate the proportion of 22 immune cells in each sample based on their gene expression profiles. Then, we quantitatively compared the infiltration of 22 immune cells between high- and low-risk groups. The estimate algorithm was also used to calculate stromal and immune scores of each sample by using R package estimate. Gene Set Enrichment Analysis (GSEA) was conducted to explore the biological difference between high- and low-risk groups. The analysis was performed by using GSEA software, and false discovery rate (FDR) < 0.25 was regarded as statistically significant according to the GSEA user guide (Subramanian et al., 2005).

There were 717 and 200 protein coding genes contained in the “GOBP inflammatory response” and “Hallmark inflammatory response” gene sets, respectively (Supplementary Tables S3 and S4). After removing the overlapped genes, a total of 841 IRRGs were obtained for further analysis. Univariate Cox regression analysis showed that there were 34 and 31 genes significantly (p < 0.001) related to OS of PDAC patients in E-MATB-6134 and PACA-AU cohorts, respectively (Supplementary Tables S5 and S6). A total of 11 intersecting genes were identified by using Venn diagram (Figure 2A). Subsequently, expression levels of these 11 genes were compared between tumor and tumor-adjacent tissues in two GEO cohorts. Taken together, eight genes were selected for LASSO analysis. Of these, ADM, DCBLD2, EREG, ITGA5, MIF, MMP14, and TREM1 were associated with poor prognosis and significantly increased in tumor, while BTG2 was related to favorable prognosis and decreased in tumor (Figures 2B–E).

A prognostic signature containing seven IRRGs was constructed based on the optimal value of penalty lambda (Supplementary Figures S1A,B). The risk score of each patient was calculated as follows: risk score = expression level of ITGA5*0.059 + expression level of TREM1*0.109 + expression level of EREG*0.132 + expression level of MIF*0.216 + expression level of ADM*0.204 + expression level of DCBLD2*0.129-expression level of BTG2*0.098 (Supplementary Figure S1C).

The Kaplan-Meier survival analysis showed that patients in the high-risk group had significantly unfavorable OS compared with those in the low-risk group (Figure 3A). The predictive efficacy of the risk score for OS was further evaluated by ROC analysis, and the area under the ROC curve (AUC) was 0.791, 0.688, and 0.678 for 1, 3, and 5-years, respectively (Figure 3B). Expression levels of the seven signature genes were significantly different between high- and low-risk groups (Figure 3C). Patients with a high-risk score had a higher death probability than those with a low-risk score (Figure 3D). The high- and low-risk groups displayed different distributed patterns in a PCA plot (Figure 3E). Correlation analysis indicated that the seven signature genes were significantly correlated with each other (Figure 3F).

We then verified the reliability and capability of the IRRGs signature in 10 PDAC cohorts. Kaplan-Meier survival analysis showed the high-risk group had a significantly worse prognosis compared to the low-risk group in all 10 cohorts (Figures 4A–J).

We conducted ROC analysis for 1, 3, and 4-years OS in the 10 validation cohorts to further validate the reliability of the signature. The AUC for 1, 3, and 4-years OS in the 10 PDAC cohorts were greater than 0.600 (Figures 5A–J). Notably, the AUC of ROC curve for 4-years OS reached 0.847, 0.829, and 0.817 in GSE71729, GSE79668, and GSE85916 cohorts, respectively. The C-index of our signature was greater than 0.600 at different survival times except for in the TCGA-PAAD cohort (Supplementary Figure S2). Additionally, compared with six previously published gene signatures, our IRRGs signature had the highest C-index in 9 of 11 PDAC cohorts.

We also explored whether our signature could be used for prediction of other important survival outcomes, including disease-free survival (DFS), disease-free interval (DFI), disease-specific survival (DSS), and progression-free interval (PFI). Kaplan-Meier survival analysis showed the patients with high-risk had a significantly worse DFS than their low-risk counterparts in the E-MTAB-6134 cohort (Figure 6A). In the PACA-AU and PACA-CA cohorts, patients with high-risk had a significantly worse DFI than those with low-risk (Figures 6B,C). In the TCGA-PAAD cohort, patients with high-risk had a significantly worse DFI and DSS than those with low-risk (Figures 6D,E). However, there were no significant differences in DFI between two groups in TCGA-PAAD cohort (Figure 6F).

In the univariate Cox analyses, high-risk was significantly correlated with poor OS (E-MTAB-6134 cohort: HR = 2.506, 95% CI = 1.842–3.411, p < 0.001; PACA-AU cohort: HR = 2.240, 95% CI = 1.597–3.143, p < 0.001; TCGA-PAAD cohort: HR = 1.651, 95% CI = 1.058–2.577, p = 0.027; PACA-CA cohort: HR = 2.080, 95% CI = 1.340–3.229, p = 0.001) (Figures 7A–D). The multivariate Cox analyses indicated that high-risk was a hazardous factor for OS after correcting for other significant clinical factors (E-MTAB-6134 cohort: HR = 2.499, 95% CI = 1.829–3.415, p < 0.001; PACA-AU cohort: HR = 2.048, 95% CI = 1.441–2.910, p < 0.001; PACA-CA cohort: HR = 1.897, 95% CI = 1.261–2.960, p = 0.005) (Figures 7E–G). ROC analyses showed that the signature had good predictive accuracy in OS of PDAC patients than other clinical variables (Figures 7H–K).

Based on the multivariate analysis results of E-MTAB-6134 cohort, we established a nomogram which could further improve survival predictive ability for PDAC patients. A total of four factors were integrated in the nomogram to predict the OS of PDAC patients (Figure 8A). The total points were calculated by adding up the corresponding points of each factor. The calibration curve suggested that the nomogram that predicted the survival rate was close to the actual situation for the 1-, 3-, and 5-year survival (Figure 8B). The time-varying AUC and C-index plots indicated that the predictive ability of the nomogram was better than that of any single factor (Figures 8C,D).

We conducted a subgroup analysis of 62 patients treated with gemcitabine in TCGA-PAAD cohort. As shown in Figure 9A, patients presented a CR had significantly lower risk score compared with those presented a PD. The survival analyses showed that patients in the high-risk group had significantly poor OS, DSS, DFI, and PFI as compared with those in the low-risk group (Figures 9B–E). By analyzing gene expression and drug sensitivity data in the GDSC database, we revealed a significant positive association between cancer cells’ risk score and gemcitabine IC50, indicating that patients with higher risk score were more likely to be resistant to gemcitabine (Figures 9F,G). We unutilized the R package pRRophetic, which had a built-in ridge regression model, to predict the drug response of patients in E-MTAB-6134 cohort based on their gene expression profiles. The estimated AUC value of each compound in each sample was obtained. As a result, a total of 354 compounds and their estimated AUC value of each patient in E-MTAB-6134 cohort were yielded and listed in Supplementary Table S7. We found that there was significantly positive correlation between risk score and AUC value of gemcitabine (Figure 9H).

We further analyze the results in Supplementary Table S6 to identify potential drugs for patients with a high-risk score. Correlation analysis between AUC value and risk score was conducted to pick compounds with negative correlation coefficient (R < −0.3). Differential analysis was then carried out between high- and low-risk groups to select compounds with lower AUC values in the high-risk group (p < 0.001 and mean_high/mean_low < 0.98). Collectively, these two analyses yielded six compounds (simvastatin, dasatinib, pluripotin, fluvastatin, BMS-536924, curcumin) (Figures 10A,B). Furthermore, there were 23 up-regulated and 34 down-regulated genes in the high-risk group as compared with the low-risk group in E-MTAB-6134 cohort (Figure 10C). The CMap analysis was then used to identify compounds of which gene expression patterns were oppositional to the expression patterns of the high-risk group (i.e., gene expression increased in the high-risk group tend to be decreased by the perturbagen of certain compounds). We found that BMS-536924 and dasatinib had CMap scores <−87 (Figure 10D), indicating that these two drugs had great potential for the treatment of high-risk patients.

We found that KRAS mutation had the most frequency in PACA-AU, PACA-CA, and E-MTAB-6134 cohorts, while TP53 mutation was the most common in TCGA-PAAD cohort (Figures 11A–C). The frequency of KRAS mutation were significantly higher in the high-risk group in all four cohorts, while the frequency of TP53 mutation was significantly higher in the high-risk group in E-MTAB-6134 and TCGA-PAAD cohorts. However, the frequency of CDKN2A mutation was significantly higher in the high-risk group only in E-MTAB-6134 cohort. We next conducted a pooled analysis to compare the overall frequencies of TP53, KRAS, and CDKN2A mutations between the high- and low-risk groups in four cohorts. We found that frequencies of TP53, KRAS, and CDKN2A mutations were significantly higher in the high-risk group compared with the low-risk group (KRAS: OR = 4.47, 95% CI = 2.71–7.39, p < 0.001; TP53: OR = 2.27, 95% CI = 1.60–3.20, p < 0.001; CDKN2A: OR = 1.95, 95% CI = 1.24–3.09, p < 0.001) (Figures 11D–F).

Based on the CIBERSORT algorithm, we compared the proportion of 22 types of immune cells between high- and low-risk groups. As shown in Figure 12A, CD8⁺ T, naïve B, plasma, and macrophages M1 cells were significantly higher in the low-risk groups, while macrophages M0, neutrophils, and macrophages M2 cells were more likely higher in the high-risk groups. Next, we used ESTIMATE algorithm to calculate stromal and immune scores of each tumor sample. We found that the immune scores were significantly negative correlation with risk scores in TCGA-PAAD, E-MTAB-6134, GSE71729, and PACA-AU cohorts, while the stromal scores were significantly positive correlation with risk scores in E-MTAB-6134 and PACA-CA cohorts (Figure 12B). We also explored the correlation of risk scores and expression of 11 well-documented immune checkpoints. We found that expression levels of immune checkpoints were significantly positive correlation with risk scores in most cohorts except for TIGIT and CTLA4 (Figure 12C). Notably, CD73 and CD276 were extremely positive correlation with risk scores in all 11 cohorts. GSEA analysis indicated that most of the 50 hallmark gene sets were upregulated in high-risk groups although there were some differences among the 11 cohorts (Supplementary Figure S3). Especially, the glycolysis and hypoxia gene sets were significantly upregulated in the high-risk groups compared with low-risk groups in 10 cohorts.