Bioinformatics analysis of PBLD as a potential biomarker in pan-cancer and colon adenocarcinoma

Introduction

Malignancies pose a critical threat to human health and are among the leading causes of global mortality (1-3). Despite diagnostic and therapeutic advances, 5-year overall survival (OS) rates remain suboptimal (4), highlighting the urgent need for novel biomarkers and therapeutic targets.

Our group’s previous multi-omics analyses integrating proteomics (5) and comparative genomics (6) revealed that phenazine biosynthesis-like protein domain-containing protein (PBLD) was significantly reduced in ulcerative colitis (UC). PBLD expression levels were negatively correlated with UC severity. Epithelial PBLD was shown to attenuate intestinal inflammatory response and enhance barrier integrity by suppressing nuclear factor kappa-light-chain-enhancer of activated B cells (NF-κB) signaling (7). Expanding on our findings, recent research established that vitexin targets the vitamin D receptor (VDR) and regulates macrophage polarization via the VDR/PBLD pathway, thereby impeding the progression from chronic colitis to colorectal cancer (8). PBLD, also known as mitogen-activated protein kinase activator with WD40 repeats (MAWD) binding protein (MAWBP), modulates transforming growth factor-beta (TGF-β) signaling (9), intestinal inflammation (7), and antiviral innate immunity (10). Elevated PBLD expression is documented in insulin resistance, folate deficiency, and hypertension (11). Previous studies have indicated tumor suppressor functions in gastric, hepatocellular, breast, and colorectal cancers (8,9,12,13), where PBLD inhibits proliferation, invasion, metastasis, and angiogenesis through the vascular endothelial growth factor (VEGF)/vascular endothelial growth factor receptor 2 (VEGFR2) pathway (14). However, previous studies of PBLD have been confined to individual cancer types or single disease models, lacking a comprehensive overview of its pan-cancer role. To bridge this gap, we designed a multi-level bioinformatic study to systematically characterize PBLD, following a logical sequence from expression to immune association, drug response, and metabolic pathway analysis.

Pan-cancer analyses enable the identification of conserved molecular features and accelerate the discovery of biomarkers. Using public databases [The Cancer Genome Atlas (TCGA) (15), Genotype Tissue Expression (GTEx)], this study comprehensively evaluated the roles of PBLD in gene expression, prognosis, genetic alterations, and interactions with tumor mutational burden (TMB), microsatellite instability (MSI), immune infiltration, drug sensitivity, and metabolic pathways. Given that PBLD deficiency promotes colitis-associated colorectal cancer (CAC), colon adenocarcinoma (COAD) can represent a biologically ideal model for investigating its tumor-related functions. This relevance, combined with the extensive, well-annotated TCGA-COAD dataset of matched tumor-normal pairs, solidified our selection of COAD for focused investigation following the pan-cancer screen. This multidimensional approach aims to delineate the potential of PBLD as a diagnostic/prognostic biomarker and a therapeutic target. We present this article in accordance with the REMARK reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1-2921/rc).

Methods

PBLD expression analysis

