Establishment and analysis of a novel miRNA prognostic risk model for bladder cancer based on TCGA database

Abstract

Bladder cancer (BCa) is the most common urological tumor, and most BCas are diagnosed at an advanced stage due to the lack of reliable diagnostic markers. In this study, transcriptomic and clinical data from BCa patients in the TCGA database were analyzed and all patients were randomly divided into a training group and a validation group. Univariate Cox regression analysis was performed to screen miRNAs significantly associated with BCa prognosis, followed by multivariate Cox regression analysis in the training group to establish a prognostic risk model for multiple miRNAs. The model was validated in the validation group and in the whole group and assessed using Norman plots. KEGG pathway and GO enrichment analyses were performed on possible target genes, and PPI networks were constructed. Finally, a prognostic model for 5-miRNA was successfully constructed. Multivariate Cox analysis showed that it had good predictive properties and could be used as a novel biomarker to predict prognosis and survival of BCa patients and may exert its biological function through IGF1.

KEYWORDS:

• MiRNA
• bladder cancer
• TCGA
• prognostic model
• IGF1

Introduction

Bladder cancer (BCa) is the genitourinary system tumor with the highest incidence in China, and its recurrence rate and mortality rate remain high (Abdollah et al. 2013; Taylor et al. 2020). Although a variety of biomarkers have been found to predict prognosis, their reliability is still controversial (Shi et al. 2019). Therefore, there is an urgent need to screen new and reliable biomarkers or therapeutic targets for the diagnosis and prognosis of BCa.

MiRNA is a class of non-coding single-stranded RNA molecules encoded by endogenous genes with a length of about 22 nucleotides (Lou et al. 2019). It regulates the expression of protein-coding genes by inhibiting translation or cutting RNA transcripts. There is growing evidence that miRNA is associated with the occurrence, progression and metastasis of BCa (Ali Syeda et al. 2020; Taheri et al. 2020). At present, a single miRNA has been used to build a prognostic model of BCa miRNA, but there are still some shortcomings, such as inaccurate results and failure to distinguish the ‘-5p’ and ‘-3p’ ends of mature miRNA (Yin et al. 2019).

In this study, we constructed, validated and evaluated a novel 5-miRNA prognostic model through the TCGA database, distinguished the ‘-5p’ and ‘-3p’ ends of miRNA, predicted the survival of patients with BCa, and functional enrichment revealed potential biological functions and signaling pathways related to BCa, which enhanced our understanding of the molecular mechanism of the model.

Materials and methods

Data download and processing

Download miRNA, mRNA sequencing data and clinical data of BCa patients from the Cancer Genome Atlas (TCGA) detabase. The expression of miRNA and mRNA in the tumor group was compared with that in the normal group by using the edgeR package of R software (version 3.6.2), and the miRNA and mRNA with corrected p-values < 0.05 and |log₂FC| > 1 were screened (Robinson et al. 2010). Perl software (version 5.31.9) was used to combine the clinical information of patients with survival time ≥ 30 days with differentially expressed miRNA expression profiles.

Construction of prognostic model

Using ‘caret’ package, the samples with complete survival information, survival time ≥ 30 days and differentially expressed miRNA expression profiles were randomly divided into two groups (training group and verification group). Univariate Cox regression analysis was used to screen prognosis-related miRNA in the training group. In order to reduce the number of miRNA with similar expression, multivariate Cox regression analysis of miRNA with p < 0.05was carried out by using packages ‘Coxph’ and ‘direction = both’ (Stel et al. 2011). From the Cox proportional hazard regression model, we use the following formula to calculate the risk score of each patient: miRNA risk score = β_miRNA1 × exp (miRNA1) + β_miRNA2 × exp (miRNA2) + … . +β_miRNAn × exp (miRNAn), where β is the regression coefficient derived from the multivariate Cox proportional hazard regression model, and exp () is the level of the expression miRNA (Sui et al. 2017). According to the median risk score, they are divided into high-risk and low-risk subgroups. Using the ‘survivalROC’ package to calculate the AUC of the 5-year dependent ROC curve to evaluate the predictive ability of the miRNA model (Heagerty et al. 2000). Kaplan–Meier curve and logarithmic rank test were used to evaluate the survival and prognosis of the patients in the training group, the verification group and the whole group.

Independent prognostic ability test of miRNA model

In the training group, the relationship between 5-miRNA model and overall survival was analyzed by univariate Cox regression and clinical variables (including age, sex and clinical stage). Further multivariate Cox regression analysis was conducted to determine whether it could be used as an independent prognostic factor. The independent prognostic factors in the results of multivariate Cox regression analysis were plotted by Norman diagram, and the 1-, 3- and 5-year survival rates were predicted, and the 5-year survival rates were analyzed by the calibration curve.

MiRNA target gene and its potential function

