Transcriptome profiles and corresponding clinicopathological information for PDAC were collected from The Cancer Genome Atlas (TCGA)-pancreatic adenocarcinoma (PAAD) and two independent Gene Expression Omnibus (GEO) cohorts (GSE62165 and GSE85916). To minimize inter-dataset variability, batch correction was implemented using the “ComBat” function within the sva R package (15, 16).

Immune cell infiltration in TCGA-PAAD, GSE62165, and GSE85916 samples was quantified using the TIME method. Weighted gene co-expression network analysis (WGCNA) was then employed to construct gene modules correlated with DC infiltration across all cohorts. Intersection of module genes via Venn diagrams yielded a PDAC-specific DC–related gene set. Based on these genes, unsupervised consensus clustering of TCGA-PAAD samples based on this gene set was performed with the ConsensusClusterPlus package. Differentially expressed genes (DEGs) between clusters were determined by the limma package, with significance thresholds of adjusted P< 0.05.

To characterize the biological significance of candidate genes, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were conducted. In addition, Gene Set Enrichment Analysis (GSEA) was performed to evaluate pathway–level enrichment (17). All analyses were performed using the clusterProfiler package.

Patients from the TCGA-PAAD dataset were first stratified into two subtypes by consensus clustering. Candidate genes were obtained from the intersection of DEGs between subtypes and those between tumor and adjacent non-tumor samples. The TCGA-PAAD cohort was randomly divided into training and validation sets using the caret package. In the training cohort, prognostic genes were screened by univariate Cox regression and subsequently refined with least absolute shrinkage and selection operator (LASSO) regression using the glmnet package (18). Risk scores were generated as a weighted sum of gene expression values, with Cox coefficients as weights. Patients were dichotomized into high- and low-risk groups according to the median risk score, and the model’s predictive performance was further confirmed in the TCGA validation set.

Single-cell RNA sequencing data from 25 PDAC patients were acquired from the GEO database (GSE242230) for comprehensive analysis. The gene expression matrix was processed using the Seurat package, with stringent normalization and filtering criteria applied: a minimum of 200 genes per cell, no more than 20% mitochondrial genes, and at least 5% ribosomal genes (19). Subsequently, nonlinear dimensionality reduction, cell clustering, and visualization were performed to delineate the expression profiles of target genes across various cell types.

Cell clustering analysis was conducted by t-distributed stochastic neighbor embedding (t-SNE) dimensionality reduction with a resolution parameter of 0.5. Clusters were visualized using the Seurat and ggplot2 packages (20), with custom color palettes applied for enhanced visualization. Cell type annotation was conducted through manual curation based on established marker genes and cluster characteristics. Gene expression visualization was performed using feature plots and dot plots to examine the expression patterns of DC markers including CLEC9A, XCR1, BATF3, CADM1, CD1C, FCER1A, CLEC10A, LILRA4, CLEC4C, and IRF7 across different clusters. Density plots for genes of interest (CD81, KCTD14, GBP1, and MYEOV) were generated using the Nebulosa package to visualize spatial expression patterns. For cell-cell communication analysis, epithelial cells were stratified based on KCTD14 expression status into KCTD14-positive and KCTD14-negative subgroups. Cell-cell communication networks were analyzed using the CellChat package with the human CellChatDB database. A subset of 3,000 randomly selected cells was used for computational efficiency. We computed interaction probabilities (cell threshold = 0; pathway threshold = 1), visualizing networks via circular plots. Pathway-specific analyses (TNF, TGFβ, SPP1, SEMA4, and MIF) included heatmaps, ligand–receptor contributions, and top interacting pairs. Interactions with P < 0.01 were deemed significant. All visualizations and statistical analyses were performed using R software with appropriate packages including Seurat, ggplot2, Nebulosa, CellChat, and ktplots.

Human pancreatic ductal epithelial cells (H6C7) and pancreatic cancer cell lines (CAPAN-1, PANC-1) were obtained from the Cell Bank of the Chinese Academy of Sciences (Shanghai, China). Cells were cultured in DMEM supplemented with 10% fetal bovine serum (FBS) and 1% penicillin-streptomycin at 37°C with 5% CO₂ (21). For KCTD14 knockdown, cells at 60–70% confluence were transfected with small interfering RNAs (siRNAs) targeting KCTD14 (si−KCTD14−1, si−KCTD14−2, si−KCTD14−3) or non−targeting control siRNA using Lipofectamine 3000 (Invitrogen). After 48 hours, cells were harvested for downstream assays.

Total RNA was extracted with TRIzol reagent (Invitrogen) and reverse−transcribed using the PrimeScript RT kit (Takara). qRT−PCR was conducted with SYBR qPCR Master Mix (Vazyme, China). Primer sequences were as follows:

KCTD14 (Forward): 5′−AGCAAACATGAACAGTAGGTTA−3′,

