We collected the transcriptome profiles of pancreatic adenocarcinoma (PAAD) available in the TCGA database (https://portal.gdc.cancer.gov/) and Genotype-Tissue Expression (GTEx) GTEx-Pancreas datasets (https://xenabrowser.net/) on July 5^th, 2020. Our study included the expression profile of 171 normal samples and 177 PAAD samples. The clinicopathological information including gender, age, tumor grade, T classification, N classification, M classification, TNM stage, follow-up time, and survival status of the patients from TCGA-PAAD was also retrieved (10). We excluded samples with follow-up time shorter than 30 days and samples with missing clinicopathological information. For validation, gene expression data and clinical data of 70 PAAD patients in GSE62452 was downloaded from Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). The flow of the study was shown in Figure S1 .

Two macrophage phenotype switch related gene sets (GSE5099_CLASSICAL_M1_VS_ALTERNATIVE_M2_MACROPHAGE_UP, GSE5099_CLASSICAL_M1_VS_ALTERNATIVE_M2_MACROPHAGE_DN) (11) were selected from the Molecular Signatures Database v7.1 (MSigDB) (12) (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp). No overlapped genes were detected and a total of 382 genes were retrieved ( Table S1 ). The expression data of the 382 MRGs was extracted. Since the data was retrieved from public datasets and we followed the respective publishing guidelines, no ethics approval was needed.

We collected the mean expression data of 382 MRGs comprising 177 PAAD and 171 non-tumor samples. The mean expression values were then normalized by log₂ transformation. The differentially expressed MRGs between tumor and normal samples were identified using the Wilcoxon signed-rank test in R (version 4.0.2, https://www.r-project.org/) with a threshold of |log(foldchange)| > 1 and a false discovery rate (FDR) < 0.05.

We built the prognostic signatures in three steps: 1) we conducted univariate Cox regression analysis of the MRGs by “survival” in R. Genes with P < 0.05 were chosen for Lasso Cox regression analysis. 2) The Lasso Cox regression analysis (13) which was implemented with “glmnet” and “survival” packages in R was utilized to remove highly correlated genes and to prevent over-fitting. 3) Risk scores were calculated as the sum of each gene’s expression levels multiplying the regression coefficient in the multivariate Cox regression model. Patients were dichotomized into high-risk group and low-risk group by median value of risk score. The Kaplan-Meier method was utilized to compare the survival outcome between high-risk group and low-risk group. Time dependent receiver operating characteristic (ROC) curves by R package “survivalROC” (14) was also used to determine the efficacy of the prognostic model.

We developed a prognostic nomogram predicting OS based on the Cox proportional hazard regression model by the “rms” package in R. A concordance index (C-index) was calculated to evaluate the performance of the nomogram.

Using R package “GSVA,” single-sample gene set enrichment analysis (ssGSEA) was conducted to quantify the activity or enrichment levels of immune cells, functions, or pathways (15, 16). The comparison of immune cell distribution between patients at high-risk group and patients at low-risk group was performed.

Besides, the relationships between the MRGs based risk score and several immune checkpoints, including PD-1, PD-L1, CTLA-4, TIM-3, and LAG-3, were analyzed.

The mutations of the four most prominent key drivers in pancreatic cancer progression, namely KRAS, TP53, CDKN2A, and SMAD4 (17), were downloaded. And the relationships between them and MRGs based risk score were gauged.

Gene set enrichment analysis (GSEA) (12) was performed in javaGSEA v. 4.1.0 based on the Molecular Signatures Database v. 7.1. C2 (curated gene sets), C5 (GO gene sets), and C6 (oncogenic signatures) were searched to identify enriched KEGG pathways, biological processes, cellular components, molecular functions, and dysregulated oncogenic signatures associated with the high-risk group (12). |NES| > 1 and FDR < 0.05 were considered statistically significant.

Human pancreatic ductal epithelium cell line HPDE6-C7 was obtained from Institute of Biochemistry and Cell Biology, Chinese Academy of Sciences (Shanghai, China). PDAC cell lines (MIA PaCa-2 and Capan-1) were acquired from the American Type Culture Collection (ATCC, Manassas, VA, USA). HPDE6-C7 was cultured in Dulbecco’s modified Eagle’s medium (DMEM) (Invitrogen, Carlsbad, CA, USA) added with 10% fetal bovine serum (FBS) (Invitrogen). MIA PaCa-2 was cultured in DMEM added with 10% FBS and 2.5% horse serum. Capan-1 was cultured in Iscove’s modified Dulbecco’s medium added with 20% FBS. All cells were cultured in a humidified 5% CO₂ incubator at a temperature of 37°C.

