A single-cell sequencing dataset (GSE155698) including 16 PDAC and 3 adjacent normal samples was obtained from the Gene Expression Omnibus (GEO). The bulk sequencing datasets were derived from The Cancer Genome Atlas (TCGA) (TCGA-PAAD), International Cancer Genome Consortium (ICGC) (ICGC-AU), GEO (GSE21501, GSE57495, GSE71729, and GSE85916), and ArrayExpress (E-MTAB-6134) databases. Both the transcriptome information and clinical information of each dataset were concurrently downloaded from the respective websites. The transcriptome data were transformed to the format of Log₂[transcripts per million (TPM) + 1]. Only PDAC tissues were included in the subsequent analysis, while other histological subtypes, such as neuroendocrine tumors, acinar cell carcinoma, and intraductal papillary mucinous neoplasms, were excluded. T-exhaust and immune checkpoint blockade (ICB) resistance signatures were downloaded from the Tumor Immune Dysfunction and Exclusion (TIDE) database.

The hazard ratio (HR) for each infiltrated immune cell against the OS of PDAC was computed with the log-rank test. The HRs of each immune cell in different bulk sequencing-based cohorts were pooled in a fixed-effects model if no robust heterogeneity was observed (I² < 50% and P > 0.05). The meta-analysis was performed using Stata 15.1, and the forest plot was depicted via GraphPad Prism 7.0.

The “Seurat” package was used to perform the single-cell sequencing analysis. The batch effect of studies was removed through regularized negative binomial regression by the “Seurat” package (24). Genes detected in <3 cells were excluded, and cells with <200 total detected genes were excluded. Afterwards, we calculated the standardized variance of each gene across different cells, and only the top 2,000 variable genes were selected for subsequent analysis. Principal component analysis (PCA) was performed to identify significant dimensions with P < 0.05 (25). Then, the t-distributed stochastic neighbor embedding (tSNE) algorithm was applied for dimensionality reduction with the 20 initial PCs and for performing cluster classification analysis across all cells (26). Non-linear dimensional reduction was also performed with the UMAP method. Then, different cell clusters were determined and annotated by the “singleR” package according to the composition patterns of the marker genes and were then manually verified and corrected with the CellMarker database (27, 28). Given that both “singleR” and CellMarker could only classify cell clusters into basic types, we also referred to previously published scRNA-seq analyses to further classify each cell cluster into more precise subtypes (14, 15, 29–31).

The enrichment scores of the hallmark genes were evaluated using ssGSEA with the R package “GSVA” (32). Hallmark genes were defined as the top 50 genes with the largest fold change (FC) in each cluster. Using ssGSEA, each sample with complete bulk sequencing data and clinical information was labeled with an enrichment score of the specific hallmark gene signature. The samples were divided into two groups based on the median enrichment score. Then, a Kaplan-Meier curve was plotted to visualize the survival difference between the two groups. The log-rank test and Gehan-Breslow-Wilcoxon test were performed to verify the statistical significance of the survival difference. A P < 0.05 was regarded as indicative of a significant difference.

We utilized the transcriptome profile of the E-MTAB-6134 cohort, which was the largest PDAC cohort with transcriptome and clinical data, to qualify and construct the co-expression network by the “WGCNA” package in R (33). Next, Pearson's correlation matrices were constructed for pairwise genes. We constructed a weighted adjacency matrix using a power function: a_mn = |c_mn|^β (c_mn = Pearson's correlation between gene m and gene n; a_mn = adjacency between gene m and gene n). The β value emphasizes strong correlations between genes and penalizes weak correlations. After choosing the appropriate β value, the adjacency matrix was transformed into a topological overlap matrix (TOM), which measures the network connectivity of a gene defined as the sum of its adjacency with all other genes for network construction. To divide genes with similar expression patterns into gene modules, average linkage hierarchical clustering was performed according to the TOM-based dissimilarity measure with a minimum size of 50 for the gene dendrogram. Then, we further calculated the dissimilarity of module eigengenes (MEs), chose a cut line for the module dendrogram and merged some modules.

