The tissue expression data and clinical information of BLCA patients were downloaded from The Cancer Genome Atlas (TCGA) database (https://www.cancer.gov/ccg/research/genome-sequencing/tcga) and the ArrayExpress database (https://www.ebi.ac.uk/biostudies/arrayexpress). The dataset underwent the following processing steps: (1) exclusion of samples lacking clinical follow-up information; (2) conversion of gene identifiers to gene symbols; (3) selection of the highest expression value in cases of multiple entries for the same gene symbol; (4) removal of samples with missing expression profile data. After preprocessing, 359 samples were obtained from TCGA-BLCA and 476 samples were obtained from the ArrayExpress database (Table 1).

We retrieved three Homo sapiens-derived BLCA tissues with comprehensive staging information from the GEO database (GSE129845) and analyzed their single-cell RNA sequencing data using the R package Seurat. Quality control filters excluded cells with fewer than 200 genes or over 15% mitochondrial genes, resulting in the retention of 13,490 cells for further analysis. Data normalization was achieved using the LogNormalize method, and batch effects were mitigated with the Harmony package. The FindVariableFeatures function identified the top 2,000 variable genes. Dimensionality reduction was performed through principal component analysis, t-SNE, and UMAP, followed by cell clustering using the FindNeighbors and FindClusters functions. Cell annotation was facilitated by R package SingleR.

Genes related to the FAM pathways (HALLMARK_FATTY_ACID_METABOLISM) were obtained from the Molecular Signature Database v7.0 (MSigDB). A total of 158 FAM-related genes were included in the analysis (Supplementary Table S1).

The TCGA expression profile data underwent filtration to exclude genes with expression levels below one in more than 50% of samples, followed by univariate Cox analysis to identify prognostic-related FAM genes using a significance threshold of p < 0.05. ConsensusClusterPlus (v1.48) was utilized for uniform clustering of TCGA samples, utilizing D2 and Euclidean distance as the clustering algorithm and measure, with parameters set as pFeature = 1, rep = 100, distance = “Spearman,” and pItem = 0.8. The limma package was used to assess molecular subtype discrepancies and conduct functional enrichment analysis, whereas DAVID was implemented to assess significantly enriched pathways, including KEGG and GO pathways, across distinct BLCA groups. Pathways were considered enriched if they met the criteria of p < 0.05 and a false discovery rate (FDR) < 0.05.

To analyze tumor immune infiltration, we performed a multi-algorithm assessment using the TCGA-BLCA cohort normalized expression matrix [log2(TPM+1)]. First, stromal and immune scores were globally evaluated via the ESTIMATE package (v1.0.13) under default parameters. Next, immune cell subtype deconvolution was conducted using CIBERSORT with the LM22 signature and 1,000 permutations, retaining only samples with a CIBERSORT output p-value <0.05. Additionally, absolute immune and stromal cell abundances were quantified using the MCP-counter package (v1.2.0) and its predefined signatures. Between-group differences were assessed via the Wilcoxon rank-sum test; pairwise correlations were examined with Spearman’s rank correlation. For all multiple comparisons, the false discovery rate (FDR) was controlled using the Benjamini–Hochberg procedure, with an adjusted p-value (FDR) < 0.05 being considered significant. To enhance robustness, final interpretations focused exclusively on immune infiltration signals that were consistently observed across all three computational methods (ESTIMATE, CIBERSORT, and MCP-counter).

Accurate prognosis is crucial for tailoring cancer treatments and improving patient outcomes. Traditional clinical factors offer limited precision, but molecular signatures can enhance prognostic accuracy by capturing tumor diversity.

Our study developed a predictive risk model using FAM-related genes. We began by extracting FAM gene expression data from a public dataset. Univariate Cox regression identified survival-associated FAM genes, which were further refined through LASSO regression to prevent overfitting. We then constructed a prognostic model using multivariate Cox regression with the selected genes. Patient risk scores were calculated based on gene expression levels and model coefficients. Patients were categorized into high- and low-risk groups using the median risk score.

We evaluated the model’s performance using ROC and Kaplan–Meier analyses. This FAM gene-based risk model aims to improve survival predictions and inform personalized treatment strategies for cancer patients.

We performed bootstrap resampling (1,000 iterations) on our training cohort to obtain more reliable performance estimates. The results showed a C-index of 0.604 (SD: 0.006), with a 95% confidence interval of 0.591–0.611. In the study, the 359 samples from the TCGA dataset were randomly divided into a training set (180 samples) and a test set (179 samples). Training and test sets were selected based on two criteria: balanced patient demographics and clinical outcomes and similar distribution of binary samples after gene expression clustering.

LASSO regression was used to identify prognostic genes and optimize the risk model. This method reduces coefficients and selects relevant variables, effectively addressing multicollinearity and promoting sparsity in the data. LASSO Cox regression analysis was performed using the glmnet R-package. The optimal model was selected through five-fold cross-validation, and the number of target genes was determined based on confidence intervals at each lambda value.

Total RNA was extracted from the T24 and 5637 cells using TRIzol reagent (Invitrogen) and reverse-transcribed into cDNA with the kit of HiScript II Q RT SuperMix for qPCR (Vazyme). The cDNA was subjected to qPCR analysis using the Universal SYBR qPCR Master Mix (Vazyme).

T24 and 5637 cells were seeded in 96-well plates at a density of 3,000 cells per well in Dulbecco’s Modified Eagle Medium (DMEM) supplemented with 10% fetal bovine serum (FBS) and 1% penicillin–streptomycin, respectively. After 48 and 72 h of incubation, 10 μL of Cell Counting Kit-8 (CCK-8, BioSharp) was added to each well, followed by incubation at 37 °C in the dark for 1 h. Cellular viability was assessed by measuring the absorbance at 450 nm. All experiments were performed in triplicate to ensure reproducibility.

T24 and 5637 cells were seeded in 6-well plates and cultured until a confluent monolayer was formed. A uniform scratch was introduced into the monolayer using a sterile 200-μL pipette tip. The cells were then washed thrice with PBS to remove detached cells and debris, followed by the addition of DMEM supplemented with 10% FBS and 1% penicillin–streptomycin. Wound healing progression was monitored and photographed at 0, 24, and 48 h using an IX71 inverted microscope (Leica Corporation). The images were analyzed using ImageJ software (v1.80), and the wound closure was calculated using the following formula: [(area at 0 h – area at 24/48 h)/area at 0 h] × 100. All the experiments were performed in triplicate to ensure reproducibility.

For the analysis of two or more continuous variables with normal distribution, t-tests or analysis of variance (ANOVA) were applied. Statistical significance was determined by an adjusted p-value <0.05.

A total of 359 patients were identified from the TCGA-BLCA database as the training and validation cohort, and 476 patients were identified from the ArrayExpress-mRNAseq-E-MTAB-4321 database as the external validation cohort. A total of 158 FAM-related genes were collected from the Molecular Signature Database v7.0 (MSigDB). The workflow of this study is shown in Figure 1.

We first obtained the expression of 158 FAM-related genes from the TCGA-BLCA expression profile data. A total of 18 genes associated with BLCA prognosis (p < 0.05) were identified through univariate Cox analysis by R (Supplementary Table S2). Based on the expression of the 18 genes, we clustered BLCA patients using non-negative matrix factorization (NMF). We chose the optimal clustering of k = 2 by synthesizing and the residuals sum of squares (RSS), and two clusters (C1 and C2) were obtained (Figures 2A–C). We then compared the clinical characteristics of the two clusters and found that C1 had a higher survival rate than C2 (Figures 2D,E). In addition, high-grade tumor was more enriched in C2 than in C1, indicating an enhanced cancer progression (Figure 2F). C1 had more early-stage tumors, including stages I and II, and less late-stage tumors, including stages III and IV, than C2, all of which indicated a better prognosis of C1 (Figure 2G).

Thorsson et al. (2018) identified six types of inter-tumoral immune stages: cluster 1 (wound healing), cluster 2 (IFN-γ dominant), cluster 3 (inflammatory), cluster 4 (lymphocyte depleted), cluster 5 (immunologically quiet), and cluster 6 (TGF-β dominant), of which cluster 3 has the best prognosis. We found that C1 was more correlated with clusters 1, 3, 4, and 5, whereas C2 was more correlated with cluster 2 (Figure 2H), indicating a better prognosis of C1, which is consist with the Kaplan–Meier (KM) curves. C1 showed a better prognosis indicated by overall survival (OS) and disease-specific survival (DSS) (Figures 2I,J).

We next evaluated the immune scores of C1 and C2 using the ESTIMATE package in R software. The result showed that C1 had a lower immune score, stromal score, and ESTIMATE score than C2 (Figure 3A). Further studies using R software packages CIBERSORT and MCPcounter indicated that C1 had higher scores for memory B cells, plasma cells, follicular helper T cells, and activated dendritic cells, whereas it had lower scores for M0 macrophage, M2 macrophages, monocytic lineage, and cytotoxic lymphocytes (Figures 3B,C). Heatmap also showed the lower immune score of C1 and different specific immune cell scores between C1 and C2 (Figure 3D).

To verify the biological characteristics between the two clusters, we next calculated the DEGs between C1 and C2 using the limma package and applied a filter with the criteria of |log2FC| > 1 and FDR < 0.01. The results showed 1,886 DEGs, with 632 upregulated and 1,254 downregulated genes based on C1 (Figures 4A,B). The full list of DEGs is provided in Supplementary Table S3. Pathways associated with steroid metabolic process, response to xenobiotic stimulus, and hormone metabolic process were enriched in C1, whereas pathways related to chemotaxis, leukocyte migration, and extracellular structure organization were downregulated in C1 based on GO functional pathway enrichment analysis (Figures 4C,D). KEGG enrichment analysis showed that pathways associated with drug metabolism, retinol metabolism, and steroid hormone biosynthesis were upregulated in C1, whereas cytokine signaling, PI3K-AKT signaling, and cytoskeleton were downregulated in C1 (Figures 4E,F).

To construct the prognostic risk model, we categorized the 359 patients into training and validation sets randomly, with the training group containing 180 patients and the validation group containing 179 patients. Based on C1 and C2, we identified 1,886 DEGs (Figure 4; Supplementary Table S3). We then identified 295 prognosis-associated genes through univariate regression Cox risk model analysis based on the survival data (with p < 0.01 as the threshold) (Supplementary Table S4). We next applied the R-package glmnet for LASSO Cox regression analysis of the 295 prognosis-associated genes and revealed that the number of coefficients of independent variables gradually increased as lambda increased (Figure 5A). Cross-validation was performed to calculate the confidence intervals under each lambda, and it revealed the optimal model at lambda = 0.07290695, with 16 genes identified (Figure 5B; Supplementary Table S5).

To identify the optimal predictive gene signature, we performed stepwise model selection based on the Akaike information criterion (AIC) using the step AIC function from the MASS package. This procedure began with a full model containing all 16 candidate genes. At each step, the variable whose removal most significantly lowered the AIC value was eliminated, thereby refining the model toward a parsimonious structure without substantially compromising the goodness-of-fit. Through this iterative process, the initial 16-gene set was reduced to a final signature comprising four genes: PATZ1, TTC6, AEBP1, and MAOA. PATZ1, TTC6, and MAOA were positively associated with BLCA prognosis, as higher expression of these genes was correlated with better survival outcomes (Figures 5C,D,F). AEBP1 was negatively associated with BLCA prognosis, as higher expression of the gene had poor survival outcomes (Figure 5E).

We next calculated the RiskScore of the 180 samples in the training cohort based on the expression of the four indicated genes using the ggRISK package. The final four-gene signature formula is as follows: FAMR = [−0.606 × log2(PATZ1+1)] + [0.0594 × log2(TTC6+1)] + [0.5085 × log2(AEBP1+1)] + [−0.0352 × log2(MAOA+1)]. The results showed that higher RiskScore represented poor prognosis (Figure 5G). Higher FAMR was associated with lower expression of PATZ1, TTC6, and MAOA and higher expression of AEBP1, which was consistent with the KM curve results (Figure 5G). Then, we evaluated the area under the curve (AUC) of prognostic prediction efficiency at 1, 3, and 5 years by R-wrapper time ROC analysis. The results showed that the AUC values of the indicated models were higher than 0.65 (Figure 5H). Furthermore, the KM analysis showed that the low FAMR group had a better prognosis than the high FAMR group (Figure 5I). These results indicated the clinical application value of FAMR.

To verify the FAMR model, we performed analysis using the 179 test patients’ data alone and the data of all 359 patients as the test cohort. Consistently, higher FAMR was correlated with poor prognosis, with lower expression of PATZ1, TTC6, and MAOA and higher expression of AEBP1 (Figures 6A,D). The AUC of prognostic prediction efficiency at 1, 3, and 5 years was higher than 0.65 in both the 179 and 359 validation cohorts (Figures 6B,E). The KM analysis results also showed that the low FAMR group was associated with a better prognosis compared with the high FAMR group (Figures 6C,F). To further validate the FAMR model using the external validation cohort, we collected data of the 476 patients from the ArrayExpress-mRNAseq-E-MTAB-4321 database and performed the FAMR analysis. The results showed that higher FAMR was associated with poor prognosis based on the external validation cohort (Figure 6G). The AUC of prognostic prediction efficiency at 1, 3, and 5 years was higher than 0.75 (Figure 6H). Consistently, the low FAMR group had a significantly better prognosis than the high FAMR group, as shown by the KM analysis (Figure 6I). Together, we verified the FAMR model using internal and external validation cohorts, and the results showed the stability and universality of the model.

We next studied the corresponding gene expression profiles of the GSEA to analyze the correlation between FAMR and biological functions. The ssGSEA scores of each sample were calculated, and the correlation between the FAMR and biological function was further analyzed. The results showed the top 10 significantly positively and negatively correlated pathways (Figures 7A,B). The positively correlated pathways included KEGG_FOCAL_ADHERION, KEGG_ECM_RECEPTOR_INTERACTION, KEGG_REGULATION_OF_ACTIN_CYTOSKELETON, and KEGG_PRION_DISEASES, among others, which were associated with chemotaxis, cytoskeleton, infections, and inflammations (Figures 7A,B). The negatively correlated pathways included KEGG_PEROXISOME, KEGG_TASTE_TRANSDUCTION, KEGG_GLYCEROPHOSPHOLIPID_METABOLISM, KEGG_RETINOL_METABOLISM, and KEGG_LINOLEIC_ACID_METABOILISM, among others, which were associated with fatty acid metabolism and metabolism-related diseases (Figures 7A,B). Furthermore, we assessed the correlations between FAMR and ImmuneScore, and the results revealed that the ImmuneScore, StromalScore, and ESTIMATEScore were positively related with FAMR (Figures 7C–E).

Additionally, we investigated whether gender, age, and clinical features such as the cancer stage affected FAMR. Therefore, KM analyses were performed in male female patients, as well as in patients aged >60 years and ≤60 years. The results showed that the low FAMR groups all had better prognosis than the high FAMR groups in the four indicated groups (Figures 8A–D), indicating the stable indicative ability of the FAMR model. We next analyzed the FAMR in the gender, age, and tumor stage groups and found that the female group had higher FAMR (Figure 8E), indicating a poor prognosis for females, which was consistent with the information that BLCA in women is often diagnosed at a higher stage with worse prognosis (Dobruch et al., 2025). The >60 group had higher FAMR (Figure 8F), which was consistent with the information that the elderly are more susceptible to BLCA (Van Hoogstraten et al., 2023). Furthermore, FAMR increased from stages II to IV, indicating worse prognosis along with cancer progress (Figure 8G). We further revealed that FAMR was correlated with the survival rate of BLCA by multivariate Cox regression analysis (Figures 8H,I). In addition to FAMR, factors such as the patient age, gender, and disease stage also exert varying degrees of influence on the prognosis of BLCA. Integrating these clinical variables with FAMR may thus serve as a more robust approach for predicting outcomes in BLCA patients.

The tumor microenvironment (TME) consists of diverse cell types, including cancer cells, stromal, and immune cells, which interact with each other and play critical roles in cancer. To deconvolute the cellular source of the four indicated FAMR-related genes, we analyzed their expression in publicly available single-cell RNA-seq datasets from three patients with BLCA. A total of six cell types were identified, encompassing epithelial cells, cancer-associated fibroblasts (CAFs), smooth muscle cells (CNN1, DES), endothelial cells (PECAM1, AQP1), myeloid cells (LYZ, MS4A7), and lymphocytes (Figures 9A–D). Epithelial cells were characterized by the high expression of KRT19, AGR2, and AQP3 (Figures 9C–E). CAFs were characterized by the high expression of COL1A1, PLAC9, C1R, and CXCL14 (Figures 9C–E). MZB1, CD79A, and CD3G were enriched in lymphocytes, indicating that this population encompassed B cells and T cells (Figures 9C–E). We next analyzed the expression of the four indicated FAMR-related genes and found that AEBP1 was mainly expressed by CAFs and some of the smooth muscle cells (Figures 9F,G). The expression of MAOA was mainly contributed by epithelial cells, indicating their cancer cell-derived origin (Figures 9F,G). In addition, PATZ1 and TTC6 were rarely expressed in those three patients within the BLCA dataset (Figures 9F,G). These data indicated that the prognosis-related genes AEBP1 and MAOA were CAF-derived and cancer cell-derived, respectively. Furthermore, we categorized 8,258 BLCA epithelial cells into four clusters based on the DEGs among different BLCA epithelial cell sub-clusters. A total of four cell types were identified, encompassing basal (FABP5 and AQP3), luminal (JUN and JUNB), immunity (UPK1A and UPK2), and cycle (STMN1 and HMGB2) (Figures 10A,B). MAOA is highly expressed in epithelial cells and is present across all four epithelial cell subpopulations (Figure 10C). Pseudo‐time analysis revealed that the four sub-clusters were distributed along an evolutionary pathway: immunity→ luminal/cycle/basal (Figures 10D–G). The expression of genes such as PATZ1, TTC6, AEBP1, and MAOA does not change with the pseudo-time progression of epithelial cells (Figure 10H).

To validate the function of the four identified genes in BLCA prognosis, we chose MAOA, the only gene predominantly expressed by epithelial tumor cells, for in vitro functional validation. We first silenced the expression of MAOA in human BLCA cell lines T24 and 5637 using short hairpin RNA (shRNA) (Figures 11A,B). Cell Counting Kit-8 (CCK-8) assays showed that knockdown of MAOA significantly promoted the proliferation of both T24 and 5637 cells (Figures 11C,D), indicating the inhibitory role of MAOA in BLCA proliferation. In addition, MAOA knockdown also resulted in enhanced migration of T24 (Figures 11E,F) and 5637 (Figures 11G,H) cells, which was supported by the wound healing assays. Collectively, these findings indicated that MAOA inhibited the proliferation and migration of BLCA cells, providing a mechanistic basis for the observed association between high MAOA expression and improved patient prognosis.