Total RNA was extracted from cells using the TRIzol reagent (Invitrogen). After quantification using NanoDrop 2000c instrument (Thermo Fisher, Waltham, MA, USA), the RNA was reverse transcribed into first-strand cDNA using the M-MLV Reverse Transcriptase (Invitrogen). qRT-PCR was carried out on StepOne Plus Real-Time PCR System (Applied Biosystems, Foster City, CA, USA) using TB Green^® Premix Ex Taq™ (Tli RNaseH Plus), ROX plus (TaKaRa, Dalian, China). GAPDH was used as an endogenous control. The quantification of RNA level expression was calculated using the 2^−ΔΔct method. The primer sequences used were as follows: for KIF23, 5’-AGTCAGCGAGAGCTAAGACAC-3’ (sense), 5’-GGTTGAGTCTGTAGCCCTCAG-3’ (antisense); for BIN1, 5’-TGAGCAGTGCGTCCAGAATTT-3’ (sense), 5’-CGATCTTGTTTGCCTCATCCC-3’ (antisense); for LAPTM4A, 5’-ATGGTGTCCATGAGTTTCAAGC-3’ (sense), 5’-CCACAGTCAGCAAAATTGCCA-3’ (antisense); for ERAP2, 5’-CACTAATGGGGAACGATTTCCTT-3’ (sense), 5’-CTGACCAAGACTTCGATCTTCTC-3’ (antisense); for ATP8B2, 5’-CGGGCTAATGACCGAGAATAC-3’ (sense), 5’-CTGCTCAAAGAGGTTGACAGG-3’ (antisense); for FAM118A, 5’-GTCGCCCATGATCTGATCCG-3’ (sense), 5’-CTCCAGGTCGTCAAACACCTC-3’ (antisense); for RGS16, 5’-ATCAGAGCTGGGCTGCGATA-3’ (sense), 5’-CAGGTCGAACGACTCTCTCC-3’ (antisense); for ELMO1, 5’-GGAGCAGGTTATGAGAGCACT-3’ (sense), 5’-GGGCGGGACTGGAAATCTTC-3’ (antisense); for RAPGEFL1, 5’-AGGGGCTGCTTCAAGAGGA-3’ (sense), 5’-CCCTGGTAAAGGGACTCGT-3’ (antisense).

All statistical analyses were conducted using R software (version 4.0.2), GraphPad Prism 7 (San Diego, CA, USA), and SPSS 20.0 (Chicago, IL, USA). The median OS was compared using the Kaplan-Meier method with log-rank test. The distributions of clinicopathological parameters between the high-risk and low-risk groups were compared using chi-square tests, T test, or one-way analysis. Two-sided P value less than 0.05 was considered as statistically significant.

In total, 198 differentially expressed MRGs meeting the criteria as |log(foldchange)| > 1 and FDR < 0.05 by the pooled analysis of TCGA-PAAD and GTEx-pancreas databases were obtained ( Figure 1A , Table S2 ).

The univariate Cox regression analysis and identified 59 MRGs that were significantly related to overall survival in TCGA-PAAD cohort ( Figure 1B ). The Lasso regression analysis excluded 40 genes that may be highly correlated with other genes ( Figures 1C, D ). The rest 19 genes (KIF23, BIN1, GBP4, LAPTM4A, ERAP2, SLC22A18AS, CDK6, KIF20B, RNFT2, TES, ATP8B2, FAM118A, ARID5A, RGS16, TNFRSF21, UBASH3B, GPRIN1, ELMO1, RAPGEFL1) were then submitted to a multivariate Cox proportional hazards model, resulting in 9 candidate genes (KIF23, BIN1, LAPTM4A, ERAP2, ATP8B2, FAM118A, RGS16, ELMO1, RAPGEFL1) as significant predictors of prognosis ( Table 1 ). The formula for the risk score of every PAAD patient was constructed: risk score = (expression value of KIF23 * 0.394) − (expression value of BIN1 * 0.098) + (expression value of LAPTM4A * 0.005) + (expression value of ERAP2 * 0.051) − (expression value of ATP8B2 * 0.082) − (expression value of FAM118A * 0.081) + (expression value of RGS16 * 0.014) − (expression value of ELOM1 * 0.112) + (expression value of RAPGEFL1 * 0.017).

The median overall survival time of patients in the high-risk group was significantly lower than that in the low-risk group ( Figure 1E , P < 0.001). Expression profile of the included 9 MRGs ( Figure 2A ), distribution of patients ( Figure 2B ), individual survival status ( Figure 2C ) stratified by risk score were shown in Figure 2 . Risk score remained to be an independent prognostic indicator in multivariate analysis (HR=3.129, 95% CI = 1.950–5.022, P < 0.001, Table 2 , Figures 3A, B ). Then we assessed the prediction efficiency of the risk score by drawing the time dependent ROC curve, which revealed that the risk score could predict the 1-year overall survival (AUC=0.765, Figure 3C ), 3-year overall survival (AUC =0.793, Figure 3D ), and 5-year overall survival (AUC =0.776, Figure 3E ) for PAAD patients effectively. And the AUC of risk score was larger than those of the other enrolled clinocopathological factors.

