Development and Validation of A Seven-microRNA Signature as a Novel Prognostic Biomarker for Bladder Cancer

doi:10.21203/rs.3.rs-70359/v1

Download PDF

Primary research

Development and Validation of A Seven-microRNA Signature as a Novel Prognostic Biomarker for Bladder Cancer

https://doi.org/10.21203/rs.3.rs-70359/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background The differential expression of miRNAs has played a significant role in bladder tumors. The aim of our study was to screen new biomarkers .

Methods Through differential analysis of bladder cancer mRNA and miRNA expression data in the TCGA, differential genes and miRNAs were screened. Furthermore, Cox univariate analysis and multifactor analysis were used to establish a prognostic prediction model . The predictive ability of the prognostic model was then verified on the patient. The action mechanism of these miRNAs was analyzed.

Results By the differential analysis and standardization of miRNA expression profiles. Differentially expressed miRNAs were screened, then all the patients were then randomly divided into train group and the test group. 23 miRNAs were revealed , then a Seven-miRNA signature prognostic biomarkers was constituting.Univariate cox regression and multivariate cox regression considering other clinical factors displayed that the seven-miRNA signature could serve as an independent prognostic factor.Target genes of these seven miRNAs were analyzed by KEGG signaling pathway and GO enrichment analysis. .

Conclusion The prognostic model constructed by seven miRNAs has possessed certain degree of sensitivity and specificity for the prediction of the survival of bladder cancer patients, which can be used as a potential new clinical marker for bladder cancer patients.

Cancer Biology

microRNA

bladder cancer

TCGA

prognosis

signature

Bladder cancer (Bca), one of the most commonly-seen malignant tumors of the urogenital system, has ranked the tenth among the most commonly-seen tumor globally, with an estimated number of 549,000 new cases and 200,000 deaths annually[1, 2]. However, most patients (70–75%) with bladder cancer were initially diagnosed as non-muscle invasive bladder cancer (NMIBC), with superficial location. About 70% of them will relapse after treatment[3]. In addition, 25–30% of patients were diagnosed as muscle invasive bladder cancer (MIBC), often accompanied by local or systemic metastases[4, 5]. In recent years, with the development of science and new technology, diagnosis and treatment of bladder cancer have made substantial progress, yet the prognosis still lags behind with poor performance. Currently, the main prognostic factors, namely the staging and pathological grading of tumor lymph node metastasis (TNM), were not sound enough to predict the clinical outcome alone [6] These clinical and pathological risk factors fail to clearly predict the risk of relapse, chemotherapy response, and overall survival (OS) of patients. Therefore, it is of great necessity to include new prognostic biomarkers into the current staging system, and these biomarkers can also be used as therapeutic targets.

miRNA functions as an endogenous small single-stranded non-coding RNA that does not contain an open reading frame (ORF) and it does not encode the protein by itself [7]. Most studies have found that it is about 22 nt in length, usually ranging between 21 and 25 nt. Its expression is characterized by specificity in tissue and time [8, 9], which can regulate its translation and stability through complementary combination with mRNA, thus playing a role of post-transcriptional regulation. Although miRNAs tend to be relatively small in number, more than 1/3 of the cell transcriptome is regulated by miRNA. miRNA targets for mRNA, which regulates gene expression at the post-transcriptional level. Through its combination with the 3'-UTR region, the following effects can be produced: during its period of translation initiation and peptide chain elongation, translation is stagnated and non-functional protein are produced with cleavage of miRNAs, and degradation of stability are caused by the removal of mRNA from poly-A tails. Because some genes regulated by miRNA are closely related to the cell cycle, they are considered to have played an essential role in tumorigenesis which may become potential proto-oncogenes or tumor suppressor genes[10]. In addition, related studies have demonstrated that miRNAs participate in cell growth, proliferation, differentiation, and apoptosis by regulating the expression of genes encoding proteins[11]. Studies have suggested that miRNAs function as potential biomarkers for early cancer detection and accurate prognosis, as well as targets for more effective treatments[12]. To date, however, only a few reports have described miRNA expression and its relationship with bladder cancer survival.

This study is based on miRNA sequencing data in the TCGA database, and a prognostic prediction model composed of multiple miRNAs is constructed to verify the prognostic predictive value of the prediction model by evaluating the differential expression of miRNAs. At last, combined with transcriptome data, the action mechanism of these miRNAs was analyzed through multi-omics analyses.

Data Download

