Single-cell transcriptomics reveal the prognostic roles of epithelial and T cells and DNA methylation-based prognostic models in pancreatic cancer

Abstract

Background

Pancreatic adenocarcinoma (PDAC) exhibits a complex microenvironment with diverse cell populations influencing patient prognosis. Single-cell RNA sequencing (scRNA-seq) was used to identify prognosis-related cell types, and DNA methylation (DNAm)-based models were developed to predict outcomes based on their cellular characteristics.

Methods

We integrated scRNA-seq, bulk data, and clinical information to identify key cell populations associated with prognosis. The TCGA dataset was used for validation, and cell composition was inferred from DNAm data. Prognostic models were constructed based on cell-type-specific DNAm markers, and genomic features were compared across risk groups. Nomograms were created to assess treatment responses in different risk levels.

Results

Epithelial and T cells were major prognostic factors. Genomic analysis showed that epithelial cells in PDAC followed a malignant trajectory. DNAm data from TCGA confirmed the association of higher epithelial and T cell proportions with worse prognosis. Prognostic models based on DNAm markers of these cells effectively predicted patient survival, especially 5-year overall survival (AUC = 0.834). High-risk group with epithelial cell model showed altered pathways (tight junctions, NOTCH, and P53 signaling), while high-risk group with T cell model had changes in glycolysis, hypoxia, and NOTCH signaling, with more KRAS or TP53 mutations. Low-risk groups in the T cell model displayed stronger antitumor immune responses. Treatment predictions and nomograms were developed for clinical use.

Conclusions

scRNA-seq and DNAm data integration enabled the creation of predictive models based on epithelial and T cell-specific methylation patterns, offering robust prognosis prediction for PDAC patients.

Introduction

Pancreatic cancer stands as a highly aggressive malignancy within the gastrointestinal tract, characterized by a bleak prognosis, and it may become the second leading cause of cancer-related deaths by 2030 [1, 2]. Pancreatic ductal adenocarcinoma (PDAC), constituting the predominant pathological type of pancreatic cancer, accounts for 95% of cases, with a mere 5% 5-year overall survival (OS) [3]. Due to its distinctive anatomical location and non-specific clinical symptoms, a significant proportion of PDAC patients receive a diagnosis at an advanced or metastatic stage, precluding the opportunity for surgical intervention in over 80% of cases [4]. Concurrently, a substantial number of patients undergoing surgery succumb to recurrence and/or metastasis within 2 years [5]. Moreover, the lack of effective treatments has resulted in minimal improvement in the prognosis and survival rates of patients over the years [6,7,8]. A critical imperative exists to elucidate the molecular mechanisms governing PDAC progression and unearth prognostic signatures for patients.

The dynamic interplay between pancreatic epithelial/cancer cells and stromal cells orchestrates the tumor microenvironment (TME), a pivotal factor in PDAC progression [9]. T cells, particularly FOXP3 + Tregs, critically shape the immunosuppressive TME in PDAC [10, 11]. The elimination of stromal cells leads to a more aggressive disease, and depleting myofibroblasts in mice has been linked to reduced overall survival [12]. Epithelial cells, serving as the origin of malignant neoplastic cells, significantly contribute to tumor development [13, 14]. Recently, single-cell RNA sequencing (scRNA-seq) has enabled high-resolution, cell-type-specific transcriptome analysis, offering insights into exact cellular and molecular changes at the single-cell level. This approach facilitates the dissection of hidden heterogeneity within cell populations [15,16,17]. Peng et al. revealed two ductal cell subpopulations in PDAC, characterized by varying malignancy levels and distinct cell markers [18]. Elyada et al. identified three types of cancer-associated fibroblasts (CAFs): myofibroblastic CAFs, inflammatory CAFs, and antigen-presenting CAFs, with their plasticity offering promising therapeutic avenues to improve patient prognosis [19]. However, it is crucial to acknowledge potential limitations stemming from the typically small sample sizes used in scRNA-seq analyses. To augment our understanding of the relationship between cell subpopulations and prognosis, it is imperative to validate scRNA-seq findings by integrating a larger patient cohort with clinical phenotypes (such as survival information). This approach can provide deeper insights into the implications of cellular subtypes for patient prognosis.

Distinctive functions and phenotypes exhibited by various cell types within tissues are, in part, regulated by epigenetic mechanisms, among which DNA methylation (DNAm) plays an indispensable role in maintaining cell identity and functional phenotypes by modulating gene expression [20,21,22,23,24]. Aberrant methylation of CpG islands and promoter regions results in the transcriptional silencing of pancreatic cancer suppressor genes and/or abnormal activation of oncogenes, thereby promoting its initiation and progression [25,26,27,28]. A prognostic model was developed to predict the survival of PDAC patients based on four DNAm-driven genes [29]. Additionally, Xiao et al. established novel prognosis-related molecular subgroups based on the DNAm signature [30]. Nevertheless, further research is imperative to elucidate whether these DNAm differences can accurately represent alterations in specific cell populations and to gain a more comprehensive understanding of the pathological processes underlying these variations. Efforts have been made to infer cell population compositions in DNAm profiles based on paired RNA data and scRNA-seq cell clusters [31,32,33]. In other words, it enables researchers to dissect the alterations in cell population underlying changes in DNAm profiles and explore their roles in prognosis.

In this study, we dissected the intratumoral heterogeneity of PDAC, identifying key prognostic cell subpopulations through scRNA-seq, notably epithelial and T cells with distinct genomic and functional characteristics. By leveraging these insights, we developed prognostic models grounded in DNAm markers specific to these cell types, providing a means to translate single-cell transcriptomic discoveries into detectable signals within bulk tissue. This approach not only deepens our understanding of PDAC’s molecular landscape but also offers a practical strategy to bridge the gap between advanced single-cell technologies and clinically actionable biomarkers, paving the way for more precise prognosis and personalized treatment strategies.

Methods

Patients, clinical information, and public datasets

Tumor and adjacent nonmalignant pancreatic tissue samples were collected from eight PDAC patients who underwent resection at Zhejiang Cancer Hospital. All participants provided informed consent, and the study was approved by the Ethics Committees of Zhejiang Provincial People’s Hospital (ZJPPHEC2023I-089) and Zhejiang Cancer Hospital (IRB-2023-1039). PDAC diagnoses were confirmed by two independent pathologists. Clinical details are provided in Table 1, and these samples were subsequently analyzed using DNA methylation microarray experiments.

Table 1 Clinical characteristics of PDAC patients

PDAC single-cell RNA sequencing datasets (GSE155698, GSE201230, GSE212966, GSE214295, and GSE156405) were retrieved from the Gene Expression Omnibus (GEO) [34]. GSE155698 includes PDAC samples from surgical (n = 6) and fine needle biopsy (n = 10) specimens, GSE201230 comprises fetal pancreas samples at 15–20 weeks gestational age. GSE212966, GSE214295, and GSE156405 include six, three, and five PDAC samples, respectively.