MEs were regarded as the major component in PCA for each gene module. We calculated the correlation between MEs and clinical traits or hallmark gene signatures to identify the relevant modules. Gene significance (GS) was defined as the log10 transformation of the P value in the linear regression between the gene expression level and clinical data. In addition, module significance (MS) was defined as the average GS for all the genes in a module. When the modules of interest were established, the core genes in a module were identified by GS > 0.2 and MS > 0.8. Specifically, in the present study, we also compared the expression levels of the selected core genes between tumor and normal tissues using GEPIA 2.0 (34), which incorporates the transcriptome data of normal pancreas tissue and thus facilitates the identification of differentially expressed genes.

The single-cell pseudotime trajectories of the scRNA-seq data were constructed using the Monocle 2 algorithm (35). This algorithm uses a machine learning technique, learning a parsimonious principal graph to reduce the given high-dimensional expression profiles to a low-dimensional space. Single cells were projected to this space and ordered into a trajectory with branch points. For data interpretation, the cells that were located in the same branch were thought to be in the same differentiation state, while cells located in different branches were thought to have different cell differentiation characteristics. Differentially expressed genes between different differentiation states were identified by the R package “limma” (36). Unsupervised consensus clustering based on the differentially expressed genes was conducted using the “ConsensusClusterPlus” package. The clustering procedure included 1,000 iterations, and 80% of the data were sampled in each iteration. The optimal number of clusters was determined by the relative change in the area under the cumulative distribution function (CDF) curves of the consensus score.

CAFs were first separated and purified from human pancreatic cancer tissues in our laboratory based on the study by Walter et al. and then subjected to immortalization treatment (37). Fresh pancreatic cancer tissue was minced into 1–3 mm³ fragments and digested with 0.25% trypsin at 37°C for 30 min. The resulting fragments were centrifuged at 600 × g for 5 min and washed once with Dulbecco's modified Eagle's medium (DMEM) containing 10% fetal bovine serum (FBS). The tissue fragments were then plated and allowed to adhere. After incubation at 37°C for several days, fibroblast outgrowth from the tissue fragments occurred. The fibroblasts were sub-cultured by trypsinization for 2–3 passages until free of epithelial cell contamination and maintained in DMEM supplemented with 10% FBS, 2% penicillin, and streptomycin (Invitrogen). The cells were grown at 37°C in a humidified atmosphere containing 5% CO₂. CAFs and the pancreatic cancer cell line SW1990 were cultured in DMEM supplemented with 10% FBS. The non-cell supernatant from SW1990 cultures was extracted through centrifugation at 800 rpm/min. Then, CAFs were treated with SW1990-derived supernatant during medium changing. TGF-beta (3 ng/ml) was added to CAF cultures as a positive control.

In past decades, many approaches have been developed to estimate or quantify the infiltration level of immune cells in tumor tissues, such as immunohistochemical staining and transcriptome-based estimation (16–18, 38, 39). We applied six algorithms, TIMER, CIBERSORT, XCELL, EPIC, QUANTISEQ, and MCP-counter, to estimate the infiltration percentage of immune cells (Supplementary Table 1). Then, we integrated the survival data of patients in seven PDAC cohorts with the infiltration level of each immune cell and performed univariate Cox regression to screen prognosis-related components (Supplementary Table 1). Given the difficulty of pooling the Cox regression-derived HR values due to the large confidence interval (CI) in this case, we performed a meta-analysis using the log-rank test-derived HR values. Six basic immune cell types (CD8+ T cells, CD4+ T cells, B cells, macrophages, neutrophils, and dendritic cells) estimated by the TIMER algorithm were selected for meta-analysis. Null or only minor heterogeneity was detected in the fixed effect model; hence, we generated the results of the meta-analysis with low bias. The results showed that neutrophil infiltration was associated with worse OS for PDAC, while CD8+ T cell, CD4+ T cell, and B cell infiltration was related to prolonged OS for PDAC with marginal statistical significance (Figure 1).