KCTD14 (reverse): 5′−GGAATGAAGGTAAGCAAAC−3′. Relative gene expression was calculated by the 2^(⁻ΔΔCt) method with β-actin as internal reference.

Cells were lysed in RIPA buffer supplemented with protease and phosphatase inhibitors (Roche). Protein concentrations were quantified by BCA assay (Pierce). Equal amounts of protein were separated on 10% SDS–PAGE gels and transferred onto PVDF membranes (Millipore). Membranes were blocked in 5% non−fat milk and incubated overnight at 4°C with primary antibodies against KCTD14, TNF−α, TNFRSF1A, TNFRSF1B, and β-actin. After washing, membranes were incubated using HRP−conjugated secondary antibodies and bands were visualized by ECL substrate (Thermo) and quantified by ImageJ (22).

Cell proliferation was assessed with the CCK-8 assay by seeding 3 × 10³ transfected cells per well in 96-well plates, with viability measured at 24, 48, 72, and 96 hours. Colony formation assays were performed by seeding 500 cells per well in six-well plates and allowing growth for 10–14 days before fixation and crystal violet staining. Colonies containing ≥50 cells were counted. Migration and invasion were examined using 8-μm pore Transwell chambers, with invasion assays performed in Matrigel-coated inserts. After incubation (24 h for migration, 48 h for invasion), cells on the lower membrane were fixed, stained, and quantified in five random fields (23, 24).

All statistical analyses were performed with R software version 4.0.1. Data were compared using Student’s t-test or Wilcoxon rank-sum test, as appropriate. Survival analyses were conducted with the Cox proportional hazards model and Kaplan–Meier method, with significance assessed by the log-rank test. Predictive accuracy was evaluated using time-dependent receiver operating characteristic (ROC) curves. For in vitro assays, all experiments were performed in three independent biological replicates, with each assay including three technical replicates unless otherwise specified. P value < 0.05 was considered statistically significant.

WGCNA was applied to three independent PDAC cohorts to identify modules most closely associated with DC infiltration. In the TCGA cohort, the turquoise module (1746 genes) showed the strongest association with DC infiltration (correlation coefficient = 0.87, P < 1e-200) ( Figures 1A, B ). Similarly, in GSE62165, the brown module (743 genes) correlated most highly with DC infiltration (correlation coefficient = 0.91, P < 1e-200) ( Figures 1C, D ). In GSE85916, the green module (1926 genes) showed a moderate but significant association with DC infiltration (correlation coefficient = 0.32, P = 8.2e-5) ( Figures 1E, F ). Taking the intersection across these three modules, we identified a set of 130 DC–related genes ( Figure 1G ). GO enrichment analysis of these 130 genes revealed significant enrichment in pathways such as leukocyte chemotaxis, immune response–regulating signaling pathway, and lymphocyte differentiation ( Figure 1H ). KEGG enrichment analysis further revealed enrichment in pathways such as cytokine-cytokine receptor interaction, chemokine signaling pathway, T cell receptor signaling pathway, and TNF signaling pathway ( Figure 1I ). Protein–protein interaction (PPI) network construction identified ITGAL, GAS7, and STAB1 as hub nodes with high connectivity ( Supplementary Figure S1 ), suggesting their central roles in DC recruitment and activation.

Using unsupervised consensus clustering based on the 130 DC–related genes, we stratified TCGA-PAAD patients into two distinct molecular clusters ( Figures 2A, B ). Principal component analysis confirmed clear separation between them ( Figures 2C, D ). The heatmap demonstrated generally higher expression of DC–related genes in cluster 2 compared to cluster 1 ( Figure 2E ). Survival analysis revealed that patients classified into cluster 1 demonstrated significantly improved overall survival (OS) compared to those in cluster 2 (P = 0.044, Figure 2F ). Differential expression analysis between the two clusters identified 1,827 DEGs (adjusted P < 0.05, |log₂FC| > 1, Figure 2G ). Subsequent GSEA indicated that these DEGs were particularly enriched in immune-associated pathways, such as IL6–JAK–STAT3 signaling, TNF-α signaling via NF-κB, and IL2–STAT5 signaling, underscoring functional divergence in immune regulation between the two clusters ( Figure 2H ).