The transcriptome data of 438 patients with bladder urothelial carcinoma were downloaded from the Cancer Genome Atlas (TCGA) official website (https://portal.gdc.cancer.gov/repository)[13]. The types of data include (transitional cell papilloma and carcinoma, epithelial tumor, adenoma and adenocarcinoma, squamous cell tumor) and miRNA expression profile data. All mature miRNA sequences (mature.fa) in Fasta format were all downloaded from the miRBase website (http://www.mirbase.org/). We combined these two sets of data in the Perl language to obtain the expression profile information of each mature miRNA and its corresponding clinical information (data format: BCR XML). Please check for specific baseline data (Table 1). They were randomly divided into two groups: train group (204 cases) and test group (204 cases).

Table 1

Summary of patient cohort information
Variables	Case	Percentage
GENDER
Male	301	73.78%
Female	107	26.22%
AGE (YEARS)
Range	37–90
Median	68	0.98%
CLINICAL STAGE
Stage I	2	0.49%
Stage II	131	32.11%
Stage III	139	34.07%
Stage IV	134	32.84%
Unknown	2	0.49%
T STAGE
T1	3	0.74%
T2	121	29.66%
T3	193	47.30%
T4	58	14.22%
Tx	1	0.25%
Unknown	32	7.83%
LYMPH NODE STATUS
N0	237	58.09%
N1	46	11.27%
N2	76	18.63%
N3	7	1.72%
Nx	36	8.82%
Unknown	6	1.47%
METASTATIC
M0	195	47.79%
M1	11	2.70%
Mx	199	48.77%

Differential Analysis Of Mrna And Mirna

We adopted the R language edgeR package to normalize and analyze the expression values of mRNA and miRNA sequencing matrices in normal and tumor samples. The screening criteria was set to be | log2FC |> 1 and P < 0.05, aiming at obtaining the expression matrix of mRNA and differentially expressed mRNAs. Afterwards, we combined them with differentially expressed and standardized miRNAs and mRNA expression profiles for subsequencing analysis. Then we used the the “pheatmap” package of the R language package to draw the visualized differential heatmap and volcano map of the differential mRNA and differential miRNA, respectively.

Model Establishment And Validation

We used the "caret" package of version 3.6.1 in R language to randomly divide the samples with complete survival information and differentially expressed miRNA expression profiles into two groups (train group and test group). Then through the survival package of R language, the Coxph function was used to conduct Cox univariate regression analysis on miRNAs in the train group. In order to reduce the number of miRNAs that expressed similarly, a stepwise multivariate Cox regression analysis was performed on miRNAs with p-values <0.005 to construct a prognostic model. In the multivariate Cox regression analysis, we took advantage of the "Coxph" and "direction = both" features of the R survival package. Therefore, based on the sum of the following results, a calculated risk score containing multiple prognostic miRNAs and an extreme coefficient of expression of each miRNA was established. In addition, we tested the hypothesis of proportional hazards in the Cox model. Based on the median risk score grouping, the model was used to evaluate the survival prognosis of each patient in the train group, the test group, and the entire group, using the Kaplan-Meier curve and Log-rank test, ie high-risk group and low-risk group. The predictive power of miRNA models was evaluated by calculating the AUC of the 3-year dependent ROC curve using the "survivalROC" package and its specificity.

Prediction of miRNA Target Genes and Construction of miRNA Regulatory Network

We downloaded miRNA prediction databases from three miRNA target gene prediction websites, including miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/), targetScan (http: // www. Targetscan.org), and miRDB (http (http://www.mirdb.org/). Besides, we used Perl language to find at least 2 target genes regulated by miRNAs in the database. In the meantime, we adopted the “VennDiagram”package of R software to draw a visualized Wayne diagram. We used the cytoscape software to draw a prepared miRNA and mRNA target gene regulatory relationship data to visualize the miRNA regulatory network diagram.

Gene Ontology Enrichment Analysis and Survival Analysis of Genes in MiRNA Regulatory Network

In order to further explore whether these miRNA target genes may be involved in the progression of bladder urothelial carcinoma, we intersected these target genes and differentially expressed genes in bladder urothelial carcinoma and extracted the target genes to be used in bladder urothelial carcinoma which is regulated by miRNAs. All these intersecting genes were analyzed with the R language "clusterProfiler" package and the (KEGG) signaling pathway. Gene ontology (GO) enrichment analysis was further conducted. "Organization.Hs.eg.db" software package set the cut-off criteria as: p adjustment < 0.05 and q value < 0.05, to draw a visualized graph. We used the survival package of R language to conduct batch survival analysis on genes in the regulatory network to visualize the graph.

Construction of Protein Interaction Network Map in MiRNA Regulatory Network and Selection of Core Genes

We used the string website (URL: http://string-db.org/cgi/input.pl) to build the protein interaction network map, set the medium parameter to 0.400, and then downloaded the network relationship map[14]. Then, through the cytoHubba plug-in in the cytoscape software, the central genes were screened based on the degree value and the core gene PPI network map was constructed.

Screening of Differentially Expressed mRNAs and miRNAs

We extracted a total of 7892 mRNAs and 856 miRNAs from the expression profile of urinary bladder urothelial carcinoma of TCGA. Through R language software, we used the edge R package to set the P value to less than 0.05 and log2FC > 1 as differentially expressed mRNA and miRNA. A total of 1111 differentially expressed genes (including 545 up-regulated genes and 566 down-regulated genes) and 473 differentially-expressed miRNAs were screened. (including 380 up-regulated genes and 93 down-regulated genes)). Visualized heatmap and volcano map were drawn by R software R package (Fig. 1)

We randomly divided the expression profile of the entire group with miRNA mature bodies (N = 408) into two groups, the train group (N = 204) and the test group (N = 204). The results showed that in the train group, we found 37 miRNAs related to the overall survival rate of bladder urothelial carcinoma patients based on the screening conditions (p < 0.05). For the reliability of the model, we improved the screening conditions (P < 0.01, log2FC > 1.5) and further screened 7 miRNAs by univariate and multivariate Cox regression analysis. We selected these seven miRNAs for further analysis (Table 2). The Kaplan-Meier method indicated that seven miRNAs include hsa-miR-33b-5p, hsa-miR-33a-5p, hsa-miR-590-3p, hsa-miR-576-5p, hsa-miR-200b-3p, hsa-miR-186-3p and hsa-miR-3934-5p, which are related to the overall survival of patients (p value < 0.05; Fig. 2).

Table 2

Univariate and multivariate Cox regression of differentially expressed miRNAs
miRNA	Univariate Cox regression		Multivariate Cox regression
miRNA	HR(95%CI)	P	HR(95%CI)	P
hsa-miR-33b-5p`	1.246(1.045ཞ1.486)	0.014	1.279(0.997ཞ1.641)	0.053
hsa-miR-33a-5p	1.216(1.005ཞ1.473)	0.045	1.196(0.934ཞ1.531)	0.156
hsa-miR-590-3p`	0.699(0.526ཞ0.928)	0.013	0.614(0.442ཞ0.853)	0.004
hsa-miR-576-5p`	0.556(0.413ཞ0.749)	0.000	0.644(0.477ཞ0.869)	0.004
hsa-miR-200b-3p`	0.876(0.785ཞ0.976)	0.017	0.775(0.691ཞ0.870)	0.000
hsa-miR-186-3p`	1.440(1.088ཞ1.906)	0.011	1.384(1.036ཞ1.850)	0.028
hsa-miR-3934-5p`	0.729(0.587ཞ0.905)	0.004	0.812(0.643ཞ1.025)	0.080

These 7 miRNAs make up the prognostic model, including hsa-miR-33b-5p, hsa-miR-33a-5p, hsa-miR-590-3p, hsa-miR-576-5p, hsa-miR-200b-3p, hsa-miR-186-3p and hsa -miR-3934-5p. Riskscore = (0.246 * hsa-miR-33b-5p` expression level) + (0.179 * hsa-miR-33a-5p` expression level) + (0.440 * hsa-miR-590-3p` expression level) + (0.641 * `hsa-miR-576-5p` expression level) + (0.254 * hsa-miR-200b-3p` expression level) + ( 0.325 * `hsa-miR-186-3p` expression level) + (0.208 * hsa-miR-3934-5p expression level). By calculating the risk score of each patient, the patients were divided into high-risk group and low-risk group based on the median value.

Roc Curve And Km Survival Analysis

Based on the model, we used the R language programming survival package and survivalROC package to draw relevant graphs. From the figure, we can observe that the area under curve(AUC) of the train group, test group, and all sample groups were 0.830, 0.722, and 0.770 respectively, all of which were greater than 0.7, indicating that this prognostic model has sensitivity and specificity in predicting the survival of patients with bladder cancer. (Fig. 3D-F). At the same time, (Fig. 3A-C) it can be seen that there was a significant difference between the high and low risk groups in the train group, the test group, and all sample groups (P < 0.05), and the survival rate decreased as the risk value increased. The train group showed that the 5-year overall survival rates for the low- and high-risk groups were 32.0% and 70.8%, respectively. The test team demonstrated that the 5-year survival rates for the overall high-risk and low-risk groups were 34.3% and 62.7%, respectively. The entire group showed that the 5-year overall survival rate of the high-risk and low-risk groups were 36% and 56.9%. We could observe that (Fig. 3H-I) it was the visualization results of the sample survival status of the train group, test group and all sample groups. From the figure, we could find that as the risk value increases, the number of death also increases. This evidence suggests that the model was sensitive and specific in predicting the survival of bladder cancer.

Independence Of Clinicopathological Factors

Univariate Cox regression analysis showed that a model constructed from 7 miRNAs was significantly related to overall patient survival (hazard ratio HR = 1.592, confidence interval 95% CI = 1.365–1.856, p = 3.19E-09; Table 3). Multivariate Cox regression analysis indicated that the comprehensive consideration of other commonly-seen clinical factors (HR = 1.605, 95% CI = 1.357–1.898, p = 3.28E-08), and other clinical factors including clinical stage, T stage, lymph node status, and distant metastasis, which still has nothing to do with overall survival rate. In addition, the patient's gender and age make it possible to become a prognostic indicator of bladder cancer in the future. At the same time, the patient's age also serves as an independent prognostic factor (HR = 1.030, 95% CI = 1.002–1.059, p = 0.037). The risk score of the model and other clinical characteristics ROC curves indicate that the risk score (0.698), clinical stage (0.637), T stage (0.623), lymph node status (0.631), and distant metastasis (0.523). These results indicate that this model still has good predictive ability (Fig. 4). In this regard, we can observe that risk score can be adopted as an independent prognostic factor for predicting the survival of bladder cancer patients.

Table 3

Univariate and multivariate Cox regression of clinical features
Clinical features	Univariate Cox regression		Multivariate Cox regression
Clinical features	HR(95%CI)	P	HR(95%CI)	P
Age (continuous variable)	1.030(1.002ཞ1.059)	0.037	1.035(1.005ཞ1.066)	0.024
Gender (male vs female)	0.618(0.353ཞ1.081)	0.091	0.551(0.311ཞ0.975)	0.041
Clinical stage (I + II vs III + IV))	1.742(1.216ཞ2.496)	0.002	1.166(0.573ཞ2.371)	0.672
T (T1 + 2 vs T3 + 4) N (N0 vs N1 + 2)	1.697(1.146ཞ2.513) 2.162(0.778ཞ6.011)	0.008 0.139	1.260(0.715ཞ2.219) 1.687(0.534ཞ5.328)	0.424 0.373
M (M0 vs M1)	1.541(1.169ཞ2.031)	0.002	1.195(0.732ཞ1.952)	0.467
Seven-miRNA signature	1.592(1.365ཞ1.856)	< 0.001	1.605(1.357ཞ1.898)	< 0.001