Interestingly, many theoretical anticancer cells, such as CD8+ T cells, did not show obvious survival relevance in the meta-analysis. We assumed that at least some of these anticancer immune cells were exhausted or experienced differentiation into protumoral cells, which deprived them of their tumor-killing capability. To investigate the alteration of stromal cells in PDAC, we reanalyzed the scRNA-seq dataset and annotated cell types according to their perturbation in the transcriptome. After carrying out the quality control procedures described in the Methods section (Supplementary Figures 1A,B), the top 2,000 variable genes were used in cell clustering (Figure 2A; Supplementary Figures 1C,D), which identified 26 clusters (Supplementary Figures 1E,F) and then classified these clusters into 17 cell types (Supplementary Figure 1G). The proportions of each cell type are shown in Figure 2B, and their clustering distribution is presented in Figures 2C–E; Supplementary Figure 1H. For some clusters, we took a conservative approach and annotated them as basic cell categories, such as T cells and myeloid cells. However, when specific cell markers were highly overexpressed in a cluster, we tended to annotate them into more precise clusters, such as NKT and CD8+ T cells. The significant gene markers are presented in Supplementary Table 2.

NKT and CD8+ T cells are theoretically capable of killing tumor cells through cytotoxic effects; however, most of these cells seem to lose their anticancer ability in the real tumor microenvironment. Through scRNA-seq, we verified that these T cells with a cytotoxic signature universally expressed exhaustion-related or ICB resistance markers (Figures 3A,B), which not only produced immune evasion among tumor cells but also caused a low response to immunotherapy, such as ICB. Then, we tested whether such signatures of exhausted cytotoxic T cells were associated with patient prognoses. As expected, no obvious association was observed between the signature of CD8+ T cells and OS in most cohorts except for GSE57495 (Supplementary Figure 2; Figure 3C). Then, we explored why the signatures of exhausted CD8+ T cells could still predict prolonged OS in some patients. Several recent studies suggested a novel mechanism by which cytotoxic cells trigger tumor cell death, as we reviewed previously (10). These studies demonstrated that cytotoxic T cells, including CD8+ T cells and NKT cells, could kill tumor cells through ferroptosis and pyroptosis induction. Then, we validated this hypothesis in GSE57495, which was the only dataset that showed a correlation between exhausted cytotoxic T cell signatures and patient prognoses. Survival analysis showed that the signature of CD8+ T cells was positively associated with prolonged OS only in samples harboring overexpression of targets for ferroptosis (SLC7A11) and pyroptosis (GSDMB, GSDMC, and GSDME) induction (Figure 3C). We continued to investigate the association between the NKT cell signature and patient OS and found that the NKT signature predicted better OS in only two datasets (E-MTAB-6134 and GSE57495). We generated a new index called the target score, which is defined as target score=HR(targets high expression)HR(targets low expression). In this equation, HR(targets high expression) reflects the influence of the NKT signature on patients' OS in PDAC samples with overexpression of specific targets for ferroptosis or pyroptosis, while HR (targets low expression) reflects the influence of the NKT signature on patients' OS in PDAC samples with downregulation of specific targets for ferroptosis or pyroptosis. A heatmap was generated to visualize the Target_score of different targets across various datasets (Figure 3D), which suggested universal survival benefits caused by high NKT signature expression in samples with high target expression. Survival analysis further demonstrated that the high expression of some ferroptosis and pyroptosis targets was a precondition for the NKT signature predicting a benefit in terms of OS (Figures 3E–H). Then, we determined the independency of the prognostic implication of the NKT signature from other infiltrated immune cells. Multivariate Cox regression showed that the cytotoxic T cell signature was associated with prolonged OS, independent of other components in the PDAC microenvironment (Supplementary Table 3).

An LTB(+)IL-7R(+)CD3(+) cell cluster was identified by scRNA-seq (cluster 0). We analyzed the top 100 upregulated genes in this cluster and showed that 51 genes were ribosome-related proteins (Figure 4A). Then, we performed survival analysis and found that the cluster 0 signatures could predict prolonged OS in PDAC in multiple datasets (Figures 4B–E).

In addition, we investigated the prognostic implications of the three states of B cells in PDAC. Overall, these signatures predict prolonged OS in only some cohorts (Supplementary Figure 3A). Notably, plasma cell infiltration predicted better OS in the three cohorts (MTAB, ICGC, and GSE57495) (Supplementary Figure 3B).