To construct a prognostic model, we intersected the 1,827 cluster–associated DEGs with the 1,943 DEGs identified between tumor and adjacent normal tissues in TCGA-PAAD, yielding 771 candidate genes ( Figure 3A ). The TCGA cohort was randomly split into equal training and validation sets. In the training set, univariate Cox regression identified 22 genes significantly associated with OS (P < 0.05, Supplementary Figure S2 ). These candidates were further narrowed down using LASSO regression and stepwise regression analysis, resulting in a four-gene prognostic signature ( Figures 3B, C ). The DCRS was calculated using the formula: Risk score = CD81 × -0.806289828 + KCTD14 × 0.739669347 + GBP1 × 0.415429681 + MYEOV × 0.17556951. Using the medium value of the risk score, patients were divided into high- and low-risk groups. In the TCGA training cohort, patients in the low-risk group demonstrated significantly better OS compared to the high-risk group (P = 0.005, Figure 3D ). The risk curves also showed higher survival probability in the low-risk group ( Figure 3E ). Similarly, in the TCGA validation cohort, the low-risk group showed significantly better OS compare to high-risk group (P = 0.006, Figure 3F ), with risk curves corroborating superior survival probability in the low-risk group ( Figure 3G ). Time-dependent ROC curves demonstrated that area under the curve (AUC) values of the DCRS for OS in the TCGA training cohort at 1-, 3-, and 5-year were 0.727, 0.694, and 0.824, respectively ( Figure 3H ). Despite a slight decrease in the AUC values of the TCGA validation cohort, the results also demonstrated relatively good predictive potential for the short- and long-term OS of PC patients (1-year:0.651; 3-year: 0.684; 5-year:0.671, Figure 3I ).

To elucidate cellular sources of the four signature genes, we analyzed single-cell RNA-seq data from PDAC tumors, clustering cells into 23 populations ( Supplementary Figure S3 ). Cell type annotation identified 13 major cell clusters, including acinar cells, B cells, DCs, endothelial cells, erythrocytes, fibroblasts, malignant cells, mast cells, monocytes, NK/T cells, plasma cells, platelets, and stellate cells ( Figure 4A ). CD81, KCTD14, and MYEOV were predominantly expressed in malignant clusters, whereas GBP1 localized primarily to DC clusters ( Figure 4B ).

In the risk model, KCTD14 carried the largest positive coefficient, suggesting the strongest prognostic weight. CD81 is a canonical exosomal tetraspanin involved in extracellular vesicles (EV) biogenesis and immune modulation across solid tumors. PDAC EV studies indicate that circulating EV features correlate with survival or treatment response, underscoring a biomarker-oriented role rather than a single-gene driver in tumor–DC crosstalk (25–27). GBP1, a prototypical IFN-induced GTPase, is enriched in myeloid/immune cells and mediates innate immune programs (28, 29). In line with this biology, our single-cell analysis localized GBP1 predominantly to DC clusters rather than malignant epithelial clusters. In contrast, MYEOV already has PDAC-specific mechanistic evidence showing it enhances SOX9 transactivity to elevate HES1 and drive tumor progression (30), which reduces incremental novelty for the present study’s mechanistic focus. Therefore, we centered subsequent analyses on KCTD14 in PDAC. The cell–cell communication analysis revealed significant negative correlations between KCDT14 ⁺ epithelial cells and DCs ( Figures 4C, D ). Pathway-specific analyses showed that the connection between KCDT14 ⁺ epithelial cells and DCs is related to a series of signaling pathways, including TNF-TNFRSF1A, SPPA-CD44, and SEMA4G-PLXNB2 ( Figure 4E ). Given that prior GSEA, GO, and KEGG results consistently highlighted TNF signaling among DC-related gene programs ( Figures 1H, I ; Figure 2H ), we prioritized TNF pathway for in-depth investigation. The results confirmed that the TNF-TNFRSF1A pathway plays a crucial role in the connection between KCDT14 ⁺ epithelial cells and DCs ( Figures 4F–I ).

KCTD14 expression and function were validated in H6C7, CAPAN-1, and PANC-1 cell lines. According to the Western Blot and qRT-PCR assays, the expression of KCTD14 was significantly up-regulated in CAPAN-1 and PANC-1 cell lines compared with H6C7 cell line ( Figures 5A, B , **P < 0.01, ***P < 0.001). Among three siRNAs, si-KCTD14–3 achieved the highest knockdown efficiency and was therefore chosen for subsequent experiments. In both CAPAN-1 and PANC-1 cell lines, the expression of KCTD14 was significantly reduced after si-KCTD14–3 transfection ( Figure 5C , **P < 0.01). Functionally, the cell viability was significantly reduced following KCTD14 knockdown in both cell lines ( Figure 5D , **P < 0.01, ***P < 0.001). Clonal formation assays demonstrated significantly fewer colonies in both cell lines after KCTD14 knockdown ( Figures 5E, F , **P < 0.01). In Transwell assays, the numbers of migrated and invaded cells significantly decreased in both cell lines following KCTD14 knockdown ( Figures 6A, B , **P < 0.01). The Western Blot results confirmed lowered protein levels of TNF-α, TNFRRSF1A, and TNFRRSF1B following KCTD14 knockdown in both cell lines ( Figure 6C ). Collectively, KCTD14 knockdown significantly reduced PC cell proliferation, colony formation, migration, and invasion, and downregulated TNF pathway–related proteins.