Prediction And Visualization Of 7 Mirna Target Genes

The seven miRNAs we selected were used to predict target genes from three databases (miRDB, miRTarBase, and Targetscan). We used the R software VennDiagram package to draw a visualized Wayne diagram to find reliable targets. (Fig. 5) In order to improve the reliability of bioinformatics, we used targets that were supported by all 3 databases, which could be used as targets for miRNAs. From: 4 we can observe that we have screened out 7 miRNA including hsa-miR-33b-5p, hsa-miR-33a-5p, hsa-miR-590-3p, hsa-miR-576-5p, hsa-miR-200b-3p, hsa-miR-186-3p, and hsa-miR-3934-5p,together with 17, 33, 169, 31, 74, 49 of hsa, and 22 overlapping bases. Seven miRNAs were predicted to have a total of 395 target genes.

Construction Of Mirna Regulatory Network

We can see from the regulatory network diagram that in order to explore whether these miRNA target genes may be involved in the progression of bladder urothelial carcinoma, 1111 differentially expressed genes obtained above (among which 545 genes were up-regulated and 566 genes were down-regulated) and 473 differentially expressed miRNAs (380 differentially expressed miRNAs up-regulated and 93 down-regulated) were used for analysis.7 up-regulated miRNAs (hsa-miR-33b-5p, hsa-miR-33a-5p, hsa-miR- 590-3p, hsa-miR-576-5p, hsa-miR-200b-3p, hsa-miR-186-3p and hsa-miR-3934-5p) and 247 down-regulated mRNAs. (The regulatory relationship between MiRNA and its target genes was shown in Fig. 6.)

