Integration of single-cell RNA and bulk RNA sequencing revealed malignant ductal cell heterogeneity and prognosis signatures in pancreatic cancer

Abstract

Introduction: Pancreatic cancer is a highly malignant tumor of the digestive system with a dismal prognosis. Despite advances in diagnosis and treatment, overall survival remains extremely low. Early diagnostic markers and an improved understanding of tumor-microenvironment interactions are essential for developing more effective therapies.

Methods: We analyzed 74 single-cell RNA sequencing (scRNA-seq) samples, performing unsupervised clustering and marker-gene expression profiling to define major cell types. Large-scale chromosomal copy-number variation (CNV) analysis distinguished malignant from non-malignant ductal cells. Non-negative matrix factorization (NMF) identified stage-associated gene modules, which were integrated with TCGA bulk-RNA data and machine-learning feature selection to pinpoint candidate prognostic genes. Two independent cohorts were used for validation. Regulatory network inference (pySCENIC) and ligand-receptor interaction analysis (CellPhoneDB) explored cross-talk between malignant cells and macrophages. Finally, in vitro knockdown of CTSV assessed its functional role in pancreatic cancer (PAC) cell proliferation and migration.

Results: Three prognosis-related genes-ANLN, NT5E, and CTSV-were selected based on their strong association with clinical stage and validated in external datasets. High expression of these genes correlated with poorer overall survival and an increased infiltration of M0 macrophages. CellPhoneDB predicted significant interactions between high-expression malignant ductal cells and M0 macrophages via CXCL14-CXCR4 and IL1RAP-PTPRF axes, with SPI1 identified as an upstream regulator of IL1RAP. In vitro CTSV knockdown significantly inhibited PAC cell proliferation and migration.

Discussion: Our integrative single-cell and bulk-RNA workflow identifies ANLN, NT5E, and CTSV as novel prognostic biomarkers in pancreatic cancer and highlights a pro-tumorigenic interaction between malignant ductal cells and macrophages. Targeting CTSV or disrupting CXCL14-CXCR4 and IL1RAP-PTPRF signaling may offer new therapeutic avenues for PAC.

Keywords: CTSV; CXCL14-CXCR4; macrophages; pancreatic cancer; tumor microenvironment.

Figure 1

Single cell transcriptional landscape of PAC. (A) Uniform Manifold Approximation and Projection (UMAP) representation of the landscape of different PCA cell types, (B) Dot‐plots for the merged scRNA‐seq data demonstrates the marker expressions in the different cell types, (C) Bar chart illustrating the cell-type proportions across different anatomical locations (e.g., pancreatic head, body, liver metastasis), (D) Bar chart comparing cell-type distributions between tumor tissues and adjacent normal tissues, (E) Bar chart depicting the relative abundance of each cell type in male vs. female patients, (F) Bar chart showing cell-type composition stratified by clinical stages (I, II, III, and IV), (G) Bar chart presenting cell-type proportions under various treatment modalities (e.g., FOLFIRINOX, gemcitabine + nab-paclitaxel, untreated).

Figure 2

Characterization of malignant ductal cell states. (A) Heatmap of inferCNV analysis, illustrating large-scale chromosomal copy number variations in malignant versus non-malignant ductal cells. (B) UMAP plot showing the distribution of malignant and non-malignant ductal cells. (C) Bar chart depicting the distribution of malignant and non-malignant ductal cells across different tumor locations. (D) GSEA analysis of malignant and non-malignant ductal cells, highlighting significant pathways. (E) CellPhoneDB analysis of cell-cell interactions between malignant and non-malignant ductal cells.

Figure 3

Characteristics of malignant ductal cells (A) UMAP plot showing the clustering of malignant ductal cells. (B-F) Bar charts illustrating the biological characteristics of each subpopulation: SAMD12+, RND3+, RPLP1+, MKI67+, and VIM+ ductal cells. (G) Distribution of different subpopulations across various tumor locations. (H) Distribution of different subpopulations across tumor stages.

Figure 4

NMF analysis of malignancy-associated module in malignant ductal cells (A, B) GSEA of SAMD12+ and MKI67+ subpopulations, respectively. (C) Heatmap showing the correlation of five NMF modules with stages. (D, E) GO and MsigDB pathway analysis of the top 200 genes in the Usage_2 module, respectively.

Figure 5

Prognostic analysis of genes associated with Usage_2 (A) Distribution of LASSO regression coefficients for survival associated AS events (left). Shrinkage parameter selection in the LASSO model using ten-fold cross-validation via minimum criteria (right). (B) KM survival curve of overall survival in TCGA-PAC. (C) KM survival curve of overall survival in ICGC. (D) KM survival curve of overall survival in GSE62452. (E) Box plot showed the ratio differentiation of 21 kinds of immune cells between TCGA-PAC samples with high-risk and low-risk groups. Significance levels were denoted as *P < 0.05. “ns” indicates no significant difference.

Figure 6

M0 Macrophages Interact with Malignant Ductal Cells (A) Bubble plots illustrate the intercellular interactions among M0 macrophages and Malignant Ductal Cells, with bubble size and color indicating the significance and strength of each ligand–receptor pair. (B) Heatmap comparing the most active TFs across different subpopulations, depicted as Z-scores.

Figure 7

In vitro experimental validation. (A-F) Representative microscopic images showing the effects of gene knockdown in Capan-2 cells. (A) Control (0h), (B) shRNA targeting CTSV A (0h), (C) shRNA targeting CTSV B (0h), (D) Control (48h), (E) shRNA targeting CTSV A (48h), and (F) shRNA targeting CTSV B (48h). (G-L) Representative microscopic images showing the effects of gene knockdown in MIAPaCa cells. (G) Control (0h), (H) shRNA targeting CTSV A (0h), (I) shRNA targeting CTSV B (0h), (J) Control (48h), (K) shRNA targeting CTSV A (48h), and (L) shRNA targeting CTSV B (48h). (M) Invasion rate of Capan-2 cells following gene knockdown with different shRNAs, as measured by a transwell assay. (N) Invasion rate of MIAPaCa cells following gene knockdown with different shRNAs, as measured by a transwell assay. Data are presented as mean ± SD from three independent experiments. **P < 0.01, ***P < 0.001 and ***P < 0.0001.