Molecular vasculogenic mimicry–Related signatures predict clinical outcomes and therapeutic responses in bladder cancer: Results from real-world cohorts

Bladder cancer (BLCA) is a heterogeneous disease, and there are many classical molecular subtypes that reflect tumor immune microenvironment (TME) heterogeneity but their clinical utility is limited and correct individual treatment and prognosis cannot be predicted based on them. To find reliable and effective biomarkers and tools for predicting patients’ clinical responses to several therapies, we developed a new systemic indicator of molecular vasculogenic mimicry (VM)–related genes mediated by molecular subtypes based on the Xiangya cohort and additional external BLCA cohorts using a random forest algorithm. A correlation was then done between the VM_Score and classical molecular subtypes, clinical outcomes, immunophenotypes, and treatment options for BLCA. With the VM_Score, it is possible to predict classical molecular subtypes, immunophenotypes, prognosis, and therapeutic potential of BLCA with high accuracy. The VM_Scores of high levels indicate a more anticancer immune response but a worse prognosis due to a more basal and inflammatory phenotype. The VM_Score was also found associated with low sensitivity to antiangiogenic and targeted therapies targeting the FGFR3, β-catenin, and PPAR-γ pathways but with high sensitivity to cancer immunotherapy, neoadjuvant chemotherapy, and radiotherapy. A number of aspects of BLCA biology were reflected in the VM_Score, providing new insights into precision medicine. Additionally, the VM_Score may be used as an indicator of pan-cancer immunotherapy response and prognosis.

1 Introduction

Bladder cancer is one of the most common urinary malignancies causing an estimated 81,180 new cases and 17,100 deaths in 2022 (Siegel et al., 2022). Surgery, chemotherapy, and targeted therapy are the effective treatments for BLCA, but most patients have a poor prognosis (Witjes et al., 2021). Due to the high immunogenicity of BLCA, the application of tumor immunotherapy, especially immune checkpoint inhibitors (ICI), has achieved certain results, but the overall efficacy is still not ideal (Witjes et al., 2021). This is because there are no reliable or effective biomarkers or tools for predicting patients’ clinical responses to these therapies, and only a minority of patients with BLCA respond to these therapies. Therefore, it is crucial to develop biomarkers for precise diagnosis and treatment.

There are three components to the tumor immune microenvironment (TME): cancer cells, immune cells, and the extracellular matrix. There can be differences in clinical responses to treatment based on the level of TME heterogeneity in patients with the same pathological stage and grade (Dagogo-Jack and Shaw, 2018). It is worth noting that high TME heterogeneity hinders the realization of BLCA precision medicine. As a result, understanding TME heterogeneity can shed light on many aspects of bladder cancer biology and lead to more effective bladder cancer treatment. The TME could be a promising pathway to precision medicine in BLCA through the development of new therapeutic response prediction markers and therapeutic target heterogeneity.

Vasculogenic mimicry (VM) has recently been detected in a number of malignant tumors that provide a novel strategy for treating them clinically. VM vessels are composed of endothelial tumor cells, periodic acid–Schiff (PAS)–positive cells, and rich external matrix components (Maniotis et al., 1999). Nutrients and oxygen-carrying red blood cells are delivered to tumors by VM. A number of mechanisms have been proposed for VM formation, which include cancer stem cells (CSCs) and epithelial–mesenchymal transformation (EMT), as well as various signaling pathways that promote VM formation, such as matrix metalloproteinases (MMPs), focal adhesion kinase (FAK), vascular endothelial (VE)–cadherin, phosphatidylinositol 3-kinase (PI3K), and hypoxia inducible factor 1a (HIF-1a) (Lu et al., 2013; Sun et al., 2017). In addition to liver cancer, ovarian cancer, gastric cancer, prostate cancer, and nasopharyngeal cancer, VM has been observed in many other cancers as well. According to a large number of clinical studies, VM is strongly associated with tumor invasiveness and poor prognosis (Seftor et al., 2012; Zhang et al., 2019). However, the significance of VM for the diagnosis, molecular typing, and treatment of BLCA remains unclear.

VM is associated with tumor microenvironments such as CSC, cancer fibroblasts, tumor-associated macrophages, and hypoxia (Luo et al., 2020). To date, the role of molecular VM-related genes in shaping TME heterogeneity remains unclear in BLCA. A correlation was assessed between the mRNA levels of these VM-related molecules and the heterogeneity of the TME, immunophenotype, clinical characteristics of the disease, and response to BLCA treatment. In order to quantify these subtypes in BLCA, the VM_Score was generated using a molecular subtype system mediated by molecular VM-related genes.

2 Material and methods

2.1 Workflow

The workflow of our research is shown in Figure 1.

FIGURE 1

FIGURE 1. Flow chart.

2.2 Data retrieval and preprocessing

2.2.1 External public cohorts

A cohort from the Cancer Genome Atlas (TCGA) was downloaded for RNA sequencing data (FPKM values) and clinical data (https://portal.gdc.cancer.gov/), and FPKM was subsequently converted to TPM value. A total of 412 BLCA samples were included in this study after genomic and clinicopathological data were filtered.

The GSE32894 dataset was downloaded from Gene Expression Omnibus (GEO), and 224 BLCA sample were included in this study. The GSE13507 and GSE48075 datasets were downloaded from GEO, 165 and 73 BLCA sample were included in this study respectively.

2.2.2 Xiangya cohort

Preprocessing and data analysis have been well described in our previous studies (Liu et al., 2021).

2.3 Determination of molecular VM-related gene sets

There is no comprehensive overview of molecular VM-related genes at present. The following search term was used in PubMed: vasculogenic mimicry to retrieve all genes associated with molecular VM. Ultimately, a total of 182 molecular VM-related genes were collected (Supplementary Table S1).

The identification of differentially expressed molecular VM-related genes (molecular VM related DEGs) and functional analysis.

To identify molecular VM-related DEGs in BLCA tissue when compared to normal tissue, the empirical Bayesian method of the limma R package was used (Supplementary Table S2). For molecular VM-related DEGs, the screening criteria were |log (fold change) |>1 and adjusted p-value <0.05. Based on the aforementioned molecular VM-related DEGs, the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) analyses were performed (Yu et al., 2012).

