Our study applied a comprehensive analysis of data publicly available online. The scRNA-seq dataset (PRJCA001063) containing 24 PDAC tumor samples and 11 control pancreases without any treatment was obtained from the Zenodo database (www.zenodo.org) (Peng et al., 2019). Bulk sequencing data in the training and validation sets were obtained from The Cancer Genome Atlas (TCGA) (TCGA- PAAD) (https://portal.gdc.cancer, updated until 07-20-2019), the Gene Expression Omnibus (GEO) (GSE71729, GSE21501) (https://www.ncbi.nlm.nih.gov/geo/query), the European Molecular Biology Laboratory (EMBL-EBI) (E-MTAB-6134) (https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-6134/) and the International Cancer Genome Consortium Data Portal (ICGC) Canada pancreatic cancer project (PACA-CA) (ICGC-CA, https://dcc.icgc.org/releases/current/Projects/PACA-CA). Mutation data in TCGA cases were downloaded from the Broad Institute TCGA Genome Data Analysis Center (http://gdac.broadinstitute.org). Transcriptomic and clinical information for each dataset was simultaneously downloaded from the respective websites when available. Only PDAC was included in the subsequent analysis, while other histological subtypes, such as Acinar Cell Carcinoma, adenosquamous carcinoma, mucinous cystadenocarcinoma, mixed ductal endocrine carcinoma were excluded.

We used the “Seurat” package deployed in R for quality control and downstream analysis of scRNA-seq data (Hao et al., 2021). The regularizing negative binomial regression was used to eliminate batch effects. Low quality cells (<200 genes/cell, <3 cells/gene and >10% mitochondrial genes) were excluded. Afterward, we calculated the standardized variance of each gene across cells to generate highly variable genes and used them for principal component analysis (PCA). “ElbowPlot” analysis and heat map visualization in “Seurat” were used to identify significant principal components (PCs). Based on PC1 to PC13, graph-based clustering was applied (res = 0.8) to identify different cell groups. Then, non-linear dimensionality reduction was performed using the “tSNE” method. Different cell clusters were identified and annotated with the “singleR” package, the CellMarker database, and previously published scRNA-seq analysis (Peng et al., 2019), (Aran et al., 2019; Zhang et al., 2019; Schlesinger et al., 2020). To identify differentially expressed genes (DEGs) between ductal cells in normal and tumor tissues, we used the “FindMarkers” function in the “Seurat” package. DEGs were filtered by |log2 (fold change) | > 0.5 and p < 0.05.

Firstly, a univariate Cox regression analysis was performed on the TCGA-PDAC cohort to determine the association between the aforementioned duct-cell-related DEGs and OS. Then, significant DEGs related to OS (p < 0.05) were included in the least absolute shrinkage and selection operator (LASSO) regression analysis to reduce multicollinearity and make the model simpler and more effective. Finally, a multivariate Cox regression analysis was conducted on the screened genes to assess the impact of each gene as an independent prognostic factor on patient survival. The risk score formula was constructed as described in the previous study (Ma et al., 2018). We used the median risk score as a cut-off value to categorize the patients in the TCGA training set into high-risk and low-risk groups. Patients in the GEO validation group were also divided into two groups by the same method. Survival differences between the two groups were assessed by the Kaplan-Meier method and log-rank test. Furthermore, we plotted the time-dependent receiver operating characteristic (ROC) curves with 1, 3, and 5 years as the defined points and calculated the corresponding area under the ROC curve (AUC) to assess and compare the predictive power of the risk model with other recently published PADC risk models. The values of AUC range from 0.5 to 1, with 1 indicating full discriminant and 0.5 indicating no discriminant.

Several clinical characteristics of PDAC patients might affect prognosis. To investigate whether the risk score was independent of relevant clinical factors (e.g., age, alcohol history, gender, chronic pancreatitis history, TNM stage, histologic grade, etc.), we performed univariate and multivariate Cox regression analyses. The significance level was set at p < 0.05.

To clarify the differences in signaling pathways and molecular mechanisms associated with gene expression profiles between high-risk and low-risk groups, we analyzed the gene expression data using GSEA analysis. Enriched gene sets with p < 0.05 and FDR <0.25 were considered statistically significant.

TMB was defined as the total number of mutations (changes) found in the DNA of cancer cells (Chan et al., 2019). The TMB of patients with mutation data in the TCGA cohort was calculated using the “tmb” function in the “MAFTOOLS” package (Mayakonda et al., 2018). Somatic mutation comparison was also completed with the same package for patients in the high- and low-risk groups. p < 0.05 was set as the level of significance.

The immune score for each patient in the TCGA cohort was calculated using the ESTIMATE algorithm. The standardized expression profile was then uploaded to the TIMER website (http://timer.cistrome.org/) to obtain an assessment of immune cell infiltration (including B cells, macrophage, myeloid dendritic cells, neutrophil, CD4⁺ T cells, and CD8⁺ T cells) (Li et al., 2020). We also applied the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm to predict the patient’s response to immunotherapy (Fu et al., 2020). At a significance level of p < 0.05, we compared immune scores, degree of immune cell infiltration, and response to immunotherapy between the high-risk and low-risk groups.

Stemness features were extracted from transcriptomic data of TCGA PDAC patients by using an innovative one-class logistic regression (OCLR) machine-learning algorithm (Malta et al., 2018). The transcriptome-based stemness index (mRNAsi) was then mapped to a range of 0–1 by linear transformation (subtracting the minimum and dividing by the maximum). To predict which compounds were likely to be effective against PDAC, we used Broad institute’s Connectivity Map (CMap) to screen drug candidates based on key genes of the risk model (Subramanian et al., 2017).

To investigate changes of ductal cells in PDAC, we reanalyzed the scRNA-seq dataset and annotated the cell types according to their transcriptomic characteristics. After performing cell quality control (as described in the methods section), we performed principal component analysis using the top 2000 variable genes (Figure 1A) and identified 13 principal components for downstream analysis via elbow plots (Figure 1B). Then these cells were divided into 33 clusters (Supplementary Figure S1) and based on marker gene expression (Supplementary Figure S1, Supplementary Tables S1, S2) these clusters were classified into acinar cells, ductal cells, endothelial cells, endocrine cells, fibroblasts, stellate cells, T cells, B cells, and macrophages (Figures 1C,D). By comparing the transcriptional data of normal and malignant pancreatic ductal cells, we identified 881 differentially expressed genes (log2FC >0.5, p < 0.05), of which 475 were upregulated and 406 were downregulated (Figure 1E, SupplementaryTable S3).

Among these DEGs, we used univariate Cox regression to initially screen out 98 genes associated with PDAC prognosis (p < 0.05) in the TCGA cohort. Then, these genes were included in LASSO regression analysis to further exclude confounding factors (Supplementary Figures S3A,B). Next, we performed a multivariate Cox regression analysis and select 10 genes to construct a prognostic prediction model (Supplementary Figure S3C, Supplementary Table S4). According to this model, we calculated the risk scores of 147 patients in the TCGA PDAC group. Using the median risk score as the cutoff point, all patients in the set were divided into a high-risk group and a low-risk group (Figure 2A). The median OS survival was 313 and 505 days for the high-risk and low-risk groups, respectively. Kaplan-Meier analysis showed that high-risk PDAC patients had significantly lower OS than low-risk PDAC patients (p < 0.001, Figure 2B). Then, the AUCs were calculated to assess the OS prediction efficiency of the risk model. It had AUCs of 0.825, 0.819, and 0.824 at 1, 3, and 5 years, indicating that this model had favorable predictive power (Figure 2C).

Furthermore, to clarify the reliability of the risk model, we divided the PDAC patients in the TCGA dataset into subgroups with different clinical characteristics (including age, gender, histologic grade, AJCC stage, and tumor location). Kaplan-Meier analysis showed that even among the different subgroups, patients in the high-risk group had significantly lower OS than those in the low-risk group (p < 0.05, Figure 3), which further demonstrated the superiority of this model.

Using two GEO datasets (GSE71729 and GSE21501), E-MTAB-6134 and ICGC-CA, we validated the reliability and stability of the ductal cell-related risk model. Risk scores were calculated in both cohorts and PDAC patients were then divided into high-risk and low-risk groups (Figures 4A,B, Supplementary Figures S5A,B). Consistent with the results of the training set, the OS of high-risk patients in both validation sets was significantly shorter than that of low-risk patients (p < 0.001, Figures 4C,D, Supplementary Figures S5C,D). In GSE71729, the median OS of high-risk patients and low-risk patients were 300 and 540 days, while they were 330 and 540 days in GSE21501, 483 and 757 days in E-MTAB-6134, and 445 and 590 days in ICGC-CA, respectively. The results of the AUC analysis also showed that the model had favorable predictive power in the validation sets. In GSE71729, the AUCs were 0.689, 0.829, and 0.954 for 1, 3, and 5 years, while they were 0.749, 0.782, and 0.788 for 1, 3, and 5 years in GSE21501, 0.713, 0.665, and 0.686 for 1, 3, and 5 years in E-MTAB-6134, and 0.675, 0.705, and 0.749 for 1, 3, and 5 years in ICGC-CA, respectively (Figures 4E,F, Supplementary Figure S6).

To further validate the predictive power of the ductal cell risk model, we compared it with three recently published pancreatic cancer risk scoring models, including the Deng model (Deng et al., 2021), the Qiu model (Qiu et al., 2020), and the Wu model (Wu et al., 2019). Notably, these three models were all built on bulk sequencing data. Few pancreatic cancer risk models based on single cell data have been reported. We compared the predictive power of each model by the area under the ROC curve (AUC). In general, the ductal cell model based on single cell RNA sequencing possessed better predictive power in all four validation datasets, especially in terms of 5-years overall survival prediction (see Supplementary Figure S6).

A univariate Cox regression analysis was performed to explore the relationship between patient clinical features and PDAC prognosis in the TCGA cohort. The results showed that only the positive nodes rate (HR, 5.453; 95% CI, 2.051 to 14.500; p = 0.001), as well as the maximum tumor dimension (HR, 1.188; 95% CI, 1.102 to1.394; p = 0.035), were significantly associated with the prognosis of TCGA-PDAC patients (Figure 5A). Then through ROC analysis, we calculated the AUCs of risk score, positive nodes rate, and maximum tumor dimension at 1, 3, and 5 years (Figures 5B–D). The results demonstrated that the predictive efficacy of the risk score was better than that of positive nodes rate and maximum tumor dimension. The risk score also acted as an independent risk factor for the PDAC prognosis in a multivariate Cox regression analysis (p < 0.001) (Figure 5E).

GSEA analysis was conducted to annotate the function of DEGs between high-risk and low-risk patient groups. 9 cancer-related gene sets were demonstrated to be significantly enriched in the high-risk patient group (Nominal p-value < 0.05, FDR <0.25), including mTORC1 signaling, MYC targets v1, MYC targets v2, G2M checkpoint, E2F targets, mitotic spindle, glycolysis, DNA repair, and unfolded protein response (Figure 6). The gene sets were found to be intimately involved in tumorigenesis, DNA repair, genome stability, and tumor nutrition and metabolism (Stine et al., 2015; Kent and Leone, 2019; Murugan, 2019; Yang et al., 2020). The result might provide clues to the potential mechanisms affecting the prognosis of PDAC patients.

The development of cancer is the result of somatic mutation and clonal selection. To explore the relationship between risk score and mutation status, we further analyzed and visualized somatic mutations between the high-risk and low-risk cohorts of TCGA. The results indicated that the top 20 genes with the highest mutation frequencies in both groups were similar and only differed in the ranking (Figures 7A,B). The genes with mutation rates above 10% included KRAS, TP53, SMAD4, CDKN2A, and TNN. Missense mutations were the most common type. In general, the high-risk cohort harbored more somatic mutations and showed a higher tumor mutation burden than the low-risk cohort (Figure 7C, p < 0.05). When comparing mutated genes in the two groups, only KRAS (82 vs. 59%, p < 0.01) and ADAMTS12 (8 vs. 0%, p < 0.05) were statistically different in mutation frequency (Figures 7D,E).

Different tumor mutation burdens might suggest different immune cell infiltration statuses. To investigate whether patients in the high-risk group with higher TMB had more immune cell infiltration, we evaluated the immune microenvironment in both groups by CIBERSORT and TIMER algorithms (Newman et al., 2015; Li et al., 2020). The results suggested that there was no statistical difference in immune scores between the two groups (Supplementary Figure S7A). In parallel with this result, the TIMER assessment also showed that only CD4⁺ T cell infiltration might differ among the six immune cell types (Supplementary Figure S7B). The TIDE algorithm was exploited to predict responsiveness to immunotherapy and there was no statistical difference between the two groups (Supplementary Figure S7C), which further suggested that the two groups might have similar immune cell infiltration statuses. Tumor stemness was considered to be one of the key factors in tumor progression and was significantly associated with patient prognosis. In our investigation, we obtained the mRNA expression-based stemness index through one-class logistic regression (OCLR) machine learning algorithm. The high-risk group with a poorer prognosis had a higher stemness index (p < 0.05) (Figure 7F).

The CMap database was used to identify potential drugs that might target DEGs between high risk and low risk patients in TCGA cohort (see Supplementary Table S5). As a result, 120 drug candidates were screened for 58 possible mechanisms of action (see Figure 8 and Supplementary Table S6). The top five drugs included AT-9283 (count = 7), sorafenib (count = 6), sunitinib (count = 6), dovitinib (count = 5), and tozasertib (count = 4). The top five possible mechanisms of action were Acetylcholine receptor antagonist (count = 25), FLT3 inhibitor (count = 14), Dopamine receptor antagonist (count = 13), Acetylcholine receptor agonist (count = 9) and Serotonin receptor antagonist (count = 8).