Identification of a dysregulated CircRNA-associated gene signature for predicting prognosis, immune landscape, and drug candidates in bladder cancer

Abstract

Increasing evidences have demonstrated that circular RNA (circRNAs) plays a an essential regulatory role in initiation, progression and immunotherapy resistance of various cancers. However, circRNAs have rarely been studied in bladder cancer (BCa). The purpose of this research is to explore new circRNAs and their potential mechanisms in BCa. A novel ceRNA-regulated network, including 87 differentially expressed circRNAs (DE-circRNAs), 126 DE-miRNAs, and 217 DE-mRNAs was constructed to better understanding the biological processes using Cytoscape 3.7.1 based on our previously high-throughput circRNA sequencing and five GEO datasets. Subsequently, five randomly selected circRNAs (upregulated circ_0001681; downregulated circ_0000643, circ_0001798, circ_0006117 and circ_0067900) in 20 pairs of BCa and paracancerous tissues were confirmed using qRT-PCR. Functional analysis results determined that 772 GO functions and 32 KEGG pathways were enriched in the ceRNA network. Ten genes (PFKFB4, EDNRA, GSN, GAS1, PAPPA, DTL, TGFBI, PRSS8, RGS1 and TCF4) were selected for signature construction among the ceRNA network. The Human Protein Atlas (HPA) expression of these genes were consistent with the above sequencing data. Notably, the model was validated in multiple external datasets (GSE13507, GSE31684, GSE48075, IMvigor210 and GSE32894). The immune-infiltration was evaluated by 7 published algorithms (i.e., TIMER, CIBERSORT, CIBERSORT-ABS, QUANTISEQ, MCPCOUNTER, XCELL and EPIC). Next, Correlations between riskscore or risk groups and clinicopathological data, overall survival, recognized immunoregulatory cells or common chemotherapeutic agents of BCa patients were performed using wilcox rank test, chi-square test, cox regression and spearman's correlation analysis; and, these results are significant. According to R package "GSVA" and "clusterProfiler", the most significantly enriched HALLMARK and KEGG pathway was separately the 'Epithelial Mesenchymal Transition' and 'Ecm Receptor Interaction' in the high- vs. low-risk group. Additionally, the functional experiments in vitro also revealed that the overexpression of has_circ_0067900 significantly impaired the proliferation, migration, and invasion capacities of BCa cells. Collectively, the results of the current study provide a novel landscape of circRNA-associated ceRNA-regulated network in BCa. The ceRNA-associated gene model which was constructed presented a high predictive performance for the prognosis, immunotherapeutic responsiveness, and chemotherapeutic sensitivity of BCa. And, has_circ_0067900 was originally proposed as tumor suppressor for patients with BCa.

Keywords: CircRNA-related prognostic signature; bladder cancer; chemotherapy; immune infiltration; immunotherapy.

Figure 1

(A) Volcano map representing the up-regulated and down-regulated circRNAs in BCa vs. para-carcinoma tissues (log2(FC) ≥2 or ≤ -2; P< 0.05). (B–D) Based on GSE11224, GSE113786 and GSE113740 datasets, top 50 up-regulated and top 50 down-regulated miRNAs in BCa tissues vs. controls were shown. (E) Venn diagrams revealed the number of overlapping up-regulated (left) or down-regulated (right) miRNAs among GSE11224, GSE113786 and GSE113740. (F, G) Based on two datasets (GSE13507 and GSE37815), top 50 up-regulated and top 50 down-regulated miRNAs in BCa tissues vs. controls were exhibited. (H) Venn diagrams represent overlapping up or down-regulated genes between GSE13507 and GSE37815 datasets. circRNA, circular RNA; miRNA, microRNA; BCa, bladder cancer; FC, fold change.

Figure 2

Construction analysis of the ceRNA network of the DE-circRNA, DE-miRNA, and DE-mRNA identified in BCa tissues compared with their adjacent non-cancerous tissues. (A) Intersection of DE-circRNA-targeting miRNAs and common DE-miRNAs of GSE11224, GSE113786 and GSE113740 datasets. (B) Conjoint analysis of the A-mentioned 126-intersection-DE-miRNA binding mRNA targets and DE-mRNA of GSE13507 and GSE37815 datasets. (C) The circRNA-miRNA-mRNA ceRNA network was established and visualized by Cytoscape 3.7.1 that based on the above data. ceRNA, competing endogenous RNA; circRNA, circular RNA; miRNA, microRNA; differentially expressed-, DE-.

Figure 3

Functional enrichment and Lasso-penalized Cox regression analyses of these DE-mRNAs in the ceRNA network. (A, B) The GO function (BP, CC and MF) and KEGG pathway annotation analyses of DE-mRNAs in the ceRNA network. (C, D) Lasso regression analysis selected prognosis–related DE-mRNAs in the ceRNA network. ceRNA, competing endogenous RNA; BP, biological process; MF, molecular function; CC, cellular components; DE-, differentially expressed-.