The XCELL algorithm calculated the infiltration percentage of HSCs in the PDAC microenvironment; however, no direct evidence indicated the existence of intratumoral HSCs in previous studies, and scRNA-seq was performed here. It is biologically plausible that some cell clusters retain parts of the molecular signatures of HSCs. For example, we found that CD34, as a classical marker of HSCs, was significantly upregulated in cluster 19, which we identified through scRNA-seq (logFC = 1.48). Furthermore, we conducted a meta-analysis to investigate whether the HSC signature was associated with the patient prognosis. The results showed that increased HSC signatures significantly predicted prolonged OS in PDAC (HR = 0.72, 95% CI 0.61–0.85) (Figure 5A). To further analyze the mechanism by which HSC signatures could predict prolonged OS in PDAC, we performed WGCNA to explore the gene modules associated with HSC signatures in the MTAB cohort, which is the largest PDAC cohort with completed transcriptome and follow-up data (Figure 5B). Interestingly, the dark turquoise module was negatively associated with HSC signatures and the OS and DFS of patients (r = −0.47, r = −0.31, and r = −0.28, respectively) (Figure 5C). According to the criteria GS > 0.2 and MS > 0.8, we identified 5 core genes associated with HSC signatures (Figures 5D,E) and patient prognoses (Figure 5F). Given that all five genes were related to unfavorable prognoses, we speculated that these genes may be differentially expressed between tumor and normal tissues. Then, we validated the prognostic implications of these five hub genes in six other cohorts, where each gene was associated with unfavorable OS in at least two validation datasets (Supplementary Figure 4). Using the integrated data of the tumor transcriptome from the TCGA database and the normal pancreas transcriptome from the Genotype-Tissue Expression (GTEx) database, we showed that LDHA, SLC2A1, and PGK1 were upregulated in tumor samples (logFC > 2, P < 0.05).

Myeloid cells are important components in the PDAC microenvironment and consist of macrophage cells, myeloid derived suppressive cells, dendritic cells, and granulocytes. Through scRNA-seq, we identified nine myeloid cell clusters (clusters 3–6, 10, 12, 22–23, and 25) and further annotated them as macrophages, mast cells, plasma dendritic cells, and other myeloid cells. Pseudotime trajectory analysis revealed seven different cell states (Figure 6A) and showed the distributions of cell states along with pseudotime flows (Figure 6B). We also mapped cell classifications to pseudotime trajectories (Figure 6C). The difference in the transcriptome between the two cell states with the longest pseudotime span was compared, wherein many antigen-presenting molecules and complement-associated genes were lost along the pseudotime flow (Figure 6D). Unsupervised consensus clustering identified two independent clusters based on the pseudotime-related differentially expressed genes in the MTAB cohort (Figures 6E–G). Notably, the OS of patients in cluster 1 was significantly better than that of patients in cluster 2 (P < 0.05) (Figure 6H). This result was marginally confirmed in the ICGC cohort (Supplementary Figure 5). Next, we compared the activity of 29 immune signatures between cluster 1 and cluster 2. The results showed that cluster 1 harbored more activated immune signatures, which may account for its survival advantage (Figures 6I,J; Supplementary Figure 4).

Fibroblasts or pancreatic satellite cells exist in normal pancreatic tissues; however, many studies have suggested that the formation of tumors alters the states and function of these cells. Through scRNA-seq analysis, we compared the mRNA expression landscape of fibroblasts and satellite cells between PDAC and normal pancreas tissues (Figure 7A). We found that multiple genes involved in stromal formation, such as COL1A2, COL3A1, COL1A1, FN1, TIMP1, DCN, and LUM, were upregulated in CAFs (logFC > 2, adjusted P < 0.05). However, the genes that were downregulated in CAFs have rarely been investigated in CAFs. Hence, we investigated whether PDAC cells could induce alterations in these genes using in vitro experiments. Given that TGF-beta is an important mediator regulating the crosstalk between CAFs and PDAC cells, we also established a positive control group using TGF-beta to mimic the activated states of CAFs. Our results showed that the relative mRNA expression of ADIRF, MT2A, MT1M, and JUNB was downregulated after TGF-beta treatment and/or tumor stimulation, while the expression level of C11orf96 was upregulated, even after tumor stimulation (Figure 7B). Gene Ontology (GO) analysis indicated that the upregulated genes in CAFs were mainly related to the extracellular matrix and structural organization, while the downregulated genes in CAFs were associated with the response to mental ions (Figure 7C). Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis also showed that the differentially expressed genes were enriched in mineral absorption and the extracellular-receptor interaction (Figure 7D).