2.4 Development and validation of VM_Score

We used the “rfsrc” function from the “randomForestSRC” package to develop a risk score based on molecular VM-related DEGs (Supplementary Table S3). The median value of the VM_Score was used as the cutoff value, and patients were divided into high- and low-risk groups based on their VM_Scores. A survival curve was plotted using the Kaplan–Meier (K-M) method; the log–rank test was implemented in the “survminer” R package; and a predictive accuracy assessment of the risk score was conducted using the timeROC function provided in the “tROC” R package. Additionally, univariate and multivariate Cox analyses were conducted on the TCGA-BLCA cohort to analyze the independent effects of gender, age, stage, and VM_Score on prognosis. An independent prognostic predictive value nomogram was plotted based on these factors. ROC and calibration curves were used to validate the nomogram.

2.5 Description of immunological features in BLCA TME

The anticancer immune response in the BLCA TME consisted of several steps: cancer antigen release and presentation (Steps 1 and 2), anticancer immune priming and activation (Step 3), immune cell trafficking (Step 4), immune cell infiltration into the TME (Step 5), T cell recognition of cancers (Step 6), and killing of cancer cells (Step 7) (Chen and Mellman, 2013). The seven steps determined the fate of tumor cells. In addition, we employed several independent algorithms to calculate the levels of tumor-infiltrating immune cells (TIICs) based on RNA-seq data, which included TIMER, TIP, quanTIseq, xCell, CIBERSORT-ABS, and MCP-counter (Newman et al., 2015; Becht et al., 2016; Li et al., 2016; Xu et al., 2018; Finotello et al., 2019; Li et al., 2020). The T-cell inflammation score (TIS), which predicts clinical response to immune checkpoint blockade (ICB), can reflect pre-existing anticancer immunity in TMEs (Ayers et al., 2017). In this study, we calculated the enrichment scores for several immunotherapy response-related pathways using ssGSEA. A total of 20 inhibitory immune checkpoints were screened and collected, which included PD-1, TIGIT, LAG3-PD-L1, and CTLA-4. In our previous research, these immunological characteristics have been well described (Hu et al., 2021a; Hu et al., 2021b).

2.6 Prediction of molecular subtypes in BLCA

A highly heterogeneous tumor, BLCA can be treated with precision medicine when individual molecular subtypes can be identified. Previous studies have clarified molecular subtype systems as a result, such as CIT, Lund, TCGA, Baylor, MDA, UNC, and consensus subtypes (Sjödahl et al., 2012; Choi et al., 2014; Damrauer et al., 2014; Rebouissou et al., 2014; Robertson et al., 2017; Mo et al., 2018; Kamoun et al., 2020). As a starting point for subtyping molecular structures, we used the packages ConsensusMIBC and BLCAsubtyping in the R language. A total of 12 BLCA heterogeneous signatures were also collected. These BLCA-specific features and molecular subtypes were correlated further with the VM_Score. As a result of the initial assignment of all samples to basal or luminal subtypes, the ROC curves were plotted to determine the accuracy of the VM_Score in predicting the molecular subtypes.

2.7 Gene set variation analysis

As an unsupervised and non-parametric method, the gene set variation analysis (GSVA) can be used to estimate activity differences of pathways or biological processes in expression of data set samples. The differences of 50 correlation pathways in the VM_Score were investigated using the “GSVA” R package (Hänzelmann et al., 2013).

2.8 Response prediction of BLCA therapy

