Expression patterns and consensus unsupervised clustering of ICD genes in pancreatic cancer

Initially, the research examined the expression pattern of ICD genes in PC samples compared to normal samples. The study revealed significant differences in gene expression between PC and normal samples. The majority of ICD genes showed higher expression levels in cancer samples than in normal samples (Fig. 1A). Subsequently, utilizing consensus unsupervised clustering on the TCGA cohort, two distinct clusters (C1 and C2) were identified (Fig. 1B). Comparative analysis of these clusters revealed that the majority of ICD genes exhibited significantly elevated expression levels in cluster C1 as opposed to cluster C2 (Fig. 1C). Therefore, cluster C1 was categorized as the ICD-high cluster, while cluster C2 was designated as ICD-low cluster. Notably, survival analysis revealed a favorable prognosis for the low cluster of ICD (Fig. 1D, E).

Fig. 1

Identification of the ICD-related clusters in pancreatic cancer. A The expression patterns of ICD genes in pancreatic cancer samples and normal pancreas tissues. B The heatmap of consensus clustering solution (k = 2) in pancreatic cancer samples. C The expression of ICD genes in two clusters. D Kaplan–Meier curve of overall survival in two clusters. E Kaplan–Meier curve of progression free survival in two clusters. F The heatmap of differently expressed genes between two clusters. G The volcano plot of differently expressed genes between two clusters. H GO analysis of highly expressed genes in the ICD-high clusters. I KEGG enrichment analysis of highly expressed genes in the ICD-high clusters. J GO analysis of highly expressed genes in the ICD-low clusters. K KEGG enrichment analysis of highly expressed genes in the ICD-low clusters. ***P < 0.001; **P < 0.01; *P < 0.05

Functional enrichment analysis of DEGs in distinct clusters of ICD

Following this, we conducted differential analysis between the low and high clusters of ICD, finding 2627 genes with varying expression levels (Fig. 1F, G). Subsequently, enrichment analyses for GO, and KEGG were carried out on these genes. GO and KEGG analyses highlighted that DEGs upregulated in ICD-high cluster were enriched in regulation of T cell activation, leukocyte proliferation, extracellular matrix organization and Cytokine-cytokine receptor interaction (Fig. 1H, I). The DEGs upregulated in ICD-low cluster were enriched in signal release, neurotransmitter secretion, and neuroactive ligand-receptor interaction (Fig. 1J, K).

Correlation of different clusters and TMB

Numerous studies have consistently demonstrated a strong correlation between TMB and patient prognosis [21,22,23]. To examine potential disparities in TMB between two clusters, we collected and organized the mutation data. The analysis revealed KRAS, TP53, SMAD4, and CDKN2A as the most frequently mutated genes in both clusters (Fig. 2A, B). Survival analysis clearly showed that patients with higher TMB had a much worse prognosis than those with lower TMB (Fig. 2C). Subsequent combined survival analysis considering TMB and clusters disclosed that patients with higher TMB and higher expression levels of ICD genes experienced the most unfavorable prognosis, whereas patients with lower TMB and lower expression levels of ICD genes showed a more favorable prognosis (Fig. 2D). Given that KRAS and TP53 were the two genes exhibiting the highest mutation frequencies in PC, we investigated the influences of mutations in these genes on patient prognosis. Our initial findings revealed a substantial decrease in the prognosis of PC patients carrying KRAS mutations, particularly those with a combination of KRAS mutations and higher expression levels of ICD genes who displayed the most unfavorable prognosis (Fig. 2E, F). Similarly, our analysis indicated an association between TP53 mutations and prognosis, with patients possessing TP53 mutations and higher expression levels of ICD genes experiencing a less favorable prognosis (Fig. 2G, H).

Fig. 2

Tumor mutational burden (TMB) analysis in the two clusters. A, B Visual waterfall plots of mutated genes in two clusters. C Kaplan–Meier curve of H-TMB and L-TMB samples. D Kaplan–Meier curve stratified by clusters and TMB status. E Kaplan–Meier curve of KRAS-Mutant and KRAS-Wild samples. F Kaplan–Meier curve stratified by clusters and KRAS mutation status. G Kaplan–Meier curve of TP53-Mutant and TP53-Wild samples. H Kaplan–Meier curve stratified by clusters and TP53 mutation status

Analysis of the TIME in clusters with low and high levels of ICD