TCGA-PAAD project, containing bulk RNA sequencing (RNA-seq) data, DNAm, and somatic mutation data, was downloaded from University of California Santa Cruz (UCSC) Xena (185 tumor and 10 normal samples) [35, 36]. Latest clinicopathological data were obtained from The Cancer Genome Atlas (TCGA) database as of October 10, 2020. ICGC-PACA-AU project, including bulk RNA-seq, DNAm, and clinical data (250 tumor samples), were sourced from the International Cancer Genome Consortium (ICGC) database [37]. Meanwhile, six bulk RNA-seq datasets comprising 659 samples were selected as external validation cohort. Five datasets were obtained from GEO, including GSE85916 (79 samples), GSE28735 (42 samples), GSE57495 (63 samples), GSE62452 (64 samples), GSE71729 (123 samples), and the ArrayExpress dataset E-MTAB-6134 containing 288 samples, was included. During the processing of expression matrices, platform annotation files were used to map probes to gene symbols, and for genes with multiple probes, the average expression was taken as the gene’s expression level. The “SVA” R package was used to remove batch effects between different datasets [38]. Both transcriptome information and clinical information were downloaded for survival analysis. Pancreatic cancer samples with survival times exceeding 30 days were included in the analysis. As publicly accessible databases, the datasets from GEO, TCGA, ICGC, and ArrayExpress are freely available for download from their respective websites.

ScRNA-seq data processing

Raw single-cell data were downloaded and processed using the “Seurat” R package (v4.3.0.1) for quality control [39]. All functions were run with default parameters unless specified otherwise. Cells were filtered out if they met the following criteria: (1) fewer than 200 or more than 8000 genes; (2) genes expressed in fewer than 3 cells; (3) mitochondrial gene content exceeded 25%.

Data were normalized using “NormalizedData” (scale factor = 10,000), and 2000 highly variable genes were identified with “FindVariableFeatures”. Gene expression levels were scaled via “ScaleData”. Batch effects were mitigated with the “Harmony” package for dataset integration across different specimens, and the “SCTransform” function was used to remove cell cycle influences.

Principal component analysis (PCA) was performed on 2000 genes using “RunPCA”, significant principal components (PCs) were identified by the “ElbowPlot”, and the top 15 PCs were chosen for further dimension-reduction analysis. Cell clustering was conducted using “FindNeighbors” and “FindClusters” functions, with Clustering trees showed that resolution set to 0.4 achieved the clearer clustering results by “clustree”. The cell atlas was visualized using Uniform Manifold Approximation and Projection (UMAP) by “RunUMAP” [40].

Cell cluster annotation

Distinct cell clusters were identified and annotated using established cell markers from the literature [18, 41,42,43,44,45]. The COSG method (v0.9.0), based on cosine similarity, was employed for accurate and efficient marker gene identification [46]. Key markers included REG1A, REG1B, PRSS1, CTRB1, CTRB2 (acinar); CD79A, MS4A1, CD79B, VPREB3 (B cell); CDH5, VWF, PLVAP, CLDN5, RAMP2 (endothelial cell); EPCAM, KRT19, KRT8, KRT18, KRT7 (epithelial cell); LUM, DCN, COL1A1, SFRP2 (fibroblast); NAMPT, FCGR3B (granulocytes); AIF1, CD68, CD14 (macrophage); CLU, CPA3, TPSAB1 (mast cell); NKG7, GZMB (NK cell); IGJ, CD79A (plasma cell); RGS5, NDUFA4L2 (stellate cell); CD3D, CD3E, CD3G, IL7R, CD2 (T cell). Expression levels of well-established marker genes and the top five differentially expressed genes for each cell populations were depicted through violin plots and dot plots, respectively. The “VlnPlot”, “FeaturePlot”, and “DotPlot” functions were utilized to present the expression levels of specific marker genes. Sub-clustering of T cells followed the aforementioned methodology.

Prognostic cell cluster identification

The “SCISSOR” R package (v2.0.0), developed by Sun et al. [47], was employed to discern bulk phenotype-associated cell subpopulations. SCISSOR quantifies the similarity between single-cell and bulk data by gauging Pearson correlations for each cell-bulk pair, followed by optimizing a regression model to associate cells to phenotypes using a sparsity penalty and graph regularization, identifying the most relevant cells associated with the phenotype with high confidence. We applied SCISSOR to integrate PDAC bulk data (ICGC-PACA-AU), single-cell data (GSE155698, GSE212966, GSE214295, and GSE156405), and survival data to identify survival-associated cells, respectively. Adhering to SCISSOR’s criteria, we set the alpha value at 0.5 to restrict cells (20% of total). A 100-fold cross-validation reliability test (P < 0.05) was performed using “reliability.test” function to confirm the robustness of the identified prognostic cells. SCISSOR classified cells into Surv⁻ (worse survival), Surv⁺ (better survival), or background (no significant association).

Differentially expressed genes (DEGs) in prognostic-related (PR) cells

DEGs in PR cells were identified using the “FindMarkers” function from the Seurat package. To reduce the impact of inherent biological functions, differential expression analysis was performed within each cell type individually. Since PR cells were unevenly distributed across most cell types, two comparison strategies were employed to ensure robust statistical analysis. PR cells were compared to background cells when PR cells dominated a cell type, while Surv⁺ and Surv⁻ cells were compared when PR cell proportions were more balanced, such as in epithelial and macrophage cells. Genes with an absolute log₂ fold change (log₂FC) ≥ 0.5 and a Bonferroni-adjusted P value < 0.05 were considered DEGs.

Single-cell CNV analysis

Single-cell copy number variation (CNV) analysis was conducted using the “inferCNV” R package (v1.15.3) [48]. Chromosomal variations were inferred by comparing gene expression intensities across genomic regions in tumor cells relative to reference immune cells. The analysis utilized a raw counts matrix, an annotation file, and a gene/chromosome position file. For epithelial cells, CNV was calculated with immune cells as the reference.

Key parameters included “denoise” default hidden Markov model (HMM) settings and a cutoff of 0.1. The Bayesian latent mixture model was applied to ascertain posterior probabilities and reduce false-positive CNV calls. A heatmap visualizing expression intensities across each chromosome was generated, with regions of gain in red and loss in blue. CNV scores were computed as the quadratic sum of CNV regions.

Modular score calculation and signal visualization

Gene expression modular scores were calculated and visualized using the Seurat R package (v4.3.0.1). Based on previously reported studies [16, 41, 49], we selected markers for epithelial-mesenchymal transition (EMT), cancer stem cells, and genes involved in malignant progression in PDAC (Table S1). The “AddModuleScore” function was employed to compute cancer stem cell, mesenchymal, and malignant scores for epithelial cells. Expression of relevant marker genes were projected onto dimension-reduction plots using “FeaturePlot” function, and the signal distribution across PR epithelial cells was visualized with “VlnPlot” function.

Cell cluster annotation in fetal pancreata

Cell-type markers and the “SingleR” R package (v2.0.0) were used to identify and annotate cell clusters within the GSE201230 dataset [50]. Based on previous studies [51, 52], we selected markers for specific pancreatic cell types: INS, GCG, PAX6 (endocrine); CPA1, PRSS1, CLPS (acinar); COL3A1, DCN, VIM (mesenchymal); TRAC, LYZ, RAC2, PTPRC (immune); VWF, ADGRL4, ANGPT2, PECAM1 (endothelial); CRYAB, CDH19, SOX10 (Schwann); HBB, HBG2, HBA1 (erythroblasts); and EPCAM, KRT19, PROX1, PDX1 (epithelial).

Pseudo-time trajectory analysis