Functional Enrichment Analysis of Related Target Genes in Bladder Urothelial Carcinoma

The top 15 terms in GO results are: Biological Processes (BP), Cellular Components (CC), and Molecular Function (MF) were demonstrated in dot plots (Fig. 7A-C). Among these three categories, BP analysis mainly includes cell junction, adhesion, migration, protease phosphorylation, regulation of enzyme activity and smooth muscle cell proliferation and differentiation. CC analysis mainly includes adhesion plaques, actin, extracellular matrix, fibrotic proteins, etc. MF analysis mainly includes metal ion transmembrane transporter activity, transcription activator activity, DNA binding, and ion channel binding. The result of the KEGG pathway of target genes associated with bladder urothelial carcinoma was 18 (Table 4), of which counts > 10 were mainly enriched in the MAPK signaling pathway, cGMP-PKG signaling pathway, the pi3k signaling pathway, the TGF-β signaling pathway, the interaction of miRNAs in tumors, and the regulation of inflammatory cell interaction. In addition, in order to provide a complex relationship between the target gene and the relative KEGG pathway, we adopted the R language software package to make a visual network diagram of the pathways (Fig. 7).

Table 4

KEGG pathways of target genes associated Bca
ID	Description	p.adjust	qvalue	Count	geneID
hsa04510	Focal adhesion	1.02E-05	8.02E-06	14	RAP1A/TLN1/LAMC1/VCL/ACTN1/ROCK2/PARVA/FLNC/ROCK1/MYLK/CAV1/CCND2/THBS1/PPP1CB
hsa05205	Proteoglycans in cancer	0.000337935	0.000264643	12	MYC/ROCK2/PRKACB/DCN/FLNC/SDC2/MMP2/ROCK1/CAV1/THBS1/ PPP1CB/CDKN1A
hsa04611	Platelet activation	0.00077537	0.000607205	9	RAP1A/PTGS1/TLN1/ROCK2/GNAQ/PRKACB/ROCK1/MYLK/PPP1CB
hsa04022	cGMP-PKG signaling pathway	0.006120072	0.004792734	9	PLN/SRF/ROCK2/MEF2A/GNAQ/MEF2D/ROCK1/MYLK/ PPP1CB
hsa05206	MicroRNAs in cancer	0.077452232	0.060654176	9	MYC/CYP1B1/SERPINB5/MCL1/ROCK1/BCL2L2/CCND2/ THBS1/ CDKN1A
hsa04921	Oxytocin signaling pathway	0.014623413	0.011451847	8	ROCK2/GNAQ/PRKACB/ROCK1/MYLK/PPP1CB/RCAN1/ CDKN1A
hsa04218	Cellular senescence	0.016462393	0.012891984	8	HIPK3/FOXO1/MYC/GADD45B/ZFP36L1/CCND2/PPP1CB/ CDKN1A
hsa04810	Regulation of actin cytoskeleton	0.047814395	0.037444275	8	GSN/VCL/ACTN1/ROCK2/ENAH/ROCK1/MYLK/PPP1CB
hsa05163	Human cytomegalovirus infection	0.06200175	0.048554637	8	MYC/CX3CL1/ROCK2/PTGER4/GNAQ/PRKACB/ROCK1/ CDKN1A
hsa04010	MAPK signaling pathway	0.13149261	0.102974126	8	RAP1A/MYC/DUSP3/SRF/PRKACB/GADD45B/FLNC/DUSP1
hsa05165	Human papillomavirus infection	0.176682825	0.138363361	8	FOXO1/LAMC1/PTGER4/PRKACB/PPP2R5A/CCND2/THBS1/CDKN1A
hsa04151	PI3K-Akt signaling pathway	0.195202228	0.152866225	8	LAMC1/MYC/PPP2R5A/SGK1/MCL1/CCND2/THBS1/CDKN1A
hsa04068	FoxO signaling pathway	0.017708049	0.013867478	7	FOXO1/FBXO32/GADD45B/SGK1/SOD2/CCND2/CDKN1A
hsa04270	Vascular smooth muscle contraction	0.017708049	0.013867478	7	CALD1/ROCK2/GNAQ/PRKACB/ROCK1/MYLK/PPP1CB
hsa05203	Viral carcinogenesis	0.077452232	0.060654176	7	GSN/SRF/ACTN1/PRKACB/IL6ST/CCND2/CDKN1A
hsa05131	Shigellosis	0.13149261	0.102974126	7	FOXO1/TLN1/VCL/ACTN1/ROCK2/BCL10/ROCK1
hsa04350	TGF-beta signaling pathway	0.017708049	0.013867478	6	SMAD7/MYC/ID4/DCN/ROCK1/THBS1
hsa04670	Leukocyte transendothelial migration	0.030279325	0.023712261	6	RAP1A/VCL/ACTN1/ROCK2/MMP2/ROCK1
hsa04261	Adrenergic signaling in cardiomyocytes	0.077452232	0.060654176	6	PLN/ACTC1/GNAQ/PRKACB/PPP2R5A/PPP1CB
hsa04390	Hippo signaling pathway	0.077452232	0.060654176	6	SMAD7/MYC/TEAD1/WWTR1/CCND2/PPP1CB
hsa04934	Cushing syndrome	0.077452232	0.060654176	6	RAP1A/GNAQ/PRKACB/PBX1/RASD1/CDKN1A
hsa04530	Tight junction	0.102825347	0.080524299	6	RAP1A/ACTN1/AMOTL2/ROCK2/PRKACB/ROCK1
hsa05202	Transcriptional misregulation in cancer	0.13149261	0.102974126	6	FOXO1/MYC/GADD45B/PBX1/CCND2/CDKN1A
hsa04062	Chemokine signaling pathway	0.13149261	0.102974126	6	RAP1A/CX3CL1/ROCK2/CXCL13/PRKACB/ROCK1
hsa05132	Salmonella infection	0.176682825	0.138363361	6	SNX9/MYC/RHOB/ROCK2/FYCO1/FLNC
hsa04024	cAMP signaling pathway	0.176682825	0.138363361	6	RAP1A/PLN/ROCK2/PRKACB/ROCK1/PPP1CB
hsa05166	Human T-cell leukemia virus 1 infection	0.178922574	0.140117348	6	TLN1/MYC/SRF/PRKACB/CCND2/CDKN1A
hsa04115	p53 signaling pathway	0.025719742	0.020141573	5	GADD45B/SERPINB5/CCND2/THBS1/CDKN1A

