Sage Journals: Discover world-class research

Abstract

OBJECTIVE:

To identify the mRNAs associated with bladder cancer (BC) recurrence.

METHODS:

The transcription profile of GSE31684 including 39 recurrent BC tumor samples and 54 non-recurrent BC tumor samples as well as transcription profile of GSE13507 including 36 recurrent BC tumor samples and 67 non-recurrent BC tumor samples were downlaoded from the Gene Expression Omnibus. Then, the differentially expressed genes (DEGs) were identified using linear models for microarray data (limma) and the intersections of DEGs from the two datasets were further screened. The weighed gene co-expression network analysis (WGCNA) was used to screen the modules related to BC recurrence. Protein-protein interaction (PPI) network analysis was used to analyze the genes interaction. Their functions were predicted by Gene Ontology and KEGG pathway enrichment. Moreover, The Comparative Toxicogenomics Database 2017 update (CTD) was used to search the BC related pathway. The univariate cox regression analysis was used to identify DEGs associated to the recurrence. Kaplan-Meier plots were used to illustrate recurrence free survival time (RFS).

RESULTS:

A total of 692 intersections DEGs were screened. WGCNA showed that 7 modules (2279 genes) were stable in both the datasets. A total of 169 intersection DEGs were mapped to the 7 modules. There existed 149 interaction relationships among 81 proteins (18 down-regulated and 63 up-regulated DEGs) in the PPI network. Two KEGG pathways including Focal adhesion and ECM-receptor interaction were enriched which were also found in the CTD. The univariate cox regression analysis showed that 3 DEGs (COL4A1, COL1A2 and COL5A1) were significant related to the prognosis. Multivariate cox regression analysis revealed that pathologic_N (N0-N1 vs N2-N3, $p=$ 0.033) were independent prognostic factors for overall survival in patients with BC.

CONCLUSION:

COL4A1, COL1A2 and COL5A1 could be associated with BC recurrence.

Keywords

Biomarker bladder cancer recurrence mRNA weighed gene co-expression network analysis

1. Introduction

Bladder cancer (BC) is the most common malignancy of the urinary system in China and the ninth most freguent cancer worldwide [1]. Transurethral resection (TUR) is generally used to manage the BC. However, a large number of patients have a 50–70% risk of recurrence after first TUR [2, 3]. Moreover, a total of 15–25% of the patients are also at risk of progression to invasive cancer [4]. No good biomarkers could be used for the prediction of BC recurrence.

With the development of molecular biology, molecular biomarkers have been proven to be a promising clinical tool for the BC prognosis. For example, by using the bioinformatics analysis of weighted gene coexpression network analysis (WGCNA), Yan et al. found four genes could be used to predict the prognosis of BC [5]; The increased expression of matriptase, which is a member of the Type II transmembrane serine protease family, could present poor prognosis of BC [6]; However, the mRNA biomarker identified for the recurrence of BC is little. Nakahara et al. found that the expression of protease activating receptor-2 (PAR-2) is positively correlated with the recurrence of non-muscle invasive bladder cancer. Therefore, it is of great need to identify the biomarkers of BC recurrence.

In this study, using the same microarray data by Riester et al. [7] and Kim et al. [8], we aimed to screen the DEGs between recurrence and non-recurrence using the linear models for microarray data (limma) [9] with the threshold of false discovery rate (FDR)-value $<$ 0.05 and $|$ log2 fold change (FC) $|$ $>$ 0.5. Additionally, the functional enrichment, WGCNA, protein-protein interaction (PPI) network analysis, The Comparative Toxicogenomics Database (CTD) were used to identify the critical genes related BC recurrence. Reportedly, analyses based on differential statistical tests could result in different outcomes [10]; Thus, we believe that some variable results might be obtained.

2. Materials and methods

2.1 Microarray data

The transcription profile of GSE31684 [7] (GP570, Affymetrix Human Genome U133 Plus 2.0 Array) including 39 recurrent BC tumor samples and 54 non-recurrent BC tumor samples from the Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/). Additionally, we also downloaded the transcription profile of GSE13507 [8] (GPL6102, Illumina human-6 v2.0 expression beadchip) including 36 recurrent BC tumor samples and 67 non-recurrent BC tumor samples from the GEO. Above all, a total of 75 recurrent BC tumor samples (recurrence) and 121 non-recurrent BC tumor samples (non-recurrence) were included in this study. The study design was shown in Fig. S1.