To fully understand the gene expression patterns and evolution of epithelial PR cell subtypes in different disease states, the “Monocle” R package was used to perform pseudo-time analysis on normal epithelial cells from the GSE201230 dataset and epithelial PR cells from the GSE155698 dataset [53]. A Monocle object was created using the “newCellDataSet” function with the parameter expressionFamily = negbinomial.size. Highly variable genes were selected for pseudo-time analysis via the “dispersionTable” function. Trajectory analysis was conducted using genes with a mean expression ≥ 0.1 and a dispersion_empirical ≥ 1 * dispersion_fit, as determined by Monocle2, to order cells along pseudo-time. Dimensionality reduction was performed using the “reduceDimension” function (reduction_method = “DDRTree”, max_components = 2), and the cell trajectory was visualized using the ‘plot_cell_trajectory’ function, displaying the minimum spanning tree. Trajectory also visualized with tree plot to represent the hierarchical relationships between cells.

Deconvolution bulk DNAm profiles based on the scRNA-seq using scDeconv

Bulk DNAm data from ICGC and TCGA cohorts were derived from tissue-level samples, limiting the ability to explore cell-type heterogeneity. To address this, we applied the scDeconv algorithm, which performs high-resolution cell-type deconvolution on bulk DNAm profiles by integrating paired bulk RNA-DNAm data with scRNA-seq data in a trans-omics approach [33]. The “scDeconv” R package was used to estimate cell-type proportions for each TCGA-PAAD patient based on scRNA-seq data (GSE155698), paired bulk RNA and DNAm data were provided. Additionally, this algorithm trained an ensemble model using DNAm features, enabling the prediction of relative cell types in other bulk DNAm datasets even without paired bulk RNA-seq data, this model further applied to the ICGC-PACA-AU dataset and our 8 PDAC samples. All parameters were set to default.

Deconvolution bulk RNA-seq data based on the scRNA-seq using BayesPrism

Bulk RNA-seq data from pancreatic cancer samples were normalized to transcripts per million (TPM) and log₂(TPM + 1) to ensure data comparability. Genes with > 50% missing values were excluded [54]. For both the TCGA and ICGC databases, bulk transcriptomic data were matched with corresponding DNA methylation data for integrated analysis by using unique sample IDs. BayesPrism algorithm (“BayesPrism” R package) was used to perform high-resolution deconvolution of bulk RNA-seq data from integrated GEO datasets (659 samples) based on the single-cell RNA-seq data (GSE155698) [55], and the cell proportion was estimated by their gene expression levels.

DNAm profiling for Zhejiang PDAC samples

Genomic DNA from Zhejiang PDAC samples was extracted using the QIAamp DNA kit (Qiagen) and quantified with a Qubit HS DNA assay (Thermo Fisher Scientific). DNA methylation analysis was performed on ≥ 0.5 µg of bisulfite-converted DNA using the EZ methylation Kit (Zymo Research), followed by amplification, and hybridization to Infinium Human MethylationEPIC BeadChip (850 K, Illumina). Arrays were scanned using the Illumina HiScan SQ scanner.

DNAm data processing

Quality control and analysis of DNAm data were conducted using the “ChAMP” R package [56]. DNAm data from TCGA and ICGC were generated using the Illumina Infinium Human Methylation 450 K BeadChip. Raw fluorescence data were initially stored in IDAT files and then converted to β-values by “champ.load” function. Methylation levels were scored on a scale of 0 to 1 (unmethylated to totally methylated). Subsequent cell proportion analyses based on DNA methylation profiles utilized the intersecting methylation sites shared between the TCGA-PAAD , ICGC-PACA-AU dataset, and the Zhejiang samples. CpG sites with > 70% missing β-values, located on sex chromosomes, or linked to single-nucleotide polymorphisms (SNPs) were excluded. Missing values were imputed using a K-nearest neighbor algorithm, and probe sites from 2 kb upstream to 200 bp downstream of the transcription start site were selected. Mean β-value across multiple CpG sites for each gene was calculated as the average methylation level. Batch effects were corrected via beta mixture quantile expansion (BMIQ) normalization and singular value decomposition (SVD).

Differentially methylated probe (DMP) analysis

DMPs between 185 tumor and 10 normal samples in TCGA-PAAD were identified using the “champ.DMP” function. P values were adjusted using the Benjamini–Hochberg method to control the false discovery rate (FDR), and the delta β-value was calculated as the absolute difference between the group mean β-value for each probe. Probes with FDR < 0.01 and |delta β-value|> 0.2 were defined as DMPs. Volcano plot visualized the DMPs distribution, and gene annotations were obtained using GRCh38 from the GENCODE project.

Consensus clustering

Unsupervised consensus clustering of the TCGA-PAAD cohort was performed using DMPs and the “ConsensusClusterPlus” R package [57]. This methodology grouped cases into distinct clusters based on signatures, iterated 1000 times using bootstraps with K-means and Euclidean distance to ensure cluster consistency. Each bootstrap used 80% of TCGA samples, with cluster numbers (k) were set from 2 to 9. The optimal number of clusters was determined by the cumulative distribution function (CDF) and Delta plot, ensuring stabilization of the CDF plot without significant increments.

Survival analysis

Kaplan–Meier (KM) survival analysis was conducted using the “survival” and “survminer” R packages. Log-rank tests were used to compare overall survival between patient groups to assess the statistical differences. We investigated the associations between estimated cell-type abundance and prognosis. PDAC patients were divided into low and high cell proportion groups using the optimal cutoff value calculated via the “surv_cutpoint” function, which is also used to determine the cutoff expression value for methylation sites-associated genes.

Prognostic model construction and validation

Prognostic risk models were constructed using 178 PDAC samples from the TCGA dataset as the training set, and validated with external set of 250 PDAC samples from the ICGC dataset. Based on a human methylome atlas [58], which offers a repository of tissue-specific biomarkers, we curated 25 pancreatic ductal cell-specific CpG sites (CpGs) and 25 T cell-specific CpGs as potential prognostic markers. Univariate Cox regression analysis (P < 0.05) identified prognostic CpGs using “survival” R package, followed by Least Absolute Shrinkage and Selection Operator (LASSO) regression with the “glmnet” R package to reduce overfitting and remove highly similar sites. Multivariate Cox regression was applied to the selected CpGs, with tenfold cross-validation calculating the confidence interval for each lambda value.

Risk scores for each patient were calculated using the formula based on the β values of CpG sites and the accompanying regression coefficients:

$${\text{Formula}}:{\text{ RiskScore}}\, = \,\left( {\beta_{{{\text{CpG1}}}} \, \times \,{\text{Coef}}_{{{\text{CpG1}}}} } \right)\, + \,\left( {\beta_{{{\text{CpG2}}}} \, \times \,{\text{Coef}}_{{{\text{CpG2}}}} } \right)\, + \cdots \, + \,\left( {\beta_{{{\text{CpGn}}}} \, \times \,{\text{Coef}}_{{{\text{CpGn}}}} } \right).$$

Patients were classified into high- and low-risk groups based on the median risk score, KM analysis and log-rank tests assessed OS between risk groups. Time-dependent receiver operating characteristic (ROC) curves were generated with “timeROC” and “survminer” R packages, and the area under the ROC curve (AUC) values were used to evaluated the prognostic model accuracy. The ICGC cohort (ICGC-PACA-AU) validated the reliability and suitability of the models.

Independent prognostic analysis and nomogram construction

Univariate and multivariate Cox regression analyses were conducted on risk scores and clinicopathological data (age, gender, TNM stage, alcohol/smoking history, histologic grade, tumor stage, and tumor grade) to identify independent prognostic factors. Nomogram models predicting 1-, 2-, and 3-year OS were constructed based on selected independent prognostic factors using the “rms” R package. Calibration curves and decision curve analysis (DCA) were used to compare predicted OS with actual outcomes, assessing predictive accuracy.

