Gaofei Zhang1
Congya Zhou- 1Department of Radiation Oncology, Shaanxi Provincial People’s Hospital, Xi’an, Shaanxi, China
- 2Department of Oncology, Shaanxi Provincial People’s Hospital, Xi’an, Shaanxi, China
Background: This study focused on epithelial cells to construct a prognostic risk model and provide targeted insights into responses to immunotherapy.
Methods: Single-cell RNA sequencing (scRNA-seq) was clustered using Uniform Manifold Approximation and Projection (UMAP) and a risk model was developed through Least Absolute Shrinkage and Selection Operator (LASSO) regression analysis. Kaplan-Meier analysis was performed to evaluate the prognosis of PAAD. The biological characteristics of LIPH were assessed using CCK-8, colony formation and Transwell assays.
Results: Eight major cell clusters were identified, revealing two developmental trajectories for malignant epithelial cells from primary to metastases. Epithelial cells were categorized into Scissor+ and Scissor- subtypes, with Scissor+ epithelial cells exhibiting more complex cellular communication with TME cells. Furthermore, we successfully developed a risk model for PAAD patients based on the Scissor findings. The prognosis for PAAD patients in the high-risk group was significantly poorer within both the TCGA and ICGC cohorts. Differences were observed in the populations of naïve B cells, CD8 T cells, M0 macrophages, and activated dendritic cells in different groups. Knockdown of LIPH significantly inhibited the growth and invasion of PAAD cells.
Conclusion: These findings underscore the significance of this risk model in predicting prognosis and immunotherapy responses, and enhancing understanding of tumor microenvironment (TME) heterogeneity in PAAD metastases.
1 Introduction
Pancreatic adenocarcinoma (PAAD) is among the deadliest human malignancies and is projected to become the third leading cause of cancer-related mortality by 2025 (Siegel et al., 2024; Bray et al., 2024; Stoop et al., 2025). The majority of PAAD patients have been diagnosed with distant metastases, and approximately 80% ultimately succumb to these metastases (Park et al., 2021). Although chemotherapies such as NALIRIFOX, a modified FOLFIRINOX, and gemcitabine combined with nanoparticle albumin-bound paclitaxel have significantly improved overall survival (OS) for metastatic PAAD patients, drug resistance typically develops within 6 months. Furthermore, the median OS for metastatic PAAD patients is less than 12 months (Von Hoff et al., 2013; Hu et al., 2024; Conroy et al., 2023). A critical factor contributing to this situation is the complexity and heterogeneity of cancer cells and the TME. Therefore, there is an urgent need for comprehensive exploration of the heterogeneity of both cancer cells and cancer-related cells within the tumor microenvironment (TME) to inform therapeutic decisions for patients with metastatic PAAD.
The TME of PAAD is a highly heterogeneous system composed of tumor cells, endothelial cells, cancer-associated fibroblasts (CAFs), immune cells, and the extracellular matrix (Han et al., 2021; Grunwald et al., 2021). The TME of PAAD plays a crucial role in tumor initiation, metastasis, metabolic reprogramming, immune evasion, and drug resistance (Raghavan et al., 2021). With advancements in high-throughput sequencing technologies, single-cell RNA sequencing (scRNA-seq) has emerged as a powerful tool for elucidating the heterogeneity of the TME in PAAD, allowing for the exploration of distinct cell subsets and the dissection of developmental trajectories among various subpopulations at the single-cell level (Fan et al., 2022; Moncada et al., 2020; Schalck et al., 2022). Lin et al. utilized scRNA-seq to identify the heterogeneity of cellular compositions in primary and metastatic PAAD tumors (Lin et al., 2020). Additionally, Zhang et al. revealed that subpopulations of CD103+PD-1+CD39+ T cells and CD73+ macrophages could promote liver metastasis in PAAD (Zhang Z. et al., 2023). PAAD originates from epithelial cells, and Kim et al. identified a novel epithelial cell cluster, demonstrating its prognostic value during cancer progression (Kim et al., 2024). However, the heterogeneity of cancer cell subpopulations, their characteristic gene expression, and their interactions with the TME in PAAD metastases remain poorly understood. Recently, the Scissor algorithm was reported to identify phenotype-associated cell subpopulations from scRNA-seq and bulk RNA-seq data (Sun et al., 2022). This advancement will aid in the comprehensive understanding of metastasis-related subclusters and facilitate the exploration of therapeutic targets for patients with metastatic PAAD.
In this research, we investigated the potential heterogeneity of primary and metastatic PAAD samples at the single-cell level by integrating scRNA-seq and bulk RNA sequencing data. Furthermore, we identified characteristic malignant epithelial cell subsets, referred to as Scissor+ epithelial cell subsets, and constructed a risk model incorporating candidate genes (TATDN1, CAV2, CLDN1, LIPH, MT1E, PSCA, and MMP3). We evaluated the prognostic and immunological implications of this risk model. Notably, we observed that LIPH expression was elevated in tumor tissues of PAAD patients compared to paired adjacent tissues. Additionally, our results indicated a significant reduction in both the proliferation and metastasis of PAAD cells following LIPH knockdown. These findings enhance our understanding of the TME heterogeneity in PAAD metastasis and provide potential prognostic and immunological signatures for PAAD patients.
2 Materials and methods
2.1 Data acquisition
ScRNA-seq data, comprising 10 primary and six metastatic PAAD samples, were downloaded from the Gene Expression Omnibus (GEO) database, specifically from series number GSE154778 (Lin et al., 2020). Additionally, bulk RNA sequencing gene expression matrices for PAAD, along with corresponding clinical information such as age, gender, and tumor stage, were retrieved from The Cancer Genome Atlas (TCGA) and the International Cancer Genome Consortium (ICGC). The RNA expression data were transformed into transcripts per million (TPM) and subsequently log-normalized. Information regarding single nucleotide mutations, including small insertions and deletions (indels), was acquired via the cBioPortal website (https://www.cbioportal.org/). Furthermore, RNA-seq data for non-small cell lung carcinoma patients treated with anti-PD-1/PD-L1 therapy were obtained from the GEO database under series number GSE135222 (Kim et al., 2020). RNA-seq data for clear cell renal cell carcinoma patients undergoing PD-1 blockade were sourced from a previously published study (Braun et al., 2020).
2.2 ScRNA-seq data processing and clustering
The scRNA-seq data were imported into the R package Seurat (Satija et al., 2015). Genes expressed in fewer than 3 cells, cells expressing fewer than 200 genes, and cells containing more than 20,000 or fewer than 1,000 unique molecular identifiers were excluded. Cells with more than 10% mitochondrial gene expression were filtered out. The R package DoubletFinder was utilized to remove potential doublets (McGinnis et al., 2019). We normalized quality-controlled data using the NormalizeData function with default parameters. Based on the natural-log transformed normalized gene expression matrix, 2,000 highly variable genes were identified using the FindVariableFeatures function with the variance stabilizing transformation (VST) method. All genes were scaled using the ScaleData function. Principal Component Analysis (PCA) was performed to reduce dimensionality using the first 2,000 highly variable genes. The first 30 dimensions were employed to cluster cells using the FindNeighbors and FindClusters functions (resolution = 0.8). The Uniform Manifold Approximation and Projection (UMAP) method was applied to further reduce dimensionality based on the top 30 principal components. Cell types were annotated according to their canonical markers.
2.3 Identifying malignant epithelial cells
Epithelial cells were extracted both copyKAT (Gao et al., 2021) and inferCNV (Patel et al., 2014) methods were used to identify malignant epithelial cells. Epithelial cell subclusters annotated as diploid by copyKAT and having relative lower copy number variations (CNVs) score calculated by inferCNV were regarded as normal epithelial cells and excluded for subsequent analysis. ScRNA-seq data processing and clustering.
2.4 Trajectory inference analysis
To elucidate the trajectory of malignant epithelial cells from primary to metastatic tissues, we first identified differentially expressed genes (DEGs) between the two groups using the FindMarkers function. These DEGs were then provided to the R package Monocle2 for dimensionality reduction via the reduce Dimension function (Qiu et al., 2017). Genes that exhibited changes along the pseudotime were measured and visualized using a heatmap, and these genes were subsequently clustered into subgroups based on their expression patterns.
2.5 Scissor analysis
To identify malignant epithelial cells associated with metastasis in scRNA-seq data, we collected clinical information from patients in TCGA-PAAD, categorizing those diagnosed with Stage IV as metastatic. The R package Scissor was employed to reference corresponding bulk assays and estimate cell subpopulations with the highest metastatic potential. Significant DEGs were identified by comparing the populations of Scissor+ and Scissor- cells using the FindMarkers function, applying a screening threshold of p < 0.05 and log2FC >0.25. Additionally, functional enrichment analysis of the DEGs was conducted using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway annotations through the R package clusterProfiler. ScRNA-seq data processing and clustering.
2.6 Cell-cell communication analysis
For comparison of the cell-cell communication between Scissor+, Scissor- malignant cells and other cell types, the R package SingleCellSignalR was used to infer the ligand-receptor (LR) interactions (Cabello-Aguilar et al., 2020). The cut off for confident predicted LR interactions was set as 0.5.
2.7 Construction and validation of prognostic Scissor+ risk score model
We first performed univariate Cox regression analysis implemented in R package survival to determine the prognostic value of upregulated DEGs in Scissor+ cells. Unfavorable genes were determined based on the criteria hazard ratio (HR) >1 and p < 0.05. Next, we performed least absolute shrinkage and selection operator (LASSO) penalized cox regression algorithm with ten-fold cross-validation implemented in R package glmnet to fit the scissor + risk score model for predicting OS of patients. We established the risk score model with the following formula: risk score = sum (normalized expression of each significant gene multiplies its coefficient). Subsequently, we used the surv_cutpoint function implement in R package survminer to determine the optimal cut-off and divided patients into High and Low risk groups. For validation, four types of patients’ survival information, namely, OS, disease specific survival (DSS), disease free interval (DFI), and progression free interval (PFI) were used to explore the potential prognostic value of risk score model.
Additionally, we performed the same scoring system to the ICGC dataset and examine the prognostic value of the model. Besides, we used R package timeROC for time-dependent receiver operating characteristic (ROC) analysis of 1 and 3-year survival, and quantifying the area under the curve.
2.8 Nomogram construction and evaluation
Univariate and multivariate Cox regression analyses, along with a nomogram incorporating patient variables such as age, gender, tumor stage, and risk score, were employed to comprehensively evaluate the 1-, 3- and 5-year survival outcomes of patients. The accuracy of the nomogram was assessed using a calibration curve based on the Hosmer-Lemeshow test. Additionally, the predictive ability and clinical benefit of the nomogram were determined through decision curve analysis.
2.9 Tumor microenvironment estimation
The immune score, stromal score, tumor purity, and ESTIMATE score of individual patients were estimated using R package estimate. Additionally, the proportion of 22 common immune cell types were predicted by conventional marker genes with the help of R package CIBERSORT. TIDE prediction score was obtained from online website (http://tide.dfci.harvard.edu/).
2.10 Somatic mutation analysis
R package maftools was utilized to evaluate somatic variant data, including frame shift insertion, frame shift deletion, missense mutation, nonsense mutation, in frame deletion, in frame insertion, splice site, and multi-hit between high and low risk score groups (https://pubmed.ncbi.nlm.nih.gov/30341162/). The tumor mutational burden (TMB) for each patient was calculated as total mutations divided total covered bases*106.
2.11 Cell culture and gene knockdown assay
PANC-1, BXPC-3, ASPC-1, CFPAC-1, CAPAN-1, CAPAN-2 and normal human pancreatic duct epithelial (HPDE6-C7) cell lines were acquired form the BeNa Culture Collection (Bejing, China). The HEK293T cells were obtained from the Cell Repository of the Chinese Academy of Sciences (Shanghai, China). All of cells were cultured at 37 °C with 5% CO2 in DMEM, RPMI 1640 or IMDM medium adding with 10% serum, 100 U/mL penicillin and 100 mg/mL streptomycin (purchased from Beyotime, China).
Gene knockdown assay which included vector construction, virus packaging and lentiviral infection was constructed as previously described (Hu et al., 2019). The target sequence for shLIPH#1 was GCCCACATATCTGGGTTTGTT, and the target sequence for shLIPH#2 was GCGTCCTATGGATGTCACATT.
2.12 qRT-PCR assay
Total RNA from different cells was isolated with Cell Total RNA Isolation Kit (FORGENE, Cat# RE-03111) according to the protocol. The extracted RNA was reversed transcription as complementary DNA (cDNA) with PrimeScript™ RT Master Mix Kit (TAKARA, Cat# RR036A). Then the RT-qPCR was performed using Bio-Rad CFX maestro system. The relative mRNA expression was calculated using the 2−ΔΔCT method. The GAPDH gene was a control gene. LIPH forward primer (5′to 3′): CAACGGGAAACCTCACCAAGAC, Reverse primer (3′to 5′): AGCCAGGTTGATCCAATCCTCC. GAPDH forward primer (5′to 3′): GTCTCCTCTGACTTCAACAGCG, reverse primer (3′to 5′): AC-CACCCTGTTGCTGTAGCCAA.
2.13 CCK-8 and colony formation assay
To assess cell viability following LIPH knockdown, treated cells were seeded into 96-well plates at a density of 2000 cells per well. Subsequently, the CCK-8 reagent was added to each well, and the cells were incubated for 2 h at room temperature. Absorbance at 450 nm was measured using a multifunction microplate reader (BioTeK) at time points of 0, 24, 48, 72, and 96 h. For the colony formation assay, treated cells were plated in 24-well plates at a density of 1,000 cells per well. After a 14-day culture period, the cells were fixed with formaldehyde (Sigma-Aldrich) and subsequently stained with crystal violet (Sigma-Aldrich) for 30 min. The cell colonies were counted after being washed with deionized water.
2.14 Transwell assay
Transwell assay was conducted to assess the invasive and migratory effect of the cells. Treated cells with serum-free medium were placed in the upper chamber, while the lower chamber contained 600 μL of medium supplemented with 10% FBS. For invasion assay, Matrigel (BD Biosciences, USA) was added to the upper chambers. After incubated with 48 h, the cells in the upper chamber were removed, and cells mi-grated to the lower chamber were fixed with paraformaldehyde and stained with crystal violet. Subsequently, the cells were photographed for analysis.
2.15 Statistical analysis
Statistical analyses were performed using R version 4.0.2 (https://www.r-project.org/). The Wilcoxon rank sum test was employed to assess differences in expression between the two groups. The Log-rank test and Kaplan-Meier curves were utilized to evaluate the survival distributions of PAAD patients. A two-tailed Student’s t-test and one-way ANOVA were conducted to assess differences among the various groups. All data analyses were performed using GraphPad Prism version 6.0, with a significance threshold set at P < 0.05. Results were visualized with R package ggplot2.
3 Results
3.1 Clustering scRNA-seq data in PAAD patients
To comprehensively explore cellular heterogeneity in primary and metastatic samples of PAAD, the scRNA-seq dataset (GSE154778) was downloaded and analyzed. Following normalization and quality control procedures, a total of 11,420 cells and 23,759 genes were obtained from six metastatic samples and 10 primary samples (Figure 1A; Supplementary Figure S1A). Subsequently, twenty-seven clusters were identified through UMAP analysis, with each cluster exhibiting its own uniquely expressed genes (Figures 1B,C). Additionally, t-distributed Stochastic Neighbor Embedding (t-SNE) analysis was also performed to identify cell clusters (Supplementary Figure S1B-E). According to the marker genes, 27 cell clusters were annotated into eight major cell types, including epithelial cells (EPCAM and KRT19), macrophage (AIF1, C1QA, CXCL2, APOE, CD68, and FCER1G), T cells (CD3D and CD3E), dendritic cells (CD1E, CD1C, FCER1A, AIF1, C1QA, APOE, CD68 and FCER1G), endothelial cells (KDR, COL1A1, COL1A2, AOL3A1 and VWF), fibroblast (CLU, LUM, ACTN2, DCN, COL1A1, COL1A2 and AOL3A1), mast cells (CLU, CD69, FCER1A, CXCL2, APOE, FCER1G and KRT19) and mixed cells (Figures 1D,E). The epithelial cells, immune cells and matrix cells were marked by EPCAM, PTPRC and DCN, respectively (Figure 1F; Supplementary Figure S1F).
Figure 1. Cell clusters in primary and metastatic samples from PAAD patients. (A) UMAP visualization of clusters from various PAAD patients. (B) UMAP visualization of PAAD cells differentiated by Seurat cluster coloration. (C) Visualization of differentially expressed genes across 27 cell clusters. (D) UMAP visualization of eight distinct cell types in PAAD patients. (E) Visualization of differentially expressed genes in these eight cell types. (F) UMAP visualization of marker genes in epithelial cells, immune cells, and matrix cell types. PAAD, Pancreatic adenocarcinoma; UMAP, Unified Flowform Approximation and Projection.
3.2 Identification of epithelial cells in PAAD
Since ductal epithelial cells are a population of cells originating from PAAD, we proceeded to further analyze the heterogeneity and potential function of epithelial cells in PAAD. Following UMAP and t-SNE analysis, four major cell types were identified across 14 clusters (Figure 2A; Supplementary Figure S2A). The aneuploidy and diploidy analysis revealed that clusters 0, 2, and 15 were classified as diploid, whereas the remaining clusters were identified as aneuploid (Figure 2B; Supplementary Figure S2B). Furthermore, inferCNV analysis was performed to assess the heterogeneity of epithelial cells. As illustrated in Figure 2C; Supplementary Figure S3A, clusters 17, 23, 11, 3, 21, and 5 exhibited a higher number of CNVs, whereas clusters 0, 2, and 15 showed lower CNVs, corroborating the findings from the aneuploidy analysis. Additionally, we analyzed the developmental trajectories of malignant epithelial cells in both primary and metastatic samples. Pseudtotime analysis, along with cell cluster and state estimation analyses, were employed to investigate the developmental trajectories of these malignant epithelial cells (Figures 2D–F). The results indicated that malignant epithelial cells progress from primary tumors to metastases, suggesting an underlying biological process driving PAAD metastasis.
Figure 2. Identification of epithelial cells in PAAD. (A) UMAP visualization of epithelial cell clusters in PAAD. (B) UMAP visualization of epithelial cell clusters with aneuploid analysis. (C) Identification of malignant epithelial cells based on inferred CNVs from various samples, revealing that non-malignant epithelial cells exhibited no significant CNVs. (D–F) The differential trajectories of epithelial cells during the progression from primary tumors to metastasis are il-lustrated according to pseudotime, cell type, and cellular state, respectively. (G) The proportion of malignant epithelial cells throughout the transition from primary tumors to metastasis. (H) CDH1 expression levels in malignant epithelial cells during the progression from primary tumors to metastasis. (I) CDH2 expression levels in malignant epithelial cells throughout the evolution from primary tumors to metastasis. PAAD, Pancreatic adenocarcinoma; UMAP, Unified Flowform Approximation and Projection; CNV, Copy number variant.
As the progression from primary tumors to metastasis occurs, the fraction of malignant epithelial cells initially decreases before subsequently increasing. Notably, two distinct developmental trajectories for malignant epithelial cells from primary tumors to metastases have been observed (Figure 2G; Supplementary Figure S3B). We analyzed the expression of the epithelial cell marker gene CDH1 and the mesenchymal cell marker gene CDH2, finding that CDH1 expression decreased while CDH2 expression increased in malignant tumor epithelial cells as the tumor progressed from primary to metastatic stages (Figures 2H,I). These findings are consistent with previous studies and support the reliability of the malignant epithelial cell subsets we identified. Overall, we identified malignant epithelial cells and characterized two developmental trajectories from primary tumors to metastasis. Differential gene analysis revealed that DCBLD2, KLK6, RBP1, LY6D, and RNF43 were expressed at higher levels in metastatic malignant epithelial cells, whereas TFF2, TFF1, CTSE, VSIG2, and CLDN18 were more abundant in primary malignant epithelial cells (Supplementary Figure S3C). Based on the expression of the top 1,000 differentially expressed genes, malignant epithelial cells can be categorized into six distinct clusters (Supplementary Figure S3D). As tumors progress from primary to metastatic stages, clusters 2, 1, and 4 exhibit a decrease, while clusters 5, 3, and 6 show an increase. GO analysis indicated that the top pathways associated with clusters 2, 1, and 4 included oxidoreduction-driven activity, electron transfer activity, molecular function inhibitor activity, viral life cycle, cytoplasmic translation, and regulation of nitric oxide metabolism. Conversely, the top pathways in clusters 5, 3, and 6 encompassed actomyosin, actin filament bundle formation, RNA splicing via transesterification reactions, cytoplasmic translation, cell-substrate junctions, and ribosomal function (Supplementary Figure S3E).
3.3 Characteristics and intercellular communication of Scissor+ epithelial cells in metastases
To identify subclusters of malignant epithelial cells associated with metastasis, we employed the Scissor algorithm to integrate scRNA-seq data with bulk RNA-seq data. As illustrated in Figure 3A and Supplementary Figure S2C, the malignant epithelial cells were categorized into Scissor+ epithelial cells, Scissor- epithelial cells, and undetermined epithelial cells. Notably, there was a higher prevalence of Scissor+ epithelial cells in primary samples compared to Scissor-epithelial cells, which were more abundant in metastatic samples (Supplementary Figure S4A). A total of 741 genes were found to be upregulated, while 793 genes were downregulated in Scissor+ epithelial cells when compared to Scissor- epithelial cells. Among the upregulated genes were KLF6, CTSE, FOSB, CLDN18, FOS, JUN, HBA2, and DPOR1, whereas downregulated genes included RBP1, LCN2, LY6D, ASPH, ADIRF, and ASPH (Figure 3B). GO functional enrichment analysis indicated that the upregulated genes in Scissor+ epithelial cells were primarily involved in cadherin binding, cell-substrate junctions, and the regulation of cell-to-cell adhesion. In contrast, the downregulated genes in Scissor+ epithelial cells were associated with structural constituents of ribosomes, ribosomal subunits, and the respiratory electron transport chain (Figure 3C). The KEGG pathway analysis revealed that the upregulated genes in Scissor+ epithelial cells were enriched in pathways related to focal adhesion, tight junctions, TNF signaling, and ECM-receptor interactions, which are characteristic of tumor metastasis (Figure 3C). Conversely, the downregulated genes were predominantly involved in prion disease, Parkinson’s disease, oxidative phosphorylation, Huntington’s disease, and chemical carcinogenesis related to reactive oxygen species (Figure 3C). Additionally, we observed that the overexpressed genes in the Scissor+ group were enriched in epithelial-mesenchymal transition (EMT) and Transforming growth factor β (TGFβ) signaling pathways (Supplementary Figure S4B).
Figure 3. Characteristics and intercellular communication of Scissor+ epithelial cells in PAAD patients. (A) A UMAP visualization highlights the epithelial cells in PAAD samples, with Scissor+ epithelial cells indicated in red and Scissor- epithelial cells in blue. (B) The volcano plot presents differential gene expression between Scissor+ and Scissor- epithelial cells, where red represents genes that are upregulated in Scissor+ epithelial cells, and blue denotes genes that are downregulated in comparison to Scissor- epithelial cells. (C) GO and KEGG analyses were performed on the upregulated and downregulated genes in Scissor+ epithelial cells relative to Scissor- epithelial cells. (D) The intercellular communication network within the TME is depicted. (E) The incoming and outgoing signaling patterns of the Scissor+ and Scissor- groups are presented. (F) CellCall analysis elucidates the distinct communication pathways between Scissor+ and Scissor- epithelial cells. PAAD, Pancreatic adenocarcinoma; TME, Tumor microenvironment; UMAP, Unified Flowform Approximation and Projection; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.
The cell-cell communication network was utilized to assess the communication probability of Scissor+ and Scissor- epithelial cells with other TME cells using CellChat. The analysis revealed that Scissor+ epithelial cells exhibited a higher level of communication with various TME cells, with the most prominent interactions occurring with mast cells, endothelial cells, dendritic cells, and macrophages (Figures 3D,E). The specific pathways through which Scissor+ and Scissor- epithelial cells interacted with other TME cells are detailed in Supplementary Figure S5. Notably, compared to Scissor- epithelial cells, Scissor+ epithelial cells engaged with mast cells, endothelial cells, dendritic cells, and macrophages via the HLA-LILRB and B2M-LILRB pathways (Figure 3F).
3.4 Construction of prognostic risk model for PAAD
The PAAD patients from the TCGA dataset were categorized into Scissor+ and Scissor- groups based on DEGs. As illustrated in Figure 4A, the Scissor+ group exhibited a poorer prognosis. To further investigate the association between Scissor+ epithelial cells and the prognosis of PAAD patients, LASSO-Cox regression analysis was employed to develop a prognostic risk model, which identified seven genes: TATDN1, CAV2, CLDN1, LIPH, MT1E, PSCA, and MMP3 (Figure 4B). Based on their coefficients, the risk score = expression level of TATDN1 * 0.0383 + expression level of CAV2 * 0.102 + expression level of CLDN1 * 0.0471 + expression level of LIPH * 0.139 + expression level of MT1E * 0.0397 + expression level of PSCA * 0.088 + expression level of MMP3 * 0.055 (Figure 4C). Supplementary Figure S6A showed the correlation analysis between risk model with seven genes. To assess the stability and reliable generalization of the risk model, PAAD patients in the TCGA and ICGC cohorts were separately classified into high-risk and low-risk groups. The high-risk group showed a higher number of deaths (Figures 4D,F). Kaplan-Meier curves indicated that OS, DSS, DFI and PFI were poorer in the high-risk group compared to the low-risk group in the TCGA cohort (Figure 4E; Supplementary Figure S6B,C). Similarly, in the ICGC cohort, PAAD patients in the high-risk group had a worse prognosis than those in the low-risk group (Figure 4G). Additionally, we utilized the ROC curve to evaluate the sensitivity and specificity of this risk model. The results revealed that the AUCs values for 1- and 3-year predictions of OS were 0.74 and 0.81, respectively, in the TCGA cohort (Figure 4H). In the ICGC cohort, the AUCs for 1- and 3-year predictions of OS were 0.56 and 0.64, respectively (Figure 4I). These findings suggest that the risk model has good predictive value in the TCGA cohort.
Figure 4. Establishment of the prognostic risk model. (A) Kaplan-Meier analysis of survival probabilities for PAAD patients in the Scissor+ and Scissor- epithelial cell groups. (B) Volcano plot illustrating differential gene expression between Scissor+ and Scissor- epithelial cells. Genes with a p-value <0.05 and an average log Fold Change >1 (TATDN1, CAV2, CLDN1, LIPH, MT1E, PSCA, and MMP3) were identified as marker genes for Scissor- epithelial cells. (C) Upper: Tuning parameters of OS-related genes were selected to cross-verify the error curve. Lower: The LASSO coefficients for seven genes in the TCGA cohort. Ten-fold cross-validation was applied to calculate the best lambda, which leads to a minimum mean cross-validated error. (D) The risk model stratified PAAD patients in the TCGA cohort into high-risk and low-risk groups. (E) Kaplan-Meier analysis of OS for PAAD patients in the TCGA cohort, comparing high-risk and low-risk groups. (F) The risk model stratified PAAD patients in the ICGC cohort into high-risk and low-risk groups. (G) Kaplan-Meier analysis of OS for PAAD patients in the ICGC cohort, comparing high-risk and low-risk groups. (H) ROC curve analysis in the TCGA cohort. (I) ROC curve analysis in the ICGC cohort. The asterisks represented the statistical P value. *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001; PAAD, pancreatic adenocarcinoma; OS, overall survival; ROC, receiver operating characteristic; LASSO, least absolute shrinkage and selection operator.
3.5 Nomogram development and validation for PAAD
We conducted both univariate and multivariate Cox regression analyses on the TCGA cohort. Our findings indicated that the risk score serves as an independent predictor of OS (Figures 5A,B). Subsequently, we developed a nomogram that incorporates this risk model along with other clinicopathological factors, including age, gender, and tumor stage, to predict the 1-, 3- and 5-year OS in the TCGA cohort (Figure 5C). The predictive performance was evaluated using decision curve analysis and calibration curves. As illustrated in Figures 5D,E, the nomogram demonstrated improved predictive accuracy for the 1-, 3- and 5-year OS of PAAD patients within the TCGA cohort.
Figure 5. Construction a nomogram model to predict the prognosis of PAAD patients. (A) OS related clinical factors identified through univariate Cox regression analysis. (B) OS related clinical factors determined by multivariate Cox regression analysis. (C) Nomogram integrating clinical factors and risk scores. (D) Decision curves comparing the nomogram and risk score alone for predicting 1-, 3- and 5-year outcomes in the PAAD cohort. (E) Calibration curve for predicting 1-, 3- and 5-year OS of PAAD patients in the TCGA cohort. PAAD, pancreatic adenocarcinoma; OS, overall survival.
3.6 Correlation of this risk model with immune microenvironment landscape
To examine the relationship between this risk model and the immune response of tumors, we conducted an analysis investigating the association between the risk model and immune score, stromal score, ESTIMATE score, and tumor purity score. As shown in Figure 6A; Supplementary Figure S7A-C, the low-risk group exhibited a higher immune score, while no significant differences were observed between the risk model and stromal score, ESTIMATE score, or tumor purity score. Additionally, the expression levels of CD274, TIDE, and PIM3 were found to be elevated in the high-risk group (Figures 6B–D). CIBERSORT was employed to analyze the proportions of immune cells, revealing significant differences in naïve B cells, CD8 T cells, M0 macrophages, and activated dendritic cells between the high and low-risk groups (Figure 6E). In addition, we performed the correlation analysis between risk model with immune cells and immune checkpoints. As shown in Supplementary Figure S7D,E, the risk model exhibited a positive association with Macrophage M0, M1, and Dendritic cells, while demonstrating a negative association with CD8+ T cells and Mast cells. Furthermore, the risk model was positively correlated with CD74, CD276, CD274, and HAVCR2, and negatively correlated with CTLA4. Overall, these results suggest that Scissor+ epithelial cells may influence metastasis through interactions with the TME.
Figure 6. Analysis of the immune microenvironment. (A) Comparison of ImmuneScore between high-risk and low-risk groups. (B) Expression levels of CD274 in high-risk versus low-risk groups. (C) Expression levels of TIDE in high-risk versus low-risk groups. (D) Expression levels of PIM3 in high-risk versus low-risk groups. (E) Box plot illustrating the distribution of immune infiltrating cells between high-risk and low-risk groups.
To evaluate the predictive effect of immunotherapy based on a risk model in tumors, we conducted an analysis of the prognostic outcomes and immunotherapy responses in patients with non-small cell lung cancer (NSCLC) and clear cell renal cell carcinoma (ccRCC) categorized by risk model. The results indicated that NSCLC patients in the high-risk group who received immunotherapy tended to experience disease progression, accompanied by unfavorable progression-free survival (PFS) outcomes (Supplementary Figure S8A-C). Additionally, the rates of progressive disease (PD) and stable disease (SD) were higher in the high-risk group, which was associated with poor OS in ccRCC patients undergoing immunotherapy (Supplementary Figure S8B-D).
Subsequently, gene mutations were analyzed to extract molecular characteristics distinguishing the high-risk and low-risk groups. KRAS, TP53, and CDKN2A exhibited the highest mutation frequencies in the high-risk group, whereas TP53, KRAS, and SMAD4 showed the highest mutation frequencies in the low-risk group (Supplementary Figure S9A,B). Additionally, the TMB was greater in the high-risk group compared to the low-risk group (Supplementary Figure S9C). These findings suggest that the mutation frequencies of the key genes differ between the high-risk and low-risk groups.
3.7 The biological effect of LIPH on PAAD cells in vitro
The risk model included seven genes (TATDN1, CAV2, CLDN1, LIPH, MT1E, PSCA, and MMP3), and it remains unclear whether LIPH influences the metastasis of PAAD cells. To investigate the biological effect of LIPH on PAAD cells in vitro, we conducted a series of experiments. As demonstrated in Figure 7A, the mRNA levels of LIPH were higher in most PAAD cell lines compared to HPDE6-C7 cells. Consequently, we selected the cell lines CFPAC-1 and CAPAN-1, which exhibited the highest expression levels of LIPH, for further study. We successfully constructed LIPH knockdown cell lines for both CFPAC-1 and CAPAN-1 (Figures 7B,C). CCK-8 assays and colony formation assays revealed that LIPH knockdown significantly inhibited the proliferation and growth of PAAD cells (Figures 7D–H). Additionally, the migration and invasion capabilities were markedly reduced in LIPH knockdown cells (Figures 7I–M). In conclusion, these findings suggest that LIPH promotes the development and invasion of PAAD cells in vitro.
Figure 7. LIPH knockdown inhibited the growth and migration of PAAD cells. (A) RT-qPCR was employed to assess the mRNA levels of LIPH in seven PAAD cell lines (PANC-1, ASPC-1, CFPAC-1, CAPAN-1, CAPAN-2, MiaPaCa-2, and BXPC-3) alongside one normal human pancreatic duct epithelial cell line (HPDE6-C7). The relative mRNA levels of LIPH in CAPAN-1 (B) and CFPAC-1 (C) cells transfected with the shLIPH virus were analyzed. The CCK-8 assay demonstrated that DLAT knockdown significantly inhibited the proliferation of CAPAN-1 (D) and CFPAC-1 (E) cells. (F) Colony formation assays were conducted to evaluate the colony forming ability of LIPH knockdown cells. The number of colonies formed was significantly reduced in CAPAN-1 (G) and CFPAC-1 (H) cells following LIPH inhibition. (I) Transwell assays were performed to assess the migratory and invasive capabilities of LIPH knockdown cells. LIPH knockdown inhibited the migration ability of CAPAN-1 (J) and CFPAC-1 (K) cells, as well as their invasion ability (L,M). The asterisks represented the statistical P value. *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001. PAAD, pancreatic adenocarcinoma.
4 Discussion
Although chemotherapy holds some value for patients with metastatic PAAD, the prognosis remains poor, underscoring the importance of elucidating the underlying molecular mechanisms of pancreatic cancer metastasis. Recent advancements in scRNA-seq have provided unique advantages in exploring tumor and TME heterogeneity, as well as in identifying novel cell subclusters within tumors, which is crucial for the development of clinically relevant therapeutic targets. This study demonstrated the heterogeneity of PAAD and identified eight distinct clusters using published scRNA-seq datasets including epithelial cells, macrophages, T cells, dendritic cells, endothelial cells, fibroblasts, mixed cells, and mast cells. Most of these clusters are involved in regulating tumor metastasis, proliferation, and responses to immunotherapy (Wei et al., 2019; Kalluri, 2016; Verneau et al., 2020). We identified a metastasis-related population of malignant epithelial cells, termed Scissor+ cells, in PAAD by integrating scRNA-seq data with bulk RNA-seq data. We observed two developmental trajectories for malignant epithelial cells transitioning from primary tumors to metastases. Furthermore, we constructed a prognostic risk model for PAAD patients based on Scissor+ epithelial cells, which effectively predicts OS, DSS, DFI and PFI in the TCGA and ICGA cohorts. Additionally, the high-risk group exhibited a significant correlation with elevated expression levels of CD274, TIDE, and PIM3, suggesting that PAAD patients in this group may benefit from immunotherapy. Finally, TATDN1, CAV2, CLDN1, LIPH, MT1E, PSCA, and MMP3 were identified as candidate genes associated with Scissor+ epithelial cells. Our study is the first to reveal that the knockdown of LIPH can inhibit the proliferation and metastasis of PAAD cells in vitro. Further investigations are necessary to elucidate the underlying molecular mechanisms.
PAAD predominantly arises from pancreatic ductal epithelial cells, and previous studies have shed light on the heterogeneity of neoplastic epithelial cells during the progression of PAAD (Grant et al., 2016). Kim et al. identified a novel pancreatic cancer cell subcluster, Ep_VGLL1, suggesting that it represents an intermediate subcluster between the basal-like and classical clusters, and confirmed its prognostic value for PAAD patients (Kim et al., 2024). Additionally, Zhang et al. found that CEACAM5+/CEACAM6+ ductal cells are associated with poor prognosis in PDAC patients through scRNA-seq (Zhang S. et al., 2023). Notably, only 22% of basal-like malignant ductal cells influenced the response to chemotherapy and OS of PAAD patients, indicating the heterogeneity present within the intra-tumoral environment and TME (Park et al., 2024). These studies collectively establish that scRNA-seq is a crucial tool for identifying tumor heterogeneity and exploring new subclusters. By analyzing scRNA-seq data, we identified 14 epithelial subclusters and detected malignant epithelial cells in PAAD samples through inferred CNV and aneuploidy analysis. We observed two developmental trajectories for malignant epithelial cells from primary tumors to metastasis. During the metastatic process, the expression of CDH1 decreased while CDH2 increased in malignant tumor epithelial cells, consistent with previous studies that implicate EMT signaling in tumor metastasis (Aiello et al., 2018; Ligorio et al., 2019).
The Scissor algorithm can identify the most phenotype-associated cell subpopulations from single-cell data (Sun et al., 2022; Guccini et al., 2021) and has been applied to lung adenocarcinoma (LUAD) to distinguish cancer cells from normal phenotypes (Lambrechts et al., 2018). Based on the high and low Scissor subtypes of dendritic cells, Cheng et al. constructed a prognostic risk model for esophageal squamous cell carcinoma (Cheng et al., 2024). Additionally, Scissor was utilized to identify epithelial cells associated with the prognostic phenotypes of hepatocellular carcinoma (HCC) patients, lymph node metastasis-related cell subpopulations in LUAD, and endothelial cells linked to the prognostic phenotypes of head and neck squamous cell carcinoma (HNSCC) (Qi and Zhang, 2023; Ji et al., 2023; Yang et al., 2024). In this study, the Scissor algorithm was employed to integrate scRNA-seq with bulk RNA-seq data. The results demonstrated that Scissor+ epithelial cells were identified as malignant epithelial cells and were positively correlated with metastasis in PAAD. A greater number of Scissor+ epithelial cells were observed in primary samples, while more Scissor-epithelial cells were found in metastatic samples. This discrepancy may be attributed to the decrease or evolution of metastasis-associated subsets during the progression of metastasis (Gerstberger et al., 2023). Furthermore, the upregulated genes in Scissor+ epithelial cells were predominantly focused on metastasis-related pathways, such as adhesion, extracellular matrix (ECM) receptors, EMT and TGFβ signaling pathways, which aligns with previous studies (Ho et al., 2020; Wang et al., 2023).
Immunotherapy often uses immune checkpoint inhibitor (ICI) drugs, such as antibodies targeting CTLA-4, PD-1, and PD-L1, which are highly effective. TMB and microsatellite instability (MSI) are indicators that help predict the likelihood of a tumor responding to ICIs. Patients with high TMB or MSI levels often had long-term OS following immunotherapy (Chalmers et al., 2017; Yamamoto et al., 2020). The TME plays a significant role in impacting the effectiveness of ICIs treatments (Samstein et al., 2019). PAAD cells interact with various types of immunosuppressive cells, including tumor-associated macrophages (TAMs), myeloid cells, CAFs and Treg cells, to establish a “cold” TME. This cold environment contributes to resistance against ICIs and results in poor prognosis for PAAD patients (Niu et al., 2024). The strong immunosuppressive microenvironment in PAAD lead to immune evasion and rapid tumor progression (Martinez-Bosch et al., 2018). Different mechanisms were involved in tumor immune evasion. Such as cancer cells decreased immune recognition by downregulating antigen presentation pathways, like the major histocompatibility complex (MHC) I proteins. Reprogramming TAMs in TME to desired phenotypes offers a promising approach for cancer immunotherapy due to its precision and low side effects (Wang et al., 2024). TGFβ signaling in TME suppresses the antitumor functions of various immune cell populations including T cells, resulting immune suppression severely limits the efficacy of ICIs (Derynck et al., 2021). Cytokines are key mediators of cell communication in TME. Some cytokines contribute to host antitumor responses (Propper and Balkwill, 2022). Through the release of cytokines and growth factors, the M2 TAMs in TME establish an immunosuppressive niche in escaping immune detection to inhibit the anti-tumor immunity and drive tumor growth (Mantovani et al., 2017). Previous study examined disulfidptosis-related genes (DRGs) in pan-cancer and indicated the DRGs linked with poor prognosis and predicted responsiveness to immunotherapy (Xu et al., 2025). Yang et al. indicated that PLRN3 mediated M2 macrophage infiltration in various cancers and thus impacted the effectiveness of immunotherapy in LUAD patients (Yang et al., 2025). TMB is a key indicator for predicting immunotherapy response in tumors and individuals with higher TMB levels also have better response rates with immunotherapy (Yu et al., 2021; Li et al., 2020). Xu et al. indicated EPHB2 expression positively correlated with TMB, MSI in pan-cancer and revealed that EPHB2 could improve the predictive effect of immunotherapy responses (Xu et al., 2024). Our analysis also confirmed a positive relationship between the risk model and high TMB, suggesting that high-risk group predict better responses to immunotherapy. In this study, the proportions of naïve B cells, CD8 T cells, M0 macrophages, and activated dendritic cells differed between Scissor+ and Scissor- epithelial cells. Compared to Scissor-epithelial cells, Scissor+ epithelial cells exhibited higher levels of cell-cell communication, particularly with dendritic cells, macrophages, and mast cells, suggesting that Scissor+ epithelial cells may influence the response to immunotherapy. The HLA-LILRB and B2M-LILRB pathways play significant roles in these interactions, consistent with the notion that HLA class I may interact with LILRB1 or LILRB2 to impair the antibody-dependent cellular phago-cytosis function of macrophages (Zeller et al., 2022). LILRB1 and LILRB2 are expressed in various tumor types and are recognized as immunosuppressive markers within the tumor TME (van der Touw et al., 2017). Inhibition of LILRB2 has been shown to enhance T cell activation (Umiker et al., 2023). Furthermore, LILRB2 can interact with HLA-G in the TME, promoting myeloid cell tumor proliferation and increasing tumor immune evasion (Carosella et al., 2021). This study identified the existence of HLA-LILRB pathways between Scissor+ epithelial cells and dendritic cells, macrophages, and mast cells in the TME, mediating the transmission of immune signal inhibition and promoting metastasis. Further investigation into the underlying molecular mechanisms is warranted.
We constructed a prognostic model and identified seven candidate upregulated genes, which included TATDN1, CAV2, CLDN1, LIPH, MT1E, PSCA, and MMP3. TATDN1 is associated with DNA nuclease activity and endodeoxyribonuclease activity; mutations in the TATDN1 gene indicate vascular invasion in HCC patients (Xu et al., 2023). In addition, six copper dependent-related genes including TATDN1were used to perform the prognostic model in PAAD and it indicated that this model can be applied into estimating the prognosis and immunological microenvironment of pancreatic cancer patients (Guan et al., 2022). CAV2 is elevated in various types of tumors and possesses oncogenic properties that promote the progression and metastasis of PAAD cells (Wang et al., 2022). CLDN1, a member of the tight junction protein family, is considered a clinical therapeutic target. It has been reported to enhance CRC stemness and chemoresistance by interacting with EPHA2 (Primeaux et al., 2023). Knockdown of CLDN1 promotes stemness, migration, and invasion of pancreatic cancer cells (Zhu et al., 2024). MT1E has been identified as a tumor suppressor, with its methylation correlating with HCC metastasis (Liu et al., 2020). Additionally, MT1E with other ten genes were used to perform the prognostic model and indicated its prognostic predictive effect for pancreatic cancer patients (Huang et al., 2022). PSCA is upregulated in most cancers, including prostate cancer, bladder cancer, and PAAD. Previous studies have developed human iPSC-derived CAR macrophages and CAR NK cells targeting PSCA, demonstrating strong anti-pancreatic cancer effects in mouse models (Shah et al., 2024). MMP3, a member of the MMP family, has been shown to induce migration of pancreatic cancer cells (Tjomsland et al., 2016). Although Zhuang et al. identified LIPH as a novel unfavorable prognostic biomarker correlated with immunosuppression in pancreatic cancer patients with RNA-seq database, the relationship between LIPH with PAAD metastasis remains unclear (Zhuang et al., 2022). In this study, we first investigated the function of LIPH in PAAD cells. The results confirmed LIPH’s pro-growth and metastatic effects in PAAD cells, thereby laying the foundation for further research into its mechanism of action in the future.
This study has several limitations. Firstly, it relies on previously published scRNA-seq and bulk RNA data, necessitating further validation of the reliability of the prognostic prediction model in future research. The risk model is validated in the TCGA and ICGC cohorts. These two cohorts lack detailed treatment data and have potential biases in survival data. The heterogeneity of PAAD patients should also be considered. Secondly, while we demonstrated that Scissor+ epithelial cells interact with dendritic cells and macrophages through the HLA-LILRB pathways, the specific molecular mechanisms underlying these interactions require further investigation. Additionally, further functional experiments will be necessary to elucidate the biological role of the LIPH gene in PAAD metastasis and to determine whether it can be targeted to enhance the effectiveness of immunotherapies and chemotherapies.
To conclude, the scRNA-seq dataset was utilized to identify the intra-tumor heterogeneity of PAAD, revealing two developmental trajectories for malignant epithelial cells from primary tumors to metastases. We constructed prognostic models based on differentially expressed genes in Scissor+ epithelial cells, demonstrating the significance of this model in predicting the prognosis and immunotherapy response of PAAD patients. Furthermore, the inhibition of LIPH was shown to reduce the metastasis of PAAD cells in vitro, suggesting its potential as a marker gene for metastasis-related malignant epithelial cells. This study significantly enhances our understanding of metastasis-specific cell subpopulations and aids in predicting prognosis and immunotherapy outcomes.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Ethics statement
Ethical approval was not required for the studies on humans in accordance with the local legislation and institutional requirements because only commercially available established cell lines were used.
Author contributions
GZ: Writing – original draft, Software, Data curation, Visualization, Funding acquisition. JY: Formal Analysis, Visualization, Funding acquisition, Writing – original draft. FZ: Visualization, Validation, Formal Analysis, Writing – original draft. FW: Visualization, Data curation, Writing – original draft. CZ: Validation, Project administration, Funding acquisition, Writing – review and editing, Formal Analysis, Writing – original draft.
Funding
The author(s) declare that financial support was received for the research and/or publication of this article. This research was funded by Key research and development projects of Shaanxi Province, grant number 2023-YBSF-329, Young Talents Fund of Association for Science and Technology in Shaanxi, grant number 20240341, Incubation Fund Program of Shaanxi Provincial People’s Hospital, grant number 2023YJY-04 and 2023YJY-49, and Scientific and Technological Talents support Program Foundation of Shaanxi Provincial People’s Hospital, grant number 2022JY-08.
Conflict of interest
The authors declare that the research was conducted with no commercial or financial relationships that could be viewed as potential conflicts of interest.
Generative AI statement
The author(s) declare that no Generative AI was used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar.2025.1646840/full#supplementary-material
References
Aiello, N. M., Maddipati, R., Norgard, R. J., Balli, D., Li, J., Yuan, S., et al. (2018). EMT subtype influences epithelial plasticity and mode of cell migration. Dev. Cell 45 (6), 681–695. doi:10.1016/j.devcel.2018.05.027
Braun, D. A., Hou, Y., Bakouny, Z., Ficial, M., Sant' Angelo, M., Forman, J., et al. (2020). Interplay of somatic alterations and immune infiltration modulates response to PD-1 blockade in advanced clear cell renal cell carcinoma. Nat. Med. 26 (6), 909–918. doi:10.1038/s41591-020-0839-y
Bray, F., Laversanne, M., Sung, H., Ferlay, J., Siegel, R. L., Soerjomataram, I., et al. (2024). Global cancer statistics 2022: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J. Clin. 74 (3), 229–263. doi:10.3322/caac.21834
Cabello-Aguilar, S., Alame, M., Kon-Sun-Tack, F., Fau, C., Lacroix, M., and Colinge, J. (2020). SingleCellSignalR: inference of intercellular networks from single-cell transcriptomics. Nucleic Acids Res. 48 (10), e55. doi:10.1093/nar/gkaa183
Carosella, E. D., Gregori, S., and Tronik-Le Roux, D. (2021). HLA-G/LILRBs: a cancer immunotherapy challenge. Trends Cancer 7 (5), 389–392. doi:10.1016/j.trecan.2021.01.004
Chalmers, Z. R., Connelly, C. F., Fabrizio, D., Gay, L., Ali, S. M., Ennis, R., et al. (2017). Analysis of 100,000 human cancer genomes reveals the landscape of tumor mutational burden. Genome Med. 9 (1), 34. doi:10.1186/s13073-017-0424-2
Cheng, M., Xiong, J., Liu, Q., Zhang, C., Li, K., Wang, X., et al. (2024). Integrating bulk and single-cell sequencing data to construct a scissor (+) dendritic cells prognostic model for predicting prognosis and immune responses in ESCC. Cancer Immunol. Immunother. 73 (6), 97. doi:10.1007/s00262-024-03683-9
Conroy, T., Pfeiffer, P., Vilgrain, V., Lamarca, A., Seufferlein, T., O'Reilly, E. M., et al. (2023). Pancreatic cancer: ESMO Clinical Practice Guideline for diagnosis, treatment and follow-up. Ann. Oncol. 34 (11), 987–1002. doi:10.1016/j.annonc.2023.08.009
Derynck, R., Turley, S. J., and Akhurst, R. J. (2021). TGFβ biology in cancer progression and immunotherapy. Nat. Rev. Clin. Oncol. 18 (1), 9–34. doi:10.1038/s41571-020-0403-1
Fan, X., Lu, P., Wang, H., Bian, S., Wu, X., Zhang, Y., et al. (2022). Integrated single-cell multiomics analysis reveals novel candidate markers for prognosis in human pancreatic ductal adenocarcinoma. Cell Discov. 8 (1), 13. doi:10.1038/s41421-021-00366-y
Gao, R., Bai, S., Henderson, Y. C., Lin, Y., Schalck, A., Yan, Y., et al. (2021). Delineating copy number and clonal substructure in human tumors from single-cell transcriptomes. Nat. Biotechnol. 39 (5), 599–608. doi:10.1038/s41587-020-00795-2
Gerstberger, S., Jiang, Q., and Ganesh, K. (2023). Metastasis. Cell 186 (8), 1564–1579. doi:10.1016/j.cell.2023.03.003
Grant, T. J., Hua, K., and Singh, A. (2016). Molecular pathogenesis of pancreatic cancer. Prog. Mol. Biol. Transl. Sci. 144, 241–275. doi:10.1016/bs.pmbts.2016.09.008
Grunwald, B. T., Devisme, A., Andrieux, G., Vyas, F., Aliar, K., McCloskey, C. W., et al. (2021). Spatially confined sub-tumor microenvironments in pancreatic cancer. Cell 184 (22), 5577–5592.e18. doi:10.1016/j.cell.2021.09.022
Guan, X., Lu, N., and Zhang, J. (2022). The combined prognostic model of copper-dependent to predict the prognosis of pancreatic cancer. Front. Genet. 13, 978988. doi:10.3389/fgene.2022.978988
Guccini, I., Revandkar, A., D'Ambrosio, M., Colucci, M., Pasquini, E., Mosole, S., et al. (2021). Senescence reprogramming by TIMP1 deficiency promotes prostate cancer metastasis. Cancer Cell 39 (1), 68–82 e9. doi:10.1016/j.ccell.2020.10.012
Han, J., DePinho, R. A., and Maitra, A. (2021). Single-cell RNA sequencing in pancreatic cancer. Nat. Rev. Gastroenterol. Hepatol. 18 (7), 451–452. doi:10.1038/s41575-021-00471-z
Ho, W. J., Jaffee, E. M., and Zheng, L. (2020). The tumour microenvironment in pancreatic cancer - clinical challenges and opportunities. Nat. Rev. Clin. Oncol. 17 (9), 527–540. doi:10.1038/s41571-020-0363-5
Hu, X., Meng, Y., Xu, L., Qiu, L., Wei, M., Su, D., et al. (2019). Cul4 E3 ubiquitin ligase regulates ovarian cancer drug resistance by targeting the antiapoptotic protein BIRC3. Cell Death Dis. 10 (2), 104. doi:10.1038/s41419-018-1200-y
Hu, Z. I., and O'Reilly, E. M. (2024). Therapeutic developments in pancreatic cancer. Nat. Rev. Gastroenterol. Hepatol. 21 (1), 7–24. doi:10.1038/s41575-023-00840-w
Huang, X., Feng, Y., Ma, D., Ding, H., Dong, G., Chen, Y., et al. (2022). The molecular, immune features, and risk score construction of intraductal papillary mucinous neoplasm patients. Front. Mol. Biosci. 9, 887887. doi:10.3389/fmolb.2022.887887
Ji, X., Wang, Z., Wang, G., Tang, L., and Han, Z. (2023). Single-cell transcriptomics reveals the drivers and therapeutic targets of lymph node metastasis in lung adenocarcinoma. Aging (Albany NY) 15 (14), 7023–7037. doi:10.18632/aging.204890
Kalluri, R. (2016). The biology and function of fibroblasts in cancer. Nat. Rev. Cancer 16 (9), 582–598. doi:10.1038/nrc.2016.73
Kim, J. Y., Choi, J. K., and Jung, H. (2020). Genome-wide methylation patterns predict clinical benefit of immunotherapy in lung cancer. Clin. Epigenetics 12 (1), 119. doi:10.1186/s13148-020-00907-4
Kim, S., Leem, G., Choi, J., Koh, Y., Lee, S., Nam, S. H., et al. (2024). Integrative analysis of spatial and single-cell transcriptome data from human pancreatic cancer reveals an intermediate cancer cell population associated with poor prognosis. Genome Med. 16 (1), 20. doi:10.1186/s13073-024-01287-7
Lambrechts, D., Wauters, E., Boeckx, B., Aibar, S., Nittner, D., Burton, O., et al. (2018). Phenotype molding of stromal cells in the lung tumor microenvironment. Nat. Med. 24 (8), 1277–1289. doi:10.1038/s41591-018-0096-5
Li, R., Han, D., Shi, J., Han, Y., Tan, P., Zhang, R., et al. (2020). Choosing tumor mutational burden wisely for immunotherapy: a hard road to explore. Biochim. Biophys. Acta Rev. Cancer 1874 (2), 188420. doi:10.1016/j.bbcan.2020.188420
Ligorio, M., Sil, S., Malagon-Lopez, J., Nieman, L. T., Misale, S., Di Pilato, M., et al. (2019). Stromal microenvironment shapes the intratumoral Architecture of pancreatic cancer. Cell 178 (1), 160–175. doi:10.1016/j.cell.2019.05.012
Lin, W., Noel, P., Borazanci, E. H., Lee, J., Amini, A., Han, I. W., et al. (2020). Single-cell transcriptome analysis of tumor and stromal compartments of pancreatic ductal adenocarcinoma primary tumors and metastatic lesions. Genome Med. 12 (1), 80. doi:10.1186/s13073-020-00776-9
Liu, Q., Lu, F., and Chen, Z. (2020). Identification of MT1E as a novel tumor suppressor in hepatocellular carcinoma. Pathol. Res. Pract. 216 (11), 153213. doi:10.1016/j.prp.2020.153213
Mantovani, A., Marchesi, F., Malesci, A., Laghi, L., and Allavena, P. (2017). Tumour-associated macrophages as treatment targets in oncology. Nat. Rev. Clin. Oncol. 14 (7), 399–416. doi:10.1038/nrclinonc.2016.217
Martinez-Bosch, N., Vinaixa, J., and Navarro, P. (2018). Immune evasion in pancreatic cancer: from mechanisms to therapy. Cancers (Basel) 10 (1), 6. doi:10.3390/cancers10010006
McGinnis, C. S., Murrow, L. M., and Gartner, Z. J. (2019). DoubletFinder: doublet detection in single-cell RNA sequencing data using artificial nearest neighbors. Cell Syst. 8 (4), 329–337. doi:10.1016/j.cels.2019.03.003
Moncada, R., Barkley, D., Wagner, F., Chiodin, M., Devlin, J. C., Baron, M., et al. (2020). Integrating microarray-based spatial transcriptomics and single-cell RNA-seq reveals tissue architecture in pancreatic ductal adenocarcinomas. Nat. Biotechnol. 38 (3), 333–342. doi:10.1038/s41587-019-0392-8
Niu, N., Shen, X., Wang, Z., Chen, Y., Weng, Y., Yu, F., et al. (2024). Tumor cell-intrinsic epigenetic dysregulation shapes cancer-associated fibroblasts heterogeneity to metabolically support pancreatic cancer. Cancer Cell 42 (5), 869–884.e9. doi:10.1016/j.ccell.2024.03.005
Park, W., Chawla, A., and O'Reilly, E. M. (2021). Pancreatic cancer: a review. JAMA 326 (9), 851–862. doi:10.1001/jama.2021.13027
Park, J. K., Jeong, H. O., Kim, H., Choi, J. H., Lee, E. M., Kim, S., et al. (2024). Single-cell transcriptome analysis reveals subtype-specific clonal evolution and microenvironmental changes in liver metastasis of pancreatic adenocarcinoma and their clinical implications. Mol. Cancer 23 (1), 87. doi:10.1186/s12943-024-02003-0
Patel, A. P., Tirosh, I., Trombetta, J. J., Shalek, A. K., Gillespie, S. M., Wakimoto, H., et al. (2014). Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science 344 (6190), 1396–1401. doi:10.1126/science.1254257
Primeaux, M., Liu, X., Gowrikumar, S., Fatima, I., Fisher, K. W., Bastola, D., et al. (2023). Claudin-1 interacts with EPHA2 to promote cancer stemness and chemoresistance in colorectal cancer. Cancer Lett. 579, 216479. doi:10.1016/j.canlet.2023.216479
Propper, D. J., and Balkwill, F. R. (2022). Harnessing cytokines and chemokines for cancer therapy. Nat. Rev. Clin. Oncol. 19 (4), 237–253. doi:10.1038/s41571-021-00588-9
Qi, W., and Zhang, Q. (2023). Insights on epithelial cells at the single-cell level in hepatocellular carcinoma prognosis and response to chemotherapy. Front. Pharmacol. 14, 1292831. doi:10.3389/fphar.2023.1292831
Qiu, X., Mao, Q., Tang, Y., Wang, L., Chawla, R., Pliner, H. A., et al. (2017). Reversed graph embedding resolves complex single-cell trajectories. Nat. Methods 14 (10), 979–982. doi:10.1038/nmeth.4402
Raghavan, S., Winter, P. S., Navia, A. W., Williams, H. L., DenAdel, A., Lowder, K. E., et al. (2021). Microenvironment drives cell state, plasticity, and drug response in pancreatic cancer. Cell 184 (25), 6119–6137.e26. doi:10.1016/j.cell.2021.11.017
Samstein, R. M., Lee, C. H., Shoushtari, A. N., Hellmann, M. D., Shen, R., Janjigian, Y. Y., et al. (2019). Tumor mutational load predicts survival after immunotherapy across multiple cancer types. Nat. Genet. 51 (2), 202–206. doi:10.1038/s41588-018-0312-8
Satija, R., Farrell, J. A., Gennert, D., Schier, A. F., and Regev, A. (2015). Spatial reconstruction of single-cell gene expression data. Nat. Biotechnol. 33 (5), 495–502. doi:10.1038/nbt.3192
Schalck, A., Sakellariou-Thompson, D., Forget, M. A., Sei, E., Hughes, T. G., Reuben, A., et al. (2022). Single-Cell sequencing reveals trajectory of tumor-infiltrating lymphocyte States in pancreatic cancer. Cancer Discov. 12 (10), 2330–2349. doi:10.1158/2159-8290.CD-21-1248
Shah, Z., Tian, L., Li, Z., Jin, L., Zhang, J., Li, Z., et al. (2024). Human anti-PSCA CAR macrophages possess potent antitumor activity against pancreatic cancer. Cell Stem Cell 31 (6), 803–817.e6. doi:10.1016/j.stem.2024.03.018
Siegel, R. L., Giaquinto, A. N., and Jemal, A. (2024). Cancer statistics, 2024. CA Cancer J. Clin. 74 (1), 12–49. doi:10.3322/caac.21820
Stoop, T. F., Javed, A. A., Oba, A., Koerkamp, B. G., Seufferlein, T., Wilmink, J. W., et al. (2025). Pancreatic cancer. Pancreat. cancer. Lancet. 405 (10485), 1182–1202. doi:10.1016/S0140-6736(25)00261-2
Sun, D., Guan, X., Moran, A. E., Wu, L. Y., Qian, D. Z., Schedin, P., et al. (2022). Identifying phenotype-associated subpopulations by integrating bulk and single-cell sequencing data. Nat. Biotechnol. 40 (4), 527–538. doi:10.1038/s41587-021-01091-3
Tjomsland, V., Pomianowska, E., Aasrum, M., Sandnes, D., Verbeke, C. S., and Gladhaug, I. P. (2016). Profile of MMP and TIMP expression in human pancreatic stellate cells: regulation by IL-1α and TGFβ and implications for migration of pancreatic cancer cells. Neoplasia 18 (7), 447–456. doi:10.1016/j.neo.2016.06.003
Umiker, B., Hashambhoy-Ramsay, Y., Smith, J., Rahman, T., Mueller, A., Davidson, R., et al. (2023). Inhibition of LILRB2 by a novel blocking antibody designed to reprogram immunosuppressive macrophages to drive T-Cell activation in tumors. Mol. Cancer Ther. 22 (4), 471–484. doi:10.1158/1535-7163.MCT-22-0351
van der Touw, W., Chen, H. M., Pan, P. Y., and Chen, S. H. (2017). LILRB receptor-mediated regulation of myeloid cell maturation and function. Cancer Immunol. Immunother. 66 (8), 1079–1087. doi:10.1007/s00262-017-2023-x
Verneau, J., Sautes-Fridman, C., and Sun, C. M. (2020). Dendritic cells in the tumor microenvironment: prognostic and theranostic impact. Semin. Immunol. 48, 101410. doi:10.1016/j.smim.2020.101410
Von Hoff, D. D., Ervin, T., Arena, F. P., Chiorean, E. G., Infante, J., Moore, M., et al. (2013). Increased survival in pancreatic cancer with nab-paclitaxel plus gemcitabine. N. Engl. J. Med. 369 (18), 1691–1703. doi:10.1056/NEJMoa1304369
Wang, Y., Wang, Y., Liu, R., Wang, C., Luo, Y., Chen, L., et al. (2022). CAV2 promotes the invasion and metastasis of head and neck squamous cell carcinomas by regulating S100 proteins. Cell Death Discov. 8 (1), 386. doi:10.1038/s41420-022-01176-1
Wang, X., Eichhorn, P. J. A., and Thiery, J. P. (2023). TGF-beta, EMT, and resistance to anti-cancer treatment. Semin. Cancer Biol. 97, 1–11. doi:10.1016/j.semcancer.2023.10.004
Wang, J., Lu, Y., Zhang, R., Cai, Z., Fan, Z., Xu, Y., et al. (2024). Modulating and imaging macrophage reprogramming for cancer immunotherapy. Phenomics 4 (4), 401–414. doi:10.1007/s43657-023-00154-6
Wei, C., Yang, C., Wang, S., Shi, D., Zhang, C., Lin, X., et al. (2019). Crosstalk between cancer cells and tumor associated macrophages is required for mesenchymal circulating tumor cell-mediated colorectal cancer metastasis. Mol. Cancer 18 (1), 64. doi:10.1186/s12943-019-0976-4
Xu, J., Zhou, Y., Dong, K., Gong, J., Xiong, W., Wang, X., et al. (2023). Gene variation profile and it's potential correlation with clinical characteristics in HBV-associated HCC patients of Sichuan Han nationality in China. Asian J. Surg. 46 (10), 4371–4377. doi:10.1016/j.asjsur.2023.02.056
Xu, S., Zheng, Y., Ye, M., Shen, T., Zhang, D., Li, Z., et al. (2024). Comprehensive pan-cancer analysis reveals EPHB2 is a novel predictive biomarker for prognosis and immunotherapy response. BMC Cancer 24 (1), 1064. doi:10.1186/s12885-024-12843-0
Xu, S., Chen, Z., Chen, X., Chu, H., Huang, X., Chen, C., et al. (2025). Interplay of disulfidptosis and the tumor microenvironment across cancers: implications for prognosis and therapeutic responses. BMC Cancer 25 (1), 1113. doi:10.1186/s12885-025-14246-1
Yamamoto, H., Watanabe, Y., Maehata, T., Imai, K., and Itoh, F. (2020). Microsatellite instability in cancer: a novel landscape for diagnostic and therapeutic approach. Arch. Toxicol. 94 (10), 3349–3357. doi:10.1007/s00204-020-02833-z
Yang, C., Cheng, X., Gao, S., and Pan, Q. (2024). Integrating bulk and single-cell data to predict the prognosis and identify the immune landscape in HNSCC. J. Cell Mol. Med. 28 (1), e18009. doi:10.1111/jcmm.18009
Yang, S., Liu, H., Zheng, Y., Chu, H., Lu, Z., Yuan, J., et al. (2025). The role of PLIN3 in prognosis and tumor-associated macrophage infiltration: a pan-cancer analysis. J. Inflamm. Res. 18, 3757–3777. doi:10.2147/JIR.S509245
Yu, G., Pang, Y., Merchant, M., Kesserwan, C., Gangalapudi, V., Abdelmaksoud, A., et al. (2021). Tumor mutation burden, expressed neoantigens and the immune microenvironment in diffuse gliomas. Cancers (Basel) 13 (23), 6092. doi:10.3390/cancers13236092
Zeller, T., Lutz, S., Munnich, I. A., Windisch, R., Hilger, P., Herold, T., et al. (2022). Dual checkpoint blockade of CD47 and LILRB1 enhances CD20 antibody-dependent phagocytosis of lymphoma cells by macrophages. Front. Immunol. 13, 929339. doi:10.3389/fimmu.2022.929339
Zhang, Z., Zhu, X. Q., Yang, F., Lai, N. N., Zhu, L., Cole, K., et al. (2023a). Single-cell mapping reveals several immune subsets associated with liver metastasis of pancreatic ductal adenocarcinoma. Med 4 (10), 728–743.e7. doi:10.1016/j.medj.2023.07.010
Zhang, S., Fang, W., Zhou, S., Zhu, D., Chen, R., Gao, X., et al. (2023b). Single cell transcriptomic analyses implicate an immunosuppressive tumor microenvironment in pancreatic cancer liver metastasis. Nat. Commun. 14 (1), 5123. doi:10.1038/s41467-023-40727-7
Zhu, L., Tang, N., Hang, H., Zhou, Y., Dong, J., Yang, Y., et al. (2024). Loss of Claudin-1 incurred by DNMT aberration promotes pancreatic cancer progression. Cancer Lett. 586, 216611. doi:10.1016/j.canlet.2024.216611
Keywords: Single-cell RNA sequencing, PAAD, Scissor, Epithelial cells, Prognosis, Tumor microenvironment
Citation: Zhang G, Yu J, Zhang F, Wang F and Zhou C (2025) Constructing Scissor+ risk model to predict prognosis and immunotherapy responses in PAAD by integrating bulk and single-cell RNA sequencing data. Front. Pharmacol. 16:1646840. doi: 10.3389/fphar.2025.1646840
Received: 14 June 2025; Accepted: 21 August 2025;
Published: 19 September 2025.
Edited by:
Chen Ling, Fudan University, ChinaReviewed by:
Chenghui Yu, Fudan University, ChinaJingkun Zhang, University of California San Francisco, United States
Copyright © 2025 Zhang, Yu, Zhang, Wang and Zhou. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Congya Zhou, c21hbnQzMkAxNjMuY29t