Construction Of Protein Interaction Network

We used the plug-in cytoHubba calculation in the cytoscape software to sort the network nodes in descending order based on the degree value, select the top ten genes as the central genes, and construct the central network map. The results indicate that 136 of the 247 target genes were filtered to the target genes for the PPI network complex, the top ten genes based on the degree value were (MYC, VCL, MMP2, ITPKB, FOXO1, SMAD7, CDKN1A, SRF, CAV1, KLF4) .(The specific genetic information were shown in Fig. 8 and Table 5)

Table 5

Identification of hub genes by cytoHubba
Node- name	MCC	DMNC	MNC	Degree	EPC	Bottle Neck	Ec Centricity	Closeness	Radiality	Betweenness	Stress	Clustering Coefficient
MYC	471	0.26055	25	30	58.693	22	0.18963	73.2	6.67437	4910.49147	15504	0.14253
VCL	223	0.21623	26	27	58.181	28	0.15802	68.45	6.48026	3446.04512	12240	0.1567
MMP2	151	0.32307	16	19	57.795	11	0.18963	64.2	6.41307	1972.27219	7728	0.21053
ITPKB	224	0.3491	14	16	56.642	29	0.18963	63.2	6.42053	1311.8855	5756	0.25833
FOXO1	296	0.40723	11	15	55.274	9	0.23704	59.16667	6.24882	1676.85358	6838	0.22857
SMAD7	84	0.3791	10	14	56.425	12	0.18963	60.28333	6.30108	1504.76416	6392	0.21978
CDKN1A	348	0.42441	12	14	55.379	6	0.23704	59.91667	6.29361	1257.4955	4962	0.31868
SRF	80	0.35915	10	12	55.853	5	0.18963	60.03333	6.33094	1118.09772	4770	0.27273
CAV1	53	0.358	9	12	54.393	7	0.15802	59.28333	6.27122	1227.38966	4312	0.22727
KLF4	180	0.39905	10	12	53.532	4	0.18963	56.28333	6.15176	648.29204	2142	0.30303

Survival Analysis

We mapped the K-M survival curve based on 247 target genes by using the R language package, and we obtained 40 genes related to survival prognosis (ACCTC1, AMOTL2, CALD1, CAV1, CLIC4, CPE, CSRP1, CXCL13, DCN, DKK3, DPYSL2, DUSP3, FAM168B, FILIP1L, FLNC, FSTL1, GSN, IL33, KLHL21, LAMC1, LIX1L, LRP1, MEF2A, MFAP4, MSRB3, MYC, MYLK, NFIC, PTGER4, RASD1, SBDS, SGCB, SLC39A14, SYNM, SYNP02, TIMP2, TMPRSS11E, VCL, WWTR1, TRIB1). (Fig. 9) From the graph of these genes, we can find that obvious differences occur in the high and low expression groups (p < 0.05) were of statistical significance.

Bladder cancer has been regarded as a highly malignant tumor that can undergo lymphatic metastasis, hematogenous metastasis, implantation metastasis, which was prone to recurrence. The recurrence rate of non-muscle invasive bladder cancer can reach as high as 50–70%. A small number of patients can progress to muscular invasive bladder cancer after relapse, which severely affects the survival prognosis of patients[15]. Therefore, it was becoming increasingly urgent for patients to find prognostic markers with high specificity and sensitivity. MiRNAs cause degradation or translational inhibition of target miRNAs through specific base pairing with target mRNAs, and regulate gene expression at the post-transcription level, thereby regulating tumor cell proliferation and apoptosis[11]. A miRNA can target one or more mRNAs, while multiple miRNAs can also simultaneously targets one mRNA [16]. Therefore, abnormally-expressed miRNAs were believed to play a critical role in tumorigenesis, which may not only be carcinogenic but also tumor suppressive [10]. In 2002, Calin reported the deletion of a small genomic region of 13q14 in patients with chronic lymphocytic leukemia, including the miR-15a and miR-16-1 genes, for the first time it demonstrated that miRNA abnormalities were closely related to tumors. In 2007, Catto Jw et al [17].. reported the miRNA expression profiles of bladder cancer for the first time. Since then, with the continuous improvement of detection ability, Chen et al[18] applied Solexa sequencing technology and found that 33 miRNA expressions were up-regulated in bladder cancer tissue compared with normal bladder tissue. The expression of 41 miRNAs was down-regulated, of which miRNA96 was identified as the most significantly up-regulated and miRNA-490-5p was most significantly down-regulated. It can be found that miRNAs have played an important role in the occurrence and development of bladder cancer, and its abnormal expression may be related to the survival prognosis of cystic cancer patients.