Clinical characteristics and relationship with Moffitt subtypes in risk groups

The distribution of clinical characteristics between the risk groups was based on the patient clinical data. We explored the association between risk groups and Moffitt subtypes, with basal-like subtypes linked to worse outcome [59]. Moffitt subtype classification of TCGA samples was performed using the Purity Independent Subtyping of Tumors (PurIST) algorithm, a clinically validated single-sample classifier for pancreatic cancer subtypes [60].

Somatic mutation analysis

Somatic mutations in TCGA PDAC samples were analyzed using mutation annotation format (MAF) data with the “Maftools” R package [61], and visualized through waterfall plots by “oncoplot” function. Tumor mutation burden (TMB) was computed by “TMB” function, and mutation frequencies and TMB were compared between risk groups.

Functional enrichment analysis and pathway analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed using the “clusterProfiler” R package (FDR < 0.05) [62] to identify enriched pathways and biological processes, applying the Hypergeometric Test. Gene Set Variation Analysis (GSVA) evaluated pathway activity variation between risk groups (FDR < 0.05) [63], utilizing the “c2.cp.kegg.v2023.1.Hs.symbols.gmt” or “h.all.v2023.1.Hs.symbols.gmt” from the Molecular Signature Database [64], and results visualized via “pheatmap” and “ggplot2” R packages. Correlations between biological functions and RiskScores (r > 0.3) were calculated using “corrplot” package. The “PROGENy” R package assessed pathway deregulation across 14 typical tumor-related pathways [65], comparing groups using the Wilcoxon test.

Immune infiltration analysis

Immune cell infiltration in the TME was estimated using the RNA-seq expression matrix of TCGA-PAAD with XCELL [66], TIMER [67], and CIBERSORT [32] algorithms, implemented via the IOBR R package [68]. Statistical differences between groups were assessed using the Wilcoxon rank-sum test. Stromal and immune cell levels were quantified using the “ESTIMATE” R package, generating StromalScore, ImmuneScore, and EstimateScore. Single-sample gene set enrichment analysis (ssGSEA) was applied to calculate immune infiltration scores for each sample, with correlations to RiskScores examined using the “corrplot” R package.

Six immune subtypes, defined by variations in macrophage/lymphocyte signatures, Th1:Th2 ratios, intratumoral heterogeneity, aneuploidy, neoantigen load, cell proliferation, and immunomodulatory gene expression, were classified using the “ImmuneSubtypesClassifier” R package [69]. Additionally, we analyzed the correlation between RiskScores and the activities of steps in the cancer-immunity cycle using data from the Tracking Tumor Immunophenotype (TIP) database [70].

Immunotherapy and chemotherapy response prediction

We selected representative immune checkpoints (ICs) from previous studies and compared their expression across groups. The Tumor Immune Dysfunction and Rejection (TIDE) online algorithm was used to assess immune escape potential for TCGA samples [71]. Using tumor transcriptional data and the Genomics of Drug Sensitivity in Cancer (GDSC) [72], we calculated half-maximal inhibitory concentration (IC50) values for common PDAC chemotherapeutic drugs via the “pRRophetic” R package [73]. Sensitivity to chemotherapy and immunotherapy between risk groups was compared using the Wilcoxon test.

Statistical analysis

R (v4.1.3) and associated packages were used for statistical analyses and data visualization. Kaplan–Meier survival analysis with log-rank tests compared OS across patient groups. LASSO and Cox regression were employed for predictive model development, while Pearson’s correlation assessed correlation coefficients. Differences between normally distributed groups were evaluated using the Student’s t test, and the nonparametric Wilcoxon rank-sum test was applied for continuous variables. The chi-square test analyzed somatic mutation frequency and immunotherapy response, while functional enrichment used the Hypergeometric test. Cox regression significance was assessed with the Wald test. P values and FDR q-values < 0.05 were considered statistically significant.

Results

Prognostic cell populations identified via scRNA-seq: epithelial and T cells

The study flowchart is shown in Fig. 1. We analyzed the scRNA-seq dataset for PDAC (GSE155698), consisting of 6 surgical and 10 fine needle biopsy specimens from untreated patients. After quality control and normalization, 43,057 cells were clustered into 25 groups and annotated to 12 cell types based on DEGs and known markers (Fig. 2A, B, D, E, Fig. S1A, B, C, D). Substantial intertumoral heterogeneity was observed across patients (Fig. 2C). To identify survival-related clusters, we integrated 250 ICGC-PACA-AU bulk samples with survival data. A total of 5,767 PR cells were identified and categorized into three types (Fig. 2F). Surv⁻ cells (4,172, 9.6%) were linked to worse survival, while Surv⁺ cells (1,595, 3.7%) were associated with better outcomes. Surv⁻ cells mainly originated from epithelial, T cell, and granulocyte populations (Fig. 2G), whereas Surv⁺ cells were mainly composed of epithelial, fibroblast, stellate, and macrophage populations (Fig. 2H). Consistent trends in the distribution of PR cell populations were observed in the GSE212966, GSE214295, and GSE156405 datasets, further supporting the robustness of the identified results (Fig. S2–S4).

Fig. 1

Workflow of the study. Step I: PR cells in PDAC were identified through scRNA-seq, and their biological functions were assessed using GO/KEGG pathway analyses. Step II: Malignancy of PR epithelial cells was assessed using InferCNV, malignant scoring, signal visualization, and pseudo-time trajectory analyses. Step III: Cell proportions for each TCGA-PAAD sample were inferred using scDeconv algorithm. Significant differences in cell proportions between TCGA-PAAD DNAm subgroups with distinct prognoses further support the association between PR cells (epithelial and T cells) and patient prognosis. Step IV: Given the critical role of epithelial and T cells in prognosis, prognostic and nomogram models were constructed and validated leveraging the corresponding cell-type-specific signatures. Step V: Relevance analysis was conducted between risk groups determined by prognostic models, covering clinical characteristics, mutation landscapes, enrichment pathways, immune cell infiltration, ICB treatment response prediction, and drug sensitivity. PR cells, prognostic-related cells; PDAC, pancreatic ductal adenocarcinoma; scRNA-seq, single-cell RNA sequencing; GO, Gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; TCGA, The Cancer Genome Atlas Program; DNAm, DNA methylation; K-M, Kaplan–Meier; ROC, receiver operating characteristic; TME, tumor microenvironment; ICB, immune-checkpoint blockade

Fig. 2

PDAC cells and distinct PR cells characterized through scRNA-seq analysis. A UMAP visualization depicting clusters within PDAC patient samples, with each color-coded area representing an individual cluster. B UMAP plot illustrating 43,057 cells and highlighting the predominant cell types in PDAC, with each color-coded region denoting a specific cell type. C Proportional distribution of each cell type across individual samples. D Identification of cell populations using top five DEGs. Dot color indicates average expression level, and dot size reflects percentage expression. E Expression patterns of well-established marker genes across the identified PDAC cell types. Each violin represents the gene expression density for a particular marker gene within the corresponding cell type, and violin color indicates the cell type. F UMAP plot showcasing PR cells identified by the SCISSOR algorithm, with red and blue dots representing Surv⁻ (poor survival) and Surv⁺ (favorable survival) cells. G, H Distribution of Surv⁻ cell populations (G) and Surv⁺ cell populations (H) across different cell types. PDAC, pancreatic ductal adenocarcinoma; PR cells, prognostic-related cells; UMAP, Uniform Manifold Approximation and Projection; DEGs, differentially expressed genes

