We retrieved 179 samples from patients with PDAC who underwent RNA sequencing from the Cancer Genome Atlas (TCGA) database (phs000178 from https://www.omicsdi.org/dataset/dbgap/phs000178), which included clinical information, transcriptome expression, and copy number variation (CNV) data, and 171 samples of normal tissues from volunteers from the Genotype-Tissue Expression (GTEx) Project. Additionally, we obtained the pancancer cohort for the gene expression profiling and interactive analyses of the cancer and normal samples from the Gene Expression Profiling Interactive Analysis (GEPIA) and GEPIA2 platforms developed by PUCH (http://gepia2.cancer-pku.cn/, Beijing, China), as well as eight datasets (GSE111627; 141017; 148673; 154763; 154778; 158356; 162708; 165399) from patients with PDAC patients from the GEO database (https://www.ncbi.nlm.nih.gov/geo/, National Institutes of Health (NIH), MD, USA). Moreover, a panel of 223 genes related to efferocytosis was explored in the Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb), a collection of annotated gene sets for use with the gene set enrichment analysis (GSEA) software (version 4.4) by the USCD (SanDiego, CA, USA).

The single-cell RNA-seq dataset from the GEO database (GSE154778, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE154778) comprises 10 samples from patients with primary PDAC and 6 samples from patients with metastatic lesions. The Seurat R package (version 4.3.3, https://cran.r-project.org/web/packages or https://www.bioconductor.org/) was used to analyze the single-cell RNA sequencing data. While, the Tumor Immune Single-Cell Hub 2 (TISCH2) database (http://tisch.compbio.cn/) provided detailed cell type annotations at the single-cell level. Cell clustering was performed via Seurat’s “FindClusters” and “FindNeighbors” functions, and the results were visualized via t-SNE mapping. Cellular annotation was conducted by referencing signature genes associated with specific cell types. Additionally, we utilized the “Add Module Score” function from the Seurat package to assess the activity levels of specific gene sets within individual cells.

We identified 39 of the total-167 genes (Supplementary Table 1) associated with efferocytosis that exhibited prognostic significance for overall survival (OS) and disease-free survival (DFS). Seven genes were subsequently found to be significantly associated with both OS and DFS. A panel of five genes was ultimately selected from these seven genes. The top ten genes were analyzed for protein‒protein interactions via the STRING database (https://cn.string-db.org/, version 12.0) to highlight their enrichment. “Consensus Cluster Plus” package (version 1.70.0, https://cran.r-project.org/web/packages or https://www.bioconductor.org/) was used for consistency clustering to categorize PDAC patients into two subgroups. Significant differences in the clinical and prognostic characteristics were observed between the two subgroups of PDAC patients across the two independent cohorts.

To explore the biological pathways and processes associated with efferocytosis-related genes (ERGs) and models (ERMs), we used the Cluster Profiler R package (version 4.10.1, https://cran.r-project.org/web/packages or https://www.bioconductor.org/) for Gene Ontology (GO, https://geneontology.org/), Kyoto Encyclopedia of Genes and Genomes (KEGG, https://www.kegg.jp/), and gene set enrichment (GSEA, https://www.gsea-msigdb.org/gsea/index.jsp) analyses. These analyses identified key functions, pathways, and molecular interactions involved in efferocytosis. Additionally, we applied gene set variation analysis (GSVA, https://www.bioconductor.org/packages/devel/bioc/html/GSVA.html) to reveal activated or suppressed pathways involved in efferocytosis.

All efferocytosis-related genes were simultaneously included as predictors in a Cox proportional-hazard model with an L1 penalty (LASSO). The procedure was run 100 times; in each run a ten-fold cross-validation identified the penalty parameter ($\lambda$) that minimized the cross-validated partial-likelihood deviance, after which genes with non-zero coefficients were recorded. Genes retained in $>$50% of the 100 iterations were defined as stable predictors with only four genes meeting this stability criterion: ASPH, EMP1, FNDC3B and LGALS3. Their mean $\beta$-coefficients across the iterations were taken as weights, yielding the final four-gene risk-score formula. Based on this signature, we constructed an ER risk scoring system, the ERGRS, via the following formula: Risk Score = ASPH $\times$ (0.0359) + EMP1 $\times$ (0.131) + FNDC3B $\times$ (0.0372) + LGALS3 $\times$ (0.3).

To establish a robust predictive framework, we implemented a random forest-based algorithm for the feature selection. Each forest was grown with ntree = 1000, nodesize = 5, bootstrap sampling (samptype = “swr”), and $\text{mtry}\approx \sqrt{p}$ (where p is the number of candidate genes), while the out-of-bag (OOB) error was tracked to confirm convergence. The permutation-based variable importance (VIMP) was extracted from every run, and the mean VIMP over 100 forests was taken as the final importance score. Genes were ranked by this average VIMP, and scores were accumulated in descending order until 90% of the total importance was reached, yielding four top-contributing genes. This process identified the four prognostic biomarkers with significant discriminatory power in PDAC.

Human pancreatic cancer cell lines (ASPC-1 and BxPC-3) were obtained from the American Type Culture Collection (ATCC, Manassas, VA, USA), and cultured in Dulbecco’s modified Eagle’s medium (DMEM) supplemented with 10% fetal bovine serum (FBS), 100 U/mL penicillin, and 100 mg/mL streptomycin (all products purchased from Thermo Fisher Scientific, Chicago, IL, USA). All cell lines were authenticated via short tandem repeat (STR) DNA fingerprinting and tested negative for mycoplasma. Small interfering RNAs (siRNAs) specifically targeting ASPH (siRNA1 5^′-TCACGTGGTTTATGGTGAT-3^′ and siRNA2 5^′-TGTGGATGATGCCAAAGTT-3^′) and a control interfering RNA were obtained from Guangzhou RiboBio Co., Ltd., Guangzhou, Guangdong, China. For the transient transfection, ASPC-1 cells were transfected with the siRNA using Lipofectamine 2000 (11668019, Invitrogen™, Waltham, MA, USA) for 12 hours following which functional assays and the subsequent experiments were conducted.

We obtained tissue samples from patients with primary PDAC patients were obtained from who had undergone duodenopancreatectomy or pancreatospenectomy at the Department of General Surgery, Beijing Chaoyang Hospital. The collection and utilization of these tissues were conducted following the protocol approved by the Ethics Committee of the Beijing Chaoyang Hospital of Capital Medical University (approval No. 2018-K-99), with informed consent obtained from all participants. Paraffin-embedded 5 µm thick sections were prepared for multiplex immunohistochemical staining. Specifically, the sections were deparaffinized and subjected to heat-induced epitope retrieval in Tris buffer (pH 8.0). A four-marker panel consisting of EREG (Epiregulin, 12048), SMA (Smooth Muscle Actin ($\alpha$-SMA), 19245), PanCK (Pan-Cytokeratin, 4545), CD206 (Mannose Receptor (MRC1), 24595) (All products are sourced from Cell Signaling Technology, Inc., located in Danvers, MA, USA.), and DAPI was employed. Horseradish Peroxidase (HRP)-conjugated anti-rabbit/mouse TSA (Cat#abs50015, Absin, Shanghai, China) was used as the secondary antibody. The visualization of each biomarker was achieved via tyramine signal amplification-conjugated fluorophores (PerkinElmer, Naperville, IL, USA). Overlapping regions of interest were analyzed per slide according to previously established protocols. For the IHC analysis, a total of eight pairs of paraffin-embedded PDAC tissues and adjacent normal tissues were collected to assess the ASPH expression. The tissue slides were processed through deparaffinization, rehydration, and peroxidase quenching. 3,3’-Diaminobenzidine Tetrahydrochloride (DAB) staining and hematoxylin counterstaining were performed, and the slides were prepared for microscopic examination as described in our previous report. The staining intensity and percentage of positive cells were scored.

Total cellular RNA was extracted using TRIzol reagent and subsequently reverse-transcribed into cDNA using Moloney murine leukemia virus (M-MLV) reverse transcriptase (28025013, Invitrogen™, Waltham, MA, USA). Quantitative PCR (q-PCR) was then carried out with SYBR Green PCR Mix (A25742, Applied Biosystems, Foster City, CA, USA) on an ABI 7500 Fast Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). The relative gene expression levels were calculated using the $\mathrm{\Delta }$Ct method: $\mathrm{\Delta }$Ct = Ct (target gene) – Ct (GAPDH). The primer sequences used were as follows: Snail: Forward Primer TCGGAAGCCTAACTACAGCGA; Reverse Primer AGATGAGCATTGGCAGCGAG; Slug: Forward Primer CGAACTGGACACACATACAGTG; Reverse Primer CTGAGGATCTCTGGTTGTGGT; Twist: Forward Primer GCCTAGAGTTGCCGACTTATG; Reverse Primer TGCGTTTCCTGTTAAGGTAGC. All primers were synthesized by Beijing Sangon Biotech (Beijing, China). q-PCR data are presented as Mean $\pm$ Standard Deviation (SD) from three independent experiments.

The cell viability and proliferation were assessed via the Cell Counting Kit-8 (CCK-8) reagent (CK04, Dojindo Molecular Technologies, Inc., Tokyo, Japan) following the manufacturer’s instructions. The treated cells were evaluated at 12, 24, 48, and 72 hours post-transfection. The absorbance at 490 nm was measured via a multimode microplate reader (BMG Labtech ClarioStar, Cary, NC, USA).

The ATP Content Assay Kit (BC0300, Solarbio, Beijing, China) was employed to quantify intracellular ATP production in accordance with the manufacturer’s instructions.

The cells were collected and washed, followed by incubation with propidium iodide (PI) staining solution (P1304MP, Invitrogen™, Waltham, MA, USA) and 10 µL of permeabilization solution for 30 minutes at room temperature in the dark. The cell cycle distribution was subsequently analyzed via flow cytometry (BD Biosciences, San Diego, CA, USA).

The results of the cell proliferation assay are presented as means $\pm$ standard deviations. The Kaplan–Meier survival curves were generated using the survminer R package (version 0.4.9, https://cran.r-project.org/web/packages or https://www.bioconductor.org/). A Wilcoxon signed-rank test was employed to compare the significant differences between the two subgroups. Differences between the groups were determined via the two-tailed Student’s t-test or the $\chi$² test with Fisher’s exact probability method via GraphPad Prism software version 7.0 (GraphPad Software, Inc., Boston, MA, USA). For multiple tests, the Bonferroni‒Holm procedure was applied. The results with *, p 0.05, **, p 0.01, ***, p 0.001, ****, p 0.0001 were considered statistically significant.

Using the RNA-seq dataset of 179 PDAC samples from the TCGA database, we acquired 167 efferocytosis-related genes (ERGs) from the literature, and employed one-way Cox regression analysis to identify 21 and 18 prognostic genes linked to the OS and DFS of patients with PDAC, respectively (Fig. 1A). A comparison of the PDAC tumor samples and 171 normal tissue samples revealed that 5 (DUSP4, ASPH, MMP14, EMP1, and EREG) of these genes, which are associated with poor survival, were upregulated in tumors (Fig. 1B). In contrast, 2 genes (JUND and EPB41L3) associated with favorable survival were also upregulated in tumors (Fig. 1C). We subsequently utilized the same database to analyze a panel composed of the differential expression patterns of these five genes between PDAC tumor and normal tissues. Our analysis revealed significant differences in the expression levels between tumor and normal samples in this panel (Fig. 1D). The Kaplan-Meier (K‒M) survival analysis revealed a correlation between the PDAC patients’ OS and gene expression with this compound (Fig. 1E). KEGG analysis revealed that these genes are involved in pathways enriched in both tumor and normal tissues, as illustrated by the molecular functions, cellular components, and biological processes depicted in Fig. 1F. For example, the transferase activity; cytokine activity; and extracellular-integrin complex and ErbB1 signalling were enriched in tumor samples, whereas GTPase activity; lysosome, exosome and VEGF; and the VEGF receptor network were enriched in normal samples. Through a protein‒protein interaction (PPI) analysis via the STRING database, we identified the top 10 genes that strongly correlated with the 5 target genes in pancreatic cancer cells, as evidenced by their high pearson correlation coefficient (PCC) scores, as shown in Fig. 1G.

Fig. 1.

Identification of 39 hub genes related to efferocytosis. (A) The significant prognostic value of OS (11 red boxes, 10 blue boxes) and DFS (13 red boxes, 5 blue boxes) upon ERGs in PDAC. (B,C) Box plots depicting upregulated (B) and downregulated (C) genes both in OS and DFS. *, p 0.05. (D) Box plots depicting the mixed RNA expression of 5 poor survival-related genes in PDAC tumor and normal tissues. *, p 0.05. (E) Kaplan-Meier survival curves illustrating that the 5 genes correlated with poor survival in PDAC patients. (F) KEGG analysis of the 5 genes involved in these pathways. (G) Top 10 genes correlated with the 5 genes in PDAC according to the PCC. (H) LGALS3 was positively correlated with the 5 genes. (I) Venn diagram showing each of the top 100 genes related to LGALS3 in PDAC tumor and normal tissues. (J) Forest plot of these genes. (K) Distribution and mutation frequency of ERGs in the TCGA-PDAC cohort. (L) The frequency of CNV alterations in ERGs, with bar height indicating mutation frequency. ERGs, efferocytosis-related genes; PDAC, pancreatic ductal adenocarcinoma; OS, overall survival; DFS, disease-free survival; KEGG, Kyoto Encyclopedia of Genes and Genomes; PCC, pearson correlation coefficient; TCGA, the Cancer Genome Atlas; CNV, copy number variation.

Given its established role as a strong prognostic marker for PDAC, as previously reported by our team and others [15, 16], LGALS3, a member of the galectin family, is positively correlated with the expression of these five signatures in both tumor and normal tissues (Fig. 1H). We further identified the top 100 genes in PDAC tumor and normal tissues that were correlated with LGALS3. As illustrated in the Venn diagram (Fig. 1I), RTN4, MAP4K4, PANX1, PLS3, FNDC3B, and RBMS1 exhibited significantly greater correlations with tumor samples; ITGB1, ADAM17, and PICALM were correlated with both tumor and normal samples; notably, only PPP3R1 was highly correlated with normal samples alone. Next, a univariate Cox regression analysis was performed to generate a significant forest plot for 14 genes, except for RTN4 (Fig. 1J). Among the 173 patients with PDAC, ER-related gene mutations were identified in 2.31% (4 patients). RBMS1 presented the highest mutation frequency, followed closely by ADAM17 and ASPH, as illustrated in Fig. 1K. Additionally, we observed varying degrees of DNA copy number variation (CNV) among the ERGs. As shown in Fig. 1L, RBMS1 displayed extensive CNV amplification, whereas some ERGs presented CNV depletion specifically.

The principal component analysis (PCA) revealed that the distribution of the aforementioned 14 genes (“14-gene panel”) in all the pancreatic cancer samples presented two subpopulations (Fig. 2A). We successfully classified the PDAC patients into two distinct types by employing unsupervised clustering of this panel. The optimal number of clusters, k = 2, was used to divide the entire cohort into Cluster C1 (consisting of 111 patients) and Cluster C2 (consisting of 33 patients) (Fig. 2B,C). The Kaplan‒Meier survival analysis demonstrated significantly improved prognosis in Cluster C2 patients compared with those in the Cluster C1 patients (p = 0.015), as shown in Fig. 2D. Notably, the RNA level of each component in the panel was weakly expressed in the Cluster C2, as shown in Fig. 2E. To further investigate the biological characteristics of these two subtypes, we conducted GSVA. The C1 subtype was significantly enriched in oncogenic pathways, including cell mitosis, protein secretion, focal adhesion, and epithelial‒mesenchymal transition (EMT). In contrast, the C2 subtype was enriched in metabolic pathways, including retinol metabolism, diabetes mellitus, and olfactory transduction signalling (Fig. 2F). The GSEA using hallmark and ontology gene sets from the Molecular Signatures Database (MSigDB) further revealed the significant enrichment of the EMT-related pathways in the C1 subgroup, whereas the C2 subgroup exhibited enrichment in the large ribosomal subunit pathways, as illustrated in Fig. 2G,H. The PDAC patients were stratified into three groups according to the “14-gene panel” expression signature, as shown in the volcano plot. The genes encoding YAP1, ITGAV, MMP7, CST6, and COL11A1 were significantly upregulated (Fig. 2I). We then conducted a KEGG functional enrichment analysis of differentially expressed genes (DEGs), which revealed that high expression of the “14-gene panel” signature is associated with the assembly and activation of the focal adhesion, phosphatidylinositol-3-kinase- protein kinase B (PI3K-AKT), mitogen-activated protein kinases (MAPK), and hippo signaling pathways. In contrast, the downregulated genes are involved in pathways related to protein and fat digestion and absorption (Fig. 2J,K).

Fig. 2.

Two types of PDAC patients were categorized by “14-gene panel” clustering. (A) The scatter plot displays the gene distribution and subgrouping of PDAC patients. (B) The TCGA-PDAC cohort was clustered into C1 and C2 molecular subgroups on the basis of 14 ERGs. (C) The clustering results were best when k = 2 from 2–5 cumulative distribution function values. (D) K‒M survival analysis of the prognosis of C1 and C2 patients. (E) ERG expression levels varied between the two clusters. (F) GSVA revealed the hallmark pathways distinguishing clusters, with red indicating promotion and blue signifying inhibition. (G,H) GSEA enrichment analysis, including hallmark gene sets (G) and ontology gene sets (H). (I) Volcano map of “14-gene panel”-associated genes in the TCGA-PDAC cohort. (J,K) KEGG analysis of upregulated (J) and downregulated (K) genes among the DEGs. K‒M, Kaplan-Meier; GSVA, gene set variation analysis; GSEA, gene set enrichment analysis; DEGs, differentially expressed genes.

It was necessary to explore the potential prognostic significance of the “14-gene panel” signature owing to its limited prognostic clustering in the TCGA cohort, as demonstrated by the receiver operating characteristic (ROC) analysis, with the area under the curve (AUC) values of 0.42, 0.34, and 0.08 at 1 year, 3 years, and 5 years, respectively (Supplementary Fig. 1A), and as further confirmed via the univariate and multivariate Cox regression analyses, which indicated that the “14-gene panel” was not a significant independent prognostic predictor for PDAC (Supplementary Fig. 1B,C). To this end, least absolute shrinkage and selection operator (LASSO) regression and multivariate Cox modelling were applied to develop ER-associated prognostic signatures. We performed tenfold cross-validation on these genes through 100 iterations, identifying the four key genes—LGALS3, EMP1, ASPH, and FNDC3B—termed the “LEAF” signature (Fig. 3A) with the best $\lambda$ (–2.96) (Fig. 3B). We further determined the selection frequency of this model, which revealed the “LEAF” signature as the most critical gene among the “14-gene panel” genes (Fig. 3C), which were scored based on their stability and contribution (Fig. 3D). These genes demonstrated promising early-stage diagnostic potential in PDAC. The PDAC patients were categorized into high-risk and low-risk groups based on the median risk score. The K‒M analysis revealed that the patients in the low-risk group had better overall survival (Fig. 3E), and the ROC analysis of the “LEAF” signature in the TCGA cohort revealed robust prognostic discrimination, with AUC values of 0.75 at 1 year, 0.89 at 3 years, and 0.95 at 5 years, indicating an excellent temporal predictive performance (Fig. 3F). We compared the high-risk and low-risk groups via the “LEAF” signature, and generated a clinical risk factor heatmap that highlighted distinct “14-gene panel” expression patterns. Notably, all the “14-gene panel” genes presented increased expression in deceased patients (Fig. 3G). We explored the high and low expression levels of the LEAFs in the TCGA database and identified the biological processes and signalling pathways associated with each group. High “LEAF” expression correlated with processes such as the mitotic spindle, Notch signaling, the G2M checkpoint, and adhesion junctions, whereas low “LEAF” expression correlated with olfactory and diabetes metabolism, as shown in the heatmap (Fig. 3H). These results aligns with the signaling pathways identified via the “14-gene panel” in Fig. 2J,K. Finally, we employed both univariate and multivariate Cox regression analyses to validate cluster grouping as a standalone prognostic indicator for PDAC (Fig. 3I,J). These data indicate that the “LEAF” signature plays a central role in the prognostic capability of the “14-gene panel” model.

Fig. 3.

Development and verification of prognostic indicators related to “LEAF” genes. (A) Coefficient selection under LASSO regression analyses with “LEAF” genes. The vertical line is drawn at the optimal $\lambda$ value determined. (B) Tenfold cross-validation for parameter tuning adjusted through LASSO regression. (C) Feature selection frequency of “14-gene panel” genes in the TCGA-PDAC cohort through 100 deep learning steps. (D) Bar chart depicting the difference in the above “14-gene panel” gene scores in terms of their stability and contribution. (E) Survival curves generated via the Kaplan‒Meier method for the high-risk and low-risk groups. (F) Time‒dependent ROC curve analysis of the TCGA-PDAC cohort after deep learning. (G) Heatmap showing the clinical risk factors associated with the “LEAF” expression pattern. (H) Heatmap of biological processes or signaling pathways in the high-risk and low-risk groups. Red symbolizes advancement, whereas blue signifies restraint. (I,J) Through the use of “LEAF” genes, univariate (I) and multivariate (J) Cox model analyses of clinicopathologic factors and gene subtypes were performed. LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristic.

To elucidate the role of efferocytosis in the prognoses of patients with PDAC, we employed random survival forest (RSF) algorithms to strategically identify the key efferocytosis genes associated with patient outcomes. We identified four prognostic biomarkers with significant discriminatory power; namely, LGALS3, EREG, ASPH, and PLS3, also called “LEAP”. Integrating the LASSO analysis results, we identified LGALS3 and ASPH as the two pivotal genes involved in efferocytosis, which are significantly correlated with the prognoses of patients with PDAC (Fig. 4A). To validate the relationship between LGALS3 and ASPH, we established a genetically modified PDAC cell line overexpressing LGALS3, as shown in Fig. 4B, in which both the RNA and protein levels of LGALS3 are low in BxPC-3 cells. RNA sequencing revealed the gene distribution across distinct groups through a Venn diagram, intersecting four dimensions: (1) highly expressed genes in LGALS3-overexpressing BxPC-3 cells versus control cells; (2) the top 100 genes correlated with LGALS3 in tumor tissue; (3) the top 100 genes correlated with LGALS3 in normal tissue; and (4) a total of 167 efferocytosis-related genes (Fig. 4C). Interestingly, LGALS3 overexpression resulted in the identification of only one common gene expressed in both the tumor and normal samples: ITGB1, as shown in Fig. 1I. Additionally, only ASPH was found at the intersection of the tumor and RNA-seq groups. Importantly, EREG was not associated with the sequencing outcomes but emerged only at the intersection of the tumor group and the efferocytosis-associated genes, implying that the EREG expression is not dependent on tumor cells but is more likely to affect other tissue cells within the tumor microenvironment, and that ASPH is specific to tumor cells.

Fig. 4.

ASPH expression promotes PDAC progression. (A) Venn diagram showing “14-gene panel” expression in tumor and normal tissues selected by the LASSO regression and random survival forest (RSF) algorithms. (B) RNA and protein levels of LGALS3 in PDAC cell lines. The red arrow pointed out the low expression of ASPH in BxPC-3. (C) Venn diagram of the top 100 genes related to LGALS3 in tumors and normal tissues; sequencing of LGALS3 genetic cells compared with vector control cells; and identification of ER genes from the Molecular Signatures Database. These genes cross four distinct groups. (D) Expression level of ASPH across cancers. (E) Pathological ASPH staining images of a male pancreatic cancer patient (HPA059303) from the Human Protein Atlas (HPA) database. Scale bar = 200 µm. (F) Positive correlation between ASPH and LGASL3. (G) Relationship between ASPH and pathological grade in the Clinical Proteomic Tumor Analysis Consortium (CPTAC) samples. (H) Specific ASPH immunoreactivity between tumor and adjacent normal tissues by IHC. Scale bar = 200 µm. (I) The positive rate and intensity of ASPH expression were calculated. (J) The proliferation rates of PDAC cells with or without ASPH knockdown. (K) The efficiency of the RNA level of ASPH. (L) Flow cytometry analysis of the proportion of G2/M phase-arresting PDAC cells treated with or without ASPH knockdown. The red arrows denote distinct phases of the cell cycle, specifically the G0/G1 phase, the S phase, and the G2/M phase. (M) The detection of ATP production using ELISA assay with three independent experiments. **, p 0.01, ***, p 0.001 and ****, p 0.0001 were considered statistically significant.

The RNA levels of ASPH are significantly greater in tumor tissues than in adjacent nontumor tissues across various solid tumor types, including liver, esophageal, and pancreatic cancer tumors (Fig. 4D). The Human Protein Atlas database (https://www.proteinatlas.org/) provides pathological images of ASPH from a male patient with pancreatic cancer (Fig. 4E). A correlation analysis revealed a direct relationship between the ASPH and LGALS3 expression levels in tumor tissues (Fig. 4F), as well as a significant association between the ASPH expression and pathological grade (Fig. 4G). An immunohistochemistry analysis was performed to analyze the ASPH localization and expression levels in eight paired PDAC tissues and adjacent non-neoplastic pancreatic tissues to better understand the ASPH protein expression in clinical cases (Fig. 4H). The results demonstrated specific ASPH immunoreactivity in all the adenocarcinoma samples, with both the positivity and staining intensity exceeding 50%. In contrast, the adjacent normal tissues exhibited minimal detectable staining, with positivity below 20% and an intensity less than 50% (Fig. 4I). Consequently, we established the pancreatic cancer cell line ASPC-1 with ASPH knockdown (Fig. 4J) to further validate the ASPH effect on the mitotic spindle function and G2/M checkpoint, as described in Fig. 3H. As shown in Fig. 4K, the ASPH silencing reduced the tumor cell proliferation. A flow cytometry analysis revealed an increased proportion of cells in the G2/M phase following the ASPH downregulation (Fig. 4L). Furthermore, an enzyme-linked immunosorbent assay (ELISA) was conducted to quantify the increase in the ATP production following the ASPH knockdown (Fig. 4M).

We further identified seven primary clusters with specific marker genes in each cell population via single-cell sequencing data from the GEO database (GSE: 154778) (Fig. 5A), based on the aforementioned exploration of EMT hallmarks via the bioinformatics analysis of the “14-gene panel” characteristics described in Fig. 2F. The uniform manifold approximation and projection (UMAP) and heatmap analyses demonstrated that the EMT population was distinctly separated from the other cell populations (Fig. 5A,B). The ER module scores highlight the EMT population, which is characterized by high scores for both “14-gene panel” genes and “LEAF” genes (Fig. 5C). The ASPH expression was significantly increased in the epithelial tumor cell (ETC) and EMT populations (Fig. 5D). To further validate the functional phenotype of ASPH, we conducted the Boyden chamber assays using ASPH-knockdown cell lines and found that the cell migration was markedly reduced following ASPH silencing (Fig. 5E). In addition, classical molecular EMT markers were detected on ASPH-silencing cells via RT-PCR. The ASPH knockdown reduced the Snail, Slug, and Twist RNA levels remarkably (Fig. 5F). These data indicate that the ASPH expression in tumor cells drives their migratory characteristics via the EMT phenotype.

Fig. 5.

ASPH is related to migration through EMT. (A) Seven cell clusters were identified in the PDAC samples via a UMAP plot. (B) The heatmap indicates 3 marker genes in each cell population specific to different cell types. (C) Violin plot showing the efferocytosis module score for the “14-gene panel” cell cluster and “LEAF” group. (D) Feature plots illustrating the expression of ASPH at the single-cell level. (E) Boyden chamber assay showing the effect of ASPH knockdown on cell migration. Scale bar = 200 µm. (F) The RNA levels of EMT markers with ASPH knockdown cells with three independent experiments. **, p 0.01 and ****, p 0.0001 were considered statistically significant. ETCs, epithelial tumor cells; TAMs, tumor-associated macrophages; EMT, epithelial‒mesenchymal transition; CAFs, cancer-associated fibroblasts; Endos, endothelial cells; TILs, tumor-infiltrating lymphocytes; MASTs, mast cells; UMAP, uniform manifold approximation and projection.

Since the EREG gene from “LEAP” was identified in the RSF analysis, as illustrated in the Venn diagram in Fig. 4C, we further compared the EREG expression levels between tumor tissues and their adjacent normal tissues across various solid tumors, including colorectal cancer, renal cell carcinoma, and pancreatic cancer tumors (Fig. 6A). Notably, we observed a significant increase in the EREG activity within tumor-associated macrophages (TAMs) compared with that in the other cell types in GSE154778 (Fig. 6B). We subsequently acquired single-cell RNA sequencing data from patients with PDAC from the CRA001160 database, adhering to the methodology outlined in the TISCH2 database. By leveraging marker genes characteristic of distinct cell types, we successfully classified the cells into twelve primary clusters: acinar cells; B cells; CD8⁺ T cells; dendritic cells (DCs); ductal cells; endocrine cells; endothelial cells; fibroblasts, malignant cells; monocytes/macrophages; plasma cells; and stellate cells. The EREG expression was significantly elevated in the monocyte/macrophage subpopulation (Fig. 6C). To elucidate the role of EREG in monocytes and macrophages, we analyzed its distribution across monocytes, M1 macrophages, and M2 macrophages. The box plot in Fig. 6D shows that the EREG expression levels were significantly higher in macrophages than those in monocytes. Additionally, as shown in Fig. 6E, the EREG expression was notably greater in M2 macrophages than that in M1 macrophages. Finally, we observed the colocalization of EREG and CD206 in pancreatic cancer tissues via multiplex immunofluorescence staining (Fig. 6F). This experimental evidence further confirms that EREG plays a role in the function of tumor-associated immune cells, particularly M2 macrophages, and facilitates a more comprehensive evaluation of the “14-gene panel” system. The workflow of this study is illustrated in Fig. 6G. This method involved the identification of 14 critical risk genes from the ER database via OS/DFS data and PCC screening, followed by deep learning analysis. The primary focus was on the role of three target genes in tumor, stromal, and immune cells within the TME: ASPH, LGALS3, and EREG.

Fig. 6.

EREG affects M2 macrophage polarization. (A) The RNA expression level of EREG across cancers. (B) Violin plot showing the EREG score for each cell type. (C) Twelve cell clusters were identified via a UMAP plot of the PDAC samples and feature plots indicating the expression of EREG at the single-cell level. (D) Box line plot of the difference in EREG expression levels between monocytes and macrophages. (E) Comparison of the EREG expression levels between M1 and M2 macrophages. (F) Multiplex immunofluorescence staining showing the colocalization of EREG and CD206 in PDAC tissues. Scale bar = 200 µm. (G) The workflow of our research. *, p 0.05 and **, p 0.01 were considered statistically significant. EREG, Epiregulin; SMA, Smooth Muscle Actin ($\alpha$-SMA); PanCK, Pan-Cytokeratin; CD206, Mannose Receptor (MRC1); PCC, Pearson correlation coefficient; GEO, Gene Expression Omnibus.