Because chemotherapy is critical for patients with end-stage BLCA, the pRRophetic software package was used to assess the efficacy of six commonly used chemotherapeutics (cisplatin, bleomycin, camptothecin, docetaxel, paclitaxel, and vinblastine). High-risk and low-risk score groups were compared in terms of 50% inhibitory concentration (IC50) of the six chemotherapeutic agents listed above. Also, GSCALite was used to analyze VM-related DEGs and chemotherapeutic drug sensitivity (http://bioinfo.life.hust.edu.cn/web/GSCALite/) (Liu et al., 2018). It is also important to consider targeted therapy, radiation therapy, and alternative treatment options. This led to the collection of several potential predictors of response to radiation therapy and targeted therapy.

2.9 Statistical analysis

In order to visualize data and perform statistical analysis, the R software (version 4.0.5) was used. In order to analyze the relationship between the continuous variables, we used Pearson’s or Spearman’s correlations. We used the t-test when the continuous variables fit the normal distribution. Otherwise, we use the Mann–Whitney U test. To assess risk score association with prognosis, univariate and multivariate Cox regression analyses were performed. With the help of the ROC curves, the accuracy of the VM_Score was calculated in terms of prognosis and molecular subtype prediction. We used p < 0.05 for all statistical tests and two-sided tests for all analyses.

3 Results

3.1 Landscape and functional analysis of VM-related DEGs in BLCA

Between BLCA and normal tissue, 66 DEGs associated with VM were screened. Among them, 28 molecular VM-related genes were highly expressed in BLCA, and 38 molecular VM-related genes had low expression (Supplementary Figures S1A, B). The GO and KEGG enrichment analysis showed that molecular VM-related DEGs were enriched in several pathways, such as, cell junction assembly, tissue remodeling, regulation of anoikis, regulation of peptidase activity, and focal adhesion. Notably, molecular VM-related DEGs were mainly associated with cell junction, adhesion, and anoikis (Supplementary Figures S1C, D). A novel focal adhesion–related gene signature was significantly correlated with tumor grade, tumor stage, immune scores, and immune infiltrate types, according to Lin et al. (2021). Anoikis-resistant mechanically stressed cancer cells exhibited enhanced motility and evaded immune surveillance by natural killer cells, according to Fanfone et al. (2022). This PPI network shows that these molecular VM-related DEGs are closely related to each other as shown in Supplementary Figure S1E.

3.2 Pan-cancer multiomics analysis of VM DEGs

In addition, we analyzed 66 molecular VM-related DEGs in pan-cancer based on multi-omics signatures. It was observed that these genes were significantly mutated in the following cancers: UCEC, SKCM, COAD, STAD, DLBC, READ, and BLCA. Among them, the SVIL gene had the highest mutation frequency in UCEC, at 76% (Supplementary Figure S2). Heterozygous amplifications and deletions are the main types of copy number variations in molecular VM-related DEGs in pan-cancer, among which the copy numbers of NTN1, SIPR1, LAMTOR5, and ZRANB2 are heterozygous deletions in most tumors. By contrast, EZH2, CDK5, AQP1, and AURKA are heterozygously amplified in most tumors (Supplementary Figure S3). It has been found that the gene copy number variation is one of the most important factors affecting the expression of molecular VM molecules. Most tumor types show a positive correlation between CNV and mRNA expression. There is a significant positive correlation between CNV and PTK2 mRNA expressions (Supplementary Figure S4). As a result, the methylation levels of molecular VM genes are negatively correlated with the levels of mRNA expression in most cancers (Supplementary Figure S5).

3.3 Development and validation of VM_Score in BLCA cohorts

We generated the VM_Score in the TCGA cohort by using the random forest algorithm (Figure 2A). Meanwhile, mRNA expression levels of these four target genes were verified in cell lines and specimens (Supplementary Figure S16). The BLCA cohort from the TCGA was divided into two risk groups based on their level of risk. Notably, we successfully validated that BLCA patients had a worse prognosis in the high-risk group (Figure 2B). More importantly, the VM_Score predicted 1-, 3-, and 5-year OS with accuracies of 0.639, 0.647, and 0.649, respectively (Figure 2C).

FIGURE 2

FIGURE 2. Development and validation of VM_Score in multiple BLCA cohorts. (A) Four best candidates were screened by random forest algorithm to determine the generation of VM_Score. (B) Kaplan–Meier analysis of OS for VM_Score in TCGA-BLCA cohort. (C) ROC curves of VM_Score for predicting OS in TCGA-BLCA cohort. (D,E) Validation of VM_Score in GSE32894. (F,G) Validation of VM_Score in Xiangya cohort.

To assess the extrapolation of risk scores, we validated the risk scores using the external cohort GSE32894, GSE13507, and GSE48075 (Figures 2D, E, and Supplementary Figure S17). In the GSE32894 cohort, patients with a higher VM_Score also had poorer survival outcomes, with prediction accuracies of 0.60, 0.67, and 0.65 at 12, 36, and 60 months, respectively. In addition to this, we constructed our own Xiangya cohort by sequencing and found that the results are consistent with the above trends and have better predictive functions. At 12, 24, and 36 months, patients with high-risk scores had poorer survival outcomes (p = 1.228e−02), with a prediction accuracy of 0.825, 0.675, and 0.905, respectively (Figures 2F, G). All these results suggest that the prediction accuracy of the VM_Score is relatively high and may be a powerful predictive tool for the BLCA operating system.

3.4 Relationship between VM_score and clinicopathological features

Figures 3A, B show that patients with higher grades and stages have higher risk scores, which is consistent with the prognostic correlation of the VM_Score. A subgroup analysis was performed by stage, sex, and age, and as expected, high-risk patients had a worse prognosis in virtually all subgroups (Supplementary Figure S6). The univariate Cox analysis also revealed that the VM_Score, stage, and age predict prognosis significantly (Figure 3C). It was also confirmed by multivariate Cox analysis that the VM_Scores were independent prognostic factors (Figure 3D). As a result, stage no longer serves as a reliable prognostic predictor in the case of patients with BLCA. The results in this study confirm that VM_Score is a reliable indicator of prognosis for BLCA patients. The purpose of this study was to improve the predictive value of the VM_Score for BLCA prognosis, we combined the VM_Score with age, tumor stage, and other factors with prognostic value in the univariate Cox regression analysis to establish a nomogram (Figure 3E). We further validated the accuracy of the nomogram in predicting BLCA prognosis using the ROC and calibration curves. The nomogram’s predictive accuracy for 1-, 3-, and 5-year OS in the TCGA-BLCA cohort was 0.715, 0.721, and 0.749, respectively (Figure 4A). This accuracy of the nomogram and its clinical significance are illustrated by the calibration curve showing that the OS predicted by the nomogram is highly consistent with the actual OS (Figure 4B). Furthermore, both validation sets of the GSE32894 and Xiangya cohort showed high accuracy of the nomograms (Figures 4C–F).

FIGURE 3

FIGURE 3. Construction of a nomogram in the TCGA-BLCA cohort. (A, B) Relationship between VM_Score and tumor grade and stage in TCGA-BLCA cohort. (C, D) Results of univariate and multivariate Cox analyses. (E) Nomogram developed based on stage, age, and VM_Score to predict overall survival at 1, 3, and 6 years.

FIGURE 4

FIGURE 4. Validation of multiple cohorts of VM_Score. (A) ROC curves of the nomogram. (B) Calibration curves of the nomogram measured using the Hosmer–Lemeshow test. (C, D) Validity of VM_Score in GSE32894. (E, F) Validity of VM_Score in Xiangya cohort.

3.5 Prediction of VM_score for immune phenotypes and clinical response of ICB in BLCA

The advent of immunotherapy has brought light to patients with advanced bladder cancer, but the efficacy of immunotherapy or neoadjuvant therapy largely depends on the TME. Most cancer immune cycles were more active in patients with high a VM_Score; for example, the release of cancer cell antigens (step 1), trafficking of immune cells to the tumor (step 4) (recruitment of CD8 T cells, CD4 T cells, macrophages, Th1 cells, NK cells, and DCs), and killing of cancer cells (Step 7) (Figure 5A). These immune cycles were also upregulated, resulting in increased levels of infiltration of the corresponding TIICs (which included CD8 T cells, CD4 T cells, NK cells, Th1 cells, macrophages, and DCs) in the BLCA TME. It was performed in six independent algorithms cross-validation (Figure 5B). The results suggest that those with a high VM_Score phenotype may be more sensitive to ICB because they show an inflammatory phenotype. Our next step was to examine the correlation between the VM_Score and predictors of ICB effectiveness. First, by analyzing the relationship between the VM_Score and TIS, we found that they were positively correlated (Figure 5C). Low VM_Score groups had lower enrichment scores for immune checkpoints (such as, PD-L1, TIGIT, and CTLA-4) and immunotherapy response genes when compared with high VM_Score groups (Figures 5D, E). In the stratified analysis of age, sex, and stage, there was a consistent trend (Supplementary Figures S7–S12).

FIGURE 5

FIGURE 5. VM_Score correlated with the tumor immune microenvironment characteristics. (A) Differences in cancer immune cycling activity between high- and low-risk groups. (B) Relationship between VM_Score and several immune cells (CD8⁺T cells, NK cells, macrophages, Th1 cells, and DCs) in six independent algorithms. (C) Relationship between VM_Score and T cell inflamed score (TIS). (D) Correlation between VM_Score and enrichment of ICB response–related pathways. (E) Relationship between VM_Score and immune checkpoints (*p < 0.05; **p < 0.01; ***p < 0.001).

All results suggested that the VM_Score was correlated with the TME and may be a potential predictor of ICB efficacy in BLCA.

3.6 Effect of VM_score for molecular subtypes and therapeutic response

In recent years, the emergence of molecular subtypes has promoted the development of precise diagnosis and treatment, and the prognosis and treatment response of BLCA can be better predicted through the refinement of molecular subtypes. After consulting the literature, it was found that there are currently seven reported molecular types of BLCA (Kamoun et al., 2020). However, their clinical application is hindered by the designation of molecular subtypes and the complex detection methods required. Initially, we compared the enrichment of several signaling pathways between high-risk and low-risk groups. In addition, the results indicated that there were significant differences between the high and low VM_Score groups with regard to biological functions (Figure 6A). A high-risk group showed a lot of signals related to complement, inflammation, and EMT, while a low-risk group showed a lot of signals related to DNA repair and xenobiotic metabolism. These results suggest that VM-Score may affect the progression of BLCA by regulating the hallmark signaling pathway. In Figures 6B, C, correlations between the VM_Score and classical molecular subtype classification are shown. According to seven molecular classifications, the high-risk group represented the basal subtype of BLCA: TCGA subtype (Robertson et al., 2017), Cartes d’Identité des Tumeurs-Curie (CIT) subtype (Rebouissou et al., 2014), MD Anderson Cancer Center (MDA) subtype (Choi et al., 2014), Lund subtype (Marzouka et al., 2018), University of North Carolina (UNC) subtype (Damrauer et al., 2014), Baylor subtype (Mo et al., 2018), and consensus subtype (Kamoun et al., 2020), which was characterized by basal differentiation, EMT differentiation, immune differentiation, interferon response, and myofibroblasts; as opposed to low VM_Scores, which are characterized by luminal differentiation, urothelial differentiation, and the Ta pathway. The ROC curves showed that the VM_Score had a high prediction accuracy for molecular subtypes, mostly exceeding 0.75 (Figure 6D).

FIGURE 6

FIGURE 6. VM_Score effectively predicted molecular subtypes and therapy response. (A) Differences in biological function between VM_Score groups. (B,C) Relationship between VM_Score and seven classical molecular subtypes. (D) Predictive accuracy of VM_Score for molecular subtypes in multiple different algorithms. (E) Relationship between VM_Score and enrichment scores of multiple therapeutic signatures. (F) Difference on the effects of six chemotherapy drugs.

A number of different treatment regimens were tested for their ability to predict their responses to VM_Score. It is possible that patients with high VM_Scores are more susceptible to EGFR-targeted therapy and radiotherapy. A significant VM_Score enrichment was found in the low VM_Score group for several immunosuppressive oncogenic pathways, such as the Wnt–β-catenin network, IDH1, PPARG network, KDM6B, and VEGFA (Figure 6E). The targeting of these oncogenic pathways may provide promising therapeutic opportunities for BLCA patients with low VM_Scores. Furthermore, chemotherapy drugs were more sensitive in patients with a high VM_Score, such as cisplatin, paclitaxel, bleomycin, camptothecin, docetaxel, and vinblastine (Figure 6F).

Overall, the VM_Score may be a cost-effective and simpler alternative to classical molecular subtypes. In addition, the VM_Score can be used to predict how well BLCA patients will respond to several treatments.

3.7 Validating roles of VM_score in Xiangya cohort and GSE32894

A further validation of the capability of the VM_Score to predict the immune status, molecular subtype, and treatment response has been performed in our cohort (Xiangya cohort). Multiple anticancer immune cycles and infiltration levels of CD8 T cells, Th1 cells, DC, NK cells, and macrophages were positively correlated with the VM_Score in six independent algorithms (Figures 7A, E). Furthermore, the VM_Score was significantly correlated with the enrichment scores of immune response–related signatures, TIS, and ICB (Figures 7B–D). Data from this study have suggested that the VM_Score is an effective tool for stratifying BLCA immunophenotypes. Additionally, the VM_Score can accurately predict classical molecular subtypes with AUCs usually around 0.90, except in the TCGA data set (AUC = 0.77) (Figures 7F, G). There is a promising correlation between VM_Score and clinical response to radiotherapy, as well as several targeted therapies. It appears that patients in the high-risk group may benefit more from radiotherapy and EGFR-targeted therapies. A targeted therapy, such as blocking the Wnt–β-catenin network or the FGFR3 network, may be better suited for patients at low risk (Figure 7H). The abovementioned results were all validated in GSE32894 (Supplementary Figure S13). In addition, we looked at the correlation between 66 VM-related DEGs and sensitivity to multiple drugs, and we found that CYLD, MYC, and ZEB1 were negatively related to most drugs, while CDH1, CLDN4, and SERPINB5 were positively associated with major drug sensitivity (Figure 8, Supplementary Figures S13, S15).

FIGURE 7

FIGURE 7. Validation of VM_Score in Xiangya cohort. (A) Relationship between VM_Score and activities of cancer immunity cycles. (B) Relationship between VM_Score and immunotherapy predicted pathways. (C) Correlations between VM_Score and several immune checkpoints. (D) Relationship between VM_Score and T cell inflammation score (TIS). (E) Relationship between VM_Score and infiltration levels of five tumor-infiltrating immune cells. (F) VM_Score accurately stratified the molecular subtypes in seven different algorithms. (G) Accuracy of VM_Score in predicting molecular subtypes in seven different algorithms. (H) Relationship between VM_Score and enrichment scores of several therapeutic signatures.

FIGURE 8

FIGURE 8. Correlation analysis between VM-related DEGs and drug sensitivity (top30 drugs in GDSC and CTRP). (A, B) Bubble chart shows the correlation analysis between VM-related genes and drug susceptibility. Red indicates positive correlation and blue indicates negative correlation. The darker the color, the higher is the correlation index. Bubble size indicates FDR.

4 Discussion

There are increasing numbers of studies showing that VM plays a critical role in the regulation of tumors; however, the role of VM in BLCA is unclear. In this study, we developed a new systemic indicator of molecular VM–related genes mediated by molecular subtypes based on the Xiangya cohort and additional external BLCA cohorts using a random forest algorithm. A correlation was then done between the VM_Score and classical molecular subtypes, clinical outcomes, immunophenotypes, and treatment options for BLCA. With the VM_Score, it is possible to predict classical molecular subtypes, immunophenotypes, prognosis, and therapeutic potential of BLCA with high accuracy. VM_Scores of high levels indicate a more anticancer immune response but a worse prognosis due to a more basal and inflammatory phenotype. the VM_Score was also associated with low sensitivity to antiangiogenic and targeted therapies that target the FGFR3, β-catenin, and PPAR-γ pathways but are associated with high sensitivity to cancer immunotherapy, neoadjuvant chemotherapy, and radiotherapy. A number of aspects of BLCA biology were reflected in the VM_Score, providing new insights into precision medicine. Additionally, the VM_Score may be used as an indicator of pan-cancer immunotherapy response and prognosis.

We collected molecular VM-related genes through keyword search, performed differential analysis in BLCA and the adjacent tissues to screen out 66 differentially expressed genes, and constructed the VM_Score by the random forest method. Then, we validated it in multiple data sets such as our own Xiangya cohort, and the results showed that the VM_Score can predict the clinical outcome, molecular subtype, and treatment response of BLCA and was related to the inflammatory TME. VM affects tumor progression and survival prognosis. Studies have found that VM was significantly associated with the survival and prognosis of various tumors (Cao et al., 2013). In this study, the KEGG and GO enrichment analyses of VM-related differential genes showed that they were mainly enriched in cell junction assembly, tissue remodeling, regulation of anoikis, regulation of peptidase activity, and focal adhesion. These pathways were closely related to tumor invasion and migration. In addition, through analyzing the differences among 50 signaling pathways between the VM_Score groups, it was found that EMT was most abundant in the high VM_Score group and DNA repair was more abundant in the low-risk group. EMT and DNA damage are involved in key steps in tumorigenesis and development (Aiello and Kang, 2019; Srinivas et al., 2019). Genomic instability promotes tumor progression (Andor et al., 2017). Through mutation analysis of VM-related DEGs, it was found that most genes have high mutation frequencies in various tumors such as UCEC, SKCM, and BLCA. The gene copy number variation was mainly heterozygous amplification and deletion. The study concluded that VM may play a significant role in the occurrence and development of tumors.

The VM_Score predicts molecular subtype and treatment response in BLCA. Researchers found that the higher the VM_Score, the higher is the tumor stage, and worse is the survival prognosis for BLCA patients. In addition, the opposite was the case with low VM_Scores. The correlation between VM-related DEGs and the sensitivity of various drugs was analyzed in CTRP and GDSC databases, and many genes were found to be closely related to most drugs. There was a greater sensitivity to chemotherapy drugs such as cisplatin, paclitaxel, bleomycin, docetaxel, camptothecin, and vinblastine, and radiotherapy and EGFR-targeted therapy in the high VM_Score group, whereas drugs that target oncogenic pathways such as Wnt–β-catenin network, PPARG network, IDH1, KDM6B, and VEGFA may be more effective for the low VM_Score group. Molecular subtype has greatly facilitated precision medicine. In recent years, a variety of molecular typing has been reported, but they have not been used in clinical practice due to factors such as complex detection methods. According to this study, the low VM_Score group mainly represented luminal types, and the high VM_Score group represented the basal type, and the prediction was more accurate. Therefore, the VM_Score is a more reliable prediction tool, which is of great help in accurate diagnosis and treatment by representing different molecular subtypes.

VM correlated with the tumor microenvironment and guided immunotherapy. Studies have shown that vasculogenic mimicry structures support the recruitment of monocytes (Tan et al., 2022). In addition, the KEGG and GO enrichment analyses of molecular VM-related DEGs found that they were significantly enriched in focal adhesion and regulation of anoikis. There was a significant association between a novel focal adhesion-related gene signature and tumor grade, tumor stage, immune scores, and immune infiltrate types, determined by Lin et al. (2021). It was discovered by Fanfone et al. (2022) that anoikis-resistant mechanically stressed cancer cells were more motile and had a better immune response against natural killer cells. The enrichment analysis of the high and low VM_Score groups showed that the complement and inflammatory responses had more abundant signals in the high VM_Score groups. In addition, according to this study, the VM_Score had a positive correlation with anti-cancer immune cycle and TIIC, suggesting that the VM_Score could be related to an inflammatory immune microenvironment, and patients with high VM_Scores had higher anti-cancer immunity. The vascular endothelial growth factor receptor-1 (VEGFR-1) is a membrane receptor which plays a crucial role in melanoma vasculogenic mimicry. VEGFR-1 mediated signaling as an effective target for reducing pre-tumor macrophage tumor invasion and improving ICI immunotherapy efficacy (Lacal et al., 2020). Myeloid-derived suppressor cells' (MDSCs) migration and differentiation are enhanced by VEGF signaling. VEGF-a knockdown in tumor cells resulted in decreased infiltration of MDSC and increased infiltration of CD8 T cells (Horikawa et al., 2017). Antiangiogenic therapy by vaccination or adoptive cell transfer a few days before immunotherapy increases the accumulation of T cells within the tumor, and thus improves the anti-cancer efficacy, when compared with vaccination alone (Fukumura et al., 2018). It is well known that ICB has achieved good results in tumor immunotherapy. In addition, the VM_Score was positively correlated with TIS and immune checkpoint scores (such as CD274, CTLA4, PDCD1, TIGIT, and LAG3), as well as the immunotherapy response–related gene signature enrichment scores, suggesting that the TMEs in the low VM_Score group had fewer immunotherapy target points and the ICB treatment effect is not ideal. By analyzing the methylation and expression levels of molecular VM-related DEGs, we found that methylation modification greatly affected the expression of VM-related DEGs. We can reverse the TME by changing the modification pattern of genes in the low VM_Score group, turning “cold tumors” into “hot tumors.” The VM_Score can predict the therapeutic response to ICB and closely relate to the TME, thus providing a new target for ICB combined treatments.

We recognize that this study has certain limitations. Firstly, this study employed bioinformatics analysis, and the accuracy of the prediction is slightly lower. Although these results have been repeatedly validated in multiple public cohorts and our own Xiangya cohort, we still have to conduct in vivo and in vitro studies on the relevant mechanisms of VM. Secondly, there is a need for prospective clinical trials to validate the clinical value of the VM_Score. Thirdly, the optimal VM_Score cutoff value was not determined.

5 Conclusion

The VM_Score developed and validated by our group can successfully predict clinical characteristics, TME characteristics, and survival outcomes of BLCA. The VM_Score was a reliable and effective biomarker and tool for predicting patients’ clinical responses to several therapies. It is possible that patients with a high VM_Score will respond more favorably to immunotherapy, chemotherapy, radiotherapy, and EGFR-targeted therapies. However, patients with low-risk scores may benefit from a number of targeted therapies, such as blocking PPARG, Wnt–β-catenin, and FGFR3.

Data availability statement

The data sets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author contributions

Conception and design: CZ, JX, JC, XZ, ZL, and PL. Provision of study materials or patients: ZX, YH, and DD. Collection and assembly of data: TY and ZL. Data analysis and interpretation: CZ and PL. Manuscript writing: all authors. Final approval of manuscript: all authors.

Funding

A grant from the National Natural Science Foundation of China supported this research [82070785, 81873626], and the Science and Technology fund project of Health and Family Planning Commission of Guizhou province (No. gzwjkj2022-100). Basic Research Project of Science and Technology Cooperation in Qiandongnan Prefecture (No. 2023-16).

Acknowledgments

The authors sincerely thank all participants in the study.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, editors, and reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar.2023.1163115/full#supplementary-material

References

Aiello, N. M., and Kang, Y. (2019). Context-dependent emt programs in cancer metastasis. J. Exp. Med. 216 (5), 1016–1026. Epub 2019/04/13. doi:10.1084/jem.20181827

PubMed Abstract | CrossRef Full Text | Google Scholar

Andor, N., Maley, C. C., and Ji, H. P. (2017). Genomic instability in cancer: Teetering on the limit of tolerance. Cancer Res. 77 (9), 2179–2185. Epub 2017/04/23. doi:10.1158/0008-5472.Can-16-1553

PubMed Abstract | CrossRef Full Text | Google Scholar

Ayers, M., Lunceford, J., Nebozhyn, M., Murphy, E., Loboda, A., Kaufman, D. R., et al. (2017). Ifn-Γ-related mrna profile predicts clinical response to Pd-1 blockade. J. Clin. Invest. 127 (8), 2930–2940. Epub 2017/06/27. doi:10.1172/jci91190

PubMed Abstract | CrossRef Full Text | Google Scholar

Becht, E., Giraldo, N. A., Lacroix, L., Buttard, B., Elarouci, N., Petitprez, F., et al. (2016). Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol. 17 (1), 218. Epub 2016/10/22. doi:10.1186/s13059-016-1070-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Cao, Z., Bao, M., Miele, L., Sarkar, F. H., Wang, Z., and Zhou, Q. (2013). Tumour vasculogenic mimicry is associated with poor prognosis of human cancer patients: A systemic review and meta-analysis. Eur. J. Cancer 49 (18), 3914–3923. Epub 2013/09/03. doi:10.1016/j.ejca.2013.07.148

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, D. S., and Mellman, I. (2013). Oncology meets immunology: The cancer-immunity cycle. Immunity 39 (1), 1–10. Epub 2013/07/31. doi:10.1016/j.immuni.2013.07.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Choi, W., Porten, S., Kim, S., Willis, D., Plimack, E. R., Hoffman-Censits, J., et al. (2014). Identification of distinct basal and luminal subtypes of muscle-invasive bladder cancer with different sensitivities to frontline chemotherapy. Cancer Cell 25 (2), 152–165. Epub 2014/02/15. doi:10.1016/j.ccr.2014.01.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Dagogo-Jack, I., and Shaw, A. T. (2018). Tumour heterogeneity and resistance to cancer therapies. Nat. Rev. Clin. Oncol. 15 (2), 81–94. Epub 2017/11/09. doi:10.1038/nrclinonc.2017.166

PubMed Abstract | CrossRef Full Text | Google Scholar

Damrauer, J. S., Hoadley, K. A., Chism, D. D., Fan, C., Tiganelli, C. J., Wobker, S. E., et al. (2014). Intrinsic subtypes of high-grade bladder cancer reflect the hallmarks of breast cancer biology. Proc. Natl. Acad. Sci. U. S. A. 111 (8), 3110–3115. Epub 2014/02/13. doi:10.1073/pnas.1318376111

PubMed Abstract | CrossRef Full Text | Google Scholar

Fanfone, D., Wu, Z., Mammi, J., Berthenet, K., Neves, D., Weber, K., et al. (2022). Confined migration promotes cancer metastasis through resistance to anoikis and increased invasiveness, 11. Epub 2022/03/09. doi:10.7554/eLife.73150Elife

CrossRef Full Text

Finotello, F., Mayer, C., Plattner, C., Laschober, G., Rieder, D., Hackl, H., et al. (2019). Molecular and pharmacological modulators of the tumor immune contexture revealed by deconvolution of rna-seq data. Genome Med. 11 (1), 34. Epub 2019/05/28. doi:10.1186/s13073-019-0638-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Fukumura, D., Kloepper, J., Amoozgar, Z., Duda, D. G., and Jain, R. K. (2018). Enhancing cancer immunotherapy using antiangiogenics: Opportunities and challenges. Nat. Rev. Clin. Oncol. 15 (5), 325–340. Epub 2018/03/07. doi:10.1038/nrclinonc.2018.29

PubMed Abstract | CrossRef Full Text | Google Scholar

Hänzelmann, S., Castelo, R., and Guinney, J. (2013). Gsva: Gene set variation analysis for microarray and rna-seq data. BMC Bioinforma. 14, 7. Epub 2013/01/18. doi:10.1186/1471-2105-14-7

CrossRef Full Text | Google Scholar

Horikawa, N., Abiko, K., Matsumura, N., Hamanishi, J., Baba, T., Yamaguchi, K., et al. (2017). Expression of vascular endothelial growth factor in ovarian cancer inhibits tumor immunity through the accumulation of myeloid-derived suppressor cells. Clin. cancer Res. official J. Am. Assoc. Cancer Res. 23 (2), 587–599. Epub 2016/07/13. doi:10.1158/1078-0432.Ccr-16-0387

CrossRef Full Text | Google Scholar

Hu, J., Othmane, B., Yu, A., Li, H., Cai, Z., Chen, X., et al. (2021). 5mc regulator-mediated molecular subtypes depict the hallmarks of the tumor microenvironment and guide precision medicine in bladder cancer. BMC Med. 19 (1), 289. Epub 2021/11/28. doi:10.1186/s12916-021-02163-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, J., Yu, A., Othmane, B., Qiu, D., Li, H., Li, C., et al. (2021). Siglec15 shapes a non-inflamed tumor microenvironment and predicts the molecular subtype in bladder cancer. Theranostics 11 (7), 3089–3108. Epub 2021/02/05. doi:10.7150/thno.53649

PubMed Abstract | CrossRef Full Text | Google Scholar

Kamoun, A., de Reyniès, A., Allory, Y., Sjödahl, G., Robertson, A. G., Seiler, R., et al. (2020). A consensus molecular classification of muscle-invasive bladder cancer. Eur. Urol. 77 (4), 420–433. Epub 2019/09/30. doi:10.1016/j.eururo.2019.09.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Lacal, P. M., Atzori, M. G., Ruffini, F., Scimeca, M., Bonanno, E., Cicconi, R., et al. (2020). Targeting the vascular endothelial growth factor receptor-1 by the monoclonal antibody D16f7 to increase the activity of immune checkpoint inhibitors against cutaneous melanoma. Pharmacol. Res. 159, 104957. Epub 2020/06/03. doi:10.1016/j.phrs.2020.104957

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, B., Severson, E., Pignon, J. C., Zhao, H., Li, T., Novak, J., et al. (2016). Comprehensive analyses of tumor immunity: Implications for cancer immunotherapy. Genome Biol. 17 (1), 174. Epub 2016/08/24. doi:10.1186/s13059-016-1028-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, T., Fu, J., Zeng, Z., Cohen, D., Li, J., Chen, Q., et al. (2020). Timer2.0 for analysis of tumor-infiltrating immune cells. Nucleic Acids Res. 48 (W1), W509–W514. Epub 2020/05/23. doi:10.1093/nar/gkaa407

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, Z., Miao, D., Xu, Q., Wang, X., and Yu, F. (2021). A novel focal adhesion related gene signature for prognostic prediction in hepatocellular carcinoma. Aging (Albany NY) 13 (7), 10724–10748. Epub 2021/04/15. doi:10.18632/aging.202871

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, C. J., Hu, F. F., Xia, M. X., Han, L., Zhang, Q., and Guo, A. Y. (2018). Gscalite: A web server for gene set cancer analysis. Bioinformatics 34 (21), 3771–3772. Epub 2018/05/24. doi:10.1093/bioinformatics/bty411

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Z., Qi, T., Li, X., Yao, Y., Othmane, B., Chen, J., et al. (2021). A novel tgf-Β risk score predicts the clinical outcomes and tumour microenvironment phenotypes in bladder cancer. Front. Immunol. 12, 791924. Epub 2022/01/04. doi:10.3389/fimmu.2021.791924

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, X. S., Sun, W., Ge, C. Y., Zhang, W. Z., and Fan, Y. Z. (2013). Contribution of the pi3k/mmps/ln-5γ2 and epha2/fak/paxillin signaling pathways to tumor growth and vasculogenic mimicry of gallbladder carcinomas. Int. J. Oncol. 42 (6), 2103–2115. Epub 2013/04/17. doi:10.3892/ijo.2013.1897

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, Q., Wang, J., Zhao, W., Peng, Z., Liu, X., Li, B., et al. (2020). Vasculogenic mimicry in carcinogenesis and clinical applications. J. Hematol. Oncol. 13 (1), 19. Epub 2020/03/15. doi:10.1186/s13045-020-00858-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Maniotis, A. J., Folberg, R., Hess, A., Seftor, E. A., Gardner, L. M., Pe'er, J., et al. (1999). Vascular Channel formation by human melanoma cells in vivo and in vitro: Vasculogenic mimicry. Am. J. Pathol. 155 (3), 739–752. Epub 1999/09/17. doi:10.1016/s0002-9440(10)65173-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Marzouka, N. A., Eriksson, P., Rovira, C., Liedberg, F., Sjödahl, G., and Höglund, M. (2018). A validation and extended description of the Lund taxonomy for urothelial carcinoma using the tcga cohort. Sci. Rep. 8 (1), 3737. Epub 2018/03/01. doi:10.1038/s41598-018-22126-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Mo, Q., Nikolos, F., Chen, F., Tramel, Z., Lee, Y. C., Hayashi, K., et al. (2018). Prognostic power of a tumor differentiation gene signature for bladder urothelial carcinomas. J. Natl. Cancer Inst. 110 (5), 448–459. Epub 2018/01/18. doi:10.1093/jnci/djx243

PubMed Abstract | CrossRef Full Text | Google Scholar

Newman, A. M., Liu, C. L., Green, M. R., Gentles, A. J., Feng, W., Xu, Y., et al. (2015). Robust enumeration of cell subsets from tissue expression profiles. Nat. Methods 12 (5), 453–457. Epub 2015/03/31. doi:10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

Rebouissou, S., Bernard-Pierrot, I., de Reyniès, A., Lepage, M. L., Krucker, C., Chapeaublanc, E., et al. (2014). Egfr as a potential therapeutic target for a subset of muscle-invasive bladder cancers presenting a basal-like phenotype. Sci. Transl. Med. 6 (244), 244ra91. 244ra91. Epub 2014/07/11. doi:10.1126/scitranslmed.3008970

PubMed Abstract | CrossRef Full Text | Google Scholar

Robertson, A. G., Kim, J., Al-Ahmadie, H., Bellmunt, J., Guo, G., Cherniack, A. D., et al. (2017). Comprehensive molecular characterization of muscle-invasive bladder cancer. Cell 171 (3), 540–556.e25. e25. Epub 2017/10/11. doi:10.1016/j.cell.2017.09.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Seftor, R. E., Hess, A. R., Seftor, E. A., Kirschmann, D. A., Hardy, K. M., Margaryan, N. V., et al. (2012). Tumor cell vasculogenic mimicry: From controversy to therapeutic promise. Am. J. Pathol. 181 (4), 1115–1125. Epub 2012/09/05. doi:10.1016/j.ajpath.2012.07.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Siegel, R. L., Miller, K. D., Fuchs, H. E., and Jemal, A. (2022). Cancer statistics, 2022. CA Cancer J. Clin. 72 (1), 7–33. Epub 2022/01/13. doi:10.3322/caac.21708

PubMed Abstract | CrossRef Full Text | Google Scholar

Sjödahl, G., Lauss, M., Lövgren, K., Chebil, G., Gudjonsson, S., Veerla, S., et al. (2012). A molecular taxonomy for urothelial carcinoma. Clin. Cancer Res. 18 (12), 3377–3386. Epub 2012/05/04. doi:10.1158/1078-0432.Ccr-12-0077-t

PubMed Abstract | CrossRef Full Text | Google Scholar

Srinivas, U. S., Tan, B. W. Q., Vellayappan, B. A., and Jeyasekharan, A. D. (2019). Ros and the DNA damage response in cancer. Redox Biol. 25, 101084. Epub 2019/01/08. doi:10.1016/j.redox.2018.101084

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, B., Zhang, D., Zhao, N., and Zhao, X. (2017). Epithelial-to-Endothelial transition and cancer stem cells: Two cornerstones of vasculogenic mimicry in malignant tumors. Oncotarget 8 (18), 30502–30510. Epub 2016/04/02. doi:10.18632/oncotarget.8461

PubMed Abstract | CrossRef Full Text | Google Scholar

Tan, L. Y., Cockshell, M. P., Moore, E., Myo Min, K. K., Ortiz, M., Johan, M. Z., et al. (2022). Vasculogenic mimicry structures in melanoma support the recruitment of monocytes. Oncoimmunology 11 (1), 2043673. Epub 2022/03/18. doi:10.1080/2162402x.2022.2043673

PubMed Abstract | CrossRef Full Text | Google Scholar

Witjes, J. A., Bruins, H. M., Cathomas, R., Compérat, E. M., Cowan, N. C., Gakis, G., et al. (2021). European association of urology guidelines on muscle-invasive and metastatic bladder cancer: Summary of the 2020 guidelines. Eur. Urol. 79 (1), 82–104. Epub 2020/05/04. doi:10.1016/j.eururo.2020.03.055

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, L., Deng, C., Pang, B., Zhang, X., Liu, W., Liao, G., et al. (2018). Tip: A web server for resolving tumor immunophenotype profiling. Cancer Res. 78 (23), 6575–6580. Epub 2018/08/30. doi:10.1158/0008-5472.Can-18-0689

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, G., Wang, L. G., Han, Y., and He, Q. Y. (2012). Clusterprofiler: An R package for comparing biological themes among gene clusters. Omics 16 (5), 284–287. Epub 2012/03/30. doi:10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, X., Zhang, J., Zhou, H., Fan, G., and Li, Q. (2019). Molecular mechanisms and anticancer therapeutic strategies in vasculogenic mimicry. J. Cancer 10 (25), 6327–6340. Epub 2019/11/28. doi:10.7150/jca.34171

PubMed Abstract | CrossRef Full Text | Google Scholar