Integrating transcriptomics, single-cell omics, and deep learning-based histopathological features to identify OLFML3 in bladder cancer

Introduction

Bladder cancer (BCa) is the most common malignancy of the urinary system and the seventh most prevalent cancer globally, with over 549,000 new cases and more than 199,000 deaths annually (1). The primary challenge in the clinical management of BCa lies in its high recurrence rate, with approximately 61% of patients experiencing recurrence within 1 year post-surgery, the highest rate among all solid tumors (2). While early-stage BCa can be treated surgically, its recurrence rate remains notably high, often reaching 60–80% (3). The 5-year survival rates for patients with local recurrence and distant metastasis are only 45% and 6%, respectively (4). Therefore, elucidating the mechanisms driving BCa recurrence and developing precise predictive tools are critical for improving patient prognosis.

Currently, international urology guidelines recommend monitoring for BCa recurrence primarily through regular upper urinary tract radiological assessments, cystoscopy, and urine cytology analysis. Although cystoscopy allows direct tumor visualization, its invasiveness may lead to urinary tract infections, causing discomfort and procedural complications for patients. Additionally, its low sensitivity potentially results in misdiagnoses of bladder malignancies in 10–40% of cases (3). While urine cytology presents a safe and cost-effective diagnostic tool, its sensitivity is suboptimal. However, it exhibits high specificity for diagnosing urothelial carcinoma (UC), of which BCa is the most prevalent malignant tumor, while upper urinary tract cancers are relatively rare (5). A significant body of research has focused on DNA methylation abnormalities and their crucial role in BCa pathogenesis and progression. Urine cytology can effectively assess patient survival by detecting the methylation status of these genes (6). Fiorentino et al. emphasized the importance of methylation analysis and the Bladder EpiCheck™ test as urine biomarkers for diagnosing non-muscle invasive bladder cancer (NMIBC) during follow-up, particularly in augmenting the efficacy of cytological and cystoscopic examinations for suspicious cases (7). Pepe et al. demonstrated that combining Bladder EpiCheck with cytology may represent the most effective method for the follow-up diagnosis of high-grade NMIBC, thereby minimizing the need for unnecessary invasive procedures (8).

However, UC exhibits significant molecular and histological heterogeneity (9), which contributes to the low sensitivity of urine cytology as a prognostic follow-up diagnostic tool for BCa recurrence. Despite the implementation of the Paris System for Reporting Urinary Cytology (TPS), a standardized evidence-based system designed to standardize reporting and enhance communication between pathologists and clinicians (5), methylation studies have demonstrated that different TPS cytological categories exhibit distinct molecular features. Further studies have confirmed that the cytological categories “suspicious for high-grade UC” (SHGUC) and “high-grade UC” (HGUC) represent distinct entities at the molecular level, warranting their continued classification as separate groups within the TPS system (10).

Additionally, the sensitivity of cytology remains poor, particularly for detecting low-grade tumors (sensitivity <50%), significantly limiting its clinical applicability (11). Despite existing clinical risk stratification tools, such as the European Organization for Research and Treatment of Cancer scoring system—which integrates various parameters (e.g., tumor stage, grade, and multifocality) to provide important prognostic references for BCa (12)—the predictive efficacy of these tools is still influenced by tumor heterogeneity, complicating the implementation of personalized precision medicine. Consequently, the development of predictive tools with higher sensitivity and accuracy is essential for implementing individualized treatment strategies and improving prognosis management.