In PR epithelial cells, Surv⁻ cells had upregulated genes related to RNA splicing, wound healing, and unfolded protein response, while downregulated genes were involved in the electron transport chain and aerobic respiration (Fig. 3A, C). KEGG analysis revealed associations with focal adhesion, extracellular matrix (ECM)-receptor interaction, and chemical carcinogenesis (Fig. 3B, D). Surv⁻ T cells showed upregulation in pathways linked to erythrocyte differentiation and NOD signaling, with downregulated genes involved in cytoplasmic translation, ribosome biogenesis, and rRNA processing (Fig. 3E, G). KEGG analysis suggested links to infection pathways (Fig. 3F, H). Functional enrichment for other cell types is detailed in Fig. S5 and S6.

Fig. 3

Functional Analysis of PR Cells. A, C Representative enriched GO terms (Biological Process) for upregulated (A) and downregulated (C) genes (|log₂FC|≥ 0.5 and adjusted P < 0.05) in Surv⁻ epithelial cells, respectively. B, D KEGG pathways associated with upregulated (B) and downregulated (D) genes (|log₂FC|≥ 0.5 and adjusted P < 0.05) in Surv⁻ epithelial cells, respectively. E, G Representative enriched GO terms (Biological Process) for upregulated (E) and downregulated (G) genes (|log₂FC|≥ 0.5 and adjusted P < 0.05) in Surv⁻ T cells, respectively. F, H KEGG pathways associated with upregulated (F) and downregulated (H) genes (|log₂FC|≥ 0.5 and adjusted P < 0.05) in Surv⁻ T cells, respectively. I–K Normalized expression of marker genes related to myofibroblasts (I), complement-related genes (J), and inflammatory cancer-associated fibroblasts (iCAFs) (K) across PR fibroblast subsets. L UMAP plot of T cells re-clustered from Fig. 1B, with each color-coded region representing a distinct cell type. M, N Expression of effector (M) and exhaustion markers (N) for PR T cells on the UMAP map. Color key from gray to pink indicates relative expression levels from low to high. PR cells, prognostic-related cells; GO, Gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; iCAFs, inflammatory cancer-associated fibroblasts

In PR fibroblasts, myofibroblastic cancer-associated fibroblasts (myCAF)-related markers were significantly elevated in Surv⁺ cells compared to background fibroblasts (Fig. 3I), indicating a myCAF-like phenotype linked to better survival [12]. Surv⁺ CAFs exhibited high expression of complement-related gene expression, including C3, C7, CFH, and CFI, potentially influencing immune regulation and responses (Fig. 3J). Conversely, iCAF-related markers expression was low in PR fibroblasts (Fig. 3K), and Surv⁺ cells showed even lower expression levels than background fibroblasts comparatively. Furthermore, CAF subtypes showed dynamic reversibility between phenotypes [74].

Recent studies reveal heterogeneity and dynamic phenotypic changes of macrophages during tumor progression [18, 75,76,77,78,79,80]. Surv⁻ macrophages (307 cells) and Surv⁺ macrophages (140 cells) represented only a small proportion (4.6%) of total macrophages (9699 cells), with inconsistent and unclear prognostic implications. Limited numbers of 284 Surv⁻ acinar, 112 mast, and 26 plasma cells were identified and excluded from further analysis.

T cell exhaustion in the immunosuppressive TME with poor PDAC prognosis

A total of 1,073 Surv⁻ T cells were identified, and all T cells clustered into eight subpopulations and annotated into three cell types (Fig. 3L; Fig. S7A–E). The proportion of PR T cells showed no significant variation across subpopulations (Fig. S7F), though Tregs exhibited a slightly higher prevalence of Surv⁻ cells, potentially suppressing effector T cells’ antitumor responses [81]. Low expression of effector markers (GZNB, GZLY; Fig. 3M) and universally high expression of exhaustion markers (EOMES, GZMK, TIGIT, KLF2; Fig. 3N) suggested T cell exhaustion and may lose their anticancer capability, linked to poor PDAC prognosis [16, 44, 45, 82].

Malignant nature of epithelial cells with higher malignancy in Surv ⁻ subtype

Surv⁻ and Surv⁺ subtypes were identified in 1095 and 747 epithelial cells, respectively, but both exhibited malignant characteristics (Fig. 4A). CNV scores were significantly elevated in both compared to reference immune cells, with Surv⁻ showing higher CNV scores (P = 5 × 10⁻⁵; Fig. 4B). Both subtypes demonstrated EMT and tumor stemness, facilitating tumor invasion and metastasis (Table S1; Fig. 4C, D). Surv⁻ cells had higher malignancy scores, along with elevated mesenchymal markers (VIM, COL1A1; P_VIM = 2 × 10^–16; P_COL1A1 = 5.7 × 10^–11; Fig. 4E, F, G) and stem cell markers (CD44, NOTCH2, NOTCH1, NES; P_CD44 = 2 × 10^–16; P_NOTCH2 = 0.00096; P_NOTCH1 = 3 × 10^–4; P_NES = 0.01; Fig. 4H, I; Fig. S1E, F). Most PR epithelial cells were classified as type II malignant ductal cells described in the previous study [18] (Fig. 4J, K). Pseudo-time analysis using PR epithelial cells and normal epithelial cells from the fetal pancreas (Fig. 4L) showed cells progressed through three states: the starting branching point (state 3) and two branches (state 1 and state 2). 15- and 18-week fetal pancreas cells were mainly at the start of trajectory, while 19- and 20-week cells were at the end of branch 2. Both Surv⁻ and Surv⁺ epithelial cells along branch 1 (Fig. 4M, N), with Surv⁻ cells at the terminus and Surv⁺ cells moving toward it (Fig. 4O). PDAC markers associated with cell proliferation and metastasis (MUC1 [83], FXYD3 [18], CEACAM6 [84], KRT19[85]) were upregulated along the late stage of branch 1 gradually, indicating tumor progression (Fig. 4P).

Fig. 4

Malignant Features of PR Epithelial Cells. A Heatmap illustrating extensive CNVs in PR epithelial cells, with red and blue colors indicating high and low CNV levels, respectively. B Distribution of CNV scores among Surv⁻ and Surv⁺ epithelial cells. C–E Malignant score (C), cancer stem cell score (D), and mesenchymal score (E) across Surv⁻ and Surv⁺ epithelial cells. F, G Comparison of expression levels of mesenchymal cell markers between Surv⁻ and Surv⁺ epithelial cells. H, I Comparison of expression levels of tumor stem cell markers between Surv⁻ and Surv⁺ epithelial cells. J, K Expression levels of type I ductal cell markers (J) and type II ductal cell markers (K) plotted on the UMAP map, respectively. Color key from gray to red indicates relative expression levels from low to high. L UMAP plot of 8,836 cells illustrating the main cell types in fetal pancreata. Each color-coded region represents a distinct cell type. M Pseudo-time analysis of PDAC Surv⁻, Surv⁺ epithelial cells, and normal epithelial cells from fetal pancreata (dark green = early state, yellow = late state). N Left panel shows pseudo-time progression, and right panel shows the differentiation branches. O Density graph of pseudo-time in PDAC Surv⁻ and Surv⁺ epithelial cells. P Expression of representative genes (MUC1, FXYD3, CEACAM6, and KRT19) along the pseudo-time trajectory plot. Color key from purple to red indicates relative expression levels from low to high. PR cells, prognostic-related cells; CNVs, copy number variations; UMAP, Uniform Manifold Approximation and Projection; PDAC, pancreatic ductal adenocarcinoma