We used the Human Protein Atlas (HPA) database (16) (version: 25.0) (https://www.proteinatlas.org/) to analyze the expression of PBLD in normal tissues. Differential expression of PBLD between tumor and adjacent normal tissues across cancer types was analyzed using the “Gene DE” module of Tumor Immune Estimation Resource version 2.0 (TIMER2.0) (17) (https://compbio.cn/timer2/), which utilizes RNA-sequencing (RNA-seq) expression data from TCGA project. Default parameters were applied for analyses. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Clinical correlations

Prognostic significance of PBLD was evaluated using Gene Expression Profiling Interactive Analysis 2 (GEPIA2) (18) (http://gepia2.cancer-pku.cn/). OS and disease-free survival (DFS) analyses were analyzed across 33 TCGA cancer types, with hazard ratios (HR) and 95% confidence intervals (CIs) derived from the Cox proportional hazards model. Survival curves were generated for each TCGA cancer type.

Gene Set Cancer Analysis (GSCA, https://guolab.wchscu.cn/GSCA/) was employed to evaluate the prognostic significance of PBLD. Progression-free survival (PFS), disease-free interval (DFI), and disease-specific survival (DSS) were analyzed using the Cox proportional hazards model, with HR and 95% CI calculated for each cancer type in TCGA dataset.

Pathological stage association analysis was performed using the “Expression DIY” module of GEPIA2. PBLD expression levels across different tumor stages were compared by one-way analysis of variance (ANOVA), with log₂(TPM+1)-transformed expression data derived from TCGA projects.

Genetic alterations of PBLD, including mutation frequency and copy-number alterations (CNAs), were analyzed using the cBioPortal (19) (https://www.cbioportal.org/). All data were derived from the TCGA PanCancer Atlas Studies dataset. The “Cancer Types Summary” module was employed to calculate the frequency of PBLD mutations and CNAs across all cancer types. Prognostic significance of PBLD genetic alterations was assessed using the “Comparison/Survival” module. Only cancer types with a minimum of three patients in the altered group were included in the survival analysis.

Molecular markers and tumor microenvironment (TME)

The correlation between PBLD expression and TMB was analyzed with the R package “TCGAplot” (20) (version 4.2.1) using the Wei Sheng Xin online platform (https://www.bioinformatics.com.cn/plot_tcgaplot_pan_cancer_gene_expression_and_tmb_correlation_radar_chart_303), and the correlation data were downloaded by searching the target gene PBLD and using Pearson’s correlation, then select “Correlation Download” to obtain the data. The same procedure was followed for MSI correlation data using the “TCGAplot” package on Wei Sheng Xin platform (https://www.bioinformatics.com.cn/plot_tcgaplot_pan_cancer_gene_expression_and_msi_correlation_radar_chart_304). Subsequently, the downloaded data were organized into an Excel file formatted with headers in the first row, including Cancer types and P value, and correlation coefficients for each cancer type listed in the second row. This file was then uploaded to the “Radar Chart” tool on the Xiantao Academic platform (https://www.xiantaozi.com/). After adjusting image parameters such as dimensions and font size as needed, the correlation radar charts were generated and visualized using the R package “ggplot2”.

Immune cell infiltration was quantified for 12 immune cell types, and its correlation with PBLD expression was analyzed by single sample Gene Set Enrichment Analysis (ssGSEA) on Xiantao Academic platform using Spearman’s method. The “immune association” module of TIMER2.0 (https://compbio.cn/timer2/) was employed to estimate cancer-associated fibroblast (CAF) infiltration through Estimating the Proportions of Immune and Cancer cells (EPIC), Microenvironment Cell Populations-counter (MCP-counter), xCell and Tumor Immune Dysfunction and Exclusion (TIDE) algorithms.

Drug sensitivity analysis of PBLD

Drug sensitivity analysis for PBLD was conducted using the “Drug” module of the GSCA database. This module incorporates half-maximal inhibitory concentration (IC₅₀) data from the Genomics of Drug Sensitivity in Cancer (GDSC) (21) and Cancer Therapeutics Response Portal (CTRP) (22). For each cancer type, we performed Pearson correlation analysis between PBLD mRNA expression and drug IC₅₀ values across all relevant cell lines. The resulting P value were adjusted for multiple testing using the false discovery rate (FDR) method. We considered gene-drug pairs with an absolute correlation coefficient (|r|) >0.1 and FDR ≤0.05 to be statistically significant. To prioritize the most robust associations, a composite score was calculated for each significant pair by multiplying |r| by −log₁₀(FDR). All pairs were subsequently ranked in descending order based on this score. The top 30 ranked gene–drug pairs are presented in the bubble plot generated by the GSCA platform.

Functional analyses

The “Most Similar Genes” module of GEPIA2 (http://gepia2.cancer-pku.cn/) was employed to identify the top 100 genes with expression patterns most similar to that of PBLD across all tumor types in the TCGA database. Subsequently, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses for these similar genes were performed using the OmicShare online platform (https://www.omicshare.com/tools/). For GO enrichment, Gene IDs were converted to Ensembl IDs and formatted to meet the platform’s specifications, with the species designated as Homo sapiens. Enrichment was examined across three GO categories: biological process (BP), cellular component (CC), and molecular function (MF). Visualization of enrichment results was performed using bar and bubble plots. KEGG enrichment analysis was conducted following the same procedure.

A co-expression gene network for PBLD was constructed using STRING (23) (version 12.0) (https://string-db.org/). The network was generated by entering “PBLD” as the protein name and selecting Homo sapiens. Parameters were adjusted as follows: active interaction sources set to “Co-expression”, a minimum required interaction score of low confidence (0.150), and a maximum of 50 interactors displayed. All other settings were kept at their defaults. The resulting protein co-expression network figure was then obtained. Enrichment analysis for the co-expressed genes was again performed using OmicShare (24). The list of co-expressed gene names from STRING was queried in the HPA database; the corresponding Ensembl IDs were retrieved from the “GENE INFORMATION” section. These IDs were compiled into an Excel sheet and uploaded to the OmicShare platform for enrichment analysis.

Correlation analyses between PBLD and its co-expressed genes were performed using the “Correlation Analysis” module of GEPIA2. For each gene pair, the analysis was conducted by selecting the two genes of interest under the “Correlation Analysis” section, adding all TCGA cancer datasets to the analysis cohort, and generating the correlation scatter plot.

Expression of PBLD in COAD

For expression analysis of unpaired samples, we used the “[Cloud] Disease vs. Non-disease” module within the “Expression Difference” category on the Xiantao Academic platform. RNA-seq data from the TCGA-COAD dataset were selected as the data parameter, and “PBLD” was entered in the molecular field under special parameters. The results were visualized using boxplots.

For expression analysis of paired samples, we employed the “[Cloud] Paired Plot” module in the “Differential Expression” category on the Xiantao Academic platform. This analysis also utilized RNA-seq data from the TCGA-COAD dataset, with “PBLD” entered in the molecular field under the special parameters. All statistical analyses and visualizations were conducted using R package, with the ggplot2 package employed for data visualization.

Prognostic analysis of PBLD in COAD

The Xiantao Academic platform’s “Clinical Significance” module was used to access the “[Cloud] Survival Curve (KM Plot)” tool within the “Prognosis” category. RNA-seq data from the TCGA-COAD dataset were used as the data parameter, with “PBLD” entered in the molecular field under special parameters. The prognostic time parameter was set to “months”. The ggplot2 and survminer packages were employed for data visualization, and the survival package was used for survival data analysis.

Immune infiltration analysis of PBLD in COAD

The “[Immune Infiltration-Cloud] Lollipop Chart” module under the “Interaction Network” category on the Xiantao Academic platform was selected. RNA-seq data from the TCGA-COAD dataset were used as the data parameter, and “PBLD” was entered as the molecular parameter. We performed the analysis with the default ssGSEA algorithm and 24 immune cell types.

Subsequently, the “[Immune Infiltration-Cloud] Correlation Scatter Plot” module under the “Interaction Network” category was selected to generate scatter plots, which visualize the correlation between PBLD expression and immune cell types identified as significantly correlated. Immune infiltration levels were quantified via the ssGSEA algorithm, employing marker genes for 24 immune cell types as defined in a prior immunity study (25). The ggplot2 package was employed for data visualization.

Enrichment analysis of PBLD

The TCGA-COAD dataset was queried via the cBioPortal database. “Mutations”, “Putative copy-number alterations from GISTIC”, and “mRNA Expression FPKM z-scores” were selected. After entering PBLD, the “Co-expression” module was accessed to download PBLD co-expressed top 100 genes for subsequent GO/KEGG enrichment analysis.

For GO enrichment analysis, the list of differentially expressed genes (DEGs) was uploaded to the “[GOKEGG] Analysis” module under the “Functional Clustering” category on the Xiantao Academic platform. In the enrichment parameters, the “Category” was set to “All GO”, with all other parameters kept at their default settings. After confirmation, the resulting table was compiled and the top five significantly enriched GO terms for each gene function category were extracted. The clusterProfiler package was employed for enrichment analysis, and org.Hs.eg.db for gene ID conversion.

Visualization was performed using the “[GOKEGG] Bubble Chart” module, also located under “Functional Clustering”. The previously saved GO analysis result was uploaded, and the top five significantly enriched GO term IDs for each gene function category were entered. Other parameters were confirmed at their default settings. The ggplot2 package was used for data visualization.

KEGG enrichment analysis was conducted in an analogous manner, with the exception that the “KEGG” category was selected in the GOKEGG enrichment analysis, and the top five significantly enriched KEGG terms were subsequently extracted.

Gene Set Enrichment Analysis (GSEA) of PBLD-related DEGs

The “[Cloud] Single Gene-Differential Analysis” module within the “Expression Difference” category on the Xiantao Academic platform was employed to download PBLD-related DEGs from TCGA-COAD for subsequent GSEA. These DEGs were then analyzed using the “[GSEA] Enrichment Analysis” module under the “Functional Clustering” category. The analysis was executed by uploading the prepared data and applying the default parameters. For visualization, the “[GSEA] Classic Visualization” module under the “Functional Clustering” category was accessed. The saved GSEA result was selected, the IDs of significantly enriched pathways were entered, and default parameters were set before confirmation to obtain the final GSEA plots.

The ggplot2 package was employed for data visualization, and enrichment analyses were performed using the clusterProfiler package, with reference gene sets obtained from the MSigDBCollections database (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp) via the msigdbr package.

Statistical analysis

All statistical analyses were performed using established methodologies. Survival analyses (OS/DFS/PFS/DSS) utilized log-rank tests with Kaplan-Meier curves generated at median expression thresholds. Pathological stage associations used one-way ANOVA. Genetic alteration prognosis log-rank testing via cBioPortal. TMB/MSI and drug sensitivity correlations employed Pearson’s correlation, with FDR adjustment where specified. Immune cell and CAF infiltration analyses applied Spearman’s rank correlation. To investigate the potential mechanisms of PBLD in pan-cancer, enrichment analyses (GO/KEGG/GSEA) were performed utilizing the hypergeometric test with FDR correction. Differential gene expression between tumor-normal pairs employed the Welch’s t-test, while paired comparisons (e.g., COAD specimens) used paired t-tests. To investigate the potential role of PBLD in COAD, enrichment analyses (GO/KEGG) were performed using the hypergeometric test with Benjamini-Hochberg correction. Statistical significance was defined as P<0.05 (or FDR ≤0.05 for pharmacogenomic analyses).

Results

PBLD downregulation predicts adverse prognosis in different cancer types

Based on datasets of the HPA, GTEx, we found that PBLD exhibited high expression in normal tissues (intestine, liver, and kidney). Moreover, based on single-cell RNA-seq, PBLD was highly expressed in enterocytes, proximal tubular cells, and hepatocytes (Figure 1A,1B). Across TCGA’s pan-cancer cohort comprising over 11,000 tumor samples spanning 33 malignancies, PBLD was significantly dysregulated in 14 cancer types, with downregulation notably observed in bladder urothelial carcinoma (BLCA), breast invasive carcinoma (BRCA), cholangiocarcinoma (CHOL), COAD, esophageal carcinoma (ESCA), kidney chromophobe (KICH), kidney renal clear cell carcinoma (KIRC), kidney renal papillary cell carcinoma (KIRP), liver hepatocellular carcinoma (LIHC), lung squamous cell carcinoma (LUSC), and rectum adenocarcinoma (READ), whereas lung adenocarcinoma (LUAD), head and neck squamous cell carcinoma (HNSC), and prostate adenocarcinoma (PRAD) displayed upregulation in tumors versus matched normal tissues (Figure 1C). Using TCGA data, we assessed the prognostic significance of PBLD expression across different tumors via GEPIA2 and GSCA. We found that low PBLD expression predicted reduced survival across multiple endpoints: OS in KIRC, brain lower grade glioma (LGG), LIHC (Figure 1D), DFS in LGG, BRCA (Figure 1E), PFS in KIRC and LGG, DFI in BRCA, and DSS in KIRC, LGG, and LIHC (Figure 1F). To the contrary, high PBLD expression predicted reduced survival across these endpoints: DFS in uveal melanoma (UVM) and adrenocortical carcinoma (ACC) (Figure 1E), PFS in UVM, and DFI in HNSC (Figure 1F), suggesting context-dependent roles.

Figure 1 The expression profiles and prognostic value of PBLD across multiple cancer types in TCGA pan cancer. (A) Consensus PBLD tissue expression based on datasets of HPA and GTEx. (B) Expression level of PBLD in different cell subpopulations through single-cell analysis. (C) The expression of PBLD in different cancers. (D) Association of PBLD expression with OS in patients with KIRC, LGG, and LIHC. Kaplan-Meier curve for OS. The red dotted line represents the 95% CI for OS in the PBLD high expression group, and the blue dotted line represents the 95% CI for OS in the PBLD low expression group. (E) Kaplan-Meier curve for DFS in patients with ACC, UVM, BRCA, and LGG. The red dotted line represents the 95% CI for DFS in the PBLD high expression group, and the blue dotted line represents the 95% CI for DFS in the PBLD low expression group. (F) Association of PBLD expression level with PFS in KIRC, LGG, and UVM; with DFI in BRCA and HNSC; with DSS in KIRC, LGG, LIHC were analyzed by GSCA. (G) Genetic alterations of PBLD in pan-cancer analysis based on various tumor types of TCGA. (H) Correlation between PBLD mutation status and OS, PFS and DSS in CESC and COAD, as analyzed by cBioPortal. (I) Correlation between PBLD expression and pathological stages of ESCA and KIRC from TCGA datasets. *, P<0.05; **, P<0.01; ***, P<0.001. CI, confidence interval; CNA, copy-number alteration; DFI, disease-free interval; DFS, disease-free survival; DSS, disease-specific survival; ESCA, esophageal carcinoma; GSCA, Gene Set Cancer Analysis; GTEx, Genotype-Tissue Expression; HPA, Human Protein Atlas; HR, hazard ratio; OS, overall survival; PFS, progression-free survival; TCGA, The Cancer Genome Atlas; TPM, transcripts per million.

Genetic alterations of PBLD in various TCGA tumor types were analyzed using cBioPortal (Figure 1G). PBLD mutations independently predicted poor OS/PFS/DSS in cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC) and COAD (Figure 1H). We found that low PBLD expression was also correlated with advanced pathological stages in ESCA and KIRC through GEPIA2 (Figure 1I).

Association of PBLD with tumor immunogenicity and immune microenvironment

Significant associations were identified between PBLD expression and TMB, demonstrating a positive correlation in CHOL, KIRP, skin cutaneous melanoma (SKCM), but negative correlations in COAD, LUAD, thyroid carcinoma (THCA), and HNSC (Figure 2A). Furthermore, PBLD expression showed a positive correlation with MSI in testicular germ cell tumors (TGCT), while a negative correlation was observed in six cancer types: diffuse large B-cell lymphoma (DLBC), COAD, pancreatic adenocarcinoma (PAAD), PRAD, HNSC and BRCA (Figure 2B). Based on a broad analysis of immunocyte infiltration patterns across cancers in TCGA data, T helper cells showed a positive correlation with PBLD expression in 20 cancer types and a negative correlation in ACC. Central memory T cells (Tcm) were positively correlated with PBLD expression in 24 cancer types and negatively correlated in ESCA. Dendritic cells (DCs) exhibited a negative correlation with PBLD in 12 cancer types. Regulatory T cells (Tregs) were positively correlated with PBLD expression in LGG, pheochromocytoma and paraganglioma (PCPG) and negatively correlated in 12 cancer types (Figure 2C). CAFs are heterogeneous stromal cells within the TME and play a critical role in tumorigenesis, progression, metastasis, and therapeutic resistance (26). Furthermore, CAFs modulate immune responses in the TME. Therefore, we used TIMER 2.0 to explore the association between CAF infiltration and PBLD expression in TCGA and found that PBLD was associated with the abundance of CAF infiltration by employing EPIC, MCP-counter, xCell and TIDE algorithms (Figure 2D). Therein, TIDE is a functional prediction tool for evaluating immunotherapy response and immune escape mechanisms, which helps reveal the immunosuppressive consequences of CAF enrichment (27). Using this algorithm, we observed positive associations between CAF infiltration and PBLD expression in CESC, LUSC, HNSC-human papillomavirus-negative (HNSC-HPV⁻), and TGCT, but negative correlations in KIRC and PAAD (Figure 2E). Therefore, PBLD may be associated with tumor immunogenicity and the immune microenvironment.

Figure 2 Association of PBLD with tumor immunogenicity and immune microenvironment in TCGA pan-cancer. (A) Radar diagram of the correlation between PBLD expression and TMB in TCGA pan-cancer. (B) Radar diagram of the correlation between PBLD expression and MSI in TCGA pan-cancer. (C) Correlation between PBLD expression and 12 types of immune cell infiltration in TCGA pan-cancer. (D) EPIC, MCP-counter, xCell and TIDE algorithms were employed to calculate the correlation between PBLD expression and CAF infiltration in TCGA pan-cancer. (E) Correlation between PBLD expression and CAF infiltration in CESC, LUSC, HNSC-HPV⁻, TGCT, KIRC, and PAAD based on TIDE algorithm. MSI, microsatellite instability; TCGA, The Cancer Genome Atlas; TGCT, testicular germ cell tumor; TIDE, Tumor Immune Dysfunction and Exclusion; TMB, tumor mutational burden; TPM, transcripts per million.

PBLD expression is a putative determinant of drug response

Comprehensive pharmacogenomic analysis through GSCA revealed robust drug sensitivity associations for PBLD expression at an FDR threshold of ≤0.05. Within the GDSC database spanning 860 cancer cell lines, elevated PBLD expression showed significant positive correlations with the IC₅₀ for 37 therapeutic agents, exemplified by SN-38, QL-VIII-58, Camptothecin, and (5Z)-7-oxozeaenol. This analysis simultaneously demonstrated negative correlations between PBLD expression and IC₅₀ for 11 compounds, such as the AKT inhibitor VIII and SGC0946 (Table S1). The top 30 ranked drugs are presented in Figure 3A. Analysis of the larger CTRP dataset that encompassed 1,001 cell lines further corroborated the predictive capacity of PBLD, identifying positive correlations between PBLD expression and IC₅₀ for 110 pharmacological agents, including alvocidib, dinaciclib, and leptomycin B. In contrast, only 4 compounds exhibited negative correlations between PBLD expression and IC₅₀, including RITA, cediranib, saracatinib, and vandetanib (Table S2). The top 30 ranked drugs are presented in Figure 3B. Collectively, these data nominate PBLD as a trans-dataset molecular predictor of responsiveness to both conventional chemotherapeutics and molecularly targeted agents.

Figure 3 Correlation between PBLD expression and drug sensitivity in pan-cancer. (A) Correlation between PBLD expression and GDSC drug sensitivity (top 30) in pan-cancer through GSCA. (B) Correlation between PBLD expression and CTRP drug sensitivity (top 30) in pan-cancer through GSCA. CTRP, Cancer Therapeutics Response Portal; FDR, false discovery rate; GDSC, Genomics of Drug Sensitivity in Cancer; GSCA, Gene Set Cancer Analysis.

Metabolic regulation underpins the function of PBLD

To further investigate the mechanism of PBLD in the BP of cancer cells, GEPIA2 was used to extract the top 100 genes with expression patterns similar to that of PBLD across all tumor types from TCGA. We found that their core co-expressed networks demonstrated significant functional convergence in metabolic regulation. Subsequent GO-BP enrichment analysis confirmed these genes were enriched in organic acid, oxoacid, alpha-amino acid and lipid metabolic process (Figure 4A). KEGG enrichment analysis indicated that PBLD-related genes were primarily associated with pathways for butyrate metabolism; valine, leucine, and isoleucine degradation; glycine, serine, and threonine metabolism; fatty acid degradation and metabolism (Figure 4B). A co-expression network for PBLD, constructed using the STRING database, identified hub nodes including ABAT, EHHADH, HGD, NR1H4, ALDH8A1, ACMSD, IYD, ADH6, A1CF, SLC22A18, HAAO, SLC2A2, RSPH3, UEVLD, MTRES1, PNPO, ENTPD5, and RBKS (Figure 4C). These 18 genes demonstrated relatively close interactions. GO and KEGG enrichment analyses of these co-expressed genes also revealed significant enrichment in metabolic processes for amino acids, organic acids, and fatty acids (Figure 4D,4E). Pearson correlation analysis revealed that PBLD expression was strongly correlated with the expression levels of ACMSD, EHHADH and NR1H4 through GEPIA2 (Figure 4F). Based on the above results, we speculate that PBLD may play an important role in tumor development, metastasis, and invasion by participating in anabolic and catabolic processes involving substances such as amino acids and fatty acids in tumor cells.

Figure 4 Functional enrichment and interaction network analyses of PBLD-related genes. (A) GO-BP enrichment analysis of the top 100 genes co-expressed with PBLD obtained by the GEPIA2; (B) KEGG enrichment analysis of the top 100 genes co-expressed with PBLD obtained by the GEPIA2; (C) Co-expression network of 18 genes co-expressed with PBLD obtained by the STRING; (D) GO-BP enrichment analysis of 18 genes co-expressed with PBLD; (E) KEGG enrichment analysis of the 18 genes co-expressed with PBLD; (F) Correlation analysis between PBLD and ACMSD, EHHADH, and NR1H4 conducted by GEPIA2 in all tumor samples from TCGA. BP, biological process; GEPIA2, Gene Expression Profiling Interactive Analysis 2; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; STRING, Search Tool for the Retrieval of Interacting Genes/Proteins; TCGA, The Cancer Genome Atlas; TPM, transcripts per million.

Epigenetic, immune, and metabolic programs mediated by PBLD in the TCGA-COAD cohort

We downloaded COAD RNA-seq data from the TCGA database and analyzed the expression level of PBLD. The results showed that PBLD expression was downregulated in COAD tissues compared to normal tissues (Figure 5A). Paired differential expression analysis further confirmed decreased PBLD expression in COAD tissues compared to matched normal tissues (Figure 5B). From the TCGA database, we obtained 477 COAD samples with complete clinical information. Using the median as the cutoff, patients were divided into a PBLD high-expression group (n=238) and a PBLD low-expression group (n=239). We subsequently analyzed the association between PBLD expression and progression-free interval (PFI) in patients with COAD. The results revealed that low PBLD expression was correlated with a shorter PFI in patients with COAD (Figure 5C).

Figure 5 Clinical significance and immune-metabolic landscape mediated by PBLD in the TCGA-COAD cohort. (A,B) PBLD expression in unpaired or paired samples from the TCGA-COAD cohort. (C) PFI curve of PBLD from TCGA database. (D) Correlations between PBLD expression and the infiltration levels of 24 immune cell types from the TCGA-COAD cohort. (E) Association of PBLD expression with NK CD56^bright and Th1 cell infiltration. (F) GO enrichment analysis of the top 100 genes positively correlated with PBLD from the cBioPortal database. (G) KEGG pathway enrichment analysis of the top 100 genes positively correlated with PBLD from the cBioPortal database. (H) GSEA of PBLD-related DEGs in the TCGA-COAD cohort identified the top five positively and negatively enriched pathways. ns, P≥0.05; *, P<0.05; **, P<0.01; ***, P<0.001. BP, biological process; CC, cellular component; COAD, colon adenocarcinoma; DEGs, differentially expressed genes; GO, Gene Ontology; GSEA, Gene Set Enrichment Analysis; HR, hazard ratio; KEGG, Kyoto Encyclopedia of Genes and Genomes; NK, natural killer; PFI, progression-free interval; TCGA, The Cancer Genome Atlas; TPM, transcripts per million.

We calculated the infiltration abundance of 24 immune cell types in COAD tissues from the TCGA-COAD cohort using the ssGSEA algorithm. Spearman correlation analysis was then performed to assess the relationship between the infiltration abundance of these immune cell types and PBLD expression (Figure 5D). In TCGA-COAD cohort, tumors with low PBLD expression exhibited enriched infiltration of NK CD56^bright cells and Th1 cells (Figure 5E).

The top 100 genes positively correlated with PBLD in the TCGA-COAD cohort (Spearman correlation, P<0.05) were obtained from the cBioPortal database for subsequent GO and KEGG enrichment analyses. Using thresholds of an adjusted P<0.05 and Q<0.25, the analysis identified seven significantly enriched BP terms, 10 CC terms, zero MF terms, and seven KEGG pathways. The top five terms with the smallest P value were selected for visualization. GO enrichment analysis revealed significant enrichment in pathways including organic acid transport, purine nucleoside bisphosphate metabolic process, ribonucleoside bisphosphate metabolic process, peroxisome, and peroxisomal matrix (Figure 5F). KEGG enrichment analysis identified pathways such as bile secretion, drug metabolism−cytochrome P450, PPAR signaling pathway, and retinol metabolism (Figure 5G). GSEA of PBLD-related DEGs revealed significant upregulation of pathways including CD22-mediated BCR regulation, FcεRI-mediated MAPK activation, FcεRI-mediated Ca²⁺ mobilization, and antigen activates B cell receptor. In contrast, pathways such as graft-versus-host disease and systemic lupus erythematosus were significantly downregulated (Figure 5H). Collectively, these interactions implied functional cross-communication spanning metabolic immunological and epigenetic dimensions.

Discussion

This study systematically delineates the pivotal role and clinical significance of PBLD in multiple malignancies through pan-cancer bioinformatics analysis, with investigation focused on COAD. Our findings establish PBLD as a putative tumor suppressor gene and prognostic biomarker. The analyses demonstrate significantly reduced PBLD expression across 14 cancer types encompassing BLCA, BRCA, CHOL, COAD, KICH, KIRC, KIRP, LIHC, LUSC, ESCA, and READ. Crucially, this downregulation independently predicts adverse prognosis in these cancers, manifesting as significantly shortened survival across multiple endpoints. These observations corroborate prior reports of PBLD as a tumor suppressor in single cancer types [gastric cancer, hepatocellular carcinoma (HCC), colorectal cancer (8,9,12,13)] while substantially extending its relevance to pan-cancer contexts. Simultaneously, elevated PBLD expression in LUAD, PRAD, and HNSC suggests cancer type-specific functionality, warranting mechanistic exploration. PBLD mutations also exhibit prognostic potential in CESC and COAD. The correlation between low expression of PBLD and progression to advanced-stages in ESCA and KIRC reinforces its tumor-suppressive role during oncogenesis. Collectively, PBLD demonstrates potential as a cross-cancer diagnostic and prognostic biomarker.

PBLD was additionally associated with tumor immunogenicity and therapy response. TMB and MSI constitute critical indicators of immunogenicity and predictors of immune checkpoint inhibitor (ICI) efficacy (28,29). Our investigation revealed that PBLD expression significantly influenced TMB/MSI status in multiple cancers. For example, in COAD, BRCA and THCA, decreased PBLD expression was associated with elevated TMB/MSI, indicating potential benefit from ICIs in these subsets of patients. Conversely, in CHOL and KIRP, low PBLD expression correlated with reduced TMB/MSI, potentially reducing immunotherapeutic responsiveness.

In accordance with these findings, PBLD expression not only predicted sensitivity to small-molecule therapeutics, but also showed associations with the TME. Pan-cancer analyses indicated significant correlations between PBLD expression and infiltration abundance of diverse immune cells, including T helper cells, Tregs, DCs, and Tcm, but showed context-specific positive/negative associations with CAF infiltration estimated by multiple algorithms. More critically, this study unveils the role of PBLD in sculpting the COAD immune infiltrate landscape. PBLD downregulation was associated with increased infiltration of NK CD56^bright cells and Th1 cells. This sophisticated regulatory network positions PBLD as a potential coordinator of immunosuppressive versus immunostimulatory states in COAD, directly influencing tumor response to immunotherapy (26,30-32). These insights provide novel perspectives for understanding and predicting TME dynamics in colorectal cancer.

Metabolic reprogramming and molecular interactions constitute fundamental mechanistic bases for the function of PBLD. To further elucidate its tumor-suppressive mechanisms in COAD, we conducted enrichment analyses of PBLD-co-expressed genes and DEGs. Pan-cancer analyses of PBLD-related and co-expressed genes (including highly correlated partners ACMSD, NR1H4, and EHHADH) alongside GO/KEGG and GSEA results in COAD collectively converged on the key involvement of PBLD in amino acid metabolism, lipid metabolism, and nucleotide metabolism—pivotal biosynthetic and catabolic processes. Metabolic reprogramming constitutes a hallmark adaptation enabling tumor cells to meet biosynthetic and bioenergetic demands for rapid proliferation (33). PBLD likely impacts cancer cell proliferation, metastasis, and therapy resistance by rewiring these metabolic pathways. Importantly, in terms of immunomodulation, PBLD-co-expressed genes and enriched pathways in COAD strongly implicate their profound engagement in immune regulation.

PBLD is implicated not only in COAD but also in stomach adenocarcinoma (STAD) and HCC, highlighting its potential pan-cancer relevance. To translate these associations into clinical insight, cross-cancer validation is the critical next step to ascertain whether its clinical value is cancer-type specific.

Conclusions

In summary, PBLD represents a highly promising novel oncological biomarker and therapeutic target. Future research should prioritize in-depth mechanistic dissection of metabolic pathways and immunomodulatory axes alongside validation of their diagnostic utility and therapeutic targeting potential in both in vitro and in vivo models.

Acknowledgments

We gratefully acknowledge Xiantao Academic platform and OmicShare platform for their support in conducting the bioinformatic analyses for this study.

Footnote

Reporting Checklist: The authors have completed the REMARK reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1-2921/rc

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1-2921/prf

Funding: This study was supported by the National Natural Science Foundation of China (No. 82470575), and Natural Science Foundation of Guangdong Province (Nos. 2024A1515012415 and 2025A1515012595).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1-2921/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Bray F, Laversanne M, Sung H, et al. Global cancer statistics 2022: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2024;74:229-63. [Crossref] [PubMed]
2. Cho WC. Contribution of oncoproteomics to cancer biomarker discovery. Mol Cancer 2007;6:25. [Crossref] [PubMed]
3. Han B, Zheng R, Zeng H, et al. Cancer incidence and mortality in China, 2022. J Natl Cancer Cent 2024;4:47-53. [Crossref] [PubMed]
4. Siegel RL, Giaquinto AN, Jemal A. Cancer statistics, 2024. CA Cancer J Clin 2024;74:12-49. [Crossref] [PubMed]
5. Zhao X, Kang B, Lu C, et al. Evaluation of p38 MAPK pathway as a molecular signature in ulcerative colitis. J Proteome Res 2011;10:2216-25. [Crossref] [PubMed]
6. Zhao X, Fan J, Zhi F, et al. Mobilization of epithelial mesenchymal transition genes distinguishes active from inactive lesional tissue in patients with ulcerative colitis. Hum Mol Genet 2015;24:4615-24. [Crossref] [PubMed]
7. Chen S, Liu H, Li Z, et al. Epithelial PBLD attenuates intestinal inflammatory response and improves intestinal barrier function by inhibiting NF-κB signaling. Cell Death Dis 2021;12:563. [Crossref] [PubMed]
8. Chen Y, Liang J, Chen S, et al. Discovery of vitexin as a novel VDR agonist that mitigates the transition from chronic intestinal inflammation to colorectal cancer. Mol Cancer 2024;23:196. [Crossref] [PubMed]
9. Li DM, Zhang J, Li WM, et al. MAWBP and MAWD inhibit proliferation and invasion in gastric cancer. World J Gastroenterol 2013;19:2781-92. [Crossref] [PubMed]
10. Chu F, Hou P, Zhu H, et al. PBLD enhances antiviral innate immunity by promoting the p53-USP4-MAVS signaling axis. Proc Natl Acad Sci U S A 2024;121:e2401174121. [Crossref] [PubMed]
11. Schvartz D, Couté Y, Brunner Y, et al. Modulation of neuronal pentraxin 1 expression in rat pancreatic β-cells submitted to chronic glucotoxic stress. Mol Cell Proteomics 2012;11:244-54. [Crossref] [PubMed]
12. Liang Y, Song X, Li Y, et al. circKDM4C suppresses tumor progression and attenuates doxorubicin resistance by regulating miR-548p/PBLD axis in breast cancer. Oncogene 2019;38:6850-66. [Crossref] [PubMed]
13. Long J, Lang ZW, Wang HG, et al. Glutamine synthetase as an early marker for hepatocellular carcinoma based on proteomic analysis of resected small hepatocellular carcinomas. Hepatobiliary Pancreat Dis Int 2010;9:296-305. [PubMed]
14. Han L, Lin X, Yan Q, et al. PBLD inhibits angiogenesis via impeding VEGF/VEGFR2-mediated microenvironmental cross-talk between HCC cells and endothelial cells. Oncogene 2022;41:1851-65. [Crossref] [PubMed]
15. Liu J, Lichtenberg T, Hoadley KA, et al. An Integrated TCGA Pan-Cancer Clinical Data Resource to Drive High-Quality Survival Outcome Analytics. Cell 2018;173:400-416.e11. [Crossref] [PubMed]
16. Pontén F, Schwenk JM, Asplund A, et al. The Human Protein Atlas as a proteomic resource for biomarker discovery. J Intern Med 2011;270:428-46. [Crossref] [PubMed]
17. Li T, Fu J, Zeng Z, et al. TIMER2.0 for analysis of tumor-infiltrating immune cells. Nucleic Acids Res 2020;48:W509-14. [Crossref] [PubMed]
18. Tang Z, Kang B, Li C, et al. GEPIA2: an enhanced web server for large-scale expression profiling and interactive analysis. Nucleic Acids Res 2019;47:W556-60. [Crossref] [PubMed]
19. Gao J, Aksoy BA, Dogrusoz U, et al. Integrative analysis of complex cancer genomics and clinical profiles using the cBioPortal. Sci Signal 2013;6:pl1. [Crossref] [PubMed]
20. Liao C, Wang X. TCGAplot: an R package for integrative pan-cancer analysis and visualization of TCGA multi-omics data. BMC Bioinformatics 2023;24:483. [Crossref] [PubMed]
21. Yang W, Soares J, Greninger P, et al. Genomics of Drug Sensitivity in Cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells. Nucleic Acids Res 2013;41:D955-61. [Crossref] [PubMed]
22. Basu A, Bodycombe NE, Cheah JH, et al. An interactive resource to identify cancer genetic and lineage dependencies targeted by small molecules. Cell 2013;154:1151-61. [Crossref] [PubMed]
23. Szklarczyk D, Kirsch R, Koutrouli M, et al. The STRING database in 2023: protein-protein association networks and functional enrichment analyses for any sequenced genome of interest. Nucleic Acids Res 2023;51:D638-46. [Crossref] [PubMed]
24. Mu H, Chen J, Huang W, et al. OmicShare tools: A zero-code interactive online platform for biological data analysis and visualization. Imeta 2024;3:e228. [Crossref] [PubMed]
25. Bindea G, Mlecnik B, Tosolini M, et al. Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity 2013;39:782-95. [Crossref] [PubMed]
26. Cords L, de Souza N, Bodenmiller B. Classifying cancer-associated fibroblasts-The good, the bad, and the target. Cancer Cell 2024;42:1480-5. [Crossref] [PubMed]
27. Jiang P, Gu S, Pan D, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med 2018;24:1550-8. [Crossref] [PubMed]
28. Holder AM, Dedeilia A, Sierra-Davidson K, et al. Defining clinically useful biomarkers of immune checkpoint inhibitors in solid tumours. Nat Rev Cancer 2024;24:498-512. [Crossref] [PubMed]
29. Luchini C, Bibeau F, Ligtenberg MJL, et al. ESMO recommendations on microsatellite instability testing for immunotherapy in cancer, and its relationship with PD-1/PD-L1 expression and tumour mutational burden: a systematic review-based approach. Ann Oncol 2019;30:1232-43. [Crossref] [PubMed]
30. de Visser KE, Joyce JA. The evolving tumor microenvironment: From cancer initiation to metastatic outgrowth. Cancer Cell 2023;41:374-403. [Crossref] [PubMed]
31. Dou A, Fang J. Heterogeneous Myeloid Cells in Tumors. Cancers (Basel) 2021;13:3772. [Crossref] [PubMed]
32. Peña-Romero AC, Orenes-Piñero E. Dual Effect of Immune Cells within Tumour Microenvironment: Pro- and Anti-Tumour Effects and Their Triggers. Cancers (Basel) 2022;14:1681. [Crossref] [PubMed]
33. Lyu S, Gildor N, Zhang Q, et al. Rewiring cancer metabolism: oncogenic signaling pathways and targeted therapeutics. Sci China Life Sci 2025;68:3281-302. [Crossref] [PubMed]