The miRNA prediction database was downloaded from three miRNA target gene prediction websites of miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/), targetScan (http://www.targetscan.org) and miRDB (http://www.mirdb.org/). The target genes of at least two miRNA models covered in the database were found by using Perl language, and the relationship between miRNA and these target genes was drawn by Venn diagram and Cytoscape software (version 3.6.1). In order to determine whether these miRNA target genes may be involved in the progression of BCa, we established the intersection of these target genes and differentially expressed genes in BCa. All the crossover genes were analyzed by using the R software packages ‘clusterProfiler’ and ‘org.Hs.eg.db’, as well as KEGG signal pathway and GO enrichment analysis (Yu et al. 2012).

PPI network construction

The PPI network is constructed by mining the relationship between the target genes from the STRING database (https://string-db.org/). The reliability of the parameters set is 0.700 (Szklarczyk et al. 2017). The top 10 central genes were determined according to Cytoscape 3.6.1 and its plug-ins (cytoHubba ranking). At the same time, Kaplan–Meier method was used to test whether the cross gene was related to survival, and the logarithmic rank test was less than 0.001.

Results

Differentially expressed miRNA and mRNA in BCa

According to the screening conditions, a total of 473 differentially expressed miRNA were found, including 379 up-regulated and 94 down-regulated, and a total of 6850 differentially expressed mRNA, of which 4478 were up-regulated and 2372 down-regulated, were shown in the volcano plot (Figure 1).

Figure 1. Volcano plot of DemRNAs (A). Volcano plot of DemiRNAs (B). Red dots represent up-regulation and green dots represent down-regulation.

Construction of 5-miRNA prognostic model

The whole group with miRNA mature expression profile (N = 160) was randomly divided into training group (N = 80) and verification group (N = 80). In the training group, univariate Cox regression analysis showed that a total of 19 miRNA were found to be related to the overall survival time of the patients (p < 0.05). Then the miRNA with similar expression was screened out by multivariate Cox regression analysis (Table 1), and a predictive miRNA model was established according to the sum of the product of each miRNA expression and its coefficients in multivariate Cox regression. The model was established as follows: prognostic model risk score = (0.3665 × hsa-miR-338-5p) + (0.3917 × hsa-miR-218-1-3p) + (−0.5233× hsa-miR-17-3p) + (0.7595 × hsa-miR-6728-5p) +(−0.1777 × hsa-miR-10a-3p).

Table 1. 5-miRNA expression.

	Single factor COX				Multivariable COX
id	HR	HR.95L	HR.95H	P value	Coef	HR	HR.95L	HR.95H	P value
hsa-miR-338-5p	1.373	1.042	1.808	0.024	0.366	1.443	1.026	2.029	0.035
hsa-miR-218-1-3p	1.524	1.118	2.079	0.008	0.392	1.480	0.995	2.199	0.053
hsa-miR-17-3p	0.607	0.370	0.995	0.048	−0.523	0.593	0.364	0.964	0.035
hsa-miR-6728-5p	1.542	1.114	2.138	0.009	0.759	2.137	1.379	3.312	0.000
hsa-miR-10a-3p	0.726	0.579	0.911	0.006	−0.178	0.837	0.661	1.060	0.140

Survival prediction of training group, verification group and whole group

According to the median risk score of the training group, each group was divided into a high-risk and low-risk subgroups. Kaplan–Meier curve showed that the overall survival rate in the high-risk subgroup was significantly lower than that in the low-risk subgroup (Figure 2(A–C)). In the training group, the verification group and the whole group, the 5-year ROC-dependent AUC for 5-miRNA prognostic model reached 0.835, 0.769 and 0.808 respectively (Figure 2(D–F)). In addition, in the three groups, the mortality rate of patients with high-risk score was significantly higher than that of patients with low-risk score (Figure 2(G–I)).

Figure 2. Survival time, ROC and survival points of patients in training group, verification group and whole group. The Kaplan–Meier curve in the training group (A), verification group (B) and the whole group (C). The ROC curve in the training group (D), verification group (E) and the whole group (F). Patients survival in the training group (G), verification group (H) and the whole group (I).

Independence analysis of 5-miRNA prognostic model

Univariate Cox regression analysis showed that the 5-miRNA prognostic risk score was significantly correlated with the overall survival rate (Figure 3(A)). Multivariate Cox regression analysis showed that 5-miRNA prognostic model was still associated with the overall survival rate considering other routine clinical factors, such as age, sex, clinical stage, T stage, lymph node status and distant metastasis (Figure 3(B)). In addition, the T stage is also an independent prognostic factor. The nomogram used to predict the OS of BCa patients is shown in Figure 4(A). Nomogram for predicting OS was based on the following two independent prognostic factors: T stage and 5-miRNA prognostic risk score. According to the sum of the distribution points of each factor in the nomogram, the higher the total score is, the worse the prognosis is.

Figure 3. Univariate and multivariate Cox regression analyses for the prognostication of overall survival in the training group. (A) univariate Cox analysis; (B) multivariate Cox analysis.

Figure 4. The nomogram used to predict the OS of BCa patients. (A) To estimate the 1-year, 3-year, and 5-year survival rate of BCa, the patient’s corresponding values were marked on each axis and the points of all variables were summed. Next, mark the sum on the total point axis and draw a straight line perpendicular to the probability axis. (B)The x-axis shows the predicted probability of 5-year survival rate of BCa, and the y-axis shows the observed actual 5-year survival rate of BCa.

Target gene of 5-miRNA prognostic model

The target genes regulated by five kinds of miRNA were predicted in the three databases of miRDB, miRTarBase and TargetScan. The overlapping target genes were identified, and a total of 3094 target genes in 5 miRNA were predicted, which were shown by the Venn map (Figure 5). In order to clarify whether the target genes of these miRNA may be involved in the progression of BCa, the mRNA predicted by the down-regulated miRNA (hsa-miR-17-3p, hsa-miR-10a-3p) in the database was interacted with the rising differentially expressed mRNA, and the mRNA predicted by the up-regulated miRNA (hsa-miR-338-5p, hsa-miR-218-1-3p, hsa-miR-6728-5p) was interacted with the down-regulated mRNA. A total of 307 target genes were included, including 134 up-regulated genes and 173 down-regulated genes. The subnetwork between 5 miRNA and its 307 target genes is shown in Figure 6.

Figure 5. Five Venn diagram of the target gene of miRNA. (A) hsa-miR-10a-3p (B) hsa-miR-338-5p (C) hsa-miR-17-3p (D) hsa-miR-218-1Mur3p and (E) hsa-miR-6728-5p.

Figure 6. MiRNA adjusts the subnetwork diagram of mRNA. Triangles represent miRNA and the circle represents mRNA. Red indicates up-regulation and green means down-regulation.

Functional enrichment analysis of BCa-related target genes

The result of GO annotation on BCa-related target genes was 979. The top 30 functions in GO results are shown in biological processes (BP), cellular components (CC) and molecular functions (MF) in the dot map (Figure 7(A–C)). The results of the KEGG pathway of target genes related to BCa are mainly concentrated in the cGMP-PKG signal pathway, and the ‘pathway-gene network’ and ‘pathway-pathway network’ are also shown in Figure 7(D–F).

Figure 7. Functional enrichment analysis of BCa-related target genes. The (A) BP, (B) CC, (C) MF, (D) KEGG signal pathway dot map shows the gene count dot map shows the gene count, (E) KEGG signal pathway cnetplot shows the ‘pathway-gene’ network the emapplot of the and (F) KEGG signal pathway shows the ‘pathway-pathway’ network.

PPI network and target genes related to survival

About 291 of the 307 target genes were filtered into the target gene PPI network, and 10 central genes (CHEK1,ESR1, CDCA8, MCM10, TTK, IGF1, WDHD1, ABL1, RRM2, CENPF) are filtered based on Cytoscape 3.6.1 and its plug-ins (cytoHubba ranking) respectively (Figure 8). Kaplan–Meier method showed that the expression of 23 genes (ADAMTS15, AKAP13, AKT3, BNC2, BOC, CHRDL1, COLEC12, DCLK1, DOK6, EBF1, ERC1, FGF2, GHR, IGF1, KCNMA1, KLHL4, MAP1A, MAP1B, MEF2A, PDZRN4, TSHZ3, UNC5C, ZFHX4) was negatively correlated with survival (Figure 9).

Figure 8. Protein–protein interaction network of target genes.

Figure 9. Kaplan–Meier method showed that the expression of 23 genes was negatively correlated with survival of BCa.

Discussion

BCa is the first disease of malignant tumor of the urinary system, which high recurrence rate and distant metastasis rate bring huge economic burden to patients and society (Leal et al. 2016). Therefore, there is an urgent need to develop new markers to reduce the burden of the disease on patients and the health care system. Current studies have shown that the expression of miRNA in urine and blood samples of patients with BCa is different from that of healthy people, which may be of diagnostic value (Blanca et al. 2017). There are several single miRNA can be used as potential prognostic indicators of BCa, such as miR-125a-5p, miR-199, miR-100 and so on (Motawi et al. 2016; Sakaguchi et al. 2018; Zhang et al. 2018). However, a large number of studies have shown that multiple miRNA models have greater advantages than a single miRNA. Prior to our study, we have established a number of prognostic models constructed by multi-miRNA, including BCa (Zhou et al. 2015; Inamoto et al. 2018). Compared with previous studies, our research updated the research methods and expanded the sample size. Most importantly, we used to distinguish between ‘-5p’ and ‘-3p’ side miRNA to verify the model.

In this study, we downloaded mature miRNA expression profiles and BC clinical information of corresponding patients from the TCGA database. All patients were randomly divided into a training group and verification group, and then five miRNA prognostic models were constructed by univariate Cox regression and stepwise multivariate Cox regression. At the same time, the 5-miRNA model is verified in the verification group and the whole group. Based on the median risk score, the Kaplan–Meier curve showed that the overall survival of the high-risk subgroup was significantly worse than that of the low-risk group. Evaluating the excessive survival of 5-miRNA signals in the three groups of patients by ROC curve showed better predictive ability. Univariate Cox regression and multivariate Cox regression analysis also indicated that the 5-miRNA model was associated with overall survival, taking into account other routine clinical factors in patients with BCa. Lu found that miR-17-3p regulates the proliferation and survival of colon cancer through targeted regulation of Par4 (Lu et al. 2018). Lazzarini found that miR-10a-3p was associated with uterine leiomyoma (Lazzarini et al. 2020). Zhou found that miR-338-5p can be used as a tumor marker of retinoblastoma and is associated with drug resistance in colorectal cancer (Xu et al. 2019; Zhou and Li 2019). Yu found that miR-218-1-3p was associated with lung adenocarcinoma (Yu et al. 2019). However, there is no report on the mechanism of miR-6728-5p in the tumor. Therefore, more experiments on miR-6728-5p are needed in the future, especiallyin BCa.

In order to further understand the regulatory mechanism of five miRNA models in BCa, the target genes of five miRNA in the model were predicted by three target gene prediction databases. At the same time, based on the study of BCa, we obtained the intersection of these miRNA target genes and differentially expressed genes from the TCGA database and analyzed the functional enrichment of these genes. GO annotation of target genes is mainly related to axonal development, axonogenesis and stem cell differentiation, synaptic membrane, postsynaptic membrane and neuronal body, metal ion transporter activity, transcriptional activator activity and DNA binding, ion channel binding. The signal pathway of the target gene is mainly concentrated in the cGMP-PKG signal pathway. Gong found that the activation of the GMP/PKG pathway contributes to the proliferation and invasion of cervical cancer cells (Gong et al. 2019). These signal pathways show their effects on the tumor in varying degrees, and these two signal pathways are only the tip of the iceberg of the target genes in the signal pathway, which suggests that the miRNA prognostic model we constructed may be involved in the regulation of tumor signalpathway.

In order to find the key nodes of the miRNA model regulating BCa, 10 hub genes (CHEK1, ESR1, CDCA8, MCM10, TTK, IGF1, WDHD1, ABL1, RRM2, CENPF) were screened according to Cytoscape 3.6.1 and its plug-ins. In addition, Kaplan–Meier method showed that 23 genes (ADAMTS15, AKAP13, AKT3, BNC2, BOC, CHRDL1, COLEC12, DCLK1, DOK6, EBF1, ERC1, FGF2, GHR, IGF1, KCNMA1, KLHL4, MAP1A, MAP1B, MEF2A, PDZRN4, TSHZ3, UNC5C, ZFHX4) were negatively correlated with survival and prognosis. Insulin-like growth factor 1 (IGF1) is not only a key gene in the PPI network but also related to the survival of patients. IGF1 is a peptide hormone that acts by binding and activating the IGF1 receptor, which is a widely distributed tyrosine kinase protein that promotes the balanced growth of circulating IGF1 among multiple tissues and organs (Maki 2010). The available verification data show that miRNA can target multiple molecules in the IGF1R pathway, including IGF1R, including miR-7, Let-7, miR-493, miR-155, miR-99, miR-145, etc. (Kumar et al. 2018). Therefore, miRNA signal may affect the survival and prognosis of patients with BCa and the progression of BCa by regulating IGF1.

The shortcoming of this study is that the model was constructed and validated using only data from the database, lacking validation from biological studies. However, the models and predictions we have developed have potential value for future research and require further validation and study.

Conclusions

To sum up, our study not only constructed a new prognostic model of the 5-miRNA model through miRNA mature expression profile but also verified and evaluated the predictive ability of the model by grouping. In addition, the potential function is inferred by predicting the target genes of the model, which deepens our understanding of the occurrence and development of the BCa. However, this study is only based on the use of bioinformatics TCGA database, lack of other databases and a large number of experiments to verify the feasibility of this prognostic model has some limitations, and further research is needed to verify this conclusion.

Data availability statement

The data that support the findings of this study are openly available in ‘figshare’ at https://doi.org/10.6084/m9.figshare.13622558.v1.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This study was supported by the Project of State Key Laboratory of Radiation Medicine and Protection, Soochow University [Grant Number GZK1202017], Gusu Medical Talent Foundation [Grant Number GSWS2020021].