PDAC level 3 expression profiles of RNA sequencing data [fragments per kilobase of transcript per million mapped reads (FPKM)], somatic mutation, copy number variation (CNV) data, and the corresponding clinical information of 143 PDAC patients were downloaded from TCGA (https://tcga-data.nci.nih.gov/tcga/). The patients who lack survival time or who had survival time of less than 30 days were removed, and only 136 PDACs were finally enrolled. Validation cohort, i.e., GSE71729, GSE102238, and GSE57495, was extracted from GEO datasets (https://www.ncbi.nlm.nih.gov/gds/). The gene expression information of normal pancreatic tissue (FPKM) was derived from Genotype-Tissue Expression (GTEx) via University of California, Santa Cruz, Xena (https://xenabrowser.net/datapages/?cohort=GTEX). All RNA sequencing data were normalized and transformed into log2(FPKM + 1); then unexpressed or extremely low-expressed genes in most of the samples (average log2(FPKM + 1) < 0.01) were filtered out. Matched clinical validation cohort, namely, Peking2020 Cohort, with 40 cases in total, originated from Cancer Center/National Clinical Research Center for Cancer/Cancer Hospital, wherein written informed consent was obtained from all patients ( Table 1 ).

The xCell online tool (https://xcell.ucsf.edu/) was used to perform enrichment analysis of tumor-infiltrating immune cell (TIIC) level, Immune Score, Stroma Score, and Microenvironment score based on gene expression data. The data analysis results of TCGA database were downloaded directly from the xCell portal website, while the GEO data needed to use the “sva” package in the R software 4.0.0 to remove the batch differences, and then we uploaded the standardized data to the xCell website for analysis. According to the median of the level of TIIC, 136 cases of PDAC patients in TCGA were divided into the high and low TIIC groups. Survival analysis in R software 4.0.0 was used to determine the cell types with prognostic value, and then the GEO datasets and Peking2020 Cohort were used for verification. Immunotherapy sensitivity was predicted through the ImmuneCellAI online perform (http://bioinfo.life.hust.edu.cn/ImmuCellAI#!/) based on gene expression data.

The risk-score of each patient was calculated based on the survival-related level of TIIC multiplied by the multivariate Cox regression coefficient (see the following formula for details), and the patients were divided into the high risk-score group and low risk-score group based on the median value of risk-score. The Kaplan–Meier (KM) survival analysis was performed on the two groups using log-rank test; receiver operating characteristic (ROC) curves of 1-year, 2-year, and 3-year survival prediction were drawn; then the area under the ROC curve (AUC) value was calculated.

where βi was the coefficient of the ith the survival-related TIIC, i.e., β(CD4⁺T_N) = −0.16, β(CLP) = 0.63. TIICi represents the level of the ith survival-related TIIC.

Differentially expressed genes (DEGs) of the high risk-score group and low risk-score group were analyzed by R software 4.0.0 “edgeR” package. The screening threshold was set to log2(fold change) >1 and false discovery rate (FDR) <0.05, and then we used “pheatmap” package and “plot” package in R software 4.0.0 to draw heatmaps and volcano maps, respectively.

“Clusterprofiler” package in the R software 4.0.0 was used to perform functional enrichment analysis on DEGs, including Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis, where GO analysis included molecular function (MF), biological process (BP), and cellular component (CC) analyses.

Through Search Tool for the Retrieval of Interacting proteins (STRING, https://string-db.org) database, the protein–protein interaction (PPI) network of DEGs, namely, PPI, was established. The comprehensive score>0.4 was used as the criterion for the existence of interaction. Then Cytoscape software was used to rebuild the PPI network, and the Cytoscape MCODE plug-in was used to find core clusters that located densely connected areas and calculated the connectivity of network nodes.

Tumor mutation burden (TMB) was defined as the total number of somatic mutations (including single-nucleotide variants (SNVs), missense mutation, and insertion–deletion (INDELs) per million bases in the coding region of exons). Somatic mutation data were analyzed using VarScan2, and then TMB was calculated as the total number of somatic mutations/all bases, and the unit was mutations/Mb. In this study, the “Maftools” package in R software 4.0.0 was used to calculate the TMB of each sample.

The copy number information of PDAC from TCGA was annotated according to the position information of grch38 genome; set the normal copy change as 0, then single copy amplification as 1, double copy even multiple copy amplification as 2, single copy deletion as −1, and double copy or multiple copy deletion as −2. Then, the chi-square test was used to select the differential copy number between the two groups in the high-risk group and low-risk groups (p < 0.05); the Kruskal test was used to determine the copy number differential genes related to gene expression; and the Gene Set Enrichment Analysis (GSEA) was used for functional analysis.

Multiplex immunofluorescence staining was performed according to a sequential multiplexed immunofluorescence protocol (17). Briefly, co-stain CD3 (Cat #ab16669, RIDD: AB_443425, Abcam), CD4 (Cat #ab133616, RIDD: AB_2750883, Abcam) and CD45RA (Cat #ab755, RIDD: AB_305970, Abcam), or CD127 (Cat #DF6362, RIDD: AB_2838326, Affinity), CD135 (Cat #DF8546, RIDD: AB_2841750, Affinity), and CCL19 (Cat #13397-1-AP, RIDD: AB_2071529, ProteinTech). Then, corresponding secondary antibodies used were fluorescein isothiocyanate–Tyramide Signal Amplification (FITC-TSA) (PPD520, Panovue) for CD3 or CCL19 detection, CY3-TSA (PPD570, Panovue) for CD4 or CD127 detection, and CY5-TSA (PPD620, Panovue) for CD45RA or CD135 detection. Nuclei were highlighted using DAPI (D9542, Sigma-Aldrich). Finally, using Phenochart software (Version 1.08, PerkinElmer) and inForm image analysis software (Version 2.4, PerkinElmer), we estimated the number of co-located cells of CD3, CD4, and CD45RA per ×100 magnification field as CD4⁺T_N level (randomly select three fields at least). Similarly, the number of co-located cells of CD127 and CD135 was calculated as CLP level. For the scoring of CCL19, we first evaluated the staining intensity of whole tumor tissue at low magnification. Samples with no staining were scored 0, weakly stained samples scored 1, mildly stained samples scored 2, and strongly stained samples scored 3. We also calculated the number of positive cells from at least three high magnification fields chosen at random as well as their mean intensities. As described above, samples with <25% positive expression were scored 1, samples within the expression range of 25%–50% scored 2, samples within the expression range of 50%–75% scored 3, and samples with expression ≥75% scored 4. The final CCL19 expression was determined by multiplying the intensity score with the positive expression value.

Statistical analysis was performed using GraphPad Prism 8.0 and R software 4.0.0. Measurement data were expressed as mean ± standard deviation, using Student’s t-test. Counting data were expressed as percentage (%), using chi-square test. Survival curve was drawn by the KM method, using log-rank test. The correlation between immune cells was evaluated by Spearman’s correlation coefficient. p < 0.05 or FDR < 0.05 was considered statistically significant in all tests.

Based on the xCell algorithm, the type and level of immune cell infiltration in 136 PDAC cases were evaluated, and a total of 36 immune cells were predicted. Then KM survival analysis showed that only CD4⁺ naïve T cell (CD4⁺T_N), CLP, and CD4⁺Th2 cell (Th2) were closely related to survival, where CD4⁺T_N was positively correlated with survival (p < 0.001, Figure 1A ), while CLP and CD4⁺Th2 cell were negatively correlated with survival (p = 0.009, Figure 1B ; p = 0.022, Supplementary Figure S1A ). Furthermore, combined with clinical factors for univariate and multivariate analyses, only CD4⁺T_N and CLP could be regarded as independent risk factors, but the p-value was only slightly significant ( Figures 1E, F and Table 2 ). However, when combining CD4⁺T_N and CLP to construct a proportional hazards model (Cox model), its clinical significance could be dramatically improved (p < 0.001, Table 2 ). Based on the individual point of risk-score calculated through the Cox model, PDAC patients were divided into the low risk-score and high risk-score groups (median set as cutoff). KM curves confirmed that the high risk-score group had much worse prognosis than the low risk-score group (p < 0.001, Figure 1C ). Additionally, the risk-score could estimate OS of patients much sensitively, as the AUC was 0.713, 0.667, and 0.609 for 1-year, 2-year, and 3-year OS, respectively, and the discrimination index (C-index) was 0.706 ( Figure 1D ). Furthermore, this result was also verified in GSE71729 dataset (p = 0.027, Supplementary Figures S1B, C ).

In succession, through the scatter plot of patient risk assessment, we found that risk-score was positively correlated with CLP level and negatively correlated with CD4⁺T_N level ( Figure 2A ). Meanwhile, with the increase of risk-score, the survival time of patients was gradually shortened, and the death events increased. Consistent with that, the risk-score was significantly higher in the short-term survival than the long-term survival ( Figure 2B , p = 0.01). Interestingly, this risk-score also happened to be negatively correlated with the Immune Score, Stroma Score, and Microenvironment score. Previous studies had shown that molecular classification of pancreatic cancer was closely related to prognosis, but its complex algorithm has also become one of the important factors limiting its clinical application. Therefore, we were trying to explore whether the risk-score of our model could well identify high-risk molecular subtypes and replace its prognostic evaluation value. According to Moffitt’s algorithm (13) for molecular classification of PDAC, we found that with the increase of the risk-score, the proportion of basal-like and activated stroma corresponding to the molecular subtype with the worst prognosis increased. Also, through ROC curve analysis, it was found that Cox model was also highly sensitive to the prediction of basal-like molecular subtype (AUC = 0.728, Figure 2C ). Furthermore, it was more sensitive to predict the joint basal-like and activated stroma molecular subtype (AUC = 0.733, Figure 2C ). In addition, we also performed molecular classification of PDAC according to Puleo’s algorithm (14), and the results were consistent with those of the former: with the increase of risk-score, the proportion of pure basal-like and stroma-activated molecular subtypes associated with the worst prognosis increased as well as high prediction sensitivity (AUC = 0.696, Figure 2C ).

Further studies suggested that risk-score could well predict the cytotoxicity of T cells in the microenvironment of PDAC and the sensitivity of PDAC patients to immune checkpoint inhibitor. The risk-score of the high cytotoxicity group and immunotherapy sensitive group was lower (p < 0.001 and p < 0.001, respectively; Figures 2D, E ). Moreover, PDAC patients in whom molecular targeted therapy is applicable had a lower risk-score value (p = 0.004, Figure 2F ). Besides that, the risk-score was positively correlated with the degree of dysplasia and tumor size; that is, the higher the degree of dysplasia, or the larger the tumor size, the higher the risk-score (p = 0.005, p = 0.045, Figures 2G, H ). Meanwhile, patients with higher risk-score usually survived with tumors, and surgery to remove the tumor totally was more difficult (p = 0.019, Figure 2I ).

In order to further verify that the Cox model we constructed could indeed assess the prognosis of PDAC patients, we selected 40 PDAC patient samples, Peking2020 cohort, and followed the above algorithm to analyze the level of intratumoral infiltration of CD4⁺T_N (CD3⁺, CD4⁺, and CD45RA⁺) and CLP (CD127⁺ and CD135⁺) to calculate the risk-score of each case by multiplex immunofluorescence staining. As expected, the results of Peking2020 cohort were consistent with the above results. Contrary to the PDAC with high risk-score, the level of CD4⁺T_N was higher in the low risk-score group ( Figure 3A ), while the level of CLP tended to be lower ( Figure 3B ). Also, compared with the short-term survival group, the level of CD4⁺T_N in the long-term survival group was also higher ( Figure 3C ), while the level of CLP and risk-score was in reverse ( Figures 3D, E ). Furthermore, we conducted KM survival analysis and found that CD4⁺T_N was positively correlated with PDAC patient survival (p < 0.001, Figure 3F ), but CLP was negatively correlated with patient survival (p < 0.001, Figure 3G ), and after the joint analysis, the significance of the negative correlation between risk-score and prognosis was further confirmed (p < 0.001, Figure 3H ). In addition, the same as the results of PDAC dataset from TCGA, the larger the tumor volume, the higher the risk-score (p < 0.001, Supplementary Figure S1D ). Interestingly, PDAC in the tail of the pancreas had a higher risk-score than the head of the pancreas (p = 0.0316, Supplementary Figure S1E ).

The above studies confirmed the role of risk-score in the prognostic evaluation of PDAC patients, so we tried to explore its potential molecular mechanism. First of all, principal component analysis (PCA) was conducted on the transcriptome data of healthy pancreas (combined with pancreas tissue from GTEx and paracancerous tissue of PDAC dataset from TCGA) and tumor tissue (PDAC dataset from TCGA). The results showed that the mRNA expression was significantly different between normal tissues and cancer tissues, while there was no significant difference between the low risk-score and high risk-score groups ( Supplementary Figures S1F, G ). However, through the analysis of heatmaps and volcano maps, we found that there were still 172 DEGs between the low risk-score and high risk-score groups, among which 162 genes were downregulated and 10 genes were upregulated in the high risk-score group ( Figures 4A, B ). Then, GO and KEGG analysis of these 172 DEGs revealed that they were mainly concentrated in immune-related functions and pathways, for example, Complement activation by classical pathway, Immunoglobulin-mediated immune response, and Cytokine–cytokine receptor interaction ( Figures 4C, D ). Furthermore, the 172 DEGs were analyzed for PPI. Then, the most important core interaction protein network was screened out through MCODE algorithm. It was found that CCL19 occupied one of the most core positions and interacted closely with other genes ( Figure 4E ). Moreover, the function of core interaction protein was mainly related to regulation of lymphocyte activation and differentiation, chemokine-mediated signaling pathway, etc. ( Figure 4F ).

Since previous studies had shown that TMB could affect tumor immune infiltration (18, 19), we further analyzed the TMB in PDAC. It was found that PDAC was a “cold” tumor with low TMB as reported in the literature (2, 3), but compared with the low risk-score, the high risk-score group still has higher TMB, such as conventional TP53 and KRAS mutation ( Figures 5A, B ), However, these mutations did not make the corresponding gene expression different at the transcriptome level, so it did not affect the expression of DEGs obtained by transcriptome analysis.

Compared with TMB, CNV between patients was more constant, and the frequency of variation was higher, which was one of the key events leading to tumor development (20–22). However, there were fewer studies on the tumor immune microenvironment. Therefore, we tried to find the main CNV that could affect the value of risk-score. By analyzing the whole-exome data of CNV, it was found that compared with the low risk-score group, the high risk-score group had more characteristic chromosome amplification and deletion ( Figures 5C, D ), such as 8q, 9p, 17, 18, 19q, and 20q chromosome amplification, and 2q, 6, 7, and 8p chromosome deletion. At the same time, a total of 1,593 differential CNVs related to expression were generated, mainly enriched in immune-related functions, such as activation and positive regulation of immune response, participating in innate immune response and antigen receptor-mediated signaling pathway ( Figure 5E ) by GSEA, where there was a total of six intersection genes with DEGs generated by the transcriptome, that is, CA9, TNNT1, FABP4, CCL21, LY86, and CCL19 ( Figure 5F ).

Compared with those in the normal pancreas (GTEx + TCGA), CA9, TNNT1, CCL21, LY86, and CCL19 were all expressed higher in tumor tissues, except for FABP4, which was expressed lower in a tumor (p < 0.001, Figures 6A – F ). In contrast with that in the low risk-score group, the expression of CA9 and TNNT1 increased in the high risk-score group, while FABP4, CCL21, LY86, and CCL19 were contrasting (p < 0.001, Figures 6G – L ). Except for FABP4, their expression was positively correlated with the change of copy number. Unlike in the low risk-score group, in the high risk-score group, CA9 showed a variety of copy number changes; that is, in addition to the normal double copy number, there was also single copy number deletion, single copy number, or even multiple copy number amplification (p = 0.023, Figure 6M ). TNNT1 and FABP4 mainly showed single copy number amplification (p = 0.002, p = 0.005, Figures 6N, O ), while CCL21, LY86, and CCL19 mainly showed single copy deletion in the high risk-score group (p < 0.001, p < 0.001, and p = 0.001, Figures 6P – R ). Furthermore, we analyzed the correlation between these six genes and risk-score, and we found that except for CA9 ( Figure 6S ), the other five genes were significantly correlated with risk-score. Among them, only the expression of TNNT1 was positively correlated with risk-score (R = 0.492, p < 0.001, Figure 6T ), FABP4, CCL21, LY86, and CCL19 were negatively correlated with risk-score (R = −0.51, R = −0.6, R = −0.616, R = −0.626, p < 0.001, Figures 6U – X ), while there was no statistical correlation between CA and risk-score.

Furthermore, we analyzed the survival of these five genes and found that only TNNT1, CCL21, and CCL19 were closely related to survival, of which TNNT1 was negatively related to survival (p = 0.026, Figure 7A ), while CCL21 (p = 0.038, Figure 7B ) and CCL19 (p = 0.019, Figure 7C ) were positively correlated with prognosis. Then univariate and multivariate analyses combined with clinical factors found that only CCL19 was positively correlated with patient survival (p = 0.003, Figure 7D ) and could play a role as an independent prognostic factor (p = 0.005, Figure 7E ). It was also verified in GSE57495 dataset concurrently (p = 0.034, Figure 7F ). Similarly, in Peking2020 cohort, we also found that CCL19 expression was significantly lower in the high risk-score group ( Figure 7G ). Compared with the long-term-survival patients, CCL19 also showed a downward trend in the short-term-survival patients group (p = 0.008, Figure 7H ). Further survival analysis showed that CCL19 was positively correlated with the prognosis (p < 0.001, Figure 7I ), consistent with the above result, which suggested that CCL19 may be an important factor affecting CD4⁺T_N and CLP infiltration.