We verified the robustness of our MRG based risk score with the patients from GSE62452 dataset. We calculated the risk score for patients and divided them into high-risk group and low-risk group in test dataset with the same formula derived from the training cohort as the TCGA-PAAD database. Patients in low-risk group had obviously better overall survival (P < 0.001, Figure 4A ). After adjusting for clinicopathological features, risk score remained to be an independent prognostic indicator in multivariate analysis (HR=2.385, 95% CI = 1.209-4.707, P = 0.012, Table 3 ). Then we assessed the prediction efficiency of the risk signature by drawing the time dependent ROC curve, which revealed that the risk score could predict the 1-year overall survival rate (AUC = 0.617, Figure 4B ), 3-year overall survival rate (AUC = 0.878, Figure 4C ), and 5-year overall survival rate (AUC = 0.761, Figure 4D ) for GSE62452 dataset effectively.

The 166 patients with complete clinical information from the TCGA-PAAD dataset were used to establish a prognostic nomogram based on the stepwise Cox regression model ( Figure 3F ). The MRGs based risk score, T classification, and N classification were finally included into the stepwise Cox regression model for the nomogram. The C-index of the nomogram was 0.723, which was higher than that of TNM stage itself as 0.535. Thus, the nomogram was regarded to be superior to the TNM stage in terms of predicting overall survival of PAAD.

Patients with high risk score had significantly higher level of APC co-inhibition (Cor=0.231, P = 0.002, Figure 5A ), MHC class I (Cor=0.190, P=0.013, Figure 5B ), parainflammation (Cor=0.167, P=0.029, Figure 5C ), and type I IFN response (Cor=0.188, P=0.014, Figure 5D ), but significantly lower level of B cell (Cor = −0.180, P = 0.019, Figure 5E ), pDCs (Cor = −0.185, P = 0.015, Figure 5F ), TIL (Cor = −0.166, P = 0.030, Figure 5G ), and type II IFN response (Cor = −0.270, P < 0.001, Figure 5H ).

Relationship between expression level of co-inhibitory immune checkpoints (PD-1, PD-L1, CTLA-4, TIM-3, LAG-3, CD8A) and MRGs based risk score were also assessed. Patients with high risk score had significantly higher level of CD274 (Cor=0.285, P < 0.001, Figure 6A ), IDO1 (Cor=0.418, P < 0.001, Figure 6C ), but significantly lower level of CD8A (Cor = −0.176, P = 0.022, Figure 6B ).

Via the analysis of the relationship between the mutations of the four “key drivers” of pancreatic cancer and risk score, we found that patients with higher risk score were enriched for KRAS mutations (P<0.001, Figure 7A ) and TP53 mutations (P < 0.001, Figure 7B ).

Relationship between MRGs prognostic index and clinicopathological features were subsequently analyzed ( Figure 8 ). Significant higher risk scores were in patients with higher tumor grade (P < 0.001). In addition, the expression level of BIN1 was significantly related with advanced M (P < 0.001). The expression level of ELMO1 was significantly related with higher tumor grade (P = 0.008). The expression level of ERAP2 was significantly related with advanced N (P = 0.027). The expression level of KIF23 was significantly related with higher tumor grade (P < 0.001). The expression level of LPTM4A was significantly related with younger age (P=0.038) and higher tumor grade (P=0.024). The expression level of RGS16 was significantly related with higher tumor grade (P=0.015).

GSEA compared the high- and low risk groups stratified by the risk score ( Figure 9 ). KEGG pathways enriched in the high-risk group included KEGG P53 signaling pathway, KEGG cell cycle, KEGG base excision repair, KEGG pancreatic cancer, and KEGG mismatch repair. KEGG pathways enriched in the low-risk group included KEGG glycine serine and threonine metabolism, KEGG chemokine signaling pathway, KEGG complement and coagulation cascades, KEGG cytokine cytokine receptor interaction, and KEGG T cell receptor signaling pathway.

The mRNA expressions of the nine MRGs (KIF23, BIN1, LAPTM4A, ERAP2, ATP8B2, FAM118A, RGS16, ELMO1, RAPGEFL1) in pancreatic ductal epithelium cell line (HPDE6-C7) and PDAC cell lines (MIA PaCa-2 and Capan-1) were assessed. Transcriptional expressions of the 9 MRGs were significantly higher in PDAC cell lines than those in HPDE6-C7 ( Figure 10 ).