We from the cancer genome atlas (TCGA) database (https://portal.gdc.cancer.gov/) to download the bladder cancer patients (including 414 samples of bladder cancer and 19 adjacent non tumor samples) expression of the spectral data, clinical information, somatic mutation. From Gene Expression Omnibus (GEO) database (https://cancergenome.nih.gov/) download GSE13507 (n = 165) and GSE32894 (n = 224) as an independent verification of the queue. The list of genes involved in glutamine metabolism was obtained from the GenesCards database.

The R package “limma” was used to perform the difference analysis with the cutoff value set to p value<0.05. Gene Ontology (GO) and Kyto Encyclopedia of Genes and Genomes (KEGGE) analyses were performed using R package clusterProfiler. The STRING database was used to analyze protein-protein interactions (PPI) (15), and Cytoscape was used to visualize PPI networks (16). We identified key modules by using the ‘MCODE’ plug-in in Cytoscape, and used seven algorithms commonly used in ‘cytohubba’ plug-in (Closeness, Degree, EPC, Radiality, Stress, MCC, MNC) to identify Hub genes. TRRUST databasewas used to predict the transcription factor (TF) of Hub gene (17).

Bladder cancer samples were clustered by R package “ConsensusClusterPlus” to identify molecular subtypes related to glutamine metabolism. The R package “Survival” was used to perform Kaplan-Meier (KM) survival curves to compare outcomes between the two clusters.

Univariate Cox regression analysis was used to screen out the genes associated with overall survival (OS) of bladder cancer, and multivariate Cox regression analysis was used to establish GMII. The Glutamine Metabolism Immunity Index (GMII) was calculated for each patient according to the following formula: Glutamine Metabolism Immunity Index (GMII) = Coef(Gene1) × Expr(Gene1) + Coef(Gene2) × Expr(Gene2) +…… Coef(Genen) × Expr(Genen). Expr(Genen) represents the expression level of a specific gene, and Coef(Genen) represents the coefficient in multivariate Cox analysis. The prognostic value of the features was verified by KM analysis and Receiver Operation Characteristic (ROC) curve, and the prognostic characteristics were verified by GSE13507. Univariate and multivariate Cox analyses were performed to determine whether the characteristics were independent risk factors. According to the clinical characteristic parameters, the correlation and stratification analysis between GMII and clinical characteristics were performed, and the nomogram was constructed to compare the consistency between predicted and actual survival rates by 1-year, 3-year and 5-year calibration maps.

The GSCALite platform was used to analyze the eight genes of GMII (18). Fourteen pairs of normal and tumor tissue samples were selected for differential expression analysis. The prognosis of eight GMII genes was analyzed in 33 pan-cancer cancers. Gene set variation analysis (GSVA) was used to score the expression of these eight genes in each pan-cancer sample, and the GSVA score was obtained. The GSVA score was used to analyze the expression, prognosis and immune cell infiltration of these genes in pan-cancer. We also performed mutational analyses of the eight genes, including single nucleotide variation, copy number variation and methylation. Through CTRP database, we analyzed the expression levels and drug sensitivity of eight genes of GMII.

Build GMII gene immunohistochemical from Human Protein Atlas (HPA) database (https://www.proteinatlas.org/) (19).

We performed GSEA analysis by R package “clusterProfiler” to evaluate the main enrichment pathways in high-GMII groups to explore the underlying biological mechanisms. Filter for | NES |> 1, nominal p value < 0.05. Sample replacement was tested 1000 times, and clustering analysis of enriched gene sets was performed using the R package “enrichplot”. Reference genomes include Hallmark, c5go, and c2kegg.

We calculated the tumor mutation burden (TMB) for each patient from somatic mutation data and compared TMB between high-GMII and low-GMII groups. The waterfall map is depicted through the R package “Maftools”, showing the mutation landscape of the high-GMII and low-GMII groups. We also performed mutually exclusive and collaborative analyses of genes with the highest frequency of mutations in the high-GMII and low-GMII groups. Finally, somatic mutations of GMII genes were identified by cBioPortal database.

CIBERSORT method used to quantify infiltrating immune cell ratio (https://cibersort.stanford.edu/) (20). The proportion of 22 immune cells (B-naive cells, B-cell memory, plasma cells, T-cell CD8, T-cell CD4 naive, T-cell follicular helper cells, T-cell CD4 memory resting, T-cell CD4 memory activation, regulatory T cells (Tregs), γδ cells, monocytes, activation) was calculated by CIBERSORT method NK cells, resting NK cells, macrophage M0, macrophage M1, macrophage M2, resting dendritic cells, activated dendritic cells, resting mast cells, activated mast cells, eosinophils, and neutrophils). Samples with P<0.05 indicated that the proportion of immune cells calculated by CIBERSORT was correct. The tumor purity, stromal score, immune score and ESTIMATE scorewere calculated for each tumor sample by R package “ESTIMATE” (21). A single sample Gene set enrichment analysis (ssGSEA) algorithm was used to assess the immune infiltration between the two groups based on 28 immune cell types. XCELL, QUANTISEQ, MCPCOUNTER, EPIC and CIBERSORT-ABS software were also used to quantify the relative proportion of immune cell infiltration.

From tumor susceptibility multiple omics (GDSC) database (https://www.cancerrxgene.org/) (22), download the cancer gene expression data of different drugs, through the calculation of R packages “pRRophetic” IC50 to assess patient response to common chemotherapy drugs.

Tumor Immune Dysfunction and Exclusion (TIDE) algorithm was used to infer the clinical response of patients to immunotherapy (23). High TIDE scores were associated with poorer immunotherapy. In addition, we extracted the IMvigor210 dataset, a group of clinical information on the treatment of urothelial carcinoma by anti-PD-L1 monoclonal antibody (atezolizumab) (24). The relationship between bladder cancer anti-PD-1 and anti-CTLA4 by Immunophenoscores(IPS) scores and GMII. The IPS score is a predictive score for a patient’s response to anti-PD-1 and anti-ctLA-4 treatments (25). was downloaded from TCIA database (https://tcia.at/home). These results were used to evaluate the predictive value of GMII for immune checkpoint therapy.

We used MOE software to screen FDA-approved drugs that bind to target proteins and perform molecular docking simulations. Protein structures of core targets were collected from the PDB database and FDA-approved drugs were collected from the zinc15 database and converted to 3D structures in MOE by energy minimization. We optimized the protonation state of the protein and the direction of hydrogen at the PH of LigX 7 and the temperature of 300K. Finally, we studied the binding mode of PPARG, SLC7A9 and GALK1 with small molecule drugs by rigid docking simulation.

Survival curves were drawn using the Kaplan-Meier method. Wilcoxon test and Kruskal-Wallis test were used for comparison between two groups and more than two groups, respectively. Correlation was assessed by Spearman correlation analysis. A P value of 0.05 or less was considered statistically significant. All statistical analyses were performed by R (version 4.1.1).

The flow chart of this study is shown in Figure S1 . We obtained 501 genes related to glutamine metabolism from the GeneCards database, and the screening criteria were that the correlation score was greater than 8 and they were protein-coding genes. Differential analysis of BC and normal bladder tissues yielded 301 genes. GO and KEGG analysis were performed to investigate the functions of these genes related to glutamine metabolism. The GO results showed that the genes related to glutamine metabolism were mainly enriched in the biological functions related to energy and metabolism. KEGG results showed that genes related to glutamine metabolism were enriched in signaling pathways such as carbon metabolism, AMPK signaling pathway and amino acid production ( Figures 1A, B ). PPI networks were analyzed by STRING database and visualized by Cytoscape to obtain three main network diagrams ( Figures 1C–E ). The common Hub genes were obtained by seven algorithms in Cytoscape software and the transcription factors of Hub genes were predicted in TRRUST database ( Figures 1F, G ).

Clustering analysis of bladder cancer patients with glutamine metabolism-related genes obtained after differential analysis showed that BC patients were best divided into two clusters, with good internal consistency and stability of each Cluster ( Figure 2A ). The general characteristics of Cluster 1 and Cluster 2 are shown in Table S1 . The survival curve showed that Cluster2 had a poor prognosis (p < 0.05) ( Figure 2B ). There were significant differences in clinical parameters such as age, gender, subtype, grade, clinical stage, M and race of bladder cancer between the two clusters. Compared with patients with C1 bladder cancer, the proportions of patients older than 60, male, non-papillary invasive subtype, high grade, high stage, M1, white and Asian in patients with C2 bladder cancer were significantly higher than those in patients with C1 bladder cancer (p < 0.05) ( Figure 2C ).

We analyzed the immune microenvironment between the two clusters. Through ESITIMATE algorithm, according to the Cluster 2 has higher ESITIMATE score, immune score, score matrix and lower purity of tumor ( Figure 2D ). The CIBERSORT algorithm showed that Cluster 2 had higher immune cell infiltration ( Figure 2E ). Patients’ response to immune checkpoint inhibitors was negatively correlated with TIDE score. We found that TIDE score of Cluster 2 patients was significantly higher than that of Cluster 1 patients (p < 0.05), indicating that Cluster 2 patients had poor effect on immune checkpoint inhibitors ( Figure 2F ).

To construct a risk model associated with glutamine metabolism and its derived GMII, eight genes with independent prognostic value were identified by univariate and multivariate Cox analyses to construct GMII ( Figure 3A ). Coefficients for each gene in GMII ( Figure 3B ). Correlation analysis between GMII and survival status showed that higher GMII was associated with higher mortality ( Figure 3C ). Patients in the high-GMII group had a worse prognosis than those in the low-GMII group (p < 0.05) ( Figure 3D ). To verify the predictive value of this GMII, GSE13507 and GSE32894 were used to verify that the mortality rate in the high-GMII group was significantly higher than that in the low-GMII group (p < 0.05) ( Figures 3E, F ). The area under the 1-year, 3-year and 5-year ROC curves were 0.747, 0.714 and 0.743, respectively ( Figure 3G ), suggesting that GMII could better predict the short-term and long-term survival status of bladder cancer patients. The general characteristics of TCGA, GSE13507 and GSE32894 patients are shown in Table S2 . Univariate and multivariate regression analyses showed that GMII was an independent risk factor ( Figures 3H, I ). Based on the results of multivariate analysis, a nomogram was constructed based on GMII, age, clinical stage, and T stage, with GMII accounting for the majority of the total score ( Figure 3J ). Calibration curves showed that the predicted and actual 1 -, 3 -, and 5-year survival rates were consistent with the reference lines ( Figure 3K ).

The expression of five of the eight GMII genes (TALDO1, AHCY, FASN, GALK1 and SLC7A9) in bladder cancer tissues was higher than that in normal tissues, while ENPP1, CYP19A1 and HSPG2 were down-regulated in tumor tissues (p < 0.05) ( Figure S2A ). Immunohistochemical data of HPA showed that FASN showed moderate staining in the cytoplasm and nucleus of urothelial cancer cells, but no staining was detected in normal bladder tissues. HSPG2 showed high staining in the cytoplasm and nucleus of normal bladder cells, but no staining was detected in urothelial carcinoma cells ( Figure S2B ). Glutamine Related Genes (GRGs) is highly expressed in most tumor tissues and is a high risk factor ( Figures S3A, B ). Gene set variation analysis (GSVA) was performed on 8 GRGs, and GSVA score was positively correlated with the expression of representative gene sets. Most tumor tissues had significantly higher GSVA scores than normal tissues ( Figure S3C ). GSVA score was significantly correlated with the prognosis of bladder cancer (p < 0.05), including overall survival (OS), progression-free survival (PFS) and disease-free survival (DSS) ( Figure S3D ). We further evaluated the correlation between GSVA score and immune cell infiltration and showed that GSVA score was significantly correlated with immune cell infiltration in bladder cancer (p < 0.05) ( Figure S3E ). We also collected the features of published prognostic models for bladder cancer and compared the GMII features with their prognostic prediction accuracy. The results showed that GMII values were superior to other models in terms of prognostic prediction ( Figure S4 ).

To further verify the clinical significance of GMII, we analyzed the differences in GMII between different clinical characteristics groups. The results showed that patients with Cluster 2, Non-papillary infiltration, lymphovascular invasion, High Grade, stage III-IV, M1 and N1-3 had higher GMII, suggesting that the higher the GMII, the more advanced the tumor ( Figure 4A ). Stratified analysis showed that GMII could significantly differentiate the prognosis of almost all clinical subgroups, with patients in the high-GMII group having a worse prognosis ( Figure S5 ). To investigate the pathways that regulate tumorigenesis in the high-GMII group, we performed GSEA analysis, and the results showed that, The high GMII group was significantly enriched in high glutamine metabolism, glutamine synthesis and decomposition (p < 0.05) ( Figure 4B ). In addition, the high-GMII group was significantly enriched in carcinogenesis, angiogenesis, and epithelial-mesenchymal transition (EMT) pathways (p < 0.05) ( Figure 4C ). Glycolysis and myosynthesis were the features with the highest NES in the high-GMII group ( Figures 4D, E ).

Previous studies have shown that glutamine metabolism and tumor microenvironment play an important role in tumor development (7, 26). The ESTIMATE algorithm found that the high GMII group had lower tumor purity ( Figure 5A ) and higher ESTIMATE score ( Figure 5B ), immune score ( Figure 5C ) and stromal score ( Figure 5D ). CIBERSORT algorithm showed that the level of CD8 T cells in the high-GMII group was significantly lower than that in the low-GMII group (p < 0.05), and the level of immunosuppressive M2-type macrophages was increased in the high-GMII group, suggesting that the high-GMII group had higher immunosuppressive activity and promoted tumor progression ( Figure 5E ). We found correlations between GMII and a variety of immune cells using different software ( Figure 5F ). The ssGSEA algorithm results showed that the high-GMII group had higher immune cell infiltration and immune-related functions and pathways than the low-GMII group ( Figure 5G ).

TMB is an important predictor of immunotherapy and chemotherapy. To further investigate the association of GMII with somatic mutations in cancer cells, we used single nucleotide variation data to investigate differences in genomic mutations between high-GMII and low-GMII groups. TP53, TTN, KMT2D, MUC16 and ARID1A were the top five genes with the highest mutation frequency in high-GMII and low-GMII populations, but the mutation frequency of each gene was different between the two groups ( Figures 6A, B ). Through mutual exclusion and cooperation analysis among mutated genes, it was found that there was gene mutation synergy between most genes, and significant mutational mutual exclusion of TP53-ARID1A and EP300-ZFHX4 was found in the high-GMII group ( Figure 6C ). Significant mutational mutual exclusion of TP53-FGFR3 and KMT2D-FGFR3 was also found in the low-GMII group ( Figure 6D ). There were no differences in TMB between the two groups or in survival curves between the high and low TMB groups ( Figure S6A, B ). After TMB combined with GMII, the prognosis of the high TMB+ low GMII group was significantly better than that of the low TMB+ high GMII group ( Figure S6C ). In addition, the mutation rate of eight genes in GMII was detected in the cBioPortal database, and it was found that the mutation rate was low ( Figure 6E ).

To investigate the potential of GMII for predicting the response to chemotherapytic therapy, We first downloaded data from the GDSC database on the response of high-GMII and low-GMII populations to common chemotherapy agents. The results showed that many common bladder cancer chemotherapy drugs had significant differences among high-GMII groups (p<0.05). Gemcitabine, as the most common chemotherapy drugs for bladder cancer, had significantly lower IC50 values in the low-GMII group than in the high-GMII group (p<0.05), suggesting that Gemcitabine may have better efficacy in the low-GMII group ( Figure 7A ). We also mapped the 3D structures of chemotherapeutic agents with differences between the two groups using the PubChem database ( Figure S7 ). The relationship between GRGs and drug sensitivity was analyzed from cellMiner database, and the histograms of the top five drugs with the highest correlation between genes and drug sensitivity were drawn. Positive correlation indicates that stronger gene expression is more sensitive to drugs, while negative correlation indicates that stronger gene expression is more resistant to drugs. Among them, up-regulation of FASN expression may lead to enhanced sensitivity of patients to most drugs ( Figure 7B ). In addition, GDSC database was used to analyze the relationship between drug sensitivity and mRNA expression of the eight GRGs used to construct GMII. Contrary to cellMiner database, positive correlation represented that gene expression was related to drug resistance, while negative correlation represented that gene expression was related to drug sensitivity. The relationship between AHCY gene expression and chemotherapy-drug sensitivity is extensive, and HSPG2 has the highest correlation with most chemotherapy-drug sensitivity, among which HSPG2 is correlated with Dasatinib sensitivity in both databases ( Figure 7C ). These results suggest that the expression changes of genes constructing GMII may be effective indicators for predicting drug response and as potential therapeutic targets.

Immune checkpoint inhibitors have provided clinical benefits in the treatment of a variety of tumors. By analyzing the correlation between genes in GMII and common immune checkpoints, we found that genes with the highest correlation with glutamine metabolism were significantly negatively correlated with PD-1/PD-L1 expression. For example, SLC7A9, a glutamine transporter, and FASN, a key enzyme in the metabolism of glutamine into fatty acids. These results indicated an inverse correlation between glutamine metabolism and PD-1/PD-L1 expression ( Figure 8A ). The expression of most immune checkpoints was different between the high-GMII group and the low-GMII group (p < 0.05) ( Figure S8 ). Higher TIDE scores were associated with poorer immune checkpoint blockade therapy and shorter survival, and higher TIDE scores in the high-GMII group suggested poorer response to immune checkpoint blockade therapy ( Figure 8B ). The results of IMvigor210 dataset showed that there was a significant difference in GMII between the immunotherapy response group and the non-response group (p<0.05), and the GMII was lower in the response group ( Figure 8C ). The scores of IPS, IPS-PD1 blocker, IPS-CTLA4 blocker, and IPS-PD1-CTLA4 co-blocker were lower in the high-GMII group, indicating that the high-GMII group had a poor effect of anti-PD1, anti-CTLA4, and anti-PD1-CTLA4 co-treatment ( Figure 8D ). These results suggest that GMII may be associated with immunotherapy in bladder cancer patients.

In order to identify the core targets of genes related to glutamine metabolism, we constructed a PPI network through the STRING database (confidence score greater than 0.4), found the most closely related subnetwork to the GMII gene and found the core gene PPARG by “MCODE” in Cytoscape software ( Figure S9 ). PPARG is located at the hub of the network and has the highest degree among all nodes. SLC7A9 is potentially used as a glutamine transporter and GAPK1 as a targeted receptor with small molecule inhibitors of kinases. We obtained X-ray structures of PPARG, SLC7A9 and GALK1 from the PDB database to screen 1379 FDA-approved drugs potentially targeting PPARG, SLC7A9 and GALK1. We used rigid docking in MOE to simulate the binding mode of small molecule drugs to PPARG, SLC7A9 and GALK1, where the interaction relationship between target proteins and candidate small molecules is shown in Table S3 . We show the top four small-molecule drugs with the highest binding ability to PPARG (Bosulif, Candesartan, Centany, and Nefazodone) ( Figures 9A–D ) and the top four small-molecule drugs with the highest binding ability to SLC7A9 (Propantheline, Naloxeg) ol, Cobicistat and Fosinopril) ( Figures 9E–H ) and the top four small-molecule drugs with the highest GALK1 binding capacity (Propantheline, Ipratropium, Cangrelor and Lopinavir) ( Figures 9I–L ). For example, Bosulif (ZINC000022448983) forms hydrogen bonds with PPARG amino acid residues Met-348, Met-364 and His-449, among which Met-348 and Met-364 act as hydrogen bond acceptor and His-449 act as hydrogen bond donor. Naloxegol (ZINC000095564694) forms hydrogen bonds with SLC7A9 non-transmembrane amino acid residues Lys-53, Ser-57, Lys-145 and Thr-434, among which Ser-57 and Thr-434 act as hydrogen bond acceptor. Lys-53 and Lys-145 act as hydrogen bond donors. Cangrelor (ZINC000085537017) formed hydrogen bond and ion bond related interaction with GAPK1 amino acid residues Arg-228 and Arg-105, among which Arg-228 and Arg-105 were hydrogen bond donors. In addition, these small molecules form van der Waals (VDW) interactions with residues around the protein receptor, which contribute to the binding energy between small molecules and PPARG.