Increasing evidence suggests that tumor-infiltrating immune cells play a crucial role in the tumor microenvironment [24]. Thus, we investigated the variances in the TIME between the two clusters. At first, we examined the ESTIMATE score, immune score, stromal score, and tumor purity in both the ICD-high and ICD-low clusters. The results showed increased stromal, immune, and ESTIMATE scores in the ICD-high cluster, along with decreased tumor purity (Fig. 3A–D). Additionally, utilizing ssGSEA, we examined potential associations between distinct clusters and immune cell infiltration (Fig. 3E). The results showed the correlation between various types of immune cells in PC (Fig. 3F). Notable variances in the presence of different types of immune cells, including B cells, CD8 T lymphocytes, and cytotoxic lymphocytes, were observed among these clusters. Specifically, these immune cells were more abundant in the ICD-high cluster, indicating an enhanced immune response in the TIME of this cluster (Fig. 3G, H). Additionally, we investigated the differences in expression of HLA molecules within these clusters and found increased levels in the ICD-high cluster, indicating potential anti-cancer immune reactions (Fig. 3I). Finally, an analysis of the differences in the main molecules of immune checkpoints among the clusters showed increased levels in the ICD-high cluster, providing possible clues for predicting the reaction to immune checkpoint inhibitors (Fig. 3J).

Fig. 3

The differences of tumor immune microenvironment between two clusters. A–D Comparison of ESTIMATE score, immune score, stromal score, and tumor purity between the two clusters. E The distributions of infiltrating immune cells between the two clusters calculated by CIBERSORT algorithm. F The correlations of infiltrating immune cells in pancreatic cancer samples calculated by CIBERSORT algorithm. G The scores of infiltrating immune cells between the two clusters calculated by MCPcounter package. H The distinction of infiltrating immune cells between the two clusters calculated by MCPcounter package. I Comparison of the expression of HLA molecules between the two clusters. J Comparison of the expression of immune checkpoint key molecules between the two clusters. ***P < 0.001; **P < 0.01; *P < 0.05

Construction of the prediction model for ICD genes in PC

We conducted a univariate COX regression analysis on the TCGA dataset to create a predictive model for genes associated with ICD. This analysis revealed that five genes, namely CASP1, IFNG, IL1R1, MYD88, and PIK3CA, were significantly correlated with patient prognosis (Fig. 4A). Subsequently, following univariate Cox regression analysis, we performed LASSO regression analysis on these five genes to construct a prognostic model (Supplementary Table S1). The results revealed that the prognostic model included four genes (CASP1, IFNG, MYD88, and PIK3CA) (Fig. 4B, C). Individual risk scores were calculated for each sample using a designated formula. Subsequently, individuals were categorized into high-risk and low-risk groups based on the median risk score. As the risk score escalated, the high-risk group demonstrated a decrease in survival rates and an increase in mortality compared to the low-risk group (Fig. 4D, E). The heatmap indicated that the expression of these four genes was increased in patients classified as high-risk, implying a notable role of these genes in the advancement of PC (Fig. 4F). Survival analysis exhibited a strong association between the risk score and patient outcomes, indicating that individuals in the high-risk category had significantly poorer outcomes compared to those in the low-risk category (Fig. 4G). To validate the feasibility of the model, the GEO dataset was employed for external validation and to assess its accuracy. The results demonstrated that the validation set showed parallel performance to the training set. With the risk score increasing, the survival rates of patients decreased (Fig. 4H, I). Analysis of gene expression showed increased levels of these four genes in patients at high risk (Fig. 4J). Additionally, patients with a higher risk level had worse outcomes in comparison to patients with a lower risk level (Fig. 4K). Collectively, these findings supported the robust predictive capacity of the model and its potential for clinical application.

Fig. 4

Construction of the risk model based on ICD genes in pancreatic cancer. A The results of the Univariate Cox regression analysis conducted in TCGA-PAAD dataset. B Cross-validation plot for the results of LASSO regression analysis. C Plots for LASSO expression coefficients of the ICD genes. D The risk score of each sample in the TCGA cohort. E The risk score and survival status of each sample in the TCGA cohort. F The expression levels of four genes in TCGA samples. G Kaplan–Meier curve of overall survival between two groups in TCGA cohort. H The risk score of each sample in the GEO cohort. I The risk score and survival status of each sample in the GEO cohort. J The expression levels of four genes in GEO samples. K Kaplan–Meier curve of overall survival between two groups in GEO cohort

Assessment of predictive capacity: risk model versus other clinical characteristics