Inference of epithelial and T cell proportions from PDAC methylation and bulk RNA-seq data validating their prognostic impact

Overall, we identified 10,406 TCGA PDAC DMPs (log₂|FC|> 0.1, adjusted P < 0.05), revealing 4163 hypermethylated and 2413 hypomethylated genes (Fig. S8A). These DMPs impacted pathways related to tumor progression (Fig. S8B–H). Methylation-based analysis divided 178 PDAC samples into two prognostic subgroups (Fig. 5A–D, Fig. S8I). Notably, Cluster 1 patients exhibited worse overall survival (P = 0.002) and higher tumor purity (P = 3.2 × 10⁻⁹), characterized by increased proportions of epithelial and T cells (P < 0.05) (Fig. 5E–G).

Fig. 5

Inferred Proportions of Epithelial Cells and T Cells in DNAm Profiles and their association with Patient Prognosis. A Unsupervised clustering (K = 2) of TCGA PDAC samples into two molecular subgroups according to DNAm profiles. B Cumulative Distribution Function (CDF) for 10,406 Differentially Methylated Positions (DMPs). C CDF Delta area curve. D Sample distribution in two subgroups displayed in T-SNE plot. E Kaplan–Meier survival plot for the two subgroups. F Comparison of tumor purity between the two subgroups. G Boxplot illustrating proportions of distinct cell types within the TME for the two subgroups. The y-axis represents the absolute cell-type proportion (0 to 1) estimated by scDeconv, and the x-axis represents different cell types. Significant changes between the subgroups were assessed using the Wilcoxon rank-sum test (*P value < 0.05, **P value < 0.01, ***P value < 0.001). H, I Overall survival comparison between high- and low-proportion groups for inferred epithelial cell (H) and T cells proportion (I). J Proportional distribution of each cell type in 8 PDAC patients we collected. K Construction of the prognostic RiskScore_Epi model using six epithelial cell-type-specific CpG sites and their corresponding coefficients through LASSO-Cox regression. L Determination of the λ value of the LASSO model via tenfold cross-validation with minimal misclassification. M, N Overall survival analysis for high- and low-risk groups with a cutoff value of the median RiskScore_Epi in the TCGA train set (M) and ICGC validation set (N), respectively. DNAm, DNA methylation; TCGA, The Cancer Genome Atlas Program; ICGC, International Cancer Genome Consortium; CDF, Cumulative Distribution Function; DMPs, differentially methylated probes; T-SNE, t-distributed stochastic neighbor embedding; TME, tumor microenvironment; PDAC, pancreatic ductal adenocarcinoma; LASSO, Least Absolute Shrinkage and Selection Operator

