The general process of this study is shown below ( Figure 1 ). First, 178 PAAD samples and 4 normal tissue samples were downloaded and collected from the The Cancer Genome Atlas (TCGA, http://portal.gdc.cancer.gov) database to obtain their quantitative gene expression data and clinical data. The TCGA-PAAD cohort was used as a training cohort. Transcriptome data from 167 normal pancreatic tissues were obtained from the Genotype-Tissue Expression (GTEx, https://gtexportal.org/) database. Matching TCGA with the data from GTEx. Next, 183 PAAD samples were downloaded from the ICGC (https://icgc.org/) database as an external validation cohort. The dataset of IRGs was obtained from the Immport (https://www.immport.org) database. The dataset of immune-related TFs was obtained from the Cistrome (https://cistrome.org) database. In addition, mutation data of PAAD and maf files of mutect versions of somatic mutations were downloaded from the TCGA database and used to analyze TP53 mutations and evaluate TMB scores. Single cell sequencing samples of 3 PAAD tissues and 1 adjacent normal tissue were obtained from GSE197177 of the Gene Expression Omnibus database (GEO, https://www.ncbi.nlm.nih.gov/geo/).

The differentially expressed genes were screened by the Wilcox test using R software on the obtained data. The pheatmap package was used to draw volcano and heat maps of the differentially expressed genes. Among the screened differentially expressed genes, differentially expressed IRGs were identified using the above method. Screening criteria were false discovery rate (FDR) < 0.05 and | log2 fold-change (FC)| > 1. COX regression analysis was performed on the differential IRGs and clinical prognostic information, which in turn filtered out prognostically relevant IRGs. In addition, 318 immune-related TFs were downloaded from Cistrome and screened for differentially expressed TFs using the same method as above. The selection criteria were the same as for the IRGs. Correlation coefficients > 0.4 and p < 0.0001 were used as thresholds for correlation analysis between differentially expressed TFs and IRGs. Cytoscape software was used to map their interaction regulatory networks.

The limma package was used to perform unsupervised clustering analysis on 178 PAAD patients from TCGA to identify different immune-related molecular subtypes. The ConsensusClusterPlus package was used to perform a consensus clustering algorithm with 1000 iterations to ensure classification stability. The survivor and survminer packages were used to plot Kaplan-Meier (K-M) survival curves. The estimate package was used to evaluate StromalScore, ImmuneScore, and ESTIMATEScore for 3 subtypes.

Discovering differentially expressed genes (DEGs) with significant prognostic value were identified using univariate Cox regression analysis. LASSO-Cox regression analysis was performed using the glmnet package to further narrow down the candidate genes (10). The prognostically relevant IRGs were selected and the TiME-score was calculated for PAAD patients (β: coefficients, Exp: gene expression level):

The TCGA-PAAD cohort was divided into high-risk and low-risk subgroups based on the median TiME-score. K-M curves and receiver operating characteristic (ROC) curves were used to assess the accuracy and specificity of the TiME-score. Moreover, the ICGC dataset was used as an external validation set to verify the predictive value of the signature.

To determine whether the TiME-score could be used as an independent predictor of prognosis, independent prognostic analyses were performed using univariate and multifactorial Cox. Patients with PAAD were divided into groups based on clinical information such as age (≤65 or >65 years), gender, pathological grade, clinical stage, T (T1-2 or T3-4), N (N0 or N1), and M (M0 or M1). On this basis, the TiME-scores of the patients were compared between groups. The beeswarm package was used to assess the correlation between prognostic characteristics and the above clinical parameters. Additionally, a nomogram was created using the rms package to accurately predict the 1-year, 2-year, and 3-year overall survival (OS) of PAAD patients. Also, the predictive performance of the nomogram was evaluated using calibration curves.

The ggplot2 package was used to analyze the differential expression of the hub genes in the TiME-score among the TCGA-PAAD cohort and the GTEx cohort. The Human Protein Atlas (HPA, https://www.proteinatlas.org/) database (11) is a database that provides extensive transcriptomic and proteomic data for specific human tissues and cells. Samples used for immunohistochemistry by the HPA do not correspond to samples in the TCGA dataset. All tissues were collected from the Uppsala Biobank and RNA samples were extracted from frozen tissue sections. The differences in protein expression of hub genes between PAAD and normal pancreatic tissues were validated using Immunohistochemistry in HPA. The Clinical Proteomic Tumor Analysis Consortium database (CPTAC, https://proteomics.cancer.gov/programs/cptac) is a research project designed to accelerate understanding of cancer through the application of large-scale proteomic and genomic analysis or proteomics. The University of ALabama at Birmingham CANcer data analysis Portal (UALCAN, https://ualcan.path.uab.edu/index.html) was used to analyze protein expression in the CPTAC database.

Single-cell RNA sequencing count matrices for all 3 PAAD tissues (GSM5910784, GSM5910787 and GSM5910789) were downloaded from GSE197177. The 3 matrices were combined into a single Seurat object using the CreateSeuratObject function (Seurat R package, version 4.3.0). Cells with >2500 or <200 genes detected in a single cell, or with UMI counts >5% from mitochondrial sources are considered low quality cells. These low-quality cells were filtered out and single-cell data were normalized. Next, the FindNeighbors function was used to identify similar cells. The FindClusters function (resolution = 0.4) was used to identify major cell clusters. Based on this, DimPlot was used to visualize the UMAP non-linear dimensionality reduction results of these cell clusters. Cell clusters were identified using the “SingleR” package and annotated based on marker genes obtained from previous studies. The accuracy of the annotation was verified using the plotScoreHeatmap function. The number and distribution of target gene expressions were marked on UMAP using the FeaturePlot function. In addition, the same analysis described above was done on a single cell sample (GSM5910786) from 1 PAAD adjacent to normal tissue in GSE197177, which served as a control.

An overview graph of TP53 mutations in PAAD was mapped using the maftools package. Using gene set enrichment analysis (GSEA) software (version: v4.2.3; www.gsea-msigdb.org/gsea/downloads.jsp) (12), determine immune-related biological process differences between TP53 samples from the TCGA-PAAD cohort (wild group n=62, mutant group n=81). An annotated gene set file (c5.bp.v7.1.symbols.gmt) was selected as the reference gene set (p<0.05). The relative proportion of all 22 tumor-infiltrating immune cells (TICs) in each sample was estimated using the CIBERSORT package (number of cycles = 1000) (13). This algorithm was run more than three times with identical results. Subsequently, the infiltration of immune cells and the correlation between immune cells in the wild and mutant groups were analyzed sequentially.

Considering that both RARRES3 and TNFSF10, the hub genes of TiME-score, are downstream genes of TP53. K-M curves and ROC curves of RARRES3 high and low expression groups were plotted with survival and survminer packages. The effects of age, gender, grade, stage, and RARRES3 expression on independent prognosis were assessed by univariate and multifactorial Cox analysis. The same method was used to analyze TNFSF10 as control and validation. In addition, the correlation between RARRES3 and clinical indicators was evaluated.

Enrichment analysis of RARRES3 was performed using the clusterProfiler package and the GSEA enrichment results of RARRES3 were plotted using the enrichplot package. This was used to assess the abundance of different immune-related functions or pathways associated with RARRES3 expression. Next, the abundance of all 22 TICs was estimated using the CIBERSORT algorithm, which in turn ranked the correlation with RARRES3 expression for each immune cell. Moreover, the correlation between RARRES3 and immune checkpoints was analyzed by the corrplot package. TNFSF10 analysis was applied in the same way, as control and validation.

The enrichplot package was used for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of the high-risk and low-risk groups. 22 TICs in each sample were estimated using the CIBERSORT package. Mutations and TMB scores were plotted for each sample using the maftools package. The correlation between TiME-score and TMB was plotted. This study also explored the association between TiME-score and the expression of key genes of immune checkpoint blockade. Also, the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm was used to evaluate the potential response of PAAD patients in different risk score groups to immunotherapy.

Half-maximal inhibitory concentrations (IC50) were estimated by the pRRophetic package. The filtering condition was P< 0.05. Integrating the results of drug sensitivity analysis and the common clinical drugs for PAAD, the hub genes of TiME-score were validated in the immune-, targeted- and chemo-therapy mice PAAD model in this study. First of all, 6-8 weeks female mice (C57BL/6), mice pancreatic cancer PANC-02 cells, Nanoparticle Albumin-Bound Paclitaxel (Abraxane, Jiangsu Hengrui Pharmaceutical Co., Ltd., Specification: 100mg, Lot No. 220310AF), gemcitabine, recombinant human vascular endothelial Inhibin Injection (Endostar, Shandong Xiangsheng Biopharmaceutical Co., Ltd., Specification: 15mg/3ml, Lot No.: 202109053), and PD1 monoclonal antibody (Jiangsu Hengrui Pharmaceutical Co., Ltd., Specification: 200mg, Lot No.: 202009043A) were prepared. Then, the cultured PAAD cell suspension (concentration: 5×107 cells/ml) of PANC-02 mice was collected and injected subcutaneously in the right axilla of mice with 0.1 ml each. The diameter of the transplanted tumor was observed, and the mice were randomly divided into 4 groups of 6 mice each when the tumor grew up to 100 mm³. Simultaneously, 4 groups of mice were started to be administered separately according to the dosing regimen ( Supplementary Table 1 ). 15 days later, immunohistochemical verification of hub genes AGT, DEFB1, GH1, IL20RB, and TRAF3 was performed in PAAD tissues of mice in each group. Supplementary Table 2 presents information about the antibodies used for immunohistochemistry.

All data analyses were based on R software (version 4.1.3). The Wilcoxon test was used to compare the differences in gene expression levels between PAAD and normal tissues. The Spearman coefficient test was used to perform correlation tests. Image-Pro Plus (version 6.0) and GraphPad Prism software (version 8.0.2) were used for immunohistochemistry analysis. The threshold of statistical significance was considered to be P < 0.05 if not explicitly mentioned.

By integrating the PAAD data of TCGA and GTEx (Normal: n=171, Tumor: n=178), 576 differentially expressed IRGs were selected ( Figures 2A, B ). Using the same method, 85 differentially expressed immune-associated TFs were obtained ( Figures 2C, D ). A total of 152 prognosis-related IRGs were obtained by COX regression analysis of PAAD immune-related differential genes and clinical prognosis data (P < 0.05), where Hazard ratio (HR) > 1 were prognostic risk factors and 0 < HR < 1 were prognostic protective factors ( Supplementary Table 3 ). Subsequently, the interaction network between prognosis-related IRGs and TFs was constructed, and most of them were found to be positively regulated ( Figure 2E ). This indicates that the immunological relevance of the 152 IRGs screened in PAAD is reliable.

Consensus cluster analysis was performed on all PAAD samples with survival information in the TCGA-PAAD cohort. Intra-group correlations were highest and inter-group correlations were lowest when the clustering variable k = 3 ( Figure 3A , Supplementary Figure 1 ). K-M survival analysis showed that cluster 1 had a much higher survival advantage than clusters 2 and 3, while cluster 3 had a slightly higher survival advantage than cluster 2 (P<0.001, Figure 3B ). Combined with the clinical indicators, it was found that the better prognosis of cluster 1 mainly corresponded to lower grade and stage, while the worse prognosis of clusters 2 and 3 mainly matched higher grade and stage ( Figure 3C ). Then, StromalScore, ImmuneScore, and ESTIMATEScore were applied to TiME for each of the three clusters ( Figures 3D-F ). The results showed that all three scoring groups of clusters 1, 2, and 3 showed a progressive relationship. Besides, the difference in programmed death-ligand 1 (PD-L1) expression between the normal and PAAD groups was not statistically significant, but there was a significant progressive trend of PD-L1 expression levels in clusters 1, 2, and 3 ( Figures 3G, H ). It indicates that this consensus clustering can better assess immune escape and guide immunotherapy than conventional classification. Interestingly, although clusters 1,2, and 3 were progressive in immune cell infiltration analysis and immune escape prediction, the survival analysis results of clusters 2 and 3 were reversed in the K-M curves. We speculate that this may be related to insufficient sample size.

To seek a better assessment method for predicting PAAD immune-related prognosis, we constructed a signature based on 152 PAAD prognosis-related IRGs. The TCGA-PAAD cohort was used as the training set and the ICGC cohort as the validation set. By using LASSO-COX regression analysis, 14 hub genes with the best prognostic value were finally identified: AGT, CXCL9, DEFB1, ERAP2, GH1, IL1R2, IL20RB, LMBR1L, MET, PLAU, RARRES3, TNFSF10, TRAF3, and TYK2. An immune-prognostic signature, TiME-score, for PAAD was constructed ( Figures 4A, B , Table 1 ). Patients from the TCGA-PAAD cohort were divided into a high-risk group (n=88) and a low-risk group (n=89) according to the median risk cutoff. Principal component analysis of the different risk groups showed satisfactory separation ( Figure 4C ). In addition, K-M survival curves showed significant differences in survival time between patients in the high-risk and low-risk groups ( Figure 4D ). Meanwhile, the area under the curve (AUC) of the ROC curve at 1, 3, and 5 years was 0.820, 0.837, and 0.946 ( Figure 4E ). This demonstrated the promising sensitivity and specificity of the TiME-score. As a validation set, the results of the analysis of the ICGC cohort further confirmed the reliability of the TiME-score ( Figures 4F-H ).

To investigate whether the prognostic value of the TiME-score was independent of other clinical factors, we conducted univariate and multifactorial COX regression analyses on TiME-score ( Figures 5A, B ). The results showed that TiME-score remained an independent prognostic factor after adjusting for relevant clinical characteristics (P<0.001). Then, the comparison of PAAD patients with different clinical indicators revealed that TiME-score was consistent with clinical malignancy, in consensus cluster grouping (P<0.001), Grade grouping (P=0.0076), T grouping (P=0.0022), and N grouping (P=0.05) ( Figures 5C-G ). Notably, the trend in the 3 clusters of consensus clustering ( Figure 5C ) was following the results of survival analysis ( Figure 3B ), which reflects the superior predictive capability of the TiME-score. Besides, to define an individualized scoring system for each patient, a nomogram combining age, gender, clinical grade, pathological stage, and TiME-score was created to predict the 1-, 2-, and 3-year survival of PAAD patients ( Figure 5H ). Specifically, the 1-, 2-, and 3-year survival rates for a given patient can be calculated from this nomogram we constructed. The points corresponding to the patient’s TiME-score, gender, age, pathological grades, and clinical stages are summed. The resulting total points correspond to the “Total points” axis at the bottom of this nomogram. A vertical line is drawn through the point. The intersection of this vertical line with the three different survival axes, Pr(futime>1), Pr(futime>2), and Pr(futime>3), corresponds to the survival rate of the patient at the three time points of 1 year, 2 years, and 3 years, respectively. The calibration curve showed that the predicted probability values of OS at 1, 2, and 3 years were similar to the actual, demonstrating the clinical value of this nomogram ( Figure 5I ).

Combining TCGA-PAAD and GTEx datasets, all 14 hub genes of TiME-score were analyzed for differential expression levels in tumor tissues and their adjacent normal pancreatic tissue samples ( Figures 6A-D , Supplementary Figure 2 ). Subsequently, the immunohistochemical results of normal and PAAD samples for each hub gene under the same antibody treatment were selected from the HLA database and analyzed for comparison ( Figure 6E ). The results showed that there were 4 genes whose immunohistochemical results were consistent with the differential analysis of transcriptome. Specifically, their immunohistochemical results in normal tissues adjacent to the cancer were AGT (Staining: Not detected, Intensity: Negative, Quantity: None), ERAP2 (Staining: Not detected, Intensity: Negative, Quantity: None), IL1R2 (Staining: Low, Intensity: Weak, Quantity: 75%-25%), and MET (Staining: Low, Intensity: Weak, Quantity: >75%). In contrast, their immunohistochemical results in PAAD tissues were AGT (Staining: Medium, Intensity: Moderate, Quantity: 75%-25%), ERAP2 (Staining: Medium, Intensity: Moderate, Quantity: > 75%), IL1R2 (Staining: High, Intensity: Strong, Quantity: 75%-25%), MET (Staining: Medium, Intensity: Moderate, Quantity: >75%). In addition, GH1 ( Supplementary Figure 2C ), PLAU ( Supplementary Figure 2F ), and TNFSF10 ( Supplementary Figure 2H ) were not significantly different immunohistochemically, although they were statistically significant in terms of expression differences. Immunohistochemical results for TRAF3 ( Supplementary Figure 2I ) and TYK2 ( Supplementary Figure 2J ) showed opposite trends to the transcriptome intergroup differences. Immunohistochemical results for CXCL9 ( Supplementary Figure 2A ), DEFB1 ( Supplementary Figure 2B ), IL20RB ( Supplementary Figure 2D ), LMBR1L ( Supplementary Figure 2E ) and RARRES3 ( Supplementary Figure 2G ) were not retrieved and need further validation. We then analyzed the proteomics of all 14 hub genes of TiME-score in the CPTAC database using UALCAN. The results showed that the differences of protein expression between normal and PAAD cells of AGT ( Supplementary Figure 3A ), ERAP2 ( Supplementary Figure 3C ), MET ( Supplementary Figure 3D ), PLAU ( Supplementary Figure 3E ), TNFSF10 ( Supplementary Figure 3F ), and TRAF3 ( Supplementary Figure 3G ) consistent with transcriptome analysis ( Supplementary Figure 3I ). DEFB1 ( Supplementary Figure 3B ) and TYK2 ( Supplementary Figure 3H ) protein expression did not match the intergroup differences between cancer and adjacent normal tissues in the transcriptome.

To further explore the expression differences of all 14 TiME-score hub genes between PAAD tissues and their adjacent normal tissues. A comparative analysis of 3 PAAD single-cell samples and 1 para-cancer sample from GSE197177 were performed. The results showed that cancer tissues were clearly separated into 10 cell clusters ( Supplementary Figure 4A ), and these cell clusters were annotated into 6 different cell types ( Figure 7A ), which were T cells, Macrophages, Epithelial cells, Tissue stem cells, NK cells, and B cells. The normal tissue was clearly separated into 13 cell clusters ( Supplementary Figure 5A ), and these cell clusters were annotated into 5 different cell populations ( Figure 7B ), named Epithelial cells, Monocytes, Macrophages, Tissue stem cells, and Endothelial cells. The accuracy of annotation was verified one by one ( Supplementary Figures 4B , 5B ). Overall, the infiltration levels of T cells, NK cells, and B cells, which are immune infiltrating cells in cancer tissues, were significantly higher than those in normal tissue samples adjacent to cancer. This suggests that the development of PAAD is accompanied by more pronounced immune cell infiltration. Furthermore, a comparison of all 14 hub genes expression counts data in single-cells of pancreatic cancer ( Supplementary Figure 4C ) and its para-cancerous tissues ( Supplementary Figure 5C ) revealed that the expression levels of RARRES3 and TNFSF10 were most significantly increased in PAAD single-cell samples compared with normal. In contrast, the expression level of DEFB1 was most significantly decreased in PAAD single-cell samples compared with normal tissues. They were consistent with the results of transcriptome analysis. Interestingly, GH1 gene expression was not detected in any cell clusters of normal tissue. This suggested that although the overall expression level of GH1 was low in both PAAD and para-cancerous tissues, the expression level of GH1 in cancerous tissues was still significantly higher than that in para-cancerous tissues, which was different from the results of our transcriptome analysis, so we verified the reliability of GH1 in subsequent animal experiments. Besides, the expression differences of the remaining 10 genes, AGT, CXCL9, ERAP2, IL1R2, IL20RB, LMBR1L, MET, PLAU, TRAF3, and TYK2, were consistent with the trend of inter-group differences in transcriptome. Subsequently, consider the value of most significantly expression in single-cell samples of PAAD cancer tissues, we labeled the expression levels and expression distribution of RARRES3 and TNFSF10 in the clustered UMAP cell cluster descending plots of cancer and para-cancer tissues. Combined with the violin plots of expression distribution, we found that RARRES3 was mainly expressed in T cells, B cells and Tissue stem cells clusters of cancer tissues ( Figure 7C ) and Epithelial cells clusters of normal tissues ( Figure 7D ). TNFSF10 was mainly expressed in Tissue stem cells clusters of cancer tissues ( Figure 7E ) while generally expressed low in the cell clusters of paraneoplastic tissue ( Figure 7F ). This suggests that the increased expression of RARRES3 is consistent with the trend of increased immune cell infiltration in tumor tissues.

RARRES3 (alternate name PLAAT4) and TNFSF10 are hub genes of TiME-score. As a co-upstream regulatory gene of RARRES3 and TNFSF10, TP53 mutation has attracted our attention for its relevance to TiME. Although the correlation between TP53 mutations and OS in PAAD patients has been well substantiated (13–15), the effect of TP53 mutation on PAAD-specific TiME has not been thoroughly investigated. We performed GSEA analysis of PAAD samples from the TP53 wild group (n=62) and mutation group (n=81) using gene expression and clinical information from TCGA. Five immune-related biological processes were selected based on the results: innate immune response activating cell surface receptor signaling pathway (NES=1.805, FDR=0.034, SIZE=115), somatic diversification of immune receptors via somatic mutation (NES=1.741, FDR=0.049, SIZE=18), interleukin 1 mediated signaling pathway (NES=1.983, FDR =0.029, SIZE=99), antigen processing and presentation of peptide antigen via MHC class I (NES=1.891, FDR=0.028, SIZE=97), and response to interleukin 7 (NES=1.876, FDR=0.025, SIZE=39) ( Figures 8A-E ). The results showed that the TP53 mutation was significantly enriched in immune-related biological processes. However, the TP53 wild group had no gene sets significantly enriched at FDR < 25%. Moreover, the mutation data analysis of the TCGA-PAAD cohort showed that missense mutations, SNPs, and C>T accounted for a large proportion of TP53 mutations, which with a high mutation rate of 63% ( Figure 8F ). Further analysis of the differences in the TiME landscape between the TP53 wild and mutant groups revealed the presence of more pronounced immunosuppression in the mutant group compared to the wild group ( Figure 8G ). The proportion of tumor-infiltrating immune cells in different subpopulations was weakly to moderately correlated ( Figure 8H ).

In this study, the prognostic value of RARRES3 was further explored. Considering that RARRES3 and TNFSF10 are both downstream genes of TP53, and TNFSF10 is a validated prognostic risk factor for PAAD, we used TNFSF10 as a comparator to evaluate the prognostic value of RARRES3. Survival analysis showed that survival time was significantly lower in the RARRES3 high-expression group than in the low-expression group in the TCGA-PAAD cohort (P=0.016) ( Figure 9A ). The ROC curves also demonstrated good sensitivity in predicting the prognosis (AUC at 1 year: 0.722, AUC at 3 years: 0.744, AUC at 5 years: 0.716) ( Figure 9B ). Furthermore, independent prognostic analysis of univariate and multifactorial showed that RARRES3 had significantly better independent prognostic value than other relevant clinical factors after integration (P<0.001) ( Figures 9C, D ). Using the same method to analyze TNFSF10, TNFSF10 exhibited a similar positive prognostic value to RARRES3 ( Figures 9E-H ). This also indirectly verified the prognostic ability of RARRES3. The clinical information of the high and low RARRES3 expression groups was additionally analyzed as a whole, and we found a significant correlation between RARRES3 and clinical grading ( Figure 9I ). Samples with high RARRES3 expression had a higher level of clinical staging than those with low expression, with a higher proportion of Grade 3 and 4. This suggests that high RARRES3 expression corresponds to worse clinical staging.

GSEA analysis was applied to RARRES3 and TNFSF10 separately, taking the top three ranked biological enrichment processes ( Figures 10A, C ). It was found that the immune-related component IMMUNOGLOBULIN COMPLEX was significantly enriched in the high expression group of both genes. Next, immune cell infiltration was analyzed for RARRES3 and TNFSF10 ( Figures 10B, D ). In particular, the high infiltration characteristics of Macrophages M1 in the high expression group attracted our additional attention. Further analysis of the correlation between all immune cells and RARRES3 expression was performed ( Figure 10E ). The results showed that Macrophages M1 was positively correlated with RARRES3 expression (R=0.2, P=0.0077) ( Figure 10F ). Additionally, RARRES3 was positively correlated with all immune checkpoint-related genes ( Figure 10G ). Then, TNFSF10 was analyzed in the same way, and the positive correlation of TNFSF10 expression levels with Macrophages M1 was significantly similar to RARRES3 ( Figures 10H-J ). The positive correlation of RARRES3 with pro-inflammatory Macrophages M1 and the negative correlation with anti-inflammatory Macrophages M2 (R=-0.004) also led to our speculation of the pro-inflammatory TiME in the RARRES3 high expression group.

The enrichment analysis of GO and KEGG showed that the hub genes of TiME-score were mainly enriched in the biological processes of immune activation, immune defense, and immune system-related diseases ( Figures 11A, B , Supplementary Table 4 ). A further assessment of the infiltration of immune cells in the high- and low-risk groups was performed ( Figure 11C ). B-cell naive phase, T-cell CD8+, T-cell CD4+ memory phase, monocytes, and macrophage M0 were found to be negatively correlated with TiME-score; while macrophage M1, activated NK cells were positively correlated with TiME-score. We also calculated TMB scores for each mutant sample in the TCGA-PAAD cohort ( Figures 11D, E ). The results showed that the high-risk group had a broader TMB distribution, and the proportion of TP53 mutations was significantly higher in the high-risk group (71%) than in the low-risk group (32%). Furthermore, TMB quantitative analysis confirmed statistically significant differences in TMB distribution between high and low risk groups (Wilcoxon test, p=8e-05; Figure 11F ). A positive correlation between TiME-score and TMB was also confirmed (R=0.33, p=7.1e-05; Figure 11G ). To further clarify the relationship between TiME-score and immunotherapy response, we also explored the correlation between immune check sites and TiME-score ( Figure 11H ). Correlation analysis showed that CD274, FAP, FEN1, LOXL2, MCM6, MSH2, MSH6, POLD3, POLE2, and TAGLN were all positively correlated with TiME-score. In addition, the TIDE scores of patients in the high-risk and low-risk groups indicated that patients with low-risk were more likely to experience immune escape and patients with high-risk were more likely to benefit from immune checkpoint inhibitor (ICI) therapy ( Figure 11I ).

Clinical treatment of PAAD with single drug is much less effective than combination therapy (16). Nano-albumin Paclitaxel in combination with Gemcitabine is the first chemotherapy choice for non-surgical PAAD patients (17). Currently, immunotherapy, represented by PD-1 monoclonal antibody, is not as effective in combination with chemotherapy alone or with other molecules (18). In contrast, some clinical therapeutic efficacy has been achieved by combining recombinant human vascular endothelial inhibitor with chemotherapy (19). Therefore, we further explored the intergroup differences in IC50 of these drugs. We selected GH1, for which immunohistochemical differences were not significant, TRAF3, for which immunohistochemical trend were opposite to transcriptomic and proteomic trends, DEFB1, for which proteomic and transcriptomic trends did not match, IL20RB, for which neither immunohistochemical nor proteomic was retrieved, and the TiME-score itself for separate analysis. AGTs with consistent transcriptomic, immunohistochemical, and proteomic trends were included in the analysis as reference. Among them, Paclitaxel had a lower IC50 in the AGT low expression group ( Figure 12A ), indicating that AGT low expressers were more sensitive to Paclitaxel. The IC50 was lower in the IL20RB ( Figure 12B ) and TRAF3 ( Figure 12C ) high expression groups, indicating that IL20RB and TRAF3 high expressers were more sensitive to Paclitaxel. For the TiME-score ( Figures 12D, E ), the IC50 of Paclitaxel was lower in the high-risk group (R = -0.55, p = 2.3e-15), indicating that those with higher scores were more sensitive to Paclitaxel. Unfortunately, there were no statistically significant group differences in GH1 and DEFB1 to Paclitaxel. There were no significant between-group differences for Gemcitabine. Recombinant human vascular endothelial inhibitors were not included in the pRRophetic package used for prediction. However, considering the value of their combination therapy mentioned above, we still included the combination regimen of Nano-albumin Paclitaxel and Gemcitabine as the chemotherapy group and included recombinant human vascular endothelial inhibitor in the targeted combination chemotherapy group. Combining the above research advances, 4 different regimen treatment cohorts of model control group, chemotherapy group (Nano-albumin Paclitaxel + Gemcitabine), recombinant human vascular endothelial inhibitor + PD-1 monoclonal antibody group and recombinant human vascular endothelial inhibitor + PD-1 monoclonal antibody + chemotherapy group were included in this study. Volume recordings of the tumors for 15 consecutive days ( Figure 12F ) showed that the tumor growth volume and rate in the recombinant human vascular endothelial inhibitor + PD-1 monoclonal antibody + chemotherapy group were significantly smaller than those in the remaining 3 groups. This suggests the value of combining immunotherapy with targeted therapy and chemotherapy for tumor suppression. Next, we performed immunohistochemical validation of GH1, TRAF3, DEFB1, and IL20RB, which are controversial above, in mice PAAD model ( Figures 12G, H ). AGT, for which transcriptomic, immunohistochemical, and proteomic trends were consistent, was also included as a reference for the feasibility and accuracy of this validation approach. The results showed that AGT with a negative coefficient (Coef. = -0.0661) was expressed at a lower level in the Nano-albumin Paclitaxel group than non-Nano-albumin Paclitaxel group. DEFB1 (Coef. = 0.0159) and IL20RB (Coef. = 0.0479) with positive coefficients had the higher expression in the Nano-albumin Paclitaxel involved chemotherapy group. This is consistent with the lower IC50 of Paclitaxel in the AGT low expression group and IL20RB high expression group. Considering the positive correlation coefficient of the hub gene in TiME-score, the expression of the gene was positively correlated with the score and those with high score were more sensitive to Nano-albumin Paclitaxel. The experiments for DEFB1 were also consistent with the predicted results. The intergroup differences for GH1 (Coef.= -3.2585) and TRAF3 (Coef.= -0.3324), although the P value was not smaller than 0.05, expression differences were consistent with the predicted trend. Moreover, considering immune cell infiltration was more pronounced in the low-risk group, DEFB1 and IL20RB with positive coefficients were significantly less expressed in the immunotherapy cohort than in the non-immunotherapy cohort. AGT with a negative coefficient had higher expression in the immunotherapy cohort. In other words, the immunotherapy cohort corresponded to a lower TiME-score, and the low-risk group had better immunotherapy outcomes and a better prognosis. Collectively, TiME-score is potentially valuable as an immune-related prognostic model for the prediction of individualized PAAD anticancer drug sensitivity.