By using the R language edgeR package for the differential analysis and standardization of miRNA expression profiles from The Cancer Genome Atlas (TCGA), 473 differentially expressed miRNAs and then all patients were randomly divided into train combination test groups. Afterwards, seven miRNA prediction models were constructed by univariate Cox regression and stepwise multivariate Cox regression in the train group. (including hsa-miR-33b-5p, hsa-miR-33a-5p, hsa-miR-590-3p, hsa-miR-576-5p, hsa-miR-200b-3p, hsa-miR-186-3p and hsa-miR-3934-5p). At the same time, the predictive ability of 7 miRNAs was verified in the test group and the entire group. According to the median risk score grouping, the Kaplan-Meier curve showed that the overall survival rate of the low-risk group was significantly lower when compared with three high-risk groups. The AUC of the train group, test group, and entire group in the ROC curve all exceeded 0.7, indicating the sensitivity and specificity of these 7 miRNA prediction models in predicting the survival of bladder cancer patients. Univariate Cox regression and multivariate Cox regression analysis also suggested that these 7 miRNAs functioned as independent prognostic factors considering the conventional clinical factors of bladder cancer patients. Most of these 7 miRNAs were reported to participate in the research progress of various kinds of tumors. Studies by X, Y et al.[19] found that MiR-33b-5p was low expressed in gastric cancer cell lines. Up-regulated miR-33b-5p may reduce the expression of HMGA2, inhibit the growth of gastric cancer cells, and make gastric cancer cells sensitive to chemotherapy drugs. Studies by H, L, et al. [20]demonstrated that the long-chain non-coding RNA HIF1A-AS2 was overexpressed in osteosarcoma tissues and cells, which could interact with miR-33b-5p and negatively regulate its expression, thereby up-regulating the expression of the target protein SIRT6 of miR-33b-5p. Studies by S, W et al[21]showed that the expression of long-chain non-coding RNA EGOT in hepatocellular carcinoma (HCC) tissues was significantly up-regulated. EGOT has a significant regulatory effect on the survival, migration and invasion of liver cancer cells. EGOT expression level was negatively correlated with miR-33a-5p expression level. It was also confirmed that EGOT can specifically bind to miR-33a-5p, thereby reducing its expression and then up-regulating the expression of histone HMGA2. Gliomas has been regarded as the most commonly-seen and the most deadly type of malignant brain tumor. Studies by Y, S, et al [22] demonstrated that miR-590-3p was down-regulated and played a tumor suppressive role in glioma cells. L, Z, et al[23] have found that long-chain non-coding RNA linc-PINT can directly regulate miR-543 and miR-576-5p and inhibit the proliferation of esophageal cancer cell lines. It was proposed that linc-pin-mir − 543 / miR-576-5p pathway could predict the prognosis of esophageal cancer and provide new targets for the treatment of esophageal cancer. Studies by L, X et al[24] indicated that the expression of miR-200b-3p and miR-200c -3p in metastatic CRPC was significantly down-regulated, and was negatively correlated with the expression level of PRKAR2B in prostate cancer tissues. At first, M, H, et al. [25]demonstrated that EGFR agonist EREG has played an important role in tamoxifen-resistant breast cancer cells by activating EGFR signaling and its downstream glycolytic genes. It was further proved that EREG was a direct target of miR-186-3p, and tamoxifen down-regulation of miR-186-3p could lead to upregulation of EREG in tamoxifen-resistant breast cancer cells [26]. Research by W, Y et al [26] found that miR-3934-5p was significantly up-regulated in neuroblastoma tissues and cell lines. Further research by the authors also found that TP53INP1 served as a direct target of miR-3934-5p Gene, miR-3934-5p can negatively regulate it

In order to predict the potential biological functions of seven-miRNA signature, target genes of these seven miRNAs were analyzed by KEGG signaling pathway and GO enrichment analysis.GO analysis of target genes mainly demonstrated that it was related with connection, adhesion, migration, protease phosphorylation and regulation of enzyme activity with cells, the proliferation and differentiation of smooth muscle cells, adhesion plaques, actin, extracellular matrix, fibrosis protein, metal ion transmembrane transporter activity, transcription activator activity, DNA binding, ion channel binding, etc. The signal pathway of target genes was mainly enriched in the interaction of MAPK signal pathway, cGMP-PKG signal pathway, pi3k signal pathway, TGF-β signal pathway, and the miRNAs in tumors to regulate the interaction of inflammatory cells. H, C et al [27]. clarified that the MAPK pathway played a crucial role in the invasion and metastasis of bladder cancer cells.S, C et al [28]have found that the upregulation of myosin light chain kinase (MYLK) and activation of cGMP-PKG signaling pathway can silence serine / arginine repeat matrix protein 2-selective splicing (SRRM2-AS), thus inhibiting nasopharyngeal carcinoma cell proliferation, colony formation, angiogenesis, blocking the cell cycle and enhancing apoptosis.Studies by Xu, X, et al [29] found that the Pi3K pathway was specifically activated, vaccine-associated kinase 1 was up-regulateed (VRK1), and then histone H2A was induced to be up-regulated, eventually promoting the progression of osteosarcoma. S, J, et al [30]clarified that the signaling pathway of transforming growth factor-β (TGF-β) was involved in many cellular processes in mature organisms and developing embryos. These processes include cell growth, cell differentiation, apoptosis, cell dynamic balance and other cell functions. It has been found through research that TGF-β1 significantly promotes autophagy and mitochondrial phagocytosis of ovarian cancer cells, and provides a new direction for targeted therapy of ovarian cancer. These signaling pathways affect various tumors to varying degrees, and these four pathways were just the tip of the iceberg of target genes involved in signaling pathways, suggesting that the miRNA prognosis model we constructed for bladder cancer may also be related to these tumor signaling pathways.

