GSE32548 (18) performed on GPL6947 was downloaded from Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/) (19, 20), which included 128 BCs with complete clinical information. In this study, we constructed co-expression network based on this data set. Another GEO data set GSE13507 (21, 22) was also downloaded from this database for validation of our findings, performed on GPL6102, which had 165 BCs with clinical information. The Cancer Genome Atlas (TCGA) database (https://genome-cancer.ucsc.edu/) (23) characterized over 20,000 primary cancer and matched normal samples spanning almost all the cancer types. Thus, we retrieved BC microarray data from this database. After excluding unqualified samples, a total of 408 BC samples with complete clinical information were used for validation in this study.

For GSE32548 and GSE13507, we first downloaded the raw expression data and then normalized and transformed the data by using R package “affy” (24). As for TCGA-bladder urothelial carcinoma (BLCA) data displayed as counts number, normalization and log2 transformation were conducted by using package “DEseq.2” (25) in R software. The flow diagram of the present study is shown in Figure 1 in detail.

We first downloaded a comprehensive list of immune-related genes (IRGs) from the ImmPort database (https://immport.niaid.nih.gov) (26), which included 2,498 IRGs. Then 1,333 genes overlapped between IRGs, and genes in GSE32548 were included for the WGCNA. Two independent methods, including gsg (goodSamplesGenes) method and sample network method, were used for outlying sample identification of the expression profile of the 1333 genes. Only samples with Z.Ku ≥ −2.5 were used to construct a co-expression network by using R package “WGCNA” (samples were regarded as outliers when Z.Ku < −2.5, which were removed from the expression profile). Soft threshold power beta (β) was chosen by using the scale-free topology criterion (27, 28). In this study, genes were classified into gene modules by using the dynamic hybrid branch cutting method. Parameters for branch splitting were set as follows: minClusterSize = 30, and deepSplit = 2. In this study, we created a cluster tree of module eigengenes (MEs) to merge some highly related modules. Modules will be merged if their pairwise correlation was larger than 0.75. Moreover, a multidimensional scaling (MDS) plot was plotted to estimate the bio-similarity of the modules.

After gene module identification, we further screen out key module related to the trait (pathological stage), which interested us most. The correlation between genes and trait was quantified by calculating gene significance (GS). Module significance (MS) was further defined based on GS. Moreover, the relationship between ME and gene expression matrix was quantified by calculating module membership (MM). Finally, the module most positively correlated with pathological stage and the module most negatively related to pathological stage were considered as key modules in the present study.

First, we screened hub genes in key modules by measuring cor.geneModuleMembership and cor.geneTraitSignificance. In this study, only genes with |cor.geneModuleMembership| > 0.80 and |cor.geneTraitSignificance| > 0.20 were regarded as hub genes in the modules. By conducting Gene Ontology (GO) (29) enrichment analysis and Kyoto encyclopedia of Genes and Genomes (KEGG) (30) pathway analysis, we could understand biological meaning behind hub genes in key modules. Both the two analyses were identified based on R package “clusterProfiler” (31). Gene sets were regarded as significantly enriched gene sets when P < 0.05, as well as KEGG signaling pathways.

Based on the Search Tool for the Retrieval of Interacting Genes (STRING) database (https://string-db.org/) (32, 33), a PPI network of hub genes in the modules was constructed with a confidence > 0.4, maximum number of interactors = 0. Based on a tool called network analyzer in Cytoscape (34), the degree of connectivity of each hub gene was further calculated. Nodes in the PPI network represented proteins and edges represented protein-protein associations. Node color and node size were changed with degree of hub gene. Edge color and width were changed with combined score. Genes with the top 10 degree of connectivity were selected for further analysis. In addition, we used the MCODE plug-in in Cytoscape software to screen out hub modules in the PPI network, by using the following criteria: degree cutoff = 2, haircut on, node score cutoff = 0.2, k-core = 2, and max. depth = 100. Gene with the top 10 MCODE score were selected for further analysis. Furthermore, we used the cytoHubba plug-in in Cytoscape software to identify hub genes in the PPI network, the top 10 hub genes were screened out by using the Betweenness algorithm, maximum clique centrality (MCC) algorithm, and stress algorithm, respectively. After finishing these steps, genes that overlapped in the five methods were regarded as hub genes in the PPI network.

To explore the prognostic value of hub genes in the PPI network, we further performed overall survival (OS) analysis by using GEPIA (http://gepia.cancer-pku.cn/) (35), an online tool based on TCGA data. Hub genes with significant P value in this analysis were considered as immune-related prognostic biomarkers in the present study. In addition, disease-free survival (DFS) analysis was also performed. Furthermore, a stage plot (I, II, III, and IV) was also drawn, the statistical significance of which was measured by one-way analysis of variance (ANOVA) test.

After the surgery, a total of 20 samples, including 10 human BC tissues and 10 adjacent normal bladder tissues, were gathered from patients at Zhongnan Hospital of Wuhan University. The samples were histopathologically confirmed by two pathologists independently. The inclusion criteria are as follows (1): the histopathological type is confirmed as bladder urothelial carcinoma (BLCA) (2), not received anti-cancer treatment before cystectomy (3), underwent radical cystectomy or partial nephrectomy (4), no history of other malignant tumors. Exclusion criteria are as follows (1): other pathological types of BC (2), metastatic BLCA or other merge tumors (3), patients who did not undergo surgery, and (4) clinical pathological data are incomplete. Each patient signed an informed consent form, and the medical ethics committee in this hospital approved the utilization of tumor tissues for the present study. The approval number for this study was 2020174 (Kelun).

To validate the prognostic value of immune-related prognostic biomarker, we performed survival analysis for hub gene based on GSE13507 and GSE32548 by R package “survival” (36). BC samples were divided into low-expression group and high-expression group in the two data sets, respectively. Then, OS using GSE32548 and cancer-specific survival (CSS) analysis using GSE13507 were performed. A P value less than 0.05 was considered significant. Furthermore, we assessed the mRNA expression levels of hub genes in BC and normal tissue by Oncomine database (https://www.oncomine.org/) (37). Two independent data sets were used in the present study. In addition, we performed quantitative real-time PCR (qRT-PCR) analysis. The expression patterns of the CD86 genes were evaluated in BCs and adjacent normal bladder tissues. The HiPure Total RNA Mini Kit and RNAiso-Plus (TAKARA, China) were used to extract total RNA from the cells and 10 pairs of bladder cancer tissue and adjacent normal bladder tissues, which collected from the Zhongnan Hospital of Wuhan University, and we used NanoDrop to quantify the RNA, which was then reverse transcribed into cDNA by ReverTra Ace qPCR RT Kit (Toyobo, China). Finally, we performed qRT-PCR analysis of cDNA with iQTM SYBR^® Green Supermix (Bio-Rad) in a final volume of 20 μl. Relative gene expression was quantified via the 2^−△△Ct approach and normalized to glyceraldehyde-3-phosphate dehydrogenase (GAPDH) expression. The primer sequences for CD86 molecule (CD86) and GAPDH were listed in Table 1 . We measured the statistical significance by conducting Student t test. Finally, we validated the translation-level expression of hub gene between normal urinary bladder tissue and bladder urothelial carcinoma tissue by using The Human Protein Atlas (HPA) database (https://www.proteinatlas.org/) (38).

Immunocytes have been proven to be independent predictors of survival in cancers, thus, in this study, we investigated the relationship between expression levels of selected IRGS and immunocytes based on tumor immune estimation resource (TIMER) (https://cistrome.shinyapps.io/timer/). TIMER was a webtool, which could estimate the abundance of immune infiltrates for six tumor-infiltrating immune cell types. We thought an immune-related prognostic biomarker is strongly related to an infiltrating level of an immune cell type when |correlation coefficient (cor) |is 0.2 or greater and a P value is less than 0.05.

In this part, we first calculated tumor purity, immune score, and stromal score for each BC sample collected from TCGA-BLCA data (n = 408) by applying estimation of stromal and immune cells in malignant tumors using expression data (ESTIMATE) algorithm (by using “estimate” package in R software) (39). Considering that immune cells played important roles in the tumor microenvironment (TME) (40, 41). Relying on ssGSEA (a R package in R software) (42, 43), the relative infiltration of 28 kinds of immunocytes was quantified. A gene list of metagenes, which contained feature genes symbolizing for each immunocyte type, was retrieved from an article in Cell Reports (44). In ssGSEA, we calculated enrichment score for each immunocyte type, which represented the relative abundance of immunocyte. Zero was the minimal meanwhile one was the maximal score. A heatmap was further composed for visualization of the relative abundance of all kinds of immune cell types. Finally, we plotted an MDS plot and constructed a Gaussian fitting model for estimation of the bio-similarity of the immune cell filtration.

In this study, we performed GSEA (45) based on TCGA-BLCA data from TCGA database to explore the potential functions of immune-related prognostic biomarkers. First, 408 BCs were divided into low expression group (n = 204), and high expression group (n = 204) by setting the median of CD86 expression as a cut line. “c2.cp.kegg.v7.0.symbols.gmt” was set as the reference gene set. In this study, we thought KEGG signaling pathways were significantly enriched when nominal P value is less than 0.05, |enrichment score (ES)| is greater than 0.6, gene size is 100 or greater, and false discovery rate (FDR) is less than 25%.

A comprehensive list of 2,498 IRGs was first downloaded from ImmPort database. Then, 1,333 genes overlapped between IRGs and genes in GSE32548 were analyzed and further included for WGCNA, the information of the 1,333 genes is shown in detail in Table S1 .

Six outliers were first removed from the expression profile of the 1,333 IRGs by using two independent methods ( Figures S1A, B ), only 122 samples (with complete clinical information) from GSE32548, which reached the standards, were included for subsequent analysis. As shown in Figures S1C–F , beta (β) = 2 (scale free R² = 0.95) was further chosen for adjacency calculation. Then, we classified genes into gene modules. Moreover, we merged modules highly correlated to each other (pairwise correlation of modules > 0.75). In total, four modules including turquoise (n = 372), yellow (n =52), blue (n = 170), and brown (n = 55) were screened out ( Figure 2A ). The rest of the IRGs, including 684 genes, showed weak correlation with other genes, which were excluded for next-step analysis (these genes were discarded to gray module). Then, we noticed that the turquoise module showed the most positive association with pathological stage (P = 5E-10, R² = 0.53), meanwhile the blue module was the most negatively associated module with pathological stage compared with others (P = 8E-10, R² = −0.52), as suggested in Figure 2B . The association between MM and GS in the two modules is shown in Figures 2C, D , separately. MM in blue module was significantly associated with GS in blue module (cor = 0.64, P = 5.7E-21), a similar trend also existed in turquoise module (cor = 0.59, P = 2.9E-36). Moreover, the MS of the two modules was significantly higher than that of any other modules ( Figure 2E ). Therefore, blue module and turquoise module were screened out and regarded as key modules in the present study. As shown in Figure 2F , the MDS plot demonstrated that each module was isolated from each other, especially the key modules.

In this study, 43 genes with |cor.geneModuleMembership| > 0.80 and |cor.geneTraitSignificance| > 0.20 were considered as hub genes in the turquoise module, meanwhile two genes that reached the same standards were regarded as hub genes in the blue module. The detailed information of these hub genes are shown in Table S4 .

As Table S2 shows, the 45 hub genes were significantly enriched in 361 biological processes (BPs). The top 10 enriched BPs were positive regulation of leukocyte activation, positive regulation of cell activation, T cell activation, positive regulation of cell-cell adhesion, positive regulation of lymphocyte activation, positive regulation of leukocyte cell-cell adhesion, regulation of leukocyte activation, leukocyte cell-cell adhesion, positive regulation of T cell activation, and regulation of leukocyte cell-cell adhesion ( Figure 3A ). Furthermore, these hub genes were significantly related to 38 KEGG signaling pathways ( Table S3 ). The top 10 KEGG pathways were rheumatoid arthritis, hematopoietic cell lineage, Th17 cell differentiation, cytokine-cytokine receptor interaction, Th1 and Th2 cell differentiation, cell adhesion molecules (CAMs), primary immunodeficiency, viral protein interaction with cytokine and cytokine receptor, T cell receptor signaling pathway, and intestinal immune network for IgA production ( Figure 3B ).

Then, we constructed a PPI network for the 45 hub genes identified before ( Figure 4A ). The degree of connectivity of each gene was calculated, which is shown in detail in Table S4 . Genes with the top 10 degree of connectivity were selected, including CD86, PTPRC, CTLA4, IL10RA, ITGB2, CCL5, TYROBP, CCL2, GZMB, and TLR8. Moreover, the MCODE plug-in revealed one important functional module in the PPI network ( Figure 4B ). Genes with the top 10 MCODE score were identified, including CD86, IL10RA, ITGB2, PTPRC, IL2RA, GZMB, IL7R, BTK, FCER1G, and TLR8. Then, 10 genes, including CD86, PTPRC, CCL5, IFI30, CTLA4, CCL2, ITGB2, TYROBP, IL10RA, and FCER1G, were screened out by using betweenness algorithm applied by cytoHubba plug-in ( Figure 4C ). Meanwhile, 10 genes, including CD86, PTPRC, IL10RA, CTLA4, CCL2, ITGAL, CCL5, IL2RB, GZMB, and IL2RA, were identified by MCC algorithm ( Figure 4D ). In addition, 10 genes, including CD86, PTPRC, CCL5, CCL2, CTLA4, ITGB2, IL10RA, TYROBP, and IFI30, were screened out by applying stress algorithm ( Figure 4E ). Finally, three genes, including CD86 molecule (CD86), IL10RA (interleukin 10 receptor subunit alpha), and protein tyrosine phosphatase receptor type C (PTPRC), overlapped among genes selected by the five methods were picked out ( Figure 4F ). We regarded the three genes as hub genes in the PPI network. We immediately performed OS analysis for the three genes, the result demonstrated that expression of CD86 could impact the survival and prognosis of patients with BC. Patients with high CD86 expression had worse OS ( Figure 5A , P = 0.049), meanwhile there was a trend that high expression of CD86 caused worse DFS compared with low expression ( Figure 5B , P = 0.2). A stage plot was shown in Figure 5C , and CD86 presented different expression levels in different stage (II, III, and IV). As for IL10RA ( Figures 5D, E ) and PTPRC ( Figures 5G, H ), there was no obvious association between their expressions and survival. Stage plots were also plotted as shown in Figures 5F, I . Thus, only CD86 was considered as immune-related prognostic biomarker in the present study.

Based on GSE32548 and GSE13507, we validated the prognostic value of CD86. The result demonstrated that BC patients with high expression of CD86 had worse OS compared with these with low expression of CD86 ( Figure 6A , P = 0.05). Also, high expression of CD86 in BC patients was significantly associated with worse CSS, as shown in Figure 6B . These results were consistent with what we found by using GEPIA. Then based on BC data from Oncomine database, we compared the expression of CD86 between BC tissues and normal bladder tissues. The result suggested that the CD86 mRNA expression was higher in BCs compared to normal tissues (P = 0.013, Figure 6C ). We further validated the mRNA expression of CD86 by using qRT-PCR analysis again. As shown in Figure 6D , CD86 was significantly up-regulated in BC tissue compared with adjacent normal bladder tissue (P = 0.0382), which was consistent with the findings based on Oncomine database. By using HPA database, we also explored the translational-level expression of CD86, but not as we imagined, there was no significant difference between the expression in normal bladder tissue and bladder cancer tissue ( Figure 6E ).

Immune infiltration was reported to be associated with survival and progression of cancers. Thus, by using TIMER (a webtool), the association between CD86 and immune infiltration level was obtained. As shown in Figure 7 , CD86 was positively associated with CD8+ T cells (cor = 0.374, P = 1.38E-13), CD4+ T cells (cor = 0.358, P = 1.76E-12), macrophages (cor = 0.258, P = 5.84E-07), neutrophils (cor = 0.713, P = 1.57E-57), and dendritic cells (cor = 0.701, P = 3.67E-55).

After calculating immune score and stromal score of each BC from TCGA-BLCA data, immune scores ranged from −1900.04 to 2903.20 meanwhile stomal scores ranged from −2496.63 to 2148.31 as the result suggested ( Table S5 ). As shown in Figure 8A , the heatmap demonstrated that CD86 expression was significantly associated with ESTIMATE score, immune score, and stromal score, positively; meanwhile, it negatively correlated to tumor purity, which was consistent with the TIMER analysis. Furthermore, CD86 high-expression samples were associated with a higher abundance of immune cell infiltration. The ssGSEA enrichment scores (after normalization) for each immune cell type in 408 BC patients are shown in detail in Table S6 . More concretely, high CD86 expression of patients was related to a higher abundance of immune cells executing anti-tumor reactivity, including activated CD4⁺ T cells, activated CD8⁺ T cells, activated dendritic cells, central memory CD4⁺ T cells, central memory CD8⁺ T cells, effector memory CD4⁺ T cells, effector memory CD8⁺ T cells, natural killer cells, natural killer T cells, and type 1 T helper cells. The abundances of immune cells, which delivered pro-tumor suppression (including macrophages, myeloid-derived suppressor cells (MDSCs), and regulatory T cells), were positively associated with CD86 expressions in patients with BC. Further analysis suggested that immune cells executing anti-tumor reactivity were positively related to immune cells delivering pro-tumor suppression within a local environment, significantly ( Figure 8B , cor = 0.9204, P < 0.001). Interestingly, these results were consistent with what we got by TIMER, which reflected that anti-tumor inflammation might facilitated the recruitment or differentiation of cells specialized for immune suppression in BC.

To explore the potential function of CD86, we performed GSEA. As shown in Table 2 , the result suggested that the high expression of CD86 was significantly correlated to nine KEGG signaling pathways, including systemic lupus erythematosus, cytokine-cytokine receptor interaction, cell adhesion molecules (CAMs), toll-like receptor signaling pathway, chemokine signaling pathway, natural killer cell-mediated cytotoxicity, T-cell receptor signaling pathway, focal adhesion, and leukocyte transendothelial migration.

Also, based on GSEA, we downloaded “DSigDBv1.0.gmt” from Drug SIGnatures Database (http://tanlab.ucdenver.edu/DSigDB/DSigDBv1.0/download.html) (46) and attempted to screen out some novel drugs for BC treatment. The result demonstrated that CD86 was associated with various drugs ( Table S7 ). Among them, suloctidil with the highest |ES| was further screened out (nominal P = 0.000, |ES| = 0.824, gene size = 137, and FDR = 6.345%), which showed powerful potential to treat BC.