Prognostic analysis based solely on scDeconv DNAm-inferred cell proportions indicated that higher epithelial and T cell proportions were linked to poor prognosis (P = 0.0036 and P = 0.034; (Fig. 5H, I). Analysis of tumor samples from our 8 PDAC patients (Zhejiang, Table 1) confirmed these findings, showing that higher epithelial and T cell proportions correlated with shorter survival, whereas higher endothelial and stellate cell proportions were linked to better prognosis (Fig. 5J). Furthermore, consistent conclusions were drawn from cell proportions inferred using the BayesPrism algorithm based on bulk RNA-seq data (659 samples), further supporting the robustness of the results (P = 0.017 and P = 0.0014; Fig. S9A, B).

Identification of epithelial and T cell-specific DNAm sites for RiskScore model construction

To develop a clinically feasible assessment, we screened prognostic methylation sites representing T cell and epithelial cell proportions, 19 specific CpG sites to epithelial cells and 14 for T cells were all significantly associated with PDAC prognosis (P < 0.05). After applying the LASSO-Cox model, six CpG sites were selected for the epithelial cell-specific model (Fig. 5K, L) and four for the T cell-specific model (Fig. S10A, B). Using β values and their corresponding risk coefficients, we calculated RiskScore_Epi and RiskScore_T following the subsequent formulas:

$$\begin{aligned} {\text{RiskScore }}\left( {{\text{Epi}}} \right)\, = & \,\left( {cg07768763\, \times \,1.254} \right)\, + \,\left( {cg11314569\, \times \, 2.050} \right) \, + \,\left( {cg05038288\, \times \, - 1.397} \right)\, \\ & + \,\left( {cg00788354\, \times \, - 2.791} \right)\, + \,\left( {cg27650540\, \times \,0.585} \right) \, + \, \left( {cg02265288\, \times \, - 1.239} \right); \\ \end{aligned}$$

$$\begin{aligned} {\text{RiskScore }}\left( T \right)\, = & \,\left( {cg16239536\, \times \,1.571} \right) \, + \,\left( {cg07545925\, \times \,1.438} \right) \, \\ & + \, \left( {cg25653947\, \times \, - 1.477} \right)\, + \,\left( {cg14150023\, \times \,0.689} \right). \\ \end{aligned}$$

Patients were stratified by the median RiskScore_Epi (− 1.104), with high-risk patients showing significantly lower survival rates (P < 0.001; Fig. 5M, 6A). The model’s AUC values for predicting 1-, 3-, and 5-year OS were 0.695, 0.785, and 0.834, respectively (Fig. 6C). In the validation cohort (ICGC-PACA-AU), survival outcomes differed significantly between high- and low-risk groups (P = 0.022; Fig. 5N, 6B), and our 8 PDAC patients with higher risk scores exhibited shorter survival time (Fig. 6G). Similarly, RiskScore_T showed predictive power (P = 0.029; Fig. S10E, H), with AUC values of 0.628, 0.733, and 0.782 for 1-, 3-, and 5-year survival (Fig. S10C). In the validation cohort, high-risk groups exhibited worse outcome (P = 0.013; Fig. S10F, G, I).

Fig. 6

Construction and Evaluation of Prognostic RiskScore_Epi Model and Nomogram Integrating Clinical Information. A, B For TCGA Train Set (A) and ICGC Validation Set (B): Top panel: Distribution plot of RiskScore_Epi in high- and low-risk groups. Middle panel: Distribution plot of survival time and survival status in high- and low-risk groups; X-axis represents patients with increasing RiskScore_Epi, and Y-axis represents survival time. Bottom panel: Beta values of the six CpG sites in high- and low-risk groups. C ROC curves and AUC analyses for 1-, 3-, and 5-year overall survival prediction of the RiskScore_Epi model in TCGA dataset. D Sankey diagram illustrating similarities and differences in sample composition before and after grouping based on two different criteria: TCGA DNAm subgroups (left) and RiskScore_Epi model (right). E, F Univariate E and Multivariate F Cox regression analysis of RiskScore_Epi and clinical information in TCGA dataset. G Heatmap of six CpG sites’ methylation levels, corresponding RiskScore_Epi, survival time, survival status, and risk groups of the 8 PDAC patients we collected. Color bar = β-value. H Nomogram for predicting 1-, 2-, and 3-year overall survival of PDAC patients based on RiskScore_Epi and prognostic factors. I Calibration curves of the nomogram at 1, 3, and 5 years. J DCA curves of the nomogram, RiskScore_Epi, and other prognostic factors. TCGA, The Cancer Genome Atlas Program; ICGC, International Cancer Genome Consortium; ROC, receiver operating characteristic; AUC, area under the curve; PDAC, pancreatic ductal adenocarcinoma; DCA, decision curve analysis

KM curves based on the expression levels of RiskScore_Epi-related genes (CUX2, GCSH) and RiskScore_T-related genes (CD3G, HMHA1) revealed significant prognostic differences, highlighting their relevance to survival outcomes (Fig. S9E–H). Over half of cluster 1 patients were in the high-risk group (Fig. 6D; Fig. S10D), and most Basal-like subtype samples were in the high-risk group, consistent with its association with patient poor prognosis, validating the robustness and predictive effectiveness of our Riskscore models (Fig. S9C, D).

Epithelial and T Cell RiskScores as independent predictors of survival

In univariate Cox analysis, RiskScore_Epi (P < 0.001), RiskScore_T, (P < 0.001), age (P = 0.008), grade (P = 0.0038), T stage (P = 0.034), M stage (P = 0.018), and N stage (P = 0.0032) were significant prognostic factors (Fig. 6E; Fig. S11A). Multivariate Cox analysis identified RiskScore_Epi (HR 5.05, P < 0.001), tumor grade (HR 1.41, P = 0.0491), and N stage (HR 2.31, P = 0.0044) as independent prognostic factors (Fig. 6F). Similarly, RiskScore_T (HR 4.59, P = 0.0042), M stage (HR 0.49, P = 0.0101), and N stage (HR 2.2, P = 0.0075) were also independent predictors (Fig. S11B).

For clinical decision-making, we constructed two prognostic nomograms based on these results. One incorporating RiskScore_Epi, tumor grade, and N stage, and the other included RiskScore_T with M and N stages (Fig. 6H; Fig. S11E). Higher total points in these nomograms were associated with worse survival. Calibration curves showed close alignment between predicted and actual survival at 1, 2, and 3 years, confirmed models’ accuracy (Fig. 6I; Fig. S11C). DCA further demonstrated superior clinical effectiveness compared to other factors (Fig. 6J; Fig. S11D).

RiskScore-based stratification reveals distinct clinical features, genomic landscape, and pathway activation

In the TCGA cohort, the high-risk group identified by RiskScore_Epi exhibited more aggressive tumor characteristics, including a higher proportion of G2/3 histological grade, T2 stage, history of alcohol and smoking, and pancreas-related diseases (Fig. S12A–M). Similar trends were observed in the RiskScore_T model, where the high-risk group comprised older patients (≥ 65 years), G2/3 grade, clinical stage II, T3 stage, and R1/2 residual tumors (Fig. S13A–M).

The high-risk group based on RiskScore_T also had significantly higher mutation frequency and TMB compared to the low-risk group (Fig. S14B; mutation frequency: 95% vs 70%, P < 0.001; Fig. 7B; P_TMB < 0.01). KRAS and TP53 mutations were more frequent in the high-risk group (Fig. 7A; KRAS: 77% vs 46%, P < 0.001; TP53: 70% vs 50%, P < 0.05). A similar pattern was seen in RiskScore_Epi (Fig. S14A, C-D).

Fig. 7

Mutation landscape, pathways, and TME between RiskScore_T model risk groups in TCGA dataset. A, B Comparative analysis of top 5 genes mutation spectrum (A) and TMB distributions (B) in the RiskScore_T high- and low-risk groups, respectively. C Heatmap showing the significant enrichment hallmark pathways of the RiskScore_T high- and low-risk groups by GSVA (P < 0.01). Color bar = enrichment score. D Boxplot illustrating proportions of distinct cell types within the TME for the RiskScore_T high- and low-risk groups. E Correlation analysis between the RiskScore and 14 tumor-related pathways. F ESTMATE analysis revealed immune and stromal infiltration differences in the RiskScore_T high- and low-risk groups. G, I Distinctions in 28 and 64 immune cell infiltration levels in the RiskScore_T high- and low-risk groups, analyzed using the ssGSEA (G) and xCell (I) algorithm, respectively. H Distribution of immune subtypes across the RiskScore_T high- and low-risk groups. (*P < 0.05, ** P < 0.01, *** P < 0.0001). TCGA, The Cancer Genome Atlas Program; TME, tumor microenvironment; TMB, tumor mutation burden; GSVA, gene set variation analysis; ssGSEA, single-sample gene set enrichment analysis

The high-risk group showed enrichment in tumor-associated pathways. In RiskScore_Epi, this included classical tumor signaling, while the low-risk group was enriched in bile acid biosynthesis (Fig. S15A; P < 0.01). For RiskScore_T, the high-risk group showed changes in tight junctions, NOTCH, and P53 signaling, while the low-risk group was enriched in immunodeficiency and type I diabetes pathways (Fig. S15B; P < 0.01). Using Hallmark gene sets, the RiskScore_Epi high-risk group showed significant enrichment in pathways like KRAS signaling, P53 pathway, EMT, E2F targets, and hypoxia (Fig. S15C; P < 0.01). RiskScore_T high-risk group was enriched in glycolysis, hypoxia, and NOTCH signaling (Fig. 7C; P < 0.01).

Correlation analysis showed positive associations between both RiskScore_Epi and RiskScore_T with pro-tumor signals and negative correlations with normal pancreatic function (correlation > 0.3; Fig. S15E-F). Pathway analysis highlighted activation of hypoxia tolerance, cell migration, and proliferation pathways (such as EGFR, Hypoxia, MAPK, WNT) in high-risk groups (Fig. S15G-H; P < 0.05). Meanwhile, RiskScore_Epi and RiskScore_T have positive correlations with these pathways (correlation > 0.3; Fig. 7E).

RiskScore stratification reveals differences in TME and immune phenotypes between high and low-risk groups

The high-risk group defined by RiskScore_Epi had significantly higher proportions of epithelial cells and T cells, while the low-risk group showed elevated endothelial cell proportions (Fig. 7D; P < 0.05). Similar trends were observed for RiskScore_T (Fig. S15D; P < 0.05). Although the RiskScore_Epi high-risk group showed lower StromalScore, ImmuneScore, and ESTIMATEScore, the differences were not significant (Fig. S16A). In contrast, the RiskScore_T high-risk group had significantly higher ImmuneScore and ESTIMATEScore (Fig. 7F; P < 0.01).

RiskScore_Epi was positively correlated with most tumor-infiltrating immune cells (TIICs) (correlation > 0.3; Fig. S16E), with the high-risk group showing enriched levels of activated CD4 T cells, gamma delta T cells, memory B cells, natural killer T cells, and neutrophils (Fig. S16C-D; P < 0.01). The low-risk group, however, exhibited diminished immune activity, indicative of immune suppression (Fig. S16C-D; P < 0.01). RiskScore_T was negatively correlated with active B cells, activated CD8 T cells, and effector memory CD4 T cells (correlation > 0.3; Fig. S16E), and its high-risk group had lower TIIC abundance (Fig. 7G, I; P < 0.05). This suggests that the low-risk group patients may benefit from antitumor immune responses, enhance immunotherapy efficacy.

Regarding immune subtypes, both RiskScore_Epi and RiskScore_T high-risk groups had higher frequencies of C1 (wound healing) and C2 (IFN-γ dominant) subtypes but a lower prevalence of C3 (inflammatory) subtype (Fig. 7H; Fig. S16B). In the tumor–immunity cycle, RiskScore_T was negatively correlated with steps like priming and activation (step 3) and recognition of cancer cells by T cells (step 6). In contrast, RiskScore_Epi was positively correlated with the release of cancer antigens (step 1) and immune cell recruitment (step 4) (Fig. S16F). These findings suggest that RiskScore models offer insights into immune status and activity of the patient’s immune response.

Therapeutic efficacy prediction for high and low-risk groups based on RiskScore stratification

Differential expression analysis of immune checkpoints revealed higher expression levels in the RiskScore_Epi high-risk group (Fig. S17A; P < 0.05), with a positive correlation between RiskScore_Epi and CD274 expression (Pearson coefficient: r = 0.21, P = 0.0063; Fig. 8C), suggesting potential benefits from IC blockade therapy. However, the RiskScore_Epi high-risk group had higher TIDE and Exclusion scores (Fig. 8D-E; P < 0.01), indicating increased immune evasion and reduced sensitivity to ICB treatment. RiskScore_Epi correlated positively with the TIDE score (r = 0.29, P = 0.00012; Fig. 8F), with fewer patients in the high-risk group predicted to respond to immunotherapy (38% vs 65%, P = 0.00001; Fig. 8G). Patients in low-risk group showed higher IC50 values for Bortezomib (P = 8e−04), Dasatinib (P = 0.00039), Imatinib (P = 0.01), Paclitaxel (P = 8.8e−05), Sorafenib (P = 0.032), Sunitinib (P = 0.036), and Vinorelbine (P = 0.0015) (Fig. 8H-N).

Fig. 8

Analysis of immune checkpoint, immunotherapy response, and drug sensitivity between different risk groups in the TCGA dataset. A Comparisons of immune checkpoint gene expression in the RiskScore_T high- and low-risk groups. B Correlation analysis between the RiskScore_T and PDCD1 expression. C Correlation analysis between the RiskScore_Epi and CD274 expression. D, E Comparisons of TIDE score (D) and T cell Exclusion (E) in the RiskScore_Epi high- and low-risk groups. F Correlation analysis between the RiskScore_Epi and TIDE score. G Bar plot showing the difference in ICB response between the RiskScore_Epi high- and low-risk groups. H–N Box plots displaying the estimated IC50 of Bortezomib (H), Dasatinib (I), Imatinib (J), Paclitaxel (K), Sorafenib (L), Sunitinib (M), and Vinorelbine (N) in the RiskScore_Epi high- and low-risk groups. (*P < 0.05, ** P < 0.01, *** P < 0.0001). TCGA, The Cancer Genome Atlas Program; TIDE, Tumor Immune Dysfunction and Exclusion; ICB, immune checkpoint blockade; IC50, half-maximal inhibitory concentration

In contrast, most ICs were upregulated in the RiskScore_T low-risk group (Fig. 8A), with PDCD1 expression showing a negative correlation with RiskScore_T (Pearson coefficient: r = − 0.23, P = 0.0025; Fig. 8B), indicating potential benefit from IC blockade therapy. Dysfunction scores were higher in the RiskScore_T low-risk group (Fig. S17C; P < 0.05), but TIDE and Exclusion scores showed no significant differences between high- and low-risk groups (Fig. S17B, D). The RiskScore_T high-risk group demonstrated increased sensitivity to Dasatinib (P = 0.025), Paclitaxel (P = 0.0011), and Erlotinib (P = 0.003) (Fig. S17F-H).

Discussion

In this study, we have delineated the PDAC cell atlas, identified prognostic cell populations, and scrutinized intracellular heterogeneity. Epithelial cells, integral to tumor development, serve as the origin of malignant neoplastic cells [13, 14]. Type I ductal cells are relatively normal and present in both non-cancerous and neoplastic tissues, while type II ductal cells, expanded in PDAC, demonstrate malignancy [18]. Despite observing two populations of epithelial cells with varying prognostic correlations, subsequent analyses affirmed their malignant nature, albeit with differing degrees of malignancy. These findings collectively imply that, in aggregate, epithelial cells are associated with a poor prognosis.

Within the PDAC TME, various subtypes of T cells with diverse functions exist. The presence of Tregs contributes to a poor prognosis by suppressing the functions of effector T cells, impeding the response to immunotherapy in patients [81, 82]. T cell exhaustion plays a pivotal role in immune dysfunction and evasion among PDAC patients. Prolonged exposure of T cells to persistent antigens and inflammation in cancer patients leads to continuous stimulation, ultimately resulting in T cell exhaustion [86]. Tregs and depletion T cells gradually increase during tumor progression [16]. Our study exclusively identifies T cell populations associated with a poor prognosis. Subsequent analyses revealed a significant upregulation of depletion markers in T cells, indicating a progressive depletion of these cells during tumor development. Furthermore, the impaired antitumor function of T cells may be exacerbated by the presence of immunosuppressive cells within the TME, including Treg cells, granulocytes, and macrophages [87, 88].

Cellular deconvolution emerges as a potent tool for analyzing DNAm in bulk tissues. TCGA-PAAD samples, stratified into two clusters based on their DNA methylation profiles, exhibit significantly different prognoses. Consequently, we opted to deconvolute these two clusters and compared the proportions of cell types between patients in these clusters, and we observe that cell populations associated with poorer prognosis are more abundant in tumors from patients with shorter survival times. This trend is further confirmed in the eight samples we collected and monitored for survival time, depicting a gradual decline in patient survival time concomitant with an increase in the proportions of T cells and epithelial cells in tumor tissues.

Prior studies identified several methylation markers for PDAC prognosis in bulk tissue [89,90,91,92]. However, tracing the impact of these markers on the functional aspects of tumor cell populations poses a significant challenge. In this study, we obtained cell-type-specific sites for T cells and epithelial cells. These markers not only denote a significant change in methylation levels but also signify an increase in specific cell types. As anticipated, prognostic models were developed based on these markers, further confirming their role as independent prognostic factors. We also filtered several clinical factors and constructed nomograms in conjunction with our prognostic models to predict survival time for each patient. Given the varied immune levels and therapeutic responses among different risk groups, these predictive models could serve as valuable tools for predicting PDAC patients who would benefit more from drug and immunotherapy. This enables the stratification of individuals and customization of treatment strategies, enhancing treatment efficacy and improving patient outcomes.

Here we attempted to transpose the gene expression signatures of prognosis-associated cell subpopulations into DNA methylation profiles. Our study represents an effort in the clinical application of scRNA-seq as we translated expensive single-cell technologies into more cost-effective DNA methylation assays. This facilitates the clinical translation of single-cell technologies and extends their utility to large clinical cohort populations.

Regrettably, the deconvolution method fell short of achieving the desired precision in prognostic cell subpopulations, impeding in-depth validation to dissect the relationship between subpopulations within a cell type and prognosis. Nonetheless, there are additional limitations in our study. Firstly, this study predominantly relied on retrospective data, necessitating further prospective studies. Additionally, the available sample size was relatively modest, warranting larger-scale studies to validate and refine our conclusions. Secondly, given the inherent nature of bioinformatics analysis, future experimental endeavors are needed to explore the molecular mechanisms and biological functions of these prognostic methylation sites in PDAC progression.

Conclusion

In this study, the integration of scRNA-seq data and DNAm data facilitated the development of predictive models based on epithelial cell and T cell-specific DNA methylation sites, enabling effective prognosis prediction for PDAC patients. Patient stratification using our models revealed distinctive prognoses and immune microenvironments, suggesting potential benefits for immunotherapy and drug therapy tailored to each patient. Notably, the construction of nomograms provides valuable tools for clinical treatment decisions. The implementation of our prognostic models offers clinicians a cost-effective and accurate method for predicting outcomes in PDAC.