Figure 4

Verification of the ceRNA-associated modeled genes expression in BCa and normal bladder tissue using the HPA database.

Figure 5

Construction and validation of a DE-mRNAs of the ceRNA network Signature. (A) The Kaplan-Meier, ROC and calibration curves of the 1-, 3-, and 5-year OS rate prediction for patients with BCa based on the risk score for the TCGA_BLCA training cohort. (B–F) Kaplan-Meier, ROC and calibration curves of the risk score for the 1-, 3-, and 5-year OS rate prediction for the testing set, including GSE13507, GSE31684, GSE48075, IMvigor210 and GSE32894. ceRNA, competing endogenous RNA; BCa, bladder cancer; DE-, differentially expressed-.

Figure 6

Identification of prognostic indicators for bladder cancer in training set. (A) The risk score curve, survival status and ceRNA-associated prognosis model gene expression heatmap were displayed from top to bottom in the TCGA_BLCA training cohort. The abscissa axis of these graphs were ranked by the risk score value. (B, C) Correlation between riskscore group and clinicopathological data of BCa patients via Wilcox rank test or Chi-square test. (D) The risk score was independent risk factors for bladder cancer. (E) Construction of a ceRNA-associated gene signature combined with clinical features nomogram. (F) Calibration curve of nomogram. (G) A multi-index ROC curve demonstrated the good discriminative abilities of the ceRNA-associated gene signature or nomograms. (H) DCA was applied to render clinical validity to the constructed gene signature or nomograms. ceRNA, competing endogenous RNA; BCa, bladder cancer; ROC, receiver operating characteristic; DCA, decision curve analysis. *P < 0.05, **P < 0.01 and ***P < 0.001.

Figure 7

Associations between ceRNA-associated gene signature and immune-cell infiltration, immune checkpoint immunotherapies. (A) There were distinct differences in many immune cells infiltration between the high- and low-risk groups. (B) The TIDE algorithm analysis and the expression of partial immune checkpoint genes was used to predict the immune response to immune checkpoint therapy in BCa patients based on riskscore. (C, D) Correlation analysis between infiltrating immune cells abundance from 7 immune-infiltration algorithm and the riskscore. ceRNA, competing endogenous RNA; TIDE, tumor immune dysfunction and exclusion; BCa, bladder cancer. *P < 0.05, **P < 0.01 and ***P < 0.001.

Figure 8

The ceRNA-associated gene signature is an important indicator of immune cell infiltration and immunotherapeutic effect. (A) The immune cell infiltration difference from 7 immune-infiltration algorithm in high-risk groups compared to low-risk groups. (B) The association between the risk score and immunotherapy outcome by wilcox test. (C) Correlations among the riskscore group and immunotherapy outcome variables. (D) Using TIDE algorithm analysis and the expression of partial immune checkpoint genes, immune checkpoint blockade therapy response based on riskscore was predicted. (E, F) Spearman association analysis between riskscore and immune cell infiltration was implemented. ceRNA, competing endogenous RNA; TIDE, tumor immune dysfunction and exclusion. *P < 0.05, **P < 0.01 and ***P < 0.001.

Figure 9

The underlying molecular mechanisms and chemotherapeutics. (A) Differences in pathway activities scored by GSVA of HALLMARK gene sets between high-risk and low-risk group. Navy blue, significantly upregulated enrichment; green, significantly downregulated enrichment; grey, no enrichment. (B) Most significant enriched 20 GSEA pathways in high- vs. low-risk group were displayed separately based on KEGG and HALLMARK gene sets. Red, significant positive enrichment; green, significant negative enrichment. (C) Wilcox group analysis and (D) spearman correlation analysis indicated that the ceRNA-associated gene model is robust to drug sensitivity of Cisplatin, Docetaxel, Paclitaxel and Vinblastine. GSVA, Geneset variation analysis; GSEA, gene set enrichment analysis; ceRNA, competing endogenous RNA.

Figure 10

circ_0067900 overexpression exerted a tumor suppressive effect in BCa cells. (A) Five circRNAs were randomly chosen to verify the circRNA-sequencing results in 20 pairs of samples by quantitative real-time PCR. (B) Lentiviral transduction efficiencies were evaluated using RT-qPCR in T24 and 253J-BV cells. (C) CCK-8 proliferation assay. (D) OE-circ_0067900 inhibited monoclonal formation ability of BCa cells. (E) Scratch assay of T24 and 253J-BV cells. Experiments were terminated after scratch for 24 h. (F) Migration and invasion were evaluated using a transwell assay without and with Matrigel, respectively. Data are showed as means ± SD; n = 3. SD, standard deviation; OE-, Overexpression; NC, negative control vector; CCK-8, cell counting kit-8; OD, optical density. *P < 0.05, **P < 0.01, ***P < 0.001 and ****P < 0.0001.