Prognostic model of 10 immune-related genes and identification of small molecule drugs in bladder urothelial carcinoma (BLCA).

BACKGROUND: We aimed to establish an immune-related gene (IRG) based signature that could provide guidance for clinical bladder cancer (BC) prognostic surveillance.
METHODS: Differentially expressed IRGs and transcription factors (TFs) between BCs and normal tissues were extracted from transcriptome data downloaded from the TCGA database. Gene Ontology (GO) function and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were carried out to identify related pathways based on differently expressed IRGs. Then, univariate Cox regression analysis was performed to investigate IRGs with prognostic values and LASSO penalized Cox regression analysis was utilized to develop the prognostic index (PI) model.
RESULTS: A total of 411 BC tissue samples and 19 normal bladder tissues in the TCGA database were enrolled in this study and 259 differentially expressed IRGs were identified. Networks between TFs and IRGs were also provided to seek the upstream regulators of differentially expressed IRGs. By means of univariate Cox regression analysis, 57 IRGs were analyzed with prognostic values and 10 IRGs were finally identified by LASSO penalized Cox regression analysis to construct the PI model. This model could significantly classified BC patients into high-risk group and low-risk group in terms of OS (P=9.923e-07) and its AUC reached 0.711. By means of univariate and multivariate COX regression analysis, this PI was proven to be a valuable independent prognostic factor (HR =1.119, 95% CI =1.066-1.175, P<0.001). CMap database analysis was also utilized to screen out 10 small molecules drugs with the potential for the treatment of BC.
CONCLUSIONS: Our study successfully provided a novel PI based on IRGs with the potential to predict the prognosis of BC and screened out 10 small molecules drugs with the potential to treat BC. Besides, networks between TFs and IRGs were also displayed to seek its upstream regulators for future researches. 2020 Translational Andrology and Urology. All rights reserved.

Introduction

Bladder cancer (BC), as one of the major causes of cancer related mortality worldwide, it attributed to nearly 549,000 new cases and 200,000 deaths in the United States, according to Global Cancer Statistics 2018 (1). Tobacco smoking and occupational exposure to certain chemical carcinogens remained the main risk factors of BC (2-4). Because of the gradual popularization of bladder B-ultrasound, transurethral cystoscopy, etc., many BC patients could be diagnosed in early time and receive timely and effective treatment, including transurethral resection, radiation therapy, and chemotherapy. However, the prognosis for BC patients was still unsatisfactory, due to frequent recurrence and metastasis (5,6).

A growing number of evidence had highlighted the important roles of immune-related genes (IRGs) in cancer initiation, progression, metastasis and so on. Therein, CD58, a cell-surface glycosylated adhesion molecule, could enhance the self-renewal of colorectal tumor-initiating cells and promote epithelial-mesenchymal transition (EMT) process through activating Wnt/β-catenin pathway in colorectal cancers (7). The overexpression of ZBTB20, a transcriptional repressor, could lead to enhancement of migration and invasion ability of gastric cancer by inhibiting IκBα to activate NF-κB (8). Besides, ZFP36L2, a RNA-binding protein, had been shown to inhibit cell proliferation in pancreatic ductal adenocarcinoma (PDAC), and it was associated with good prognosis of PDAC patients (9). IRGs also played very important roles in BC. Martínez et al. demonstrated that BMP-4, secreted by BC cells, induced monocyte/macrophage polarization toward an M2 phenotype, and promoted tumor progression (10). Ramakrishnan et al. found that inhibition of EZH2 in muscle-invasive BC with KDM6A and SWI/SNF family member mutations could activate a natural killer (NK) cell-based immune response to drive tumor differentiation and death in BC cells and xenografts (11). Cheah et al. elucidated the importance of CD14 in BC, and showed that high levels of CD14 resulted in the increase of inflammation mediators and drove tumor proliferation to promote tumor growth (12). Although a number of previous studies had proposed the roles of IRGs in BC, there were few focusing on their importance to systematically predict overall survival (OS) of BC patients.