First of all, PCA was employed to assess the risk signature’s grouping capability. By combining whole-genome expression profiles with 34 ICD genes and a risk model, the analysis successfully divided samples into two distinct risk groups, demonstrating the risk model’s ability to effectively stratify samples (Fig. 5A–C). We next compared the prognostic prediction capabilities of the risk model with other clinicopathological characteristics through ROC curve analysis. The analysis validated the risk score as the most effective predictor among the variables studied (Fig. 5D–G). Moreover, these results were corroborated through the use of the c-index (Fig. 5H). Cox regression was used to evaluate the model’s individual influence on the outcome of PC patients. The results confirmed the risk model as a significant and independent prognostic factor, exerting a substantial influence on PC patient outcomes (Fig. 5I, J). The development of a nomogram ultimately offered a tangible and numerical approach for forecasting outcomes in individuals with PC (Fig. 5K). Calibration curves were used to assess the dependability and consistency of the nomogram, revealing a strong correlation between the projected and observed survival rates (Fig. 5L). These findings collectively indicated that the risk model developed in this study exhibited stable and reliable prognostic prediction capabilities for patients with PC.

Fig. 5

Validation of the stability of the risk model. A Principal component analysis of all genes. B Principal component analysis of ICD genes. C Principal component analysis of risk score. D The ROC curves of the risk model in the TCGA cohort. E–G The ROC curves of 1-, 3-, and 5-year overall survival for the risk score and other clinical characteristics. H The c-index for the risk score and other clinical characteristics. I The results of Univariate Cox regression analysis in TCGA cohort based on risk score and other clinical characteristics. J The results of Multivariate Cox regression analysis in TCGA cohort based on risk score and other clinical characteristics. K The nomogram established by risk score and other clinical characteristics in the TCGA cohort. L 1-, 3-, 5-year nomogram calibration curves. ***P < 0.001; **P < 0.01; *P < 0.05

Assessment of the validity and predictive significance of genes within the model