In recent years, the rapid advancement of multi-omics technologies has provided unprecedented insights into the mechanisms underlying BCa recurrence. This study innovatively integrated multi-omics data analysis methods and, for the first time, revealed the key role of OLFML3 in BCa recurrence within 1 year post-surgery. Based on this finding, we developed a ResNet50-based model to accurately locate tumor areas in whole-slide images (WSIs). By extracting deep learning features that more accurately reflect the biological characteristics of the tumor, we constructed a predictive model that integrates pathological information with deep learning features. This model can accurately assess the pathological features of patients with BCa and perform static risk stratification for high‑risk individuals, providing a basis for personalized interventions. We present this article in accordance with the TRIPOD reporting checklist (available at https://tau.amegroups.com/article/view/10.21037/tau-2025-365/rc).

Methods

Data acquisition and preprocessing

First, RNA-seq data and corresponding clinical information for patients with BCa were obtained from The Cancer Genome Atlas (TCGA; https://www.cancer.gov/tcga) bladder urothelial carcinoma (BLCA) project. Based on clinical follow-up data, patients were stratified into two groups: those with recurrence within 1 year and those without recurrence. The raw RNA-seq data underwent background correction, normalization, and quality control (QC), excluding genes with low expression or high variability to ensure data quality and reliability for subsequent analyses.

Weighted gene co-expression network analysis (WGCNA) based on TCGA data: screening gene modules positively correlated with 1-year BCa recurrence

We initially transposed the expression matrix of the TCGA dataset, treating gene expression values as rows and samples as columns for subsequent analysis. To reduce computational load and eliminate low-variance genes, we selected the top 25% of genes based on variance. Specifically, we calculated the variance of each gene across all samples and retained those with variance exceeding the upper quartile. This approach ensured that genes with higher expression variation were preserved, thereby minimizing the influence of low-variance genes on the results.

Thereafter, we utilized the goodSamplesGenes function from the WGCNA package for gene and sample QC. This function checks for missing values and identifies abnormal genes, removing any that do not satisfy established quality criteria. Genes and samples with missing or abnormal values were automatically excluded, resulting in a dataset that adhered to the quality standards. To identify potential outliers, we performed sample clustering. A sample dendrogram was constructed using the hclust function, and an appropriate cutoff value of 95 was selected based on clustering results. Outlier samples identified in the dendrogram were excluded from further analysis. After data preprocessing, we incorporated the phenotype data, ensuring that they matched the gene expression data by removing any samples excluded in the previous filtering steps. These phenotype data provided the necessary context for subsequent module-phenotype association analysis.

Subsequently, we selected an appropriate soft threshold using the pickSoftThreshold function, which evaluates the impact of various soft thresholds (ranging from 1 to 30) on scale-free network topology fit and determines the optimal value. In this study, we selected a soft threshold of 4 to ensure that the constructed gene co-expression network maintained a scale-free topology.

Thereafter, we constructed the gene co-expression network using the blockwiseModules function, which employs dynamic tree-cutting algorithms to identify multiple gene modules. To ensure stability, we set a minimum module size of 50 and a merge cut height of 0.6 to combine similar modules. The deepSplit parameter was adjusted to enhance sensitivity in module division, facilitating the identification of numerous modules. After constructing the modules, we calculated the module eigengenes for each module and performed Pearson correlation analysis with the phenotype data. Modules that significantly positively correlated with 1-year BCa recurrence (P<0.05) were selected for further investigation regarding their association with clinical recurrence.

Finally, a heatmap was used to visualize module-phenotype correlation.

Univariate Cox regression analysis of red module (MEred) genes

To explore the relationship between MEred genes and patient survival, we initially extracted pre-processed clinical and expression data, which included survival time, event status, and gene expression values. Using the survival, survminer (version 0.5.0), and forestplot (version 3.1.6) packages in R, we analyzed each candidate gene within MEred using univariate Cox regression models (coxph function). For each gene variable (excluding “time” and “event” in the dataframe), we calculated regression coefficients, hazard ratios with 95% confidence intervals, and extracted corresponding P values. Genes fulfilling the predefined significance threshold (P<0.001) were selected as survival-associated candidates.

Prognostic gene screening using least absolute shrinkage and selection operator (LASSO)-penalized Cox regression models

To further screen the prognosis-associated genes initially identified in the univariate analysis, we employed the LASSO-penalized Cox regression model. First, we converted the gene expression data into a numeric matrix and constructed survival objects based on overall survival (OS) time and follow-up status (censor). Subsequently, we used the cv.glmnet function from the glmnet package (version 4.1.8) to determine the optimal penalty parameter lambda through cross-validation, where lambda.min corresponds to the minimum mean squared error, while lambda.1se follows the one-standard-error rule. With a maximum of 1,000 iterations, we constructed the LASSO-Cox model and generated coefficient path plots to display variable coefficient variations under different lambda values. Finally, candidate genes with significant prognostic value were selected based on non-zero coefficients at the optimal lambda, serving as key targets for subsequent studies.

Construction of an eight-gene prognostic risk score model based on LASSO-Cox regression and validation across multiple Gene Expression Omnibus (GEO) datasets

The LASSO-penalized Cox regression model identified eight prognosis-associated genes with corresponding regression coefficients. We calculated each patient’s risk score by summing the products of these eight genes’ expression levels and their respective coefficients. Using survival time (OStime.y) and status (OSstatus.y), we determined the optimal risk score cutoff using the surv_cutpoint function from the survminer package, stratifying patients into high- and low-risk groups. Kaplan-Meier survival curves were generated, and between-group survival differences were assessed using the log-rank test. Survival curves with corresponding P values were visualized using the ggsurvplot function (version 0.5.0), and grouping results and figures were saved for subsequent analysis.

Following the establishment of the eight-gene prognostic model in the TCGA cohort, we further validated its predictive efficacy in two independent GEO (https://www.ncbi.nlm.nih.gov/geo) datasets (GSE13507 and GSE31684).

Independent survival analysis of eight genes and association analysis with clinical variables

After screening the eight prognosis-associated genes using the LASSO-penalized Cox regression model, we constructed individual survival risk plots for each gene to analyze their relationships with patient outcomes. To further investigate associations between these genes and clinical variables, we generated survival violin plots for the eight genes across clinical parameters.

Differential gene expression between high and low tumor stromal BCa subtypes

Single-cell RNA sequencing (scRNA-seq) and data processing

scRNA-seq was performed on BCa tissue samples obtained from six patients within the GD2H cohort. This cohort comprised three patients with low tumor stromal content (low-stromal BCa) and three with high tumor stromal content (high-stromal BCa). All raw and processed scRNA-seq data generated in this study have been deposited in the National Genomics Data Center, Beijing Institute of Genomics, Chinese Academy of Sciences, under accession number GSA-Human: HRA009938. Raw sequencing data underwent comprehensive preprocessing, with initial QC applied to filter out low-quality cells. Cells were retained based on thresholds for the number of detected genes, total RNA content (counts), and the percentage of mitochondrial gene expression. Subsequently, data normalization and batch effect correction were performed using the Seurat package (v5.1.0). Principal component analysis was conducted on highly variable genes, followed by graph-based clustering to partition cells into distinct subpopulations.

Cell type annotation

Cell clusters were annotated into major lineages and specific cell types based on canonical marker gene expression patterns. The primary marker genes used for annotation were as follows: immune lineage: B cells: VPREB3; plasma cells: CD79A; T cells: CD3D, PTPN22; natural killer cells: GNLY; myeloid lineage: dendritic cells: LAMP3; macrophages: APOC1; monocytes: FCN1; neutrophils: EREG; epithelial lineage: epithelial cells: CENPA, SOX4, ITGA2, EPCAM; fibroblast lineage: fibroblasts: COL1A1; inflammatory cancer-associated fibroblasts: PDGFRA; and myofibroblasts: RGS5.

Differential gene expression analysis

Following cell type annotation, we conducted differential gene expression analysis between the high- and low-stromal BCa patient groups. Differentially expressed genes (DEGs) were identified using standard comparative analysis methods within the Seurat framework. To visualize the expression patterns of key DEGs and OLFML3, we generated bar plots depicting its average expression per cell type/group, along with uniform manifold approximation and projection (UMAP) plots highlighting its expression across cellular subpopulations.

Deep learning feature-based prediction of OLFML3 expression in BCa using whole-slide hematoxylin and eosin (H&E) images

WSIs of patients with BCa were selected from the TCGA database based on the following inclusion criteria: a pathological diagnosis of bladder UC, available RNA-seq data, complete clinical-pathological records, and high-quality H&E-stained slides. The exclusion criteria entailed the absence of identifiable lesions, poor image quality, missing feature data, and a history of preoperative treatments. For external validation, H&E-stained slides and matched RNA-seq data from patients with BCa at two centers—Shanghai Tenth People’s Hospital (STPH) and Guangdong Second Provincial People’s Hospital (GD2H)—were also collected. These patients met the same inclusion and exclusion criteria, and data collection occurred between September 2017 and May 2024. RNA-seq data from TCGA, STPH, and GD2H were merged, and batch correction was applied to the “fragments per kilobase of exon per million mapped reads” expression values of OLFML3 using ComBat. Based on the corrected data, patients were categorized into low and high OLFML3 expression groups using the median expression value as the cutoff. Considering that TCGA provided the largest dataset, its WSIs were randomly divided into training and internal validation sets at a 7:3 ratio, while those from STPH and GD2H served as the external validation set. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The study was approved by the respective committees at STPH (Approval No. 24KT68) and GD2H (Approval No. 2024-KY-KZ-128-01). Informed consent was collected from all participants. A summary of the clinical features of the patients is provided in Table S1.

To develop a deep learning model for distinguishing tumor and non-tumor regions in BCa H&E WSIs, 70 WSIs from the TCGA dataset were randomly selected for model training, with 30 WSIs designated for internal validation. Additionally, 10 WSIs from each of the STPH and GD2H datasets were utilized for external validation. All selected WSIs were manually annotated for tumor regions by experienced pathologists using QuPath software (v0.3.2). OpenSlide software (v4.0.0) was employed to segment the annotated WSIs into 224×224-pixel image patches at 20× magnification. Edge detection algorithms from OpenCV (v4.10.0) with a threshold of 0.02 were applied to automatically identify and exclude image patches predominantly containing blank backgrounds. To ensure color consistency across WSIs from diverse sources, the Reinhard color normalization method was applied to the remaining tissue-containing patches.

Using the PyTorch 2.4.1 framework, a transfer learning model based on the ResNet50 architecture was constructed to differentiate between tumor and normal tissue patches. The model was trained over 50 epochs using the stochastic gradient descent optimizer with the following key parameters: an initial learning rate of 0.1, a momentum of 0.9, a weight decay of 0.001, and a cross-entropy loss function. After training on the TCGA training set, the model’s performance was evaluated on the TCGA internal validation set as well as the external validation sets from STPH and GD2H.

Additionally, a pre-trained RetCCL model was employed for deep feature extraction. RetCCL is a variant of ResNet50 optimized through contrastive learning on 22,000 WSIs, specifically designed to extract discriminative histopathological image features (13). To focus on tumor-specific features, tumor image patches were initially selected using the trained ResNet50 model. By concentrating on tumor-specific regions, the ResNet50 model mitigates the influence of normal tissue, a common source of noise in histopathological image analysis, thereby enhancing the overall feature extraction process for downstream machine learning models (14). The final deep learning features were subsequently extracted from the 2,048-dimension global average pooling layer of the RetCCL model.

For each patient (i.e., each WSI), the 2,048-dimension feature vector of all tumor image patches was summarized using seven statistical metrics: mean, median, standard deviation (std), range (max − min), first quartile (Q1), third quartile (Q3), and interquartile range (Q3 − Q1). This aggregation yielded a comprehensive feature vector comprising 14,336 features (2,048 dimensions × 7 statistics) for each WSI.

In the model development phase, all 14,336 WSI-level features underwent Z-score standardization. The mean and standard deviation necessary for standardization were retained and consistently applied to both the internal and external validation sets. Feature selection was performed using minimal redundancy maximum relevance analysis, which identified the top 20 most predictive and least redundant features.

Based on the selected features, a random forest (RF) model was trained and validated, using 70% of the TCGA patient data for training. During hyperparameter tuning, a grid search was conducted to optimize key hyperparameters. The parameter space was defined as follows:

param_grid = {
'n_estimators': range(10, 301), # Number of trees, from 10 to 300, step size 1
'max_depth': list(range(1, 11)), # Integers from 1 to 10
'min_samples_split': [2, 5, 10], # Minimum number of samples required to split an internal node
'min_samples_leaf': [1, 2, 4], # Minimum number of samples required to be at a leaf node
'bootstrap': [True, False] # Whether to use bootstrap sampling
}

The optimal hyperparameters were determined via grid search. Model performance was assessed using the internal validation set from TCGA (30% of TCGA patient data) as well as the external validation sets from STPH and GD2H, with performance evaluated across several metrics: area under the curve (AUC), accuracy, sensitivity, specificity, positive predictive value, and negative predictive value, to assess the model’s ability to detect OLFML3 expression.

Statistical analysis

Statistical analyses were performed using R software (version 4.4.2), with the ggplot2 package (version 3.5.1) employed to visualize expression differences. A P value of less than 0.05 was considered statistically significant.

Results

WGCNA analysis: association between gene modules and recurrence status

To investigate the relationship between gene expression modules and clinical characteristics (e.g., recurrence status) in patients with BCa, we employed WGCNA to identify recurrence-associated gene modules. First, the association between modules and recurrence was evaluated by ascertaining the correlation between module eigengenes and recurrence status. Figure 1A presents a heatmap illustrating sample phenotypes with BCa recurrence or non-recurrence within 1 year, where recurrence and non-recurrence samples are indicated by a colored bar.

Figure 1 Transcriptomics workflow diagram. (A) Sample dendrogram and trait heatmap: this panel displays hierarchical clustering of samples and visualizes recurrence status (recurrence vs. no recurrence) as clinical traits aligned below the dendrogram. (B) Soft threshold selection and gene connectivity analysis: this section displays scale independence and mean connectivity across soft thresholds. A threshold power value of 4 was selected to balance scale-free topology and gene connectivity stability. (C) Gene clustering dendrogram: hierarchical clustering of genes within the co-expression network demonstrates module distribution. (D) Module-clinical trait relationships: color gradients indicate correlation strength between modules and recurrence/non-recurrence status.

Figure 1B demonstrates the soft-thresholding selection process. A soft threshold power (power =4) was selected to balance scale-free topology (fit index ≥0.9) and gene connectivity stability, thereby ensuring the robustness of co-expression network construction. Figure 1C displays the hierarchical clustering dendrogram generated using the selected soft threshold (power =4), confirming the appropriate threshold cutoff and revealing sample clustering patterns.

Figure 1D illustrates the correlation between modules and recurrence status. MEred exhibited a significantly positive correlation with recurrence status (r=0.16, P=0.006), whereas the blue module (MEblue) displayed a significantly negative correlation (r=−0.24, P=3e−5). Other modules (e.g., MEbrown, MEgreen, etc.) demonstrated weaker and non-significant correlations. This analysis indicates that specific gene modules may play critical roles in BCa recurrence, providing valuable insights for further mechanistic investigations.

Univariate Cox regression analysis: screening recurrence-associated genes from MEred

To identify genes associated with BCa recurrence within MEred, we performed univariate Cox regression analysis. This analysis revealed that 24 genes significantly correlated with BCa recurrence, yielding P values less than 0.001. These genes, including OLFML3, LRRC10B, RBP1, GLCE, TCF4, LEF1, TPD52L1, CNTN1, and PDE5A, among others, are detailed in Figure 2 and Table S2. These genes may serve as important biomarkers for predicting BCa recurrence.

Figure 2 Forest plot of HRs for the 24 genes identified through univariable Cox regression analysis. HR, hazard ratio.

Screening of prognostic candidate genes via LASSO-penalized Cox regression model

Further analysis employing a LASSO-penalized Cox regression model identified eight genes with significant prognostic value: OLFML3, GLCE, TCF4, TPD52L1, PAM, SIX1, TSPAN5, and TMEM158, along with their corresponding regression coefficients, as outlined in Table S3. These genes represent a potential set of prognostic markers for BCa.

Kaplan-Meier survival curve analysis in TCGA and GEO databases

Kaplan-Meier survival analysis was performed to assess the prognostic implications of the identified genes within both the TCGA and GEO datasets. In the TCGA dataset, the high-risk group, characterized by elevated risk scores, demonstrated a significantly lower survival probability than the low-risk group (P<0.0001), as presented in Figure 3A. A similar trend was observed in the GEO dataset (GSE13507), where the high-risk group exhibited reduced survival (P=0.007), as illustrated in Figure 3B. Furthermore, a third independent GEO cohort (GSE31684) corroborated these findings, revealing statistically significant survival differences between the risk groups (P=0.047), as depicted in Figure 3C. These results affirm the reliability and robustness of our risk score model across diverse cohorts.

Figure 3 Survival analysis results based on the LASSO-Cox regression model, validated across different datasets. (A) In the TCGA dataset, there is a significant survival difference between the high-risk and low-risk groups. (B) The model’s efficacy was validated in the GEO database (GSE13507), showing a significant difference between the two groups. (C) Similar survival differences were further validated in another GEO dataset (GSE31684). GEO, Gene Expression Omnibus; LASSO, least absolute shrinkage and selection operator; OS, overall survival; TCGA, The Cancer Genome Atlas.

Expression differences and survival analysis of eight candidate genes across BCa clinical features

The expression and prognostic significance of the eight candidate genes were assessed using Kaplan-Meier survival curves for the high- and low-risk BCa groups. For example, survival analysis for the GLCE gene revealed a marginal survival difference between the high- and low-expression groups (P=0.07), which did not achieve statistical significance (Figure 4A). In contrast, OLFML3 exhibited a significant survival difference between the high- and low-expression groups (P=0.003), suggesting its clinical relevance (Figure 4B). Similarly, PAM demonstrated a statistically significant survival difference (P=0.02), further supporting its prognostic potential (Figure 4C). Nevertheless, the SIX1 gene did not display any significant survival difference (P=0.91), indicating its limited prognostic relevance (Figure 4D). Among the other genes, TCF4 and TMEM158 exhibited significant survival curve divergences (P=0.01 and P=0.03, respectively), suggesting their potential as predictive biomarkers (Figure 5A,5B). Additionally, TPD52L1 demonstrated a statistically significant survival distinction (P=0.02), underscoring its clinical applicability in prognosis evaluation (Figure 5C). Although TSPAN5 displayed observable trends, the survival difference was not statistically significant (P=0.43), suggesting limited utility (Figure 5D).

Figure 4 Survival analysis results of genes associated with BCa prognosis. (A) Survival analysis of the GLCE gene, showing its association with BCa prognosis (P=0.07). (B) Survival analysis of the OLFML3 gene, demonstrating its association with BCa prognosis (P=0.003). (C) Survival analysis of the PAM gene, assessing its association with BCa prognosis (P=0.02). (D) Survival analysis of the SIX1 gene, evaluating its association with BCa prognosis (P=0.91). BCa, bladder cancer; OS, overall survival.

Figure 5 Presents survival analysis results of genes associated with BCa prognosis. (A) Survival analysis of the TCF4 gene, showing its association with BCa prognosis (P=0.01). (B) Survival analysis of the TMEM158 gene, demonstrating its association with BCa prognosis (P=0.03). (C) Survival analysis of the TPD52L1 gene, assessing its association with BCa prognosis (P=0.02). (D) Survival analysis of the TSPAN5 gene, evaluating its association with BCa prognosis (P=0.43). BCa, bladder cancer; OS, overall survival.

Differential gene expression analysis between high and low tumor stromal BCa subtypes

Cell type annotation using the UMAP algorithm was conducted to visualize the spatial distribution of cellular populations in high and low tumor stromal BCa samples, as illustrated in Figure 6A. This analysis revealed distinct populations of T cells, epithelial cells, fibroblasts, macrophages, monocytes, and plasma cells, establishing a foundation for investigating cellular heterogeneity within the tumor microenvironment. The expression of the eight candidate genes was further analyzed in high and low tumor stromal BCa samples. Seven genes, namely, TSPAN5, OLFML3, TPD52L1, PAM, TMEM158, GLCE, and SIX1, exhibited statistically significant differences in expression (P<0.05), as shown in https://cdn.amegroups.cn/static/public/tau-2025-365-1.xlsx. Notably, OLFML3 demonstrated a striking difference in expression between high and low tumor stromal groups, with 65.4% expression in high stromal groups compared with 38.4% in low stromal groups, suggesting its potential as a molecular marker for distinguishing tumor stromal content (Figure 6B). Further analysis of OLFML3 expression across various cell types revealed significantly higher expression levels in macrophages and fibroblasts than in other cell types, highlighting its critical role in the tumor microenvironment (Figure 6C). This was substantiated by a spatial expression profile of OLFML3 across tumor tissue, where highly expressed cells were predominantly localized in macrophage and fibroblast regions (Figure 6D).

Figure 6 Single-cell workflow diagram. (A) UMAP plot displaying the distribution of distinct cell types (e.g., T cells, epithelial cells, fibroblasts, macrophages, plasma cells, etc.) between bladder cancers with HTS and LTS. (B) Comparative analysis demonstrating differential gene expression profiles between the HTS and LTS groups. (C) Absolute expression levels of the OLFML3 gene across distinct cell types, including T cells, epithelial cells, fibroblasts, macrophages, monocytes, and plasma cells. (D) Spatial distribution of OLFML3 expression in the single-cell map, annotated with cell-type localization and depicting expression intensity within cellular populations. FC, fold change; HTS, high tumor stroma; LTS, low tumor stroma; UMAP, uniform manifold approximation and projection.

Gene expression profiles of eight genes in BCa: association with invasiveness, histological grade, and tumor stage

Using the TCGA database, the expression of the eight candidate genes was analyzed in relation to clinical features such as invasiveness, histological grade, and tumor stage in BCa. For instance, genes such as OLFML3, SIX1, and TMEM158 were significantly upregulated in invasive BCa, with OLFML3 exhibiting the most pronounced differential expression, as depicted in Figure 7A. Likewise, OLFML3, TCF4, TPD52L1, PAM, SIX1, TSPAN5, and TMEM158 were significantly overexpressed in high-grade BCa, with OLFML3 displaying the most significant upregulation (Figure 7B). Expression differences across various tumor stages (T stages) revealed that genes such as OLFML3, GLCE, TCF4, TPD52L1, PAM, and TMEM158 were significantly upregulated in the T2, T3, and/or T4 stages. Nonetheless, only OLFML3 exhibited high expression across all stages (T2, T3, and T4), as shown in Figure 8A. Additionally, OLFML3 demonstrated significantly higher expression in the N1 stage, with a notable, statistically significant differential expression between N0 and N1 (P=0.004), further supporting its role in BCa progression (Figure 8B).

Figure 7 Violin plots of gene expression in relation to clinical variables. (A) Differential expression of eight genes in BCa subtypes (invasive vs. non-invasive). (B) Differential expression of eight genes across histological grades (high-grade vs. low-grade). BCa, bladder cancer.

Figure 8 Violin plots of gene expression in relation to clinical variables. (A) Differential expression of eight genes across tumor T stages (T1, T2, T3, T4). (B) Differential expression of eight genes across node N stages.

Successful prediction of OLFML3 high/low expression using an RF model

Finally, we developed an RF model based on pathological data to predict high and low OLFML3 expression levels in clinically diagnosed patients with BCa. Figure 9A illustrates the workflow for developing the ResNet50 model, employed to differentiate tumor and normal tissue sections in patients with BCa. This process encompassed model training, prediction, and feature extraction. Optimal parameters were determined via grid search, resulting in the following settings: n_estimators =91, max_depth =5, min_samples_split =10, min_samples_leaf =1, and bootstrap = TRUE. The model exhibited excellent performance, yielding AUC values of 0.86, 0.87, 0.82, and 0.89 for the training, internal validation, external validation A1, and external validation B1 sets, respectively (Figure 9B). The accuracy metrics for the respective sets were 0.77, 0.81, 0.73, and 0.77 (Figure 9C). Additionally, prediction probability waterfall plots across four datasets (training, validation, external1, and external2) evidently demonstrated the model’s ability to discriminate between low and high OLFML3 expression groups (Figure 9D). This underscores the potential clinical diagnostic utility of the RF model in predicting OLFML3 expression levels and stratifying patients accordingly.

Figure 9 Deep learning model construction workflow for pathological features using a RF approach. (A) Workflow of the deep learning model for processing H&E-stained whole-slide images (20×), including image preprocessing, feature extraction, and model training. (B) ROC curve of the model evaluating OLFML3 expression prediction performance in BCa across the training set, internal validation set, and external validation sets A1 and B1. (C) Performance metrics of the model (including AUC, accuracy, sensitivity, specificity, PPV, and NPV) in the training set, internal validation set, and external validation sets A1 and B1. (D) Visualization of predicted probabilities for individual patients generated by the model. AUC, area under the curve; BCa, bladder cancer; IQR, interquartile range; NPV, negative predictive value; PPV, positive predictive value; RF, random forest; ROC, receiver operating characteristic; SD, standard deviation; WSI, whole-slide image.

Discussion

Cancer development is characterized by the evasion of immune surveillance by cancer cells, a phenomenon known as the immunoediting process, alongside the accumulation of genomic instability. These factors are intricately linked to the acquisition of malignant phenotypes and subsequent metastatic spread, serving as critical drivers of advanced tumorigenesis. In BCa, genomic instability is particularly correlated with the risk of recurrence. Therefore, exploring genomic-level biomarkers for BCa recurrence is imperative. For example, Russo et al. identified the loss of the Y chromosome as a potential early indicator of BCa and a biomarker for predicting pan-genomic instability (15). To investigate genomic features associated with BCa recurrence, we initially selected patients with BCa from the TCGA public database and categorized them into two groups: those with recurrence within 1 year and those without recurrence after 1 year (Figure 1). Subsequently, we employed WGCNA to cluster genes into modules based on similar expression patterns across samples, correlating these modules with clinical traits (16). This analysis revealed gene modules significantly associated with 1-year BCa recurrence, with MEred exhibiting the strongest positive correlation with recurrence. In the subsequent univariate Cox regression analysis, we examined the prognostic potential of 24 genes related to the survival time of patients with BCa (Figure 2). The LASSO-penalized Cox regression model demonstrated strong predictive performance for recurrence risk. For instance, a previous study indicated that a penalized Cox regression model outperformed the LASSO group and RF models in predicting recurrent membranous nephropathy (17). In this study (Figure 3), we similarly applied the LASSO-penalized Cox regression model to identify eight genes significantly associated with BCa prognosis (e.g., OLFML3, GLCE, TCF4, etc.) and constructed a risk score model based on these genes [risk score calculation formula: LYSscore = Σ(coefficient × gene expression value)] (18). Validation in the TCGA dataset revealed a significant survival difference between high- and low-risk patient groups, confirming the model’s efficacy in predicting BCa prognosis. Moreover, the model exhibited robust predictive performance in an independent GEO dataset, validating its broad applicability and reliability.

We conducted independent survival analyses of the eight prognosis-related genes identified using the LASSO-penalized Cox regression model. Using Kaplan-Meier survival curve analysis, we assessed the relationship between the expression level of each gene and the survival of patients with BCa, establishing the independent prognostic values of these genes. Notably, genes such as OLFML3, PAM, TCF4, TMEM158, and TPD52L1 yielded independent prognostic values in BCa (Figures 4,5). To further investigate the role of specific genes in assessing BCa recurrence, we collected specimens from six groups of patients with BCa at Guangdong Provincial Second People’s Hospital for single-cell sequencing. We compared the gene expression differences between high and low tumor stromal BCa samples. As shown in Figure 6B, the results indicate that OLFML3 was highly expressed in high tumor stromal BCa samples, yielding significant expression differences compared with that in low tumor stromal samples. Furthermore, as shown in Figure 6C,6D, OLFML3 expression was notably elevated in fibroblasts, and cancer-associated fibroblasts (CAFs) exhibited a strong correlation with BCa recurrence. Fibroblasts significantly contribute to the tumor microenvironment by secreting cytokines, matrix metalloproteinases, and other molecules, thus promoting tumor microenvironment remodeling, which supports BCa recurrence and progression. These findings suggest that OLFML3 potentially participates in the BCa recurrence process by modulating CAF function, positioning it as a potential biomarker for BCa recurrence and offering novel avenues for future prevention and treatment research (19). To further assess the correlation of OLFML3 with other clinical factors (e.g., tumor stage, grade, etc.), we observed a significant increase in OLFML3 expression in invasive BCa compared with that in non-invasive BCa (Figure 7A). Moreover, a notable survival disparity exists between invasive and non-invasive BCa, with approximately 70% of patients with BCa presenting as NMIBC, which has a favorable prognosis and a 5-year OS rate approximating 90%. In contrast, invasive BCa is associated with a lower survival rate, with several patients diagnosed with muscle-invasive BCa (MIBC) progressing to distant metastatic disease, indicating more severe clinical features (20). As illustrated in Figure 7B, OLFML3 expression significantly increased in high-grade BCa compared with that in low-grade BCa. High-grade BCa is typically associated with a higher incidence of recurrence and progression (21). As depicted in Figure 8A, OLFML3 expression gradually increased across stages T2, T3, and T4. BCa is categorized into NMIBC (stages Tis, Ta, and T1) and MIBC (≥ stage T2). BCa beyond stage T2 indicates infiltration of the bladder muscle layer, suggesting disease progression and malignancy. Clinically, radical cystectomy combined with pelvic lymph node dissection is the preferred treatment for MIBC. The correlation between OLFML3 expression and the T stage of BCa underscores its potential role in the malignant progression of the disease, particularly as elevated OLFML3 expression in high-stage (≥ stage T2) BCa may closely relate to tumor invasiveness and recurrence (22). In summary, incorporating OLFML3 gene expression levels into routine clinical practice may enhance the early identification of patients with BCa at high risk of recurrence following transurethral resection of bladder tumor or radical cystectomy. Such early identification allows for the implementation of targeted interventions and preventive strategies—such as the consideration of more aggressive adjuvant therapies—while enabling enhanced imaging surveillance frequency and intensity for high-risk groups. It further supports the implementation of closer follow-up intervals for high-risk patients to achieve earlier detection of recurrence or progression. This biomarker-based risk stratification paradigm, exemplified by OLFML3, mirrors the importance of utilizing indicators like the systemic inflammation response index in clinical risk prediction, facilitating the early identification of high-risk patients to inform personalized monitoring and therapeutic decisions, ultimately improving patient prognosis (23).

Digital pathology (DP) addresses diagnostic challenges by digitizing slides, facilitating remote consultations and standardized diagnostic protocols (24). When combined with deep learning, DP can automate image analysis, enhance tumor detection capabilities, and improve diagnostic reproducibility (25). This integration signifies a paradigm shift that could bridge healthcare gaps, ensuring that high-quality diagnostic services are accessible to all. This approach is particularly crucial for early tumor diagnosis, prognostic evaluation, and the personalized development of treatment plans, displaying great promise in the field of urological tumors. For example, one study explored pathology-based deep learning features for predicting basal and luminal subtypes in BCa, demonstrating the effectiveness of machine learning models in this context (14). Another study developed an artificial intelligence-based diagnostic model for detecting lymph node metastases in BCa, showcasing significant clinical potential for improving the accuracy and efficiency of pathologists’ work (26). As depicted in Figure 9, we constructed an RF model using deep learning pathology features to predict OLFML3 expression levels in patients with BCa. This innovative application highlights the potential of deep learning pathology features in predicting OLFML3 expression and provides a valuable tool for the clinical prediction of BCa recurrence, augmenting personalized management and prognostic evaluation for patients with BCa.

Limitations

We must acknowledge certain limitations. First, the association between OLFML3 and BCa recurrence necessitates validation using larger clinical samples. Although external validation was conducted using the GEO dataset, the BCa single-cell and bladder pathology databases from GD2H, as well as the BCa pathology database from STPH, most of the patient data appeared to derive from Chinese cohorts. Consequently, the applicability and generalizability of the findings across diverse global populations should be approached with caution. Specifically, factors such as patient race, geographical location, and lifestyle habits might have influenced the results, warranting further exploration and validation in future research. Second, while OLFML3 may be associated with BCa recurrence, its specific mechanism requires further experimental validation. Future studies should investigate its role in BCa recurrence through functional studies and animal model research.

Conclusions

In this study, we analyzed patients with BCa from the TCGA database using WGCNA, univariate analysis, and a risk scoring system based on the LASSO-penalized Cox regression model. This analysis identified eight prognosis-related genes, including OLFML3, GLCE, and TCF4. Kaplan-Meier survival analysis demonstrated that these genes significantly differentiate between high- and low-risk groups, with statistically significant survival differences observed in both the TCGA and GEO datasets. We subsequently subjected each of the eight genes to Kaplan-Meier survival analysis. Through single-cell analysis, we found that OLFML3 may play a crucial role in BCa recurrence and the tumor microenvironment, particularly owing to its elevated expression in CAFs, which are associated with recurrence. Moreover, OLFML3 expression was correlated with clinical pathological features of BCa, such as tumor stage, invasiveness, and pathological grade. The integration of DP and deep learning enhances the early diagnosis and prognostic evaluation of BCa. Our RF model successfully predicted OLFML3 expression levels, further validating the potential of deep learning in predicting BCa recurrence.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tau.amegroups.com/article/view/10.21037/tau-2025-365/rc

Data Sharing Statement: Available at https://tau.amegroups.com/article/view/10.21037/tau-2025-365/dss

Peer Review File: Available at https://tau.amegroups.com/article/view/10.21037/tau-2025-365/prf

Funding: This work was supported by the National Natural Science Foundation of China (Grant No. 82302304); the Doctoral Workstation Foundation of Guangdong Second Provincial General Hospital, China (Grant No. 2022BSGZ011); the Elevate Engineering Foundation of Guangdong Second Provincial General Hospital, China (Grant No. TJGC2022009); the 2024 Hospital Pazhou Scholar Program Research Project (Grant No. PZ2024-001); and the Science and Technology Program of Guangzhou, China (Grant Nos. 202201020582, 2023A03J0254 and 2024A04J4159).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tau.amegroups.com/article/view/10.21037/tau-2025-365/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. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The study was approved by the respective committees at STPH (Approval No. 24KT68) and GD2H (Approval No. 2024-KY-KZ-128-01). Informed consent was collected from all participants.

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. Borgna V, Lobos-González L, Guevara F, et al. Targeting antisense mitochondrial noncoding RNAs induces bladder cancer cell death and inhibition of tumor growth through reduction of survival and invasion factors. J Cancer 2020;11:1780-91. [Crossref] [PubMed]
2. Kiss B, van den Berg NS, Ertsey R, et al. CD47-Targeted Near-Infrared Photoimmunotherapy for Human Bladder Cancer. Clin Cancer Res 2019;25:3561-71. [Crossref] [PubMed]
3. Shi ZD, Han XX, Song ZJ, et al. Integrative multi-omics analysis depicts the methylome and hydroxymethylome in recurrent bladder cancers and identifies biomarkers for predicting PD-L1 expression. Biomark Res 2023;11:47. [Crossref] [PubMed]
4. Lu Y, Liu P, Wen W, et al. Cross-species comparison of orthologous gene expression in human bladder cancer and carcinogen-induced rodent models. Am J Transl Res 2010;3:8-27. [PubMed]
5. Nikas IP, Seide S, Proctor T, et al. The Paris System for Reporting Urinary Cytology: A Meta-Analysis. J Pers Med 2022;12:170. [Crossref] [PubMed]
6. Zhang S, Zhang J, Zhang Q, et al. Identification of Prognostic Biomarkers for Bladder Cancer Based on DNA Methylation Profile. Front Cell Dev Biol 2022;9:817086. [Crossref] [PubMed]
7. Fiorentino V, Pizzimenti C, Franchina M, et al. Bladder Epicheck Test: A Novel Tool to Support Urothelial Carcinoma Diagnosis in Urine Samples. Int J Mol Sci 2023;24:12489. [Crossref] [PubMed]
8. Pepe L, Fiorentino V, Pizzimenti C, et al. The Simultaneous Use of Bladder Epicheck(®) and Urinary Cytology Can Improve the Sensitivity and Specificity of Diagnostic Follow-Up of Urothelial Lesions: Up-to-Date Data from a Multi-Institutional Cohort. Diseases 2024;12:219. [Crossref] [PubMed]
9. Lobo A, Collins K, Kaushal S, et al. Advances, recognition, and interpretation of molecular heterogeneity among conventional and subtype histology of urothelial carcinoma (UC): a survey among urologic pathologists and comprehensive review of the literature. Histopathology 2024;85:748-59. [Crossref] [PubMed]
10. Pierconti F, Martini M, Cenci T, et al. Methylation study of the Paris system for reporting urinary (TPS) categories. J Clin Pathol 2021;74:102-5. [Crossref] [PubMed]
11. Dudley JC, Schroers-Martin J, Lazzareschi DV, et al. Detection and Surveillance of Bladder Cancer Using Urine Tumor DNA. Cancer Discov 2019;9:500-9. [Crossref] [PubMed]
12. Sylvester RJ, van der Meijden AP, Oosterlinck W, et al. Predicting recurrence and progression in individual patients with stage Ta T1 bladder cancer using EORTC risk tables: a combined analysis of 2596 patients from seven EORTC trials. Eur Urol 2006;49:466-5; discussion 475-7. [Crossref] [PubMed]
13. Wang X, Du Y, Yang S, et al. RetCCL: Clustering-guided contrastive learning for whole-slide image retrieval. Med Image Anal 2023;83:102645. [Crossref] [PubMed]
14. Zheng Z, Dai F, Liu J, et al. Pathology-based deep learning features for predicting basal and luminal subtypes in bladder cancer. BMC Cancer 2025;25:310. [Crossref] [PubMed]
15. Russo P, Bizzarri FP, Filomena GB, et al. Relationship Between Loss of Y Chromosome and Urologic Cancers: New Future Perspectives. Cancers (Basel) 2024;16:3766. [Crossref] [PubMed]
16. Tian Z, He W, Tang J, et al. Identification of Important Modules and Biomarkers in Breast Cancer Based on WGCNA. Onco Targets Ther 2020;13:6805-17. [Crossref] [PubMed]
17. Chung EYM, Blazek K, Teixeira-Pinto A, et al. Predictive Models for Recurrent Membranous Nephropathy After Kidney Transplantation. Transplant Direct 2022;8:e1357. [Crossref] [PubMed]
18. Song D, Zhao L, Zhao G, et al. Identification and validation of eight lysosomes-related genes signatures and correlation with immune cell infiltration in lung adenocarcinoma. Cancer Cell Int 2023;23:322. [Crossref] [PubMed]
19. Liang T, Tao T, Wu K, et al. Cancer-Associated Fibroblast-Induced Remodeling of Tumor Microenvironment in Recurrent Bladder Cancer. Adv Sci (Weinh) 2023;10:e2303230. [Crossref] [PubMed]
20. Funt SA, Rosenberg JE. Systemic, perioperative management of muscle-invasive bladder cancer and future horizons. Nat Rev Clin Oncol 2017;14:221-34. [Crossref] [PubMed]
21. Yun SJ, Kim SK, Kim WJ. How do we manage high-grade T1 bladder cancer? Conservative or aggressive therapy? Investig Clin Urol 2016;57:S44-51. [Crossref] [PubMed]
22. Su H, Jiang H, Tao T, et al. Hope and challenge: Precision medicine in bladder cancer. Cancer Med 2019;8:1806-16. [Crossref] [PubMed]
23. Russo P, Foschi N, Palermo G, et al. SIRI as a biomarker for bladder neoplasm: Utilizing decision curve analysis to evaluate clinical net benefit. Urol Oncol 2025;43:393.e1-8. [Crossref] [PubMed]
24. Janowczyk A, Zlobec I, Walker C, et al. Swiss digital pathology recommendations: results from a Delphi process conducted by the Swiss Digital Pathology Consortium of the Swiss Society of Pathology. Virchows Arch 2024;485:13-30. [Crossref] [PubMed]
25. Yoo JW, Koo KC, Chung BH, et al. Deep learning diagnostics for bladder tumor identification and grade prediction using RGB method. Sci Rep 2022;12:17699. [Crossref] [PubMed]
26. Wu S, Hong G, Xu A, et al. Artificial intelligence-based model for lymph node metastases detection on whole slide images in bladder cancer: a retrospective, multicentre, diagnostic study. Lancet Oncol 2023;24:360-70. [Crossref] [PubMed]