As a result, we aimed to establish a signature based on IRGs that could predict OS for BC patients in this study. We utilized the transcriptome data downloaded from The Cancer Genome Atlas (TCGA) to develop a prognostic index (PI) by using univariate Cox regression analysis and least absolute shrinkage and selection operator (LASSO) regression analysis. To evaluate the clinical values of this PI, we analyzed whether this signature was associated with the survival outcome of BC patients and clinicopathological factors. We present the following article in accordance with the TRIPOD reporting checklist (available at http://dx.doi.org/10.21037/tau-20-696).

Methods

Data acquisition

Transcriptome profiling data and related clinical information of bladder urothelial carcinoma (BLAC) were downloaded and extracted from The Cancer Genome Atlas (TCGA) Data Portal (https://tcga-data.nci.nih.gov/tcga/; accessed August 2019). Overall, we enrolled a total of 430 BLAC cases containing 411 BC tissue samples and 19 normal bladder tissues into this study. All procedures performed in this study were in accordance with the Declaration of Helsinki (as revised in 2013).

Identification and enrichment analysis of differently expressed IRGs

We applied “Lima” package in R statistical software to estimate differently expressed IRGs between BC and normal bladder tissues. IRGs exhibiting at least 2-fold changes with an adjusted P value less than 0.05 were identified as the remarkably differentially expressed IRGs. We also performed gene functional enrichment analyses, including Gene Ontology (GO) functional annotations and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis, to analyze the biological processes, cellular components, molecular functions and signaling pathways of targeted genes by use of “clusterProfiler” R package.

Identification of differently expressed transcription factors (TFs) and construction of the correlation network between TFs and IRGs

We applied the Cistrome database (http://www.cistrome.org/) to predict TF targets and extract enhancer profiles in cancers. The prediction was from integrative analysis of TCGA expression profiles and public ChIP-seq profiles. Then we identified differently expressed TFs by use of “Lima” package in R statistical software between BLAC and normal bladder tissues. Correlation test between differently expressed TFs and IRGs was conducted by R programming language. Moreover, we set correlation coefficient at least 0.4 with a P value less than 0.01 as the remarkably correlated.

Construction of the PI model based on IRGs

We performed univariate Cox regression analyses to select IRGs significantly associated with OS of BC patients. LASSO-penalized Cox regression model was used to identify the most significantly survival related IRGs. The formula used for calculating the PI is as follow β1 × gene1 expression + β2 ×gene2 expression + ····· + βn × gene expression, where β corresponded to the correlation coefficient.

Evaluation of the PI model

We calculated the risk scores of each patient, and set the median value as the cut-off point to divide the BC patients enrolled into a high-risk group and a low-risk group. Kaplan-Meier plots were generated to analyze the performance of PI model in predicting survival outcomes. The receiver operating characteristic curves (ROC) were created, and the area under the curve (AUC) values was calculated to assess the model specificity and sensitivity. The risk score distribution of patients in different risk groups, and the heatmap of finally confirmed IRGs were also displayed. A prognostic nomogram was also generated to analyze the relationship between individual predictors and survival rates in BC patients. To further determine whether the PI model can be used as an independent prognostic factor, univariate and multivariate cox regression analysis were performed to calculate the HR values for the PI and other clinicopathologic parameters including age, sex, stage, race, grade as well as T, M, N stage. Furthermore, we also conducted a clinical correlation analysis to analyze the correlation between PI and clinicopathologic parameters described above.

Identification of candidate small molecule drugs

Connectivity map (CMap) facilitated researchers to quickly identify molecule drugs highly correlated with diseases and discover its possible mechanism (13). Connectivity scores ranged from −1 to 1 and they were utilized to estimate how closely a compound was connected to the query signature. Negative score indicated that the query signature could be repressed by a drug, while a positive score could be promoted by a drug.

Statistical analysis

We applied R software (version 3.4.1, https://www.r-project.org/) to complete date statistical analyses in this article. It was considered statistically significant when the difference met a joint satisfaction of FDR <0.05 and fold changes ≥2, To establish a prognostic model, univariate and LASSO-penalized COX regression analyses were performed to assess the relationship between IRGs expression and survival data. The receiver operating characteristic curves were created by the “survival ROC” package of R and AUC values were also calculated by this package too. Univariate and multivariate cox regression analyses were to further determine whether our established PI model could be used as an independent prognostic factor. All statistical tests were two-sided and P<0.05 was regarded to be statistically significant.

Results

Differentially expressed IRGs

The overall flow diagram for this study was displayed in . RNA-seq data of 411 BC tissue samples and 19 normal bladder tissues were downloaded from TCGA. A total of 259 differentially expressed IRGs were extracted with their expression values. The heatmap of differentially expressed IRGs between the BC tissues and normal bladder tissues were demonstrated in . According to the criteria of a FDR <0.05 and |log2(fold change)| >1, we finally obtained 119 up-expressed and 140 down-expressed IRGs ().

Figure 1

The flow diagram for this study

Figure 2

Differentially expressed immune-related genes (IRGs) between bladder cancer (BC) tissues and normal bladder tissues. (A) Heatmap of differentially expressed IRGs expression levels; (B) the volcano plot for differentially expressed IRGs from the TCGA data portal. Red indicates high expression and green low expression. Black shows those genes showed no difference between BC and normal bladder tissues; (C) the histogram of the numbers of the up- and down-regulated IRGs.

Functional enrichment analysis

On the basis of the RNA-seq data downloaded from TCGA, we carried out GO term enrichment analysis and KEGG pathway analysis to study the biological function of the differentially expressed IRGs. The GO terms function and KEGG pathway enrichment of these genes were summarized in .

Table 1

GO and KEGG analysis of differentially expressed immune-related genes

Category	ID	Term	Gene ratio
Biological process	GO:0032103	Positive regulation of response to external stimulus	41/252
Biological process	GO:0050920	Regulation of chemotaxis	35/252
Biological process	GO:0030335	Positive regulation of cell migration	45/252
Biological process	GO:0060326	Cell chemotaxis	35/252
Biological process	GO:0050900	Leukocyte migration	43/252
Biological process	GO:0050921	Positive regulation of chemotaxis	25/252
Biological process	GO:0048660	Regulation of smooth muscle cell proliferation	25/252
Biological process	GO:0048659	Smooth muscle cell proliferation	25/252
Biological process	GO:0030595	Leukocyte chemotaxis	28/252
Biological process	GO:0002685	Regulation of leukocyte migration	26/252
Cellular component	GO:0043235	Receptor complex	24/254
Cellular component	GO:0009897	External side of plasma membrane	15/254
Cellular component	GO:0005788	Endoplasmic reticulum lumen	15/254
Cellular component	GO:0060205	Cytoplasmic vesicle lumen	16/254
Cellular component	GO:0031983	Vesicle lumen	16/254
Cellular component	GO:0031012	Extracellular matrix	19/254
Cellular component	GO:0034774	Secretory granule lumen	15/254
Cellular component	GO:0002116	Semaphorin receptor complex	3/254
Cellular component	GO:0030670	Phagocytic vesicle membrane	6/254
Cellular component	GO:0005884	Actin filament	7/254
Molecular function	GO:0048018	Receptor ligand activity	89/252
Molecular function	GO:0008083	Growth factor activity	39/252
Molecular function	GO:0005125	Cytokine activity	42/252
Molecular function	GO:0005126	Cytokine receptor binding	43/252
Molecular function	GO:0005179	Hormone activity	24/252
Molecular function	GO:0001664	G-protein coupled receptor binding	29/252
Molecular function	GO:0008009	Chemokine activity	14/252
Molecular function	GO:0042379	Chemokine receptor binding	15/252
Molecular function	GO:0070851	Growth factor receptor binding	19/252
Molecular function	GO:0005539	Glycosaminoglycan binding	22/252
KEGG pathway	hsa04060	Cytokine-cytokine receptor interaction	49/201
KEGG pathway	hsa04061	Viral protein interaction with cytokine and cytokine receptor	21/201
KEGG pathway	hsa04080	Neuroactive ligand-receptor interaction	33/201
KEGG pathway	hsa04010	MAPK signaling pathway	30/201
KEGG pathway	hsa04657	IL-17 signaling pathway	17/201
KEGG pathway	hsa05163	Human cytomegalovirus infection	23/201
KEGG pathway	hsa04062	Chemokine signaling pathway	21/201
KEGG pathway	hsa04668	TNF signaling pathway	16/201
KEGG pathway	hsa05167	Kaposi sarcoma-associated herpesvirus infection	20/201
KEGG pathway	hsa04650	Natural killer cell mediated cytotoxicity	16/201

As to GO term enrichment analysis, for biological process, the differentially expressed IRGs enriched in positive regulation of response to external stimulus, regulation of chemotaxis, positive regulation of cell migration, cell chemotaxis, leukocyte migration, positive regulation of chemotaxis, regulation of smooth muscle cell proliferation, smooth muscle cell proliferation, leukocyte chemotaxis, regulation of leukocyte migration; for cellular component, the differentially expressed IRGs enriched in receptor complex, external side of plasma membrane, endoplasmic reticulum lumen, cytoplasmic vesicle lumen, vesicle lumen, extracellular matrix, secretory granule lumen, semaphorin receptor complex, phagocytic vesicle membrane, actin filament; For molecular function, the differentially expressed IRGs enriched in receptor ligand activity, growth factor activity, cytokine activity, cytokine receptor binding, hormone activity, G-protein coupled receptor binding, chemokine activity, chemokine receptor binding, growth factor receptor binding, glycosaminoglycan binding ().

Figure 3

Function enrichment analysis of the differentially expressed IRGs between BC tissues and normal bladder tissues. (A) The bubble plot of enriched GO terms. Greed circles showed the biological process, red corresponded to the cellular component, and blue indicated the molecular function category. (B) KEGG analysis of differentially expressed IRGs. The outer circle shows a scatter plot for each term of the logFC of the assigned genes. Red circles display up-regulation, and blue ones down-regulation. (C) The heatmap of the relationship between IRGs and pathways.

With respect to KEGG pathway analysis, the result showed that the differentially expressed IRGs mainly enriched in cytokine-cytokine receptor interaction, viral protein interaction with cytokine and cytokine receptor, neuroactive ligand-receptor interaction, MAPK signaling pathway, IL-17 signaling pathway, and so on (). The Z-score of enriched pathways less than zero indicated that the related pathways were more likely to be decreased. The heatmap of the relationship between IRGs and pathways was also displayed ().

Network analysis of TF-IRG interaction

To explore the interaction between transcription factors (TFs) and IRGs, we extracted 77 differentially expressed TFs by analyzing the RNA-seq data described above. The heatmap of TFs with markedly different expression between the BC tissues and normal bladder tissues were shown in with 41 up-expressed and 36 down-expressed TFs (). The TF-IRG interaction network was detailed in .

Figure 4

Network analysis of TF-IRG interaction. (A) Heatmap of differentially expressed TFs; (B) Volcano map of differentially expressed TFs; (C) network shows the relationships between TFs and IRGs.

Construction of PI model

Univariate Cox regression analysis was then performed to investigate the prognostic values of the 259 differentially expressed IRGs. As shown in and , a total of 57 IRGs were found to be significantly associated with OS in BC patients (P<0.05). Next, LASSO regression analysis was performed to develop the PI model with 10 relevant IRGs (). Their specific coefficients of the 10 genes are shown in . The PI was calculated by use of the formula described above: Risk score = (−0.000824 × TAP1 expression) + (0.003970 ×RBP7 expression) + (−0.000031 × STAT1 expression) + (0.005626 × PDGFRA expression) + (0.004692 ×AHNAK expression) + (0.000290 × OLR1 expression) + (0.009473 × RAC3 expression) + (0.001882 × EDNRA expression) + (0.078175 × IGF1 expression) + (-0.009884 ×SH3BP2 expression).

Figure 5

Construction of PI model. (A) 57 IRGs with significant prognostic values after univariate cox regression analysis; (B) LASSO coefficient profiles of these 57 IRGs; (C) tuning the penalty parameter in LASSO using 10-fold cross validation.

Table S1

Original file for forest plot of 57 IRGs with significant prognostic values after univariate cox regression analysis

ID	HR	HR.95L	HR.95H	P value
CALR	1.000832	1.00005	1.001615	0.037021
TAP1	0.994956	0.991441	0.998484	0.005106
TAP2	0.975999	0.956671	0.995716	0.017283
THBS1	1.00428	1.001331	1.007239	0.004428
CXCL10	0.997415	0.995293	0.999542	0.017235
CXCL12	1.012369	1.004198	1.020606	0.002946
ZC3HAV1L	1.121497	1.060489	1.186014	5.87E-05
MMP9	1.000238	1.000062	1.000413	0.007942
FABP6	0.986059	0.973539	0.99874	0.031285
PAEP	1.03988	1.011951	1.068581	0.004875
RBP7	1.01316	1.006578	1.019786	8.45E-05
TFRC	1.00385	1.000564	1.007147	0.021634
IFIH1	0.972317	0.948467	0.996766	0.026721
ADIPOQ	1.086546	1.040249	1.134902	0.000187
STAT1	0.995682	0.992268	0.999109	0.013571
ISG15	0.999186	0.998378	0.999995	0.048524
ELN	1.015566	1.00412	1.027143	0.007564
CACYBP	1.01984	1.003647	1.036293	0.016138
BST2	0.998986	0.998026	0.999947	0.038713
PDGFRA	1.043412	1.01624	1.071311	0.001596
AHNAK	1.00646	1.004475	1.00845	1.64E-10
APOBEC3H	0.865297	0.758099	0.987653	0.032027
KCNH2	1.027955	1.015025	1.041049	1.96E-05
PTX3	1.007883	1.000932	1.014883	0.026157
IRF9	0.863686	0.77308	0.96491	0.009551
ANXA6	1.007813	1.00062	1.015058	0.033205
OAS1	0.986895	0.978456	0.995407	0.002608
OLR1	1.0073	1.001886	1.012744	0.008165
RAC3	1.025184	1.01452	1.035959	3.13E-06
NFATC1	1.084531	1.004895	1.170479	0.037028
NFATC4	1.046833	1.001545	1.094168	0.042519
SLIT2	1.136333	1.023269	1.261889	0.016842
EDNRA	1.084633	1.03594	1.135614	0.000527
CMTM8	1.029703	1.005087	1.054922	0.01774
IGF1	1.33133	1.17812	1.504465	4.48E-06
IL34	1.037426	1.010843	1.064708	0.005533
KITLG	1.022423	1.00301	1.042213	0.023378
MDK	0.998664	0.99734	0.999989	0.048131
NAMPT	1.009265	1.001335	1.017258	0.021938
NTF3	0.573751	0.332651	0.989595	0.045763
OGN	1.035817	1.000941	1.071909	0.044031
PDGFD	1.073897	1.028337	1.121476	0.001267
PGF	1.03297	1.016504	1.049702	7.60E-05
PPY	1.016415	1.007357	1.025554	0.000364
SPP1	1.000133	1.000022	1.000244	0.018862
TGFB3	1.030459	1.005894	1.055625	0.014797
ADRB2	0.899625	0.812037	0.996659	0.042971
AGTR1	1.117795	1.017286	1.228234	0.020532
ANGPTL1	1.024018	1.004578	1.043833	0.015221
IL17RD	1.063009	1.012719	1.115797	0.013471
IL17RE	1.041394	1.009157	1.074662	0.011468
NR3C2	1.175863	1.018721	1.357245	0.026872
NRP2	1.040675	1.008015	1.074394	0.014259
OXTR	1.030442	1.004925	1.056607	0.019081
PTGER3	1.312566	1.118447	1.540377	0.000866
TGFBR2	1.014858	1.004466	1.025358	0.004978
TNFRSF25	0.945021	0.898439	0.994017	0.028333
SH3BP2	0.899187	0.832921	0.970725	0.006515

Table 2

Coefficients of these 10 key prognostic immune-related genes (IRGs)

Gene	Coef
TAP1	−0.000823994
RBP7	0.003969584
STAT1	−3.11E-05
PDGFRA	0.005626014
AHNAK	0.004692076
OLR1	0.000290185
RAC3	0.009472789
EDNRA	0.001881597
IGF1	0.078174754
SH3BP2	−0.009883798

IRGs, immune-related genes; Coef, coefficients.

Evaluation of our established PI model

According to the formula described above, we calculated the risk scores of each patient to divide the BC patients into high-risk groups and low-risk groups, based on the median value as the cut-off point. As shown in , with the increase of the risk scores, the patients in high-risk group increased, and the number of dead persons grew. The heatmap of the 10 key prognostic genes expression profiles in the TCGA dataset was also demonstrated. To identify the performance of PI in predicting the clinical outcome of BC patients, the Kaplan-Meier plots were carried out to analyze the different survival time between the high- and low-risk groups. The results showed that patients in the high-risk group had a shorter OS time than patients in the low-risk group (P=9.923e-07, ). The ROC curve analysis was performed to evaluate how well the PI predicts the prognoses of BC patients. The AUC for the PI was 0.711, demonstrating the important value of the PI to predict survival of BC patients (). For providing a quantitative approach to predict cancer patient survival, we assembled a nomogram, integrating all the 10 key prognostic genes, which could help estimate 1-, 3-, and 5-year survival probabilities ().

Figure 6

Evaluation of this PI model. (A) The risk score distribution of patients in the BC patients; the number of dead persons grew with increased risk scores; the heatmap of the 10 key prognostic IRGs expression profiles in the TCGA dataset. (B) Kaplan-Meier plot represents that patients in the high-risk group had significantly shorter overall survival time than those in the low-risk group (P=9.923e-07). (C) The ROC curve of this PI model with its AUC of 0.711. (D) Prognostic nomogram for BC patients, integrating 10 key prognostic IRGs and our established PI model.

We next performed univariate and multivariate Cox regression analysis to evaluate prognostic value of the PI and other clinical parameters including age, gender, race, tumor stage, T, M, N. The hazard ratios (HRs) for PI in the univariate Cox regression analyses was HR =1.131, 95% confidence interval (CI) =1.080–1.185 (P<0.001, and ), and in the multivariate Cox regression analyses was HR =1.119, 95% CI =1.066–1.175 (P<0.001, and ). As displayed in , it detailed the AUC of all prognostic factors (age, gender, race, tumor stage, T, M, N and risk score) and our established PI model had the highest AUC value. We also integrated both the PI and the clinicopathologic risk factors mentioned above to assemble a nomogram, helping estimate 1-, 3-, and 5-year survival probabilities (). The univariate and multivariate Cox regression analysis of the PI and OS-associated clinical features were shown in .

Figure 7

Univariate/multivariate Cox regression analysis, multi-ROC analyses and prognostic nomogram; (A) Univariate Cox regression analysis of our established PI model and seven clinical features; (B) multivariate Cox regression analysis of our established PI model and seven clinical features; (C) Multi-ROC analyses of our established PI model and seven clinicopathologic risk factors; (D) prognostic nomogram for BC patients, integrating our established PI model and seven clinicopathologic risk factors.

Table S2

Original file for forest plot of univariate Cox regression analysis of our established PI model and seven clinical features

ID	HR	HR.95L	HR.95H	P value
Age	1.037102	1.018346	1.056203	9.14E-05
Gender	0.795505	0.547983	1.154833	0.228979
Race	1.203104	0.858987	1.685078	0.282066
Stage	1.879957	1.490326	2.371453	9.98E-08
T	1.731316	1.348288	2.223155	1.69E-05
M	1.21501	1.022681	1.443509	0.026756
N	1.238577	1.099637	1.395072	0.000424
Risk score	1.131252	1.079995	1.184941	1.86E-07

Table S3

Original file for forest plot of multivariate Cox regression analysis of our established PI model and seven clinical features

ID	HR	HR.95L	HR.95H	P value
Age	1.039269	1.019976	1.058927	5.61E-05
Gender	0.780127	0.532055	1.143863	0.203512
Race	0.860632	0.605305	1.223659	0.403232
Stage	1.397834	1.037331	1.883622	0.02775
T	1.445184	1.05297	1.983492	0.022639
M	1.177367	0.982265	1.411221	0.077331
N	1.197347	1.026196	1.397043	0.022106
Risk score	1.119068	1.065739	1.175066	6.31E-06

Table 3

Univariate and multivariate Cox regression analysis of OS for BLCA patients based on our established PI model and clinical features

Variables	Univariate Cox regression analysis					Multivariate Cox regression analysis
	HR	HR.95L	HR.95H	P value		HR	HR.95L	HR.95H	P value
Age	1.037102	1.018346	1.056203	9.14E-05		1.039269	1.019976	1.058927	5.61E-05
Gender	0.795505	0.547983	1.154833	0.228979		0.780127	0.532055	1.143863	0.203512
Race	1.203104	0.858987	1.685078	0.282066		0.860632	0.605305	1.223659	0.403232
Stage	1.879957	1.490326	2.371453	9.98E-08		1.397834	1.037331	1.883622	0.02775
T	1.731316	1.348288	2.223155	1.69E-05		1.445184	1.05297	1.983492	0.022639
M	1.21501	1.022681	1.443509	0.026756		1.177367	0.982265	1.411221	0.077331
N	1.238577	1.099637	1.395072	0.000424		1.197347	1.026196	1.397043	0.022106
Risk score	1.131252	1.079995	1.184941	1.86E-07		1.119068	1.065739	1.175066	6.31E-06

OS, overall survival; BLCA, bladder urothelial carcinoma; PI, prognostic index; HR, hazard ratio.

Clinical correlation analysis of the 10 genes and the PI

We also performed clinical correlation analysis to study the relationships between the 10 prognostic IRGs, the PI model and the clinicopathologic risk factors. As detailed in , we found that the PI was positively correlated with race, grade, M and N period (P<0.05), and negatively correlated with stage and T period (P<0.05).

Table 4

Clinical correlation analysis between the 10 prognostic IRGs, our established PI and clinical features

ID	Age	Gender	Race	Grade	Stage	T	M	N
TAP1	0.766 (0.444)	0.044 (0.965)	21.259 (2.419e-05)	3.834 (7.88e-04)	1.405 (0.162)	1.066 (0.288)	12.221 (0.002)	1.683 (0.431)
RBP7	−0.63 (0.529)	1.835 (0.070)	10.61 (0.005)	6.223 (2.682e-09)	−3.894 (1.184e-04)	−3.947 (9.727e-05)	35.454 (2.001e-08)	12.519 (0.002)
STAT1	0.874 (0.383)	0.684 (0.495)	23.347 (8.518e-06)	2.142 (0.045)	0.334 (0.739)	0.018 (0.986)	6.09 (0.048)	0.416 (0.812)
PDGFRA	0.209 (0.835)	1.988 (0.049)	9.6 (0.008)	5.512 (1.445e-06)	−4.07 (5.901e-05)	−4.442 (1.2e-05)	18.352 (1.035e-04)	12.091 (0.002)
AHNAK	1.059 (0.291)	0.659 (0.511)	18.639 (8.964e-05)	1.949 (0.063)	−2.967 (0.003)	−2.607 (0.010)	3.282 (0.194)	10.461 (0.005)
OLR1	−0.07 (0.944)	−0.077 (0.938)	27.569 (1.032e-06)	2.747 (0.010)	−2.238 (0.026)	−2.635 (0.009)	10.368 (0.006)	5.318 (0.070)
RAC3	−1.518 (0.130)	0.404 (0.687)	1.732 (0.421)	5.217 (1.806e-06)	−2.417 (0.016)	−1.662 (0.098)	9.431 (0.009)	1.697 (0.428)
EDNRA	−0.562 (0.574)	1.336 (0.184)	23.088 (9.694e-06)	9.051 (2.364e-12)	−4.686 (4.405e-06)	−4.396 (1.536e-05)	25.497 (2.906e-06)	12.358 (0.002)
IGF1	−1.466 (0.144)	1.589 (0.116)	20.381 (3.753e-05)	4.635 (6.224e-06)	−4.023 (7.308e-05)	−4.052 (6.622e-05)	22.39 (1.374e-05)	8.091 (0.017)
SH3BP2	1.106 (0.270)	−3.395 (8.436e-04)	5.688 (0.058)	−1.061 (0.302)	3.191 (0.002)	1.964 (0.051)	8.717 (0.013)	6.324 (0.042)
Risk score	0.766 (0.445)	0.135 (0.893)	21.692 (1.948e-05)	2.614 (0.009)	−2.376 (0.018)	−2.137 (0.034)	15.52 (4.264e-04)	14.002 (9.108e-04)

IRG, immune-related genes; PI, prognostic index.

Related small molecule drugs screening

We screened out candidate small molecule drugs of BLAC based on differently expressed IRGs by performing CMap database analysis. By analyzing consistent differently expressed probesets between BLAC samples and adjacent normal samples, we identified the related small molecule drugs with highly significant correlations. 10 small molecule drugs were screened out by the number of instances (n>2) and P value (<0.05). They were shown in . Among these 10 small molecule drugs, ketorolac, fluorocurarine, AR-A014418, lycorine, mebeverine, hydroquinine, and cefamandole showed negative correlation, and dexpanthenol, W-13, and benzbromarone showed positive correlation. The 7 negatively related small molecules have the potential to treat BLAC, so as the antagonists of the 3 positively related.

Table 5

Results of CMap analysis

Rank	cmap name	Mean	n	Enrichment	P	Specificity	Percent non-null
1	Ketorolac	−0.487	4	−0.644	0.03811	0.0567	75
2	Fluorocurarine	−0.439	4	−0.706	0.0154	0.0052	75
3	AR-A014418	−0.417	3	−0.823	0.01114	0.0296	66
4	Lycorine	−0.406	5	−0.755	0.00162	0.08	60
5	Mebeverine	−0.328	4	−0.79	0.00402	0.012	50
6	Hydroquinine	−0.313	4	−0.64	0.04032	0.0426	50
7	Cefamandole	−0.31	4	−0.629	0.04645	0.1083	50
8	Dexpanthenol	0.314	4	0.628	0.0479	0.0073	50
9	W-13	0.642	2	0.943	0.00606	0	100
10	Benzbromarone	0.723	3	0.964	0.00006	0	100

Discussion

As one of the most lethal cancers worldwide (1,14,15), BC has been paid much attention to all the time. In scientific studies, increasing efforts have been put to research the specific mechanisms of BC initiation, progression and prognosis, for the diagnosis, treatment and prognostic surveillance of BC. Many kinds of urinary and blood biomarkers (16-19) were researched for the detection of BC, some of which showed high sensitivity and specificity and may have clinical application value in the future. As reported, long noncoding RNA LNMAT1 and BLACAT2 could promote BC lymphatic metastasis (20,21). And circular RNAs like FNDC3B, ITCH, BCRC-3 and ACVR2A were reported to play important roles in suppressing BC progression (22-25). However, due to high recurrence rate, the long-term efficacy of BC treatment was always unsatisfactory in actual clinical work. Previous research results were not well translated into clinical applications. Hence, it was imperative to develop novel approach to evaluate the prognosis of BC patients, to provide guidance for clinical BC surveillance. In the current study, we concentrated on developing IRGs-based prognostic signature to investigate whether IRG could function as prognostic biomarkers in BC, by analyzing transcriptome data downloaded from TCGA.

In the present study, we first extracted differently-expressed IRGs between BC tissue samples and normal bladder tissues. To identify differently-expressed IRGs with prognostic values, we performed univariate Cox regression analysis and confirmed 57 IRGs. Then, LASSO penalized Cox regression analysis was conducted and 10 IRGs (TAP1, RBP7, STAT1, PDGFRA, AHNAK, OLR1, RAC3, EDNRA, IGF1 and SH3BP2) were identified to develop the PI. Therein, TAP1, a component of the antigen-presenting machinery, was involved in tumor immune escape. It was reported that down-regulation of TAP1 could elicit immune evasion and predict poor prognosis (26,27). Elmasry et al. speculated that in colon cancer, RBP7 functioned as a promotor of tumor migration and invasion through EMT, and indicated poor cancer specific survival (28). As one member of signal transducer and activator of transcription (STAT) family proteins, STAT1 played important roles in tumor microenvironment (29). Meyer Zu Horste et al. demonstrated that STAT1 participated in the process of TH17 differentiation into TH1 (30). PDGFRA mutations were common in gastric, omental and mesenteric tumor (31), and much effort had been made to promote treatments targeting PDGFRA (32-34). AHNAK is a giant structural nucleoprotein involved in diverse biological process including muscle membrane repair, blood-brain barrier, and cell proliferation and migration (35-38). As to OLR1, a lectin-like scavenger receptor, it was mainly expressed in vascular cells and vasculature-rich organs (39), and could function as inflammation promotor in atherogenesis (40). Coincidently, RAC3 could promote ox-LDL induced endothelial dysfunction (41). RAC3 also had correlation with autophagy (41-43), cell migration and invasion (44-46). Much evidence revealed the importance of EDNRA in cardiovascular and cerebrovascular diseases (47-49). As a research hotspot, IGF1 had been receiving much attention and its roles in endocrine disorders, atherosclerosis and cancer promoted development of disease diagnosis and treatment (50-52). SH3BP2 was an adapter protein involved in adaptive and innate immune response signaling, and mainly participated in inflammatory bone loss (53-55). Subsequently, univariate and multivariate Cox regression analysis were carried out to shed light on that this PI model could be an independent predictor. We also analyzed the relationship between this PI model and some clinicopathologic risk factors to confirm their clinical relevance.

Besides, in order to preliminarily explore the relationship between the 57 differentially expressed IRGs and TFs, we constructed a network for them. We first extracted 77 differentially expressed TFs by analyzing the RNA-seq data described above. Then the TF-IRG interaction network was conducted and it could provide a basis for further research. On the one hand, this network analysis was to explore the interactions between TFs and IRGs. On the other hand, it also provided guidance and assistance for further basic researches to seek the upstream regulators of differently expressed IRGs. By analyzing CMap database, we screened out several small molecule drugs with potential therapeutic efficacy to identify candidate small molecule of BLAC. Ketorolac, as a member of non-steroidal anti-inflammatory drugs (NSAIDs) family, was mainly used for inflammation and pain treatment after surgery. In a retrospective analysis of two series of 827 and 1,007 breast cancer patients, Desmedt et al. found that the intra-operative administration of ketorolac was statistically significantly associated with a reduction of distant recurrences (56). AR-A014418, a glycogen synthase kinase-3 inhibitor, could inhibit cell proliferation and induce cell apoptosis in cancers including pancreatic cancer (57) and neuroblastoma (58,59). It also could sensitize pancreatic cancer cells to gemcitabine (60). Wang et al. demonstrated that Lycorine, an alkaloid extracted from Amaryllidaceae genera, could induce apoptosis of BC T24 cells (61). It was reported that anti-calmodulin agent W-13 had the ability to suppress breast cancer cell growth (62) and colony formation (63). As to benzbromarone, it was found to have cytotoxic effects in in human hepatocarcinoma FLC4 cells (64). To our best knowledge, there have no cancer treatment researches of the fluorocurarine, mebeverine, hydroquinine, cefamandole, and dexpanthenol. Overall, our study provided the feasibility of some potential biomarkers or molecular targets in the treatment of BC.

In recent decades, much attention has been paid to researches concentrating on immune in cancer biological and clinical research. T-regulatory (Treg) cells were found to be enriched in BC, and CD4+CD25+ T cells were elevated in patient blood compared with healthy individuals (65). Treg cells were showed to exert the ability to suppress essential antitumor responses (66). Related researches demonstrated that Th17 cells were negatively correlated with Treg cells in BC (67). IL-2 is primarily produced by T-lymphocytes. It was reported that IL-2 could promote the conversion of Treg cells to TH17 cells (68-70). BCG immunotherapy was an important component of BC treatment. Intravesical BCG instillation activated urothelial cells and antigen-presenting cells (APCs), leading to the production of cytokines and chemokines, to further recruit immune cells (71). Cytokines and chemokines like IL-2, IL-6, IL-8, GM-CSF, CC-chemokine ligand 2 (CCL2), and CCL3 could be found in the urine of patients who has underwent BCG treatment (68,72,73). It was reported that IL-2 had correlation with BC recurrence, and along with IL-6, the two cytokines were suggested to be used to monitor BCG treatment effect (74).

The strength of this article was that it was the first time for us to perform a systematic analysis of the roles of immunity in BLCA with a robust statistical approach. A novel PI was successfully established and carefully evaluated. Moreover, networks between IRGs and TFs were also constructed for future basic research. Our study has some limitations. Firstly, we studied only IRGs and other genes associated with OS of BC were not enrolled. Secondly, this study is retrospective and our results should be further validated in prospective investigations.

Conclusions

Taken together, our study successfully developed a novel PI model to predict the OS of BC and it was carefully evaluated. Moreover, networks between TFs and IRGs were also displayed to seek its upstream regulators for future researches. Besides, we also screened out 10 small molecules drugs with the potential to treat BC, based on differently expressed IRGs. Our study highlights the importance of IRGs in clinical application. Although our established PI model stills needs to be further validated by more prospective studies, it was promising to assist clinicians in clinical BC surveillance.

The article’s supplementary files as