In order to find the key nodes that regulate the prognosis model of bladder cancer miRNA, the network nodes were sorted in descending order based on the degree value according to Cytoscape and its cytoHubba, and the top ten genes were selected as the central genes (MYC, VCL, MMP2, ITPKB, FOXO1, SMAD7, CDKN1A, SRF, CAV1, KLF4). In addition, Kaplan-Meier showed that the expression of 39 genes (ACCTC1, AMOTL2, CALD1, CAV1, CLIC4, CPE, CSRP1, CXCL13, DCN, DKK3, DPYSL2, DUSP3, FAM168B, FILIP1L, FLNC, FSTL1, GSN, IL33, KLHL21, LAMC1, LIX1L, LRP1, MEF2A, MFAP4, MSRB3, MYC, MYLK, NFIC, PTGER4, RASD1, SBDS, SGCB, SLC39A14, SYNM, SYNP02, TIMP2, TMPRSS11E, VCL, WWTR1) was negatively correlated with survival prognosis, while the high expression of TRIB1 level was positively correlated with the survival prognosis. Among them, MYC, VCL, and CAV1 were not only the central genes of the miRNA prognosis model, but their expression was also negatively correlated with the survival prognosis. MYC genes were a group of oncogenes which were discovered earlier, including C-MYC, N-MYC, L-LYC. The MYC gene family and its products can promote cell proliferation, immortalization, dedifferentiation and transformation, etc, which has played a critical role in the formation of many types of tumors[31]. Related studies have shown that MYC gene was affected by DNA amplification in some types of tumors. MYC gene was amplified at a higher frequency in small cell lung cancer and was also amplified in other types of tumors such as breast cancer [32, 33].The VCL product was a highly conserved protein involved in cell proliferation, migration, and adhesion. Studies have shown that high expression degree of VCL may lead to significantly lower overall survival rates in patients with gastric cancer. Studies by S, F et al. have shown that Caveolin-1 (CAV1) functioned as the key component of caveolae protein. The CAV1 polymorphism affected mRNA expression through the effect of mutant alleles and has been proven to affect tumorigenesis, including bladder cancer, colon cancer, liver cancer, stomach cancer, breast cancer and lung cancer [34]. Therefore, miRNA may affect the progress and prognosis of bladder cancer by regulating MYC, VCL, and CAV1.

In summary, our study not only constructed a new miRNA model related to the prognosis of bladder cancer through the expression profile of miRNA, but also verified and evaluated the predictive ability of the model, and found that it could be used as an independent prognostic factor for bladder cancer. In addition, inferring its potential functions by predicting the target genes of the model has enhanced our understanding of the development of bladder cancer. However, this was only a study of the TCGA database based on bioinformatics. It is hoped that other databases and a large number of experiments will be able to verify the feasibility of this prognostic model and provide reliable prediction and treatment targets for bladder cancer patients.

Declarations

Acknowledgements

None.

Funding

Not applicable.

Conflict of interests

All authors declare no conflict of interests.

Data Availability Statement

All data generated or analyzed during the present study was downloaded from TCGA database and miRBase website.

Author Contributions

Study concepts: Ning Xu; Yong Wei.

Study design: Jin-Xing Lv;Ning Xu.

Data acquisition: Shao-Hao Chen;Peng-Fei Zhuang.

Quality control of data and algorithms: Xue-Yi Xue.

Data analysis and interpretation: Ye-Hui Chen; Qing-Shui Zheng.

Statistical analysis: Zhi-Bin Ke.

Manuscript preparation: Jin-Xing Lv; Fei Lin; Yun-Zhi Lin; Zhi-Bin Ke.

Manuscript editing: Jin-Xing Lv; Fei Lin.

Manuscript review: Ning Xu; Yong Wei.

Consent for publication

Not applicable.

Ethics approval and consent to participate

Not applicable.