2.2 The dataset preprocessing and DEGs analysis

The R3.4.1 oligo (version: 3.6, http://www.bioconduc tor.org/packages/release/bioc/html/oligo.html) [11] and Limma package in R language (Version 3.34.0, https:// bioconductor.org/packages/release/bioc/html/limma. html) [9] was employed to perform the background correction and the normalization. Then, we used the limma package in R language to screen the DEGs between recurrence and non-recurrence. The FDR-value $<$ 0.05 and $|$ log2 FC $|$ $>$ 0.5 were used as the cutoff criteria for DEGs. The intersections of DEGs from the GSE31684 and GSE13507 were screened. The pheatmap package in R (version 1.0.8, https://cran.r-project.org/package =pheatmap) [12] was employed to conduct the bidirectional hierarchical clustering.

2.3 WGCNA

WGCNA (version 1.61, https://cran.r-project.org/web /packages/WGCNA/index.html) [13] was used to discover gene clusters (or modules) in high-throughput data that were highly correlated with recurrence. The GSE13507 was used as the training dataset, and the GSE31684 was used as the validation dataset. Briefly, By setting a series of soft-thresholding power values, the correlation coefficient and the average connection degree of the connection degrees k and p(k) under each power value were calculated. The cutHeight $=$ 0.99 and the gene number $\geqslant$ 80 within a module were used as the threshold. Then, the intersection DEGs were mapped to the validated modules to screen the module related DEGs.

Table 1
The top ten up-regulated and the down-regulated DEGs between recurrence and non-recurrence bladder cancer in the two datasets

DEGs	Genes	FDR	logFC	Genes	FDR	logFC
	GSE31684			GSE13507
Up-regulated	KRT20	0.0106956	1.92	SCUBE2	0.0017478	2.	315
	EFEMP1	0.0027141	1.78	MPDZ	0.0330861	2.	22
	IGFBP5	0.0006684	1.68	GULP1	0.0224745	2.	1425
	HLA-DRB4	0.0045324	1.63	TCEA3	0.0136794	2.	1175
	COMP	0.0048096	1.5	PLN	0.0007131	2.	0075
	PTGIS	0.0059337	1.47	NRP2	0.0078915	1.	9575
	SFRP2	0.0341571	1.46	ITGBL1	0.0025401	1.	8575
	SULF1	0.0158658	1.44	CNN1	0.0106611	1.	75
	POSTN	0.016968	1.44	TGFB2	0.0207231	1.	7125
	PRSS2	0.0006639	1.42	CTHRC1	0.0046896	1.	6475
Down-regulated	AKR1B10	0.0003297	$-$ 2.16	DKK1	0.0059748	$-$ 2.	0775
	CLCA2	0.0096522	$-$ 1.51	ADGRV1	0.0011859	$-$ 2.	0725
	GJB6	0.0269205	$-$ 1.35	PEG3	0.0118275	$-$ 1.	8675
	CYP24A1	0.0184545	$-$ 1.32	SPRR3	0.0189708	$-$ 1.	835
	SPRR1B	0.0307749	$-$ 1.3	TSPAN1	0.0091401	$-$ 1.	8175
	S100A2	0.0318657	$-$ 1.25	LGALS4	0.0147888	$-$ 1.	7725
	S100A9	0.0219258	$-$ 1.21	CEL	0.0108606	$-$ 1.	6975
	MMP1	0.027594	$-$ 1.18	SCIN	0.0159102	$-$ 1.	6725
	S100A14	0.0142584	$-$ 1.17	SPINK1	0.0478206	$-$ 1.	6275
	KRTDAP	0.0259278	$-$ 1.17	TMEM116	0.0003351	$-$ 1.	5575

DEGs, differentially expressed genes; FC, fold change.

2.4 PPI network construction

The STRING (Search Tool for the Retrieval of Interacting Genes/Proteins, version 10.0, http://string-db. org/) online software [14] was used to search interrelationships of module related DEGs. Then the PPI network was constructed and visualized using the Cytoscape (version 3.6.1, http://www.cytoscape.org/) [15]. Furthermore, DAVID (The Database for Annotation, Visualization and Integrated Discovery, version 6.8, https://david.ncifcrf.gov/) was used to perform the GO and KEGG pathway enrichment for the DEGs in the PPI network with the threshold of $p$ -value $<$ 0.05 [16].

2.5 BC related DEGs screening

The CTD (http://ctd.mdibl.org/) [17] were used to search the BC related pathway with keyword of bladder cancer. Then the searched pathways were further filterd based on the enriched pathways of the genes in the PPI network. Finally, the BC related DEGs were screened.

2.6 Recurrence related DEGs screening

The univariate cox regression analysis of survival package (Version 2.41-1, http://bioconductor.org/packa ges/survivalr/) in R language [18] was used to identify DEGs associated to the recurrence in the dataset of GSE31684. In addition, using the median of DEGs (associated to the recurrence) expression values as the cut-off, the samples in the GSE31684 dataset were divided into high expression group (above the median) and low expression group (under the median). Then, Kaplan-Meier plots were used to illustrate recurrence free survival time (RFS). The multivariate cox regression analysis for these three DEGs were conducted.

We also did univariate and multivariate cox regression analysis using backward selection to test the independent signifcance of different clinical factors; the $p$ value threshold was 0.1 ( $p>$ 0.1) for removing non-signifcant variables from the analysis, and marginally signifcant variables (0.05 $<p<$ 0.1) remained in the final cox model. Covariates included radiotherapy (yes vs no), pathologic_M (M0 vs M1), pathologic_N (N0-N1 vs N2-N3) and so on.

The genes expression profile (Illumina HiSeq 2000 RNA Sequencing) of 336 BC tumor samples with the complete clinical information (recurrence or not) were downloaded from the the cancer genome atlas (TCGA, https://gdc-portal.nci.nih.gov/) as the validation dataset.

Figure 1.

A, the volcano plot of DEGs between recurrence and non-recurrence; left for GSE31684; right for GSE13507; The two red vertical dotted lines represent $|$ log ${}_{2}$ FC $|$ $>$ 0.5; the red Horizontal dotted line represents FDR $<$ 0.05; Pink and blue dots indicated significantly up-regulated and down-regulated DEGs, respectively.

Figure 2.

A, the venn showed the intersections of DEGs between the two datasets; B, the two-way clustering of differentially expressed genes (DEGs); horizontal axis, the samples; vertical axis, DEGs between recurrence and non-recurrence; left for GSE31684; right for GSE13507; color key, the log2FC of DEGs.

Figure 3.

Protein-protein interaction network analysis; the colors of the node indicated the log2FC; the size of the nodes represent the degree.

3. Results

3.1 DEGs analysis

Based on the microarray data analysis between recurrence and non-recurrence by Limma, we obtained 1496 DEGs (362 down-regulated and 1134 up-regulated DEGs) in GSE31684 and 1570 DEGs (850 down-regulated and 720 up-regulated DEGs) in GSE13507. Table 1 lists the top 10 markedly up-regulated and downregulated DEGs in the two datasets. We used volcano plots for the visualization and assessment of the variation (or reproducibility) of genes expression between recurrence and non-recurrence for GSE31684 (Fig. 1A) and GSE13507 (Fig. 1B). A total of 692 intersections DEGs (221 up-regulated and 420 down-regulated DEGs) were identified in these two datasets (Fig. 2A). Furthermore, the bidirectional hierarchical clustering revealed that DEGs could distinguish these two groups (Fig. 2B).

3.2 Modules and genes screened by WGCNA

According to the method, the power value when the correlation coefficient squared value of connection degree k and p(k) for the first time to reach 0.9 were selected, that was, power $=$ 12; under this power parameter, the mean connectivity degree of the geness was 1, which conformed to the small-world property in a scale-free network. Based on clustering and dynamic pruning, 3855 highly correlated genes were clustered into 10 modules in the training dataset GSE13507 (Fig. S2A). Among which, 7 modules (2279 genes) including modules black, blue, brown, green, pink, turquoise and yellow were validated using the validation dataset GSE31684 (Table S1). A total of 169 intersection DEGs were mapped to the 7 modules.

3.3 PPI network analysis

In the PPI network (Fig. 3), there existed 149 interaction relationships among 81 proteins (18 down-regulated and 63 up-regulated DEGs) from 169 proteins. The 2 KEGG pathways including Focal adhesion and ECM-receptor interaction were enriched with the DEGs in the PPI network (Table 2).

Table 2
The GO and KEGG pathway enrichment for the DEGs in the PPI network

Category	Term	Count	PValue
Biology process	GO: 0007155 $\sim$ cell adhesion	18	4.24E-08
	GO: 0022610 $\sim$ biological adhesion	18	4.33E-08
	GO: 0001568 $\sim$ blood vessel development	9	2.99E-05
	GO: 0001944 $\sim$ vasculature development	9	3.55E-05
	GO: 0022604 $\sim$ regulation of cell morphogenesis	6	5.04E-04
	GO: 0008283 $\sim$ cell proliferation	9	0.0014835
	GO: 0019318 $\sim$ hexose metabolic process	6	0.0027646
	GO: 0010035 $\sim$ response to inorganic substance	6	0.0036614
	GO: 0048514 $\sim$ blood vessel morphogenesis	6	0.0041389
	GO: 0005996 $\sim$ monosaccharide metabolic process	6	0.0051278
	GO: 0051050 $\sim$ positive regulation of transport	6	0.0052253
	GO: 0007010 $\sim$ cytoskeleton organization	8	0.0061376
	GO: 0007167 $\sim$ enzyme linked receptor protein signaling pathway	7	0.0073337
	GO: 0006928 $\sim$ cell motion	8	0.0096284
	GO: 0010033 $\sim$ response to organic substance	10	0.0097112
	GO: 0009725 $\sim$ response to hormone stimulus	7	0.010191
	GO: 0009719 $\sim$ response to endogenous stimulus	7	0.0159407
	GO: 0043933 $\sim$ macromolecular complex subunit organization	9	0.0254293
	GO: 0009628 $\sim$ response to abiotic stimulus	6	0.0374444
	GO: 0065003 $\sim$ macromolecular complex assembly	8	0.0488364
KEGG pathway	hsa04510: Focal adhesion	7	5.22E-04
	hsa04512: ECM-receptor interaction	5	9.10E-04

3.4 BC related target gene screening

In the CTD database, 210 KEGG pathways were searched, among which vascular smooth muscle contraction and ECM-receptor interaction were intersected with the enriched KEGG pathway for the DEGs (including 7 genes: COL4A1, VEGFA, COL6A3, COL1A2, COL6A1, COL5A1 and MYL9) in the PPI network.

Table 3
Clinical factor statistics and cox regression analysis

Clinical characteristics	Uni-variables cox			Multi-variables cox
	HR	95% CI	$p$	HR	95% CI	$p$
Age (years, mean $\pm$ sd)	1.013	0.992–1.034	2.33E-01	–	–	–
Gender (male/female)	0.976	0.596–1.598	9.22E-01	–	–	–
Pathologic_M (M0/M1)	5.813	2.002–16.88	2.46E-04	1.755	0.379–8.130	4.72E-01
Pathologic_N (N0-N1/N2-N3)	1.719	1.354–2.184	3.12E-04	2.332	1.069–5.085	3.33E-02
Pathologic_T (T1-T2/T3-T4)	1.652	1.191–2.291	1.92E-03	3.267	0.499–21.35	2.16E-01
Pathologic_stage (I-II/III-IV)	1.812	1.380–2.379	1.58E-03	0.7638	0.115–5.071	7.80E-01
Pathologic_grade (high/low)	4.065	0.564–29.28	1.31E-01	–	–	–
Radiotherapy (yes/no/-)	1.67E+00	0.768–3.619	1.91E-01	–	–	–

Note: TCGA, The Cancer Genome Atlas; HR, hazard ratio; $p$ , $p$ value.

Figure 4.

Kaplan-Meier curves of recurrence free survival according to low-expression or high-expression categories; (A) For the dataset of GSE31684; (B) For the dataset of TCGA.

3.5 Selection of Candidate DEGs

The univariate cox regression analysis obtained 3 DEGs (including COL4A1, COL1A2 and COL5A1) which were significant related to the prognosis.

Using the median of the expression of the 3 DEGs (COL5A1, COL1A2 and COL4A1) as the cut-off point, each individual was then classified into low expression ( $n=$ 46) or high expression ( $n=$ 47) categories in the GSE31684 dataset. The high expression patients had poorer RFS (hazard ratio [HR], 2.041; 95% CI, 1.058 to 3.936; $P=$ 0.0298) than did the low expression patients in the training cohort for COL5A1; hazard ratio [HR], 1.876; 95% CI, 1.281 to 3.591; $P=$ 0.0453 for COL1A2 and hazard ratio [HR], 1.872; 95% CI, 0.983 to 3.564; $P=$ 0.0236 for COL4A1. The validation analyses were performed in an internal validation cohort. A total of 168 high expression patients had poorer RFS (hazard ratio [HR], 1.961; 95% CI, 1.255 to 3.065; $P=$ 0.00259; Fig. 2B) than did the 168 low expression patients in the validation cohort for COL5A1; hazard ratio [HR], 1.893; 95% CI, 1.218 to 2.943; $P=$ 0.00394 for COL1A2 and hazard ratio [HR], 1.585; 95% CI, 1.018 to 2.470; $P=$ 0.0398 for COL4A1 (Fig. 4).

Figure 5.

Histogram of three important genes expression in GSE31684 (A) and TCGA (B) between recurrence and non-recurrence of bladder cancer.

Additionally, the expression of the 3 DEGs (COL4A1, COL1A2 and COL5A1) was analyzed in the datasets of TCGA and GSE31684 and all the expressions of the three DEGs were significantly higher in recurrence than non- recurrence (Fig. 5).

The univariate and multivariate analysis of overall survival by clinical factors in the training dataset was shown in Table 3. Using univariate analysis, overall survival times were revealed to be significantly associated with pathologic_M (M0 vs M1, $p=$ 0.000246), pathologic_N (N0-N1 vs N2-N3, $p=$ 0.000312), Pathologic_T (T1-T2 vs T3-T4, $p=$ 0.000192) and pathologic_stage (I-II vs III-IV, $p=$ 0.00158). The multivariate analysis revealed that pathologic_N (N0-N1 vs N2-N3, $p=$ 0.033) was independent prognostic factors for overall survival in patients with BC (Table 3).

The multivariate analysis on the three DEGs showed that in the dataset of GSE31684, COL4A1 was the independent prognostic factor for overall survival in BC (HR $=$ 1.22689, 95% CI $=$ 1.0197 to 1.637, $p=$ 0.0164) while in the TCGA dataset, it is not significant ( $p=$ 0.0729) (Table 4). For the other two genes, they were not the independent prognostic factor for overall survival in BC.

Table 4

Multivariate analysis for COL5A1, COL4A1, COL1A2

	Gene	Coef	$p$	HR	95% CI
GSE31684 dataset	COL5A1	$-$ 0.03013	0.875	0.97032	0.6661	–1.413
	COL4A1	0.20448	0.0164	1.22689	1.0197	–1.637
	COL1A2	0.13795	0.697	1.14792	0.573	–2.3
TCGA dataset	COL5A1	$-$ 0.01545	0.906	0.98466	0.7615	–1.273
	COL4A1	0.10004	0.0729	1.10521	0.9277	–1.946
	COL1A2	0.00363	0.99	1.00364	0.5837	–1.726

Note: TCGA, The Cancer Genome Atlas; coef, the coefficients, HR, hazard ratio; $p$ , $p$ value, CI, confidence interval.

4. Discusion

In this study, by using bioinformatics methods, we used genes microarray data from GEO and RNA sequencing data TCGA to develop and validate a novel prognostic tool based on 3 genes (COL4A1, COL1A2 and COL5A1). Our results indicated that COL4A1, COL1A2 and COL5A1 identified in this study could be used to categorize patients into low expression and high expression who had significantly different RFS.

COL4A1, COL1A2 and COL5A1 belong to the collagen family which are the most abundant components of extracellular matrix (ECM). They have been discovered to play an important role in the structural integrity and tensile strength of human tissues and organs [19]. As for cancer, collagens could regulate the physical and the biochemical properties of the tumor microenvironment and modulate cancer cell polarity, migration and signaling [20].

COL4A1 was significantly up-regulated in the recurrence than non-recurrence (log ${}_{2}$ FC $=$ 0.811, FDR $=$ 0.0038085). COL4A1 (Collagen type IV alpha 1) is the subunit of the type IV collagen, which is the most abundant and essential in basement membranes (BMs). The BM acts as a critical component regulating tumor cell behavior [21]. Former studies have illustrated the importance of COL4A1 in angiogenesis and tumor progression [22]. Moreover, COL4A1 was found to be produced in cancer cells and promote tumor budding at the invasion front in human urothelial carcinoma of the bladder [23].

Additionally, COL1A2 was significantly up-regulated in recurrence than non-recurrence (log ${}_{2}$ FC $=$ 0.564, FDR $=$ 0.031287). The COL1A2 (collagen type 1 $\alpha$ 2) gene beonging to the type I collagen whoses remodeling in the ECM microenvironment accompanies stromal invasion and epithelial tumorigenesis in other cancer types [24]. And it was discovered that there is a positive association of high COL1A2 mRNA expression to progression and decreased overall survival in non-muscle invasive bladder cancer [25].

At last, the COL5A1 was also significantly up-regulated in the recurrence than non-recurrence (log ${}_{2}$ FC $=$ 0.901, FDR $=$ 0.0331431). The COL5A1 is in smaller fibrillar collagen in mammals. The high expression of COL5A1 indicated poor prognosis of the BC patients [26]. COL5A1 was associated with PI3K-Akt and TNF signaling pathways in BC cells and the down-regulated of COL5A1 could constrain the proliferative and migrated abilities of BC cells [27].

The KEGG pathway enrichment showed that COL4A1, COL1A2 and COL5A1 were significantly enriched in Focal adhesion and ECM-receptor interaction. The ECM-receptor interaction pathway and focal adhesion pathway are associated with the regulation of heterogeneity for muscle-invasive BC [28]. And ECM-receptor interaction pathway and focal adhesion pathway are found to be related to the progression of muscle-invasive BC [29].

Our study had several limitations. First, the samples from the TCGA were muscle-invasive BC while the datasets from the GEO were non muscle-invasive BC. It is not ideal to use the TCGA dataset as the validation dataset and the further study was still needed. Secondly, since the study was to screen the biomarkers for BC, so the blood samples which is the noninvasive were needed to validate the results.

In conclusion, COL4A1, COL1A2 and COL5A1 involved in the pathway of focal adhesion and ECM-receptor interaction could be associated with BC recurrence.

Footnotes

Conflict of interest

None to report.

Supplementary data

Table S1

Module validation statistics

ID	Color	Z-score
Module 1	Black	35.070
Module 2	Blue	10.571
Module 3	Brown	20.446
Module 4	Green	7.855
Module 5	Grey	2.599
Module 6	Magenta	$-$ 0.415
Module 7	Pink	5.159
Module 8	Red	0.409
Module 9	Turquoise	21.160
Module 10	Yellow	46.270

Note: Z-score represented the stability of the modules in the two datasets, Z-score $>$ 5 means stable.

Figure S1.

The schematic representation of the workflow.

Figure S2.

WGCNA of genes from the two datasets; (A) Hierarchical cluster tree showing co-expression modules identified by WGCNA (GSE31684 and GSE13507). Each leaf in the tree represents one gene. The major tree branches constitute 10 modules, labeled with different colors; (B) Module-trait association. Each row corresponds to a module, labeled with a color as in (A); The column corresponds to a clinical trait-recurrence; The color of each cell at the row – column intersection indicates the correlation coefficient between the module and the trait; A $p$ value was also listed in the brackets under the coefficient.

References

Chen

Zheng

Baade

P.D.

Zhang

Zeng

Bray

Jemal

X.Q.

and He

, Cancer statistics in China, 2015, CA Cancer J Clin 66 (2016), 115–132.

Bohle

, A single immediate postoperative instillation of chemotherapy decreases the risk of recurrence in patients with stage Ta T1 bladder cancer: A meta-analysis of published results of randomized clinical trials, Int Braz J Urol 30 (2004), 353–354.

Sylvester

R.J.

Oosterlinck

Holmang

Sydes

M.R.

Birtle

Gudjonsson

De Nunzio

Okamura

Kaasinen

Solsona

Ali-El-Dein

Tatar

C.A.

Inman

B.A.

N’Dow

Oddens

J.R.

and Babjuk

, Systematic review and individual patient data meta-analysis of randomized trials comparing a single immediate instillation of chemotherapy after transurethral resection with transurethral resection alone in patients with stage pTa-pT1 urothelial carcinoma of the bladder: Which patients benefit from the instillation? Eur Urol 69 (2016), 231–244.

Frantzi

Latosinska

Fluhe

Hupe

M.C.

Critselis

Kramer

M.W.

Merseburger

A.S.

Mischak

and Vlahou

, Developing proteomic biomarkers for bladder cancer: Towards clinical application, Nat Rev Urol 12 (2015), 317–330.

Yan

Guo

Z.X.

Liu

X.P.

Feng

Y.J.

Zhao

Y.J.

Liu

T.Z.

and Li

, Four novel biomarkers for bladder cancer identified by weighted gene coexpression network analysis, J Cell Physiol (2019).

Yamasaki

Mukai

Nagai

Nakahara

Fujii

Terada

Ohno

Sato

Toda

Kataoka

and Kamoto

, Matriptase-induced phosphorylation of MET is significantly associated with poor prognosis in invasive bladder cancer; an immunohistochemical analysis, Int J Mol Sci 19 (2018).

Riester

Taylor

J.M.

Feifer

Koppie

Rosenberg

J.E.

Downey

R.J.

Bochner

B.H.

and Michor

, Combination of a novel gene expression signature with a clinical nomogram improves the prediction of survival in high-risk bladder cancer, Clin Cancer Res 18 (2012), 1323–1333.

Kim

W.J.

Kim

E.J.

Kim

S.K.

Kim

Y.J.

Y.S.

Jeong

Kim

M.J.

Yun

S.J.

Lee

K.M.

Moon

S.K.

Lee

S.C.

Cha

E.J.

and Bae

S.C.

, Predictive value of progression-related gene classifier in primary non-muscle invasive bladder cancer, Mol Cancer 9 (2010), 3.

Ritchie

M.E.

Phipson

Law

C.W.

Shi

and Smyth

G.K.

, limma powers differential expression analyses for RNA-sequencing and microarray studies, Nucleic Acids Res 43 (2015), e47.

10.

Afsari

Geman

and Fertig

E.J.

, Learning dysregulated pathways in cancers from differential variability analysis, Cancer Inform 13 (2014), 61–67.

11.

Parrish

R.S.

and Spencer

H.J.

, 3rd, Effect of normalization on significance testing for oligonucleotide microarrays, J Biopharm Stat 14 (2004), 575–589.

12.

Wang

Cao

Zeng

Wang

Cheng

Zhu

Nian

and Wang

, RNA-seq analyses of multiple meristems of soybean: Novel and alternative transcripts, evolutionary and functional implications, BMC Plant Biol 14 (2014), 169.

13.

Langfelder

and Horvath

, WGCNA: An R package for weighted correlation network analysis, BMC Bioinformatics 9 (2008), 559.

14.

Szklarczyk

Morris

J.H.

Cook

Kuhn

Wyder

Simonovic

Santos

Doncheva

N.T.

Roth

Bork

Jensen

L.J.

and von Mering

, The STRING database in 2017: Quality-controlled protein-protein association networks, made broadly accessible, Nucleic Acids Res 45 (2017), D362–D368.

15.

Shannon

Markiel

Ozier

Baliga

N.S.

Wang

J.T.

Ramage

Amin

Schwikowski

and Ideker

, Cytoscape: A software environment for integrated models of biomolecular interaction networks, Genome Res 13 (2003), 2498–2504.

16.

da Huang

Sherman

B.T.

and Lempicki

R.A.

, Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources, Nat Protoc 4 (2009), 44–57.

17.

Davis

A.P.

Grondin

C.J.

Johnson

R.J.

Sciaky

McMorran

Wiegers

T.C.

and Mattingly

C.J.

, The comparative toxicogenomics database: Update 2019, Nucleic Acids Res 47 (2019), D948–D954.

18.

Wang

Hang

Zou

and Mao

J.H.

, A novel gene expression-based prognostic scoring system to predict survival in gastric cancer, Oncotarget 7 (2016), 55343–55351.

19.

van der Rest

and Garrone

, Collagen family of proteins, FASEB J 5 (1991), 2814–2823.

20.

Gilkes

D.M.

Semenza

G.L.

and Wirtz

, Hypoxia and the extracellular matrix: drivers of tumour metastasis, Nat Rev Cancer 14 (2014), 430–439.

21.

Kalluri

, Basement membranes: structure, assembly and role in tumour angiogenesis, Nat Rev Cancer 3 (2003), 422–433.

22.

Zhang

Teng

Liu

and Zhang

, Gene expression profile analyze the molecular mechanism of CXCR7 regulating papillary thyroid carcinoma growth and metastasis, J Exp Clin Cancer Res 34 (2015), 16.

23.

Miyake

Hori

Morizawa

Tatsumi

Toritsuka

Ohnishi

Shimada

Furuya

Khadka

V.S.

Deng

Ohnishi

Iida

Gotoh

Nakai

Inoue

Anai

Torimoto

Aoki

Tanaka

Konishi

and Fujimoto

, Collagen type IV alpha 1 (COL4A1) and collagen type XIII alpha 1 (COL13A1) produced in cancer cells promote tumor budding at the invasion front in human urothelial carcinoma of the bladder, Oncotarget 8 (2017), 36099–36114.

24.

Egeblad

Rasch

M.G.

and Weaver

V.M.

, Dynamic interplay between the collagen scaffold and tumor evolution, Curr Opin Cell Biol 22 (2010), 697–706.

25.

Brooks

Krasnow

P.L.

Lee

Y.C.

Xiao

Kurtova

Lerner

Godoy

Jian

Castro

Chen

Rowley

Ittmann

and Chan

K.S.

, Positive association of collagen type I with non-muscle invasive bladder cancer progression, Oncotarget 7 (2016), 82609–82619.

26.

Chen

and Yan

, Bladder cancer stage-associated hub genes revealed by WGCNA co-expression network analysis, Hereditas 156 (2019), 7.

27.

Jin

Yao

Xie

Cai

Zhang

and Li

, Melittin constrains the expression of identified key genes associated with bladder cancer, J Immunol Res 2018 (2018), 5038172.

28.

Liu

P.F.

Cao

Y.W.

Jiang

H.P.

Wang

Y.H.

Yang

X.C.

Wang

X.S.

and Niu

H.T.

, Heterogeneity research in muscle-invasive bladder cancer based on differential protein expression analysis, Med Oncol 31 (2014), 21.

29.

Bhat

Heinzel

Mayer

Perco

Muhlberger

Husi

Merseburger

A.S.

Zoidakis

Vlahou

Schanstra

J.P.

Mischak

and Jankowski

, Protein interactome of muscle invasive bladder cancer, PLoS One 10 (2015), e0116404.

Identifying the mRNAs associated with Bladder cancer recurrence

Abstract

OBJECTIVE:

METHODS:

RESULTS:

CONCLUSION:

Keywords

1. Introduction

2. Materials and methods

2.1 Microarray data

2.2 The dataset preprocessing and DEGs analysis

2.3 WGCNA

Table 1 The top ten up-regulated and the down-regulated DEGs between recurrence and non-recurrence bladder cancer in the two datasets

2.5 BC related DEGs screening

2.6 Recurrence related DEGs screening

3.1 DEGs analysis

3.2 Modules and genes screened by WGCNA

3.3 PPI network analysis

Table 2 The GO and KEGG pathway enrichment for the DEGs in the PPI network

Table 3 Clinical factor statistics and cox regression analysis

Footnotes

Conflict of interest

Supplementary data

References

Table 1
The top ten up-regulated and the down-regulated DEGs between recurrence and non-recurrence bladder cancer in the two datasets

Table 2
The GO and KEGG pathway enrichment for the DEGs in the PPI network

Table 3
Clinical factor statistics and cox regression analysis