Based on the preceding results, a prognostic model for PC was constructed utilizing ICD genes. The model showed robust predictive ability and acted as an independent prognostic indicator for PC patients, including four distinct genes. In order to investigate the functions of these four genes in PC more comprehensively, their expression profiles were initially downloaded from the GEPIA (http://gepia.cancer-pku.cn/). The findings revealed elevated expression levels of CASP1, MYD88, and PIK3CA in PC tissues (Fig. 6A–D). Subsequently, the levels of expression of these genes in cancer cells were examined utilizing single-cell datasets. The findings indicated that in addition to their expression in immune cells, these genes demonstrated elevated expression in cancer cells (Fig. 6E–I). Furthermore, leveraging the TCGA dataset, we explored the influence of these four genes’ expression on the prognoses of patients with PC. The findings revealed that elevated expression of these four genes collectively exerted an adverse effect on patient prognosis (Fig. 6J–M). The discovery of these genes suggested that they were crucial in the advancement of PC.

Fig. 6

Validation of the expression and prognostic value of genes in the model. A–D Comparison of expression levels of four genes in tumor and normal samples. E The cell types revealed by the scRNA-seq dataset GSE111672. F–I The expression levels of four genes in different cell types. J–M Kaplan–Meier curves of four genes in the model. ***P < 0.001; **P < 0.01; *P < 0.05

Validation of the risk model in our own PC samples

To enhance the model’s accuracy, we validated it using our PC samples. Initially, we assessed the expression levels of the four genes in PC and adjacent normal tissues through qRT-PCR. Our findings revealed significant upregulation of CASP1, MYD88, and PIK3CA in PC compared to the adjacent normal samples, aligning with earlier results and highlighting the pivotal roles of these genes in PC progression (Fig. 7A–D). Subsequently, the expression levels of these four genes were integrated using the established risk score formula, enabling the categorization of patients into high-risk and low-risk groups based on their respective risk scores. Survival analysis results indicated that individuals in the high-risk group had significantly worse outcomes than those in the low-risk group, consistent with previous research findings (Fig. 7E). Furthermore, the two groups were compared in terms of patient stage, as well as CA19-9 and CEA levels. The examination showed that there was a greater percentage of patients in advanced stages III and IV within the group at high risk, as well as increased levels of CA19-9 and CEA (Fig. 7F–H). To confirm the expression levels of these genes in PC, we conducted IHC staining on both tumor tissues and adjacent normal tissues from PC patients, demonstrating elevated expression of these genes within the tumors (Fig. 7I). In summary, the risk model exhibited remarkable accuracy and robustness, demonstrating its efficacy not only within the database but also within our patient cohort.

Fig. 7

validation of the risk model in our own samples. A–D Comparison of expression levels of four genes between normal and tumor samples. E Kaplan–Meier curve of overall survival between two groups in our own samples. F The differences of tumor stage between two groups in our own samples. G, H The differences of tumor markers (CA19-9 and CEA) between two groups in our own samples. I Representative CASP-1, MYD88, and PIK3CA IHC staining images in normal and tumor samples. Scale bar, 25 μm. ***P < 0.001; **P < 0.01; *P < 0.05

MYD88 promoted PC cells proliferation and migration in vitro

Studies have demonstrated the significant involvement of MYD88 in tumorigenesis. However, its specific role in PC remains uncertain. Therefore, we investigated its function in the forthcoming experiments. We evaluated the impact of MYD88 on PC progression by examining the proliferation and migration of PC cells. Subsequent to the transfection of siRNA targeting MYD88, the efficiency of knockdown was validated at the mRNA and protein levels (Fig. 8A, B). In addition, two cell lines were selected for exogenous overexpression (Fig. 8C). The CCK-8 assay indicated that reducing MYD88 expression effectively suppressed in vitro cell proliferation (Fig. 8D). On the contrary, the overexpression of MYD88 led to the enhancement of their proliferation (Fig. 8E). Similarly, the colony-formation assay produced consistent results (Fig. 8F–H). Moreover, we explored the effect of MYD88 on the migration rate of PC cells. The results showed a notable decrease in the migratory capacity of PC cells following MYD88 knockdown, as well as a significant increase in their migratory ability following MYD88 overexpression (Fig. 8I–K). The above data showed that MYD88 had a significant impact on the growth and migration of PC cells.

Fig. 8

The effects of MYD88 on PC cells in vitro. A The mRNA knockdown efficiency of MYD88 in MIA PaCa-2 and Patu8988. B The knockdown efficiency of MYD88 in MIA PaCa-2 and Patu8988 validated by Western blot assay. C The overexpression efficiency of MYD88 in SW-1990 and PDC 0034 validated by Western blot assay. D The effect of MYD88 knockdown on the proliferation of PC cells demonstrated through the CCK-8 assay. E The effect of MYD88 overexpression on the proliferation of PC cells demonstrated through the CCK-8 assay. F–H The effects of MYD88 knockdown or overexpression on the proliferation of PC cells demonstrated through the colony-formation assay. I–K The effects of MYD88 knockdown or expression on the migration of PC cells demonstrated through the transwell assay. ***P < 0.001; **P < 0.01; *P < 0.05

MYD88 genetic inhibition suppressed the growth PC cells in vivo

To verify the influence of MYD88 genetic inhibition on the proliferation of PC cells in vivo, we established the subcutaneous xenograft model. The result showed that mice injected with MYD88 knockdown cells displayed a reduced tumor burden compared to the control group (Fig. 9A). In terms of tumor volume, the mice injected with NC cells exhibited larger tumors compared to the mice injected with shMYD88 cells (Fig. 9B, C). Additionally, immunohistochemical results also showed that the ki-67 positivity rate in the shMYD88 group was significantly lower than that in the NC group (Fig. 9D). Subsequently, we established the orthotopic xenograft model to corroborate the findings. Consistent with the subcutaneous tumor results, the tumor burden in the mice injected with sgMyd88 cells was significantly smaller than that in the mice injected with sgLacZ cells (Fig. 9E-G). The immunohistochemical results indicated that the ki-67and Myd88 positivity rates in the sgMyd88 group was significantly lower than that in the sgLacZ group (Fig. 9H). All these results indicated that MYD88 knockout significantly inhibits the proliferation of PC cells in vivo.

Fig. 9

PC cell proliferation is suppressed by genetic inhibition of MYD88 in vivo. A The protein level of MYD88 in shNC and shMYD88 Patu8988 cells validated by Western blot assay. B, C Growth of subcutaneous xenografts from shMYD88 cells. D Representative MYD88, ki-67 IHC and H&E staining images showed the role of MYD88 in vivo. E The protein level of Myd88 in sgLacZ and sgMyd88 KPC1199 cells validated by Western blot assay. F, G Growth of orthotopic xenografts from sgMyd88 cells. H Representative Myd88, ki-67 IHC and H&E staining images showed the role of MYD88 in vivo. ***P < 0.001; **P < 0.01; *P < 0.05