Feng D, et al. Comparison of perioperative and oncologic outcomes between robot-assisted and laparoscopic radical cystectomy for bladder cancer: a systematic review and updated meta-analysis. International urology and nephrology. 2020.
Lin YZ, et al. Bioinformatics Analysis of the Expression of Key Long Intergenic Non-Protein Coding RNA Genes in Bladder Cancer. Med Sci Monit. 2020;26:e920504.
Lin N, et al. Risk factors for upper tract urothelial recurrence following local excision of bladder cancer. Cancer medicine. 2018;7(8):4098–103.
Nargund VH, et al. Management of non-muscle-invasive (superficial) bladder cancer. Semin Oncol. 2012;39(5):559–72.
Xu N, et al. ARPC4 promotes bladder cancer cell invasion and is associated with lymph node metastasis. J Cell Biochem. 2020;121(1):231–43.
Solomon JP, Hansel DE. Prognostic factors in urothelial carcinoma of the bladder: histologic and molecular correlates. Adv Anat Pathol. 2015;22(2):102–12.
Wang G, et al. Up-regulation of microRNA in bladder tumor tissue is not common. Int Urol Nephrol. 2010;42(1):95–102.
Fabbri M, et al. MicroRNAs Cancer journal (Sudbury Mass). 2008;14(1):1–6.
Carthew RW, Sontheimer EJ. Origins and Mechanisms of miRNAs and siRNAs. Cell. 2009;136(4):642–55.
Sassen S, et al. MicroRNA: implications for cancer. Virchows Arch. 2008;452(1):1–10.
Catto JW, et al. MicroRNA in prostate, bladder, and kidney cancer: a systematic review. Eur Urol. 2011;59(5):671–81.
Liu N, et al. Prognostic value of a microRNA signature in nasopharyngeal carcinoma: a microRNA expression analysis. The Lancet Oncology. 2012;13(6):633–41.
Xu N, et al. Identification of survival-associated alternative splicing events and signatures in adrenocortical carcinoma based on TCGA SpliceSeq data. Aging. 2020;12(6):4996–5009.
Liang YC, et al. LIMK1 depletion enhances fasudil-dependent inhibition of urethral fibroblast proliferation and migration. J Cell Biochem. 2019;120(8):12977–88.
Ambros V. The functions of animal microRNAs. Nature. 2004;431(7006):350–5.
Hanlon K, et al. Investigating the targets of MIR-15a and MIR-16-1 in patients with chronic lymphocytic leukemia (CLL). PLoS One. 2009;4(9):e7169.
Catto JW, et al. Distinct microRNA alterations characterize high- and low-grade bladder cancer. Cancer Res. 2009;69(21):8472–81.
Chen YH, et al. Characterization of microRNAs expression profiling in one group of Chinese urothelial cell carcinoma identified by Solexa sequencing. Urol Oncol. 2013;31(2):219–27.
Yang X, et al. MiR-33b-5p sensitizes gastric cancer cells to chemotherapy drugs via inhibiting HMGA2 expression. J Drug Target. 2017;25(7):653–60.
Lin H, et al. Long noncoding RNA HIF1A-AS2 facilitates cell survival and migration by sponging miR-33b-5p to modulate SIRT6 expression in osteosarcoma. Biochemistry and cell biology = Biochimie et biologie cellulaire. 2019:1–9.
Wu S, et al. Long Non-Coding RNA EGOT Promotes the Malignant Phenotypes of Hepatocellular Carcinoma Cells and Increases the Expression of HMGA2 via Down-Regulating miR-33a-5p. OncoTargets therapy. 2019;12:11623–35.
Song Y, et al. Inhibition of the aberrant A1CF-FAM224A-miR-590-3p-ZNF143 positive feedback loop attenuated malignant biological behaviors of glioma cells. J Exp Clin Cancer Res. 2019;38(1):248.
Zhang L, et al. Linc-PINT acted as a tumor suppressor by sponging miR-543 and miR-576-5p in esophageal cancer. Journal of cellular biochemistry. 2019;120(12):19345–57.
Xia L, et al. Transcriptional regulation of PRKAR2B by miR-200b-3p/200c-3p and XBP1 in human prostate cancer. Biomed Pharmacother. 2020;124:109863.
He M, et al. The miR-186-3p/EREG axis orchestrates tamoxifen resistance and aerobic glycolysis in breast cancer cells. Oncogene. 2019;38(28):5551–65.
Ye W, et al. Downregulation of microRNA-3934-5p induces apoptosis and inhibits the proliferation of neuroblastoma cells by targeting TP53INP1. Exp Ther Med. 2019;18(5):3729–36.
Chao H, et al. RAB14 activates MAPK signaling to promote bladder tumorigenesis. Carcinogenesis. 2019;40(11):1341–51.
Chen S, et al. Silencing of long noncoding RNA SRRM2-AS exerts suppressive effects on angiogenesis in nasopharyngeal carcinoma via activating MYLK-mediated cGMP-PKG signaling pathway. Journal of cellular physiology. 2019.
Xu X, Yu H. Ras-PI3K pathway promotes osteosarcoma progression via regulating VRK1-mediated H2A phosphorylation at threonine 120. Artif Cells Nanomed Biotechnol. 2019;47(1):4274–83.
Jin S, et al. TGF-beta1 fucosylation enhances the autophagy and mitophagy via PI3K/Akt and Ras-Raf-MEK-ERK in ovarian carcinoma. Biochem Biophys Res Commun. 2020;524(4):970–76.
Venkateswaran N, Conacci-Sorrell M. MYC leads the way. Small GTPases. 2020;11(2):86–94.
Katsuta E, et al. High MYC mRNA Expression Is More Clinically Relevant than MYC DNA Amplification in Triple-Negative Breast Cancer. Int J Mol Sci. 2019;21(1).
Zhang M, et al. Vinculin promotes gastric cancer proliferation and migration and predicts poor prognosis in patients with gastric cancer. Journal of cellular biochemistry. 2019;120(8):14107–15.
Fan S, et al. CAV1 polymorphisms rs1049334, rs1049337, rs7804372 might be the potential risk in tumorigenicity of urinary cancer: A systematic review and meta-analysis. Pathol Res Pract. 2019;215(1):151–58.

Download PDF

Version 1

posted

You are reading this latest preprint version

Development and Validation of A Seven-microRNA Signature as a Novel Prognostic Biomarker for Bladder Cancer

Status:

Version 1

Abstract

Figures

Introduction

Materials And Methods

Data Download

Differential Analysis Of Mrna And Mirna

Model Establishment And Validation

Results

Screening of Differentially Expressed mRNAs and miRNAs

Roc Curve And Km Survival Analysis

Independence Of Clinicopathological Factors

Prediction And Visualization Of 7 Mirna Target Genes

Construction Of Mirna Regulatory Network

Construction Of Protein Interaction Network

Survival Analysis

Discussion

Conclusion

Declarations

References

Status:

Version 1