A novel prognostic biomarker for muscle invasive bladder urothelial carcinoma based on 11 DNA methylation signature

ABSTRACT

Muscle-invasive bladder urothelial carcinoma (MIBC) is a highly invasive cancer, which leads to prevalent recurrence and poor prognosis. Exploring the association of DNA methylation and the prognosis of MIBC will thus be of important value in clinical management and treatment. Bumphunter method and adaptive lasso regression were used to explore the relationship between different methylation regions (DMRs) and the prognosis of MIBC. Next, we constructed a risk prognosis model and validated this model. Moreover, the performance of this risk model was examined by using time-dependent receiver operating characteristic curve (ROC). We identified 58,449 different methylation sites and 490 different methylation regions. Among them, 11 DMRs were associated with the prognosis of MIBC through rigorous screening. Through the linear combination of 11 DMRs, a putative marker was developed, which can distinguish the survival risk in both the training dataset (HR = 2.58, 95% CI = (1.64, 4.05)) and the verification dataset (HR = 2.77, 95% CI = (1.25, 6.15)). Relatively high predictive values were observed from this model for training dataset (AUC = 0.791) and verification dataset (AUC = 0.668). Stratified analysis showed that the association was independent of gender. A nomogram was additionally generated to predict 5-year survival probability containing risk score and pathological stage. Its performance was evaluated by applying calibration curve. The methylation signature risk model based on 11 DMRs may be a reliable prognostic signature for MIBC, which provides new insights into development of individualized therapy for MIBC.

KEYWORDS:

• Muscle invasive bladder urothelial carcinoma (MIBC)
• methylation
• nomogram
• prognostic
• risk model

Introduction

Bladder cancer is a highly prevalent cancer arising from urothelial epithelium. 95% of bladder cancer are urothelial carcinoma, also called transitional cell carcinoma.¹ On average, 430,000 incident bladder cancer cases and 165,000 deaths are reported approximately each year worldwide.² Among them, about 25% of new cases are diagnosed as muscle invasive bladder urothelial carcinoma (MIBC), and 30–40% of nonmuscle invasive bladder urothelial carcinoma (NMIBC) will gradually infiltrate the muscular wall of the bladder to progress to MIBC.³ The pathological grade of MIBC is usually higher and is more prone to distal metastasis, which brings about an unfavorable 5-year survival of being 30–50%.⁴

In clinical practice, radical cystectomy of bladder in conjunction with neo-adjuvant therapy such as cisplatin is still the main treatment for MIBC.⁵ But the response rate to cisplatin-based therapy is below 50%, and 63% of patients are forced to adjust doses because of toxicity from chemotherapy.^6,7 As a result, the survival rate of MIBC has not improved in the past 30 years.⁸ In consideration of the highly invasive nature and the genetic variability of MIBC, identifying important risk markers associated with MIBC prognosis at the molecular level will thus be of great value in exploring new individualized therapies and improving survival rates.

DNA methylation is a kind of epigenetic modification of human genome, which is related to several key processes, for example, genomic imprinting, repression of transposable elements, carcinogenesis, and so on.^9,10 Association of DNA methylation and development of cancer has been reported in previous studies.^11–13 However, few studies have examined the role of a single DNA methylation site as prognostic and metastatic markers for MIBC among limited sample size, which hindered the generalization of findings to some extent.^14,15 While existing evidences support that some functional changes are more related to the genomic region than to a single CpGs (5ʹ-C-phosphate-G-3ʹ) site.¹⁶ Therefore, more attentions are suggested to be given to the abnormal methylation of gene region.

In this study, we used bumphunter algorithm to analyze the level 3 methylation data of 419 samples from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/) in order to screen differential methylation regions (DMRs) between tumor and normal tissues. Among them, we identified those DMRs related to overall survival, constructed a risk score model based on 11 DMRs in training dataset, and verified the model in validation dataset. Further, a time-dependent receiver operating characteristic (ROC) curve was applied to evaluate the predictive value of the 5-year survival rate of the model.

Material and methods

Data collection

We downloaded level 3 methylation data of 421 muscle invasive bladder urothelial carcinoma (MIBC) samples (Illumina Infinium Human DNA Methylation 450 platform, beta values) and the corresponding clinical information (GDC) (Bethesda, MD, USA). Among them, 400 were tumor samples and 21 were adjacent normal tissue samples. Meanwhile, batch information for methylation data were downloaded from MD Anderson Cancer Center.

Data preprocessing

The TCGA level 3 methylation data had been converted to a β value, representing the degree of methylation at the CpGs site. First, we used the ChIPseeker¹⁷ R package to re-annotate the location of all probes based on the GRCh38 annotation file and the promoter region defined as upstream 2.0 kb to downstream 0.5kb of the transcriptional start site. Further, we deleted the XY chromosome probes, the SNP probes, the non-CpG site probes, the multiple matching probes, and the nonpromoter region probes. Meanwhile, we eliminated probes with missing values greater than 5% and applied the ChAMP pipeline to do the missing value filling, quality control, and BMIQ¹⁸ normalization.¹⁹ Subsequently, we used ComBat²⁰ to correct for potential confounders, such as time and sampling location. Finally, 399 tumor samples and 20 normal samples with a total of 169,730 probes were included in the follow-up analysis. Also, we categorized the age into two groups (<60 and ≥60), divided pathological Stage into Stage II group and Stage III and IV group, and limited the follow-up time to 5 years (60 months) to investigate their effects on the prognosis.

Screening differential methylation regions

In order to reduce the heteroscedasticity, the β value was transformed into M value according to the formula $\mathrm{M}=\mathrm{l}\mathrm{o}\mathrm{g}2\left(\frac{\beta }{1-\beta }\right)$.²¹ First, the empirical Bayes method in limma R package was applied to screen the differentially methylated probes (DMPs) in both tumor and normal tissues.²² Second, we identified tumor and normal differential methylation regions (DMRs) using bumphunter method,²³ in which clusters were defined by neighboring CpGs site of at an interval at most 500bp and 1,000 bootstrap samples were used to generate a null distribution of regions. Meanwhile, we used the 95% quantile as the cutoff value for the candidate regions and each candidate region had at least 7 CpGs probes. In addition, the false positive rate was controlled by using the family-wise error rate (FWER) in consideration for its conservation compared with the BH and FDR methods. Third, we annotated all DMRs gene coordinates using wANNOVA.²⁴ Forth, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were carried out to understand the underlying pathways of DMRs, using “clusterProfiler” R package²⁵ based on a threshold of P < .05.

Model development and validation

The methylation level of each DMR was expressed by the average β value of the probes contained in the DMR. MIBC patients were randomly divided into a training dataset (n = 319) and a validation dataset (n = 80). First, we used univariate cox regression to screen DMRs in training dataset and found 20 DMRs that were correlated with overall survival. Second, adaptive lasso penalty cox regression (ALasso)²⁶ which avoided the problem of over fitting was performed to filter the important DMRs where the λ value was obtained by a 10-fold cross validation. Third, the risk model was established according to the formula $\mathrm{R}\mathrm{i}\mathrm{s}\mathrm{k}={\beta }_{1}DMR{s}_1+{\beta }_{2}DMR{s}_2+{\beta }_{3}DMR{S}_3+\dots .+{\beta }_{n}DMR{s}_n$ (β is the ALasso cox regression coefficient and DMRs is the methylation level of DMR). Forth, patients were divided into high and low risk groups using the optimal cutoff value of the risk score, which was determined by the X‐tile software (version 3.6.1, Yale University, CT). Meanwhile, we validated the risk model in the validation dataset using the same cutoff value. Fifth, we conducted a time-dependent receiver operating characteristic (ROC) curve to evaluate the predictive value of the five-year survival rate of the model using the survival ROC package.²⁷ In addition, total dataset was separated by gender and Z test was used to compare the predictive value of the risk model between genders. Finally, we plotted nomogram based on the model-based risk and the participants’ clinical characteristics. Kaplan-Meier and log-rank method were used to test the survival differences of different risk groups.

Results

Patient characteristics

A total of 399 patients were included in this study after preprocessing and they were randomly divided into training datasets (n = 319) and validation datasets (n = 80). The demographic and clinical characteristics of both datasets are listed in Table 1. The median survival time were 35 months in training dataset, 22 months in validation dataset, and 33 months for all patients combined, respectively. Age, gender, race and pathological grade were comparable between training dataset and validation dataset (P > .05). However, the distribution of participants were significantly different by pathological stage between two datasets, with more high stages being in validation dataset (P = .013).

Table 1. Characteristics of training dataset and validation dataset

		Training dataset		Validation dataset		Total
Characteristics		No	%	No	%	No	%
All		319	100.0	80	100.0	399	100.0
Age
	<60	67	21.0	19	23.8	86	21.6
	≥60	252	79.0	61	76.2	313	78.4
Gender
	Male	228	71.5	64	80.0	292	73.2
	Female	91	28.5	16	20.0	107	26.8
Race
	White	253	79.3	65	81.3	318	79.7
	Asian	38	11.9	4	5.0	42	10.5
	Black OR African American	15	4.7	8	10.0	23	5.8
	Unknown	13	4.1	3	3.7	16	4.0
Stage ^a
	Stage II	110	34.5	16	20.0	126	31.6
	Stages III and IV	209	65.5	64	80.0	273	68.4
Grade
	Low Grade	19	6.0	1	1.3	20	5.0
	High Grade	300	94.0	79	98.7	379	95.0
Survival
	Death	129	40.4	40	50.0	169	42.4
	Alive	190	59.6	40	50.0	230	57.6
	MST ^b	35		21.9		33.1

a. Stage = Pathological Stage.

b. MST = Median Survival Time (Months)

Identification of differential methylation regions

In this study, 58,449 differentially methylated probes (DMPs) were identified from 169,730 probes according to the criteria of “P (Benjamini adjusted P value) < 0.05.” Bumphunter algorithm was performed to analyze these DMPs in 399 MIBC samples and 20 normal samples. A total of 490 DMRs had been found based on the criteria of FWER being less than 0.05 and the number of probes being more than 7. KEGG pathway analysis showed that these DMRs were mainly involved in cAMP signaling pathway, Calcium signaling pathway, and transcriptional misregulation in cancer. (Figure 1).

Figure 1. Enrichment of KEGG paths with significant P values. The x-axis represents the number of DMRs involved in the pathway

Risk model construction

We first used univariate Cox regression to evaluate the association of DMRs with overall survival in the training dataset, which generated a total of 23 differential DMRs (Supplementary table). Next, we used 10-fold cross-validation adaptive lasso cox regression to screen out stable DMRs combinations from 23 DMRs, from which 11 DMRs were identified and annotated, namely, ARHGDIB, SHANK2, PROM1, RASGRF2, SEC31B, PXDNL, TRIM27, GPR75-ASB3, GALR1, PP12613, and NRN1 (Supplementary table). Furthermore, we applied the ALasso cox regression coefficient to build a risk model as followed,

$Risk=0.2907\times DMRs1-0.7774\times DMRs2-0.8036\times DMRs3+0.8632\times DMRs4-1.8386\times DMRs5-1.6511\times DMRs6+0.2448\times DMRs7-0.69\times DMRs8+1.3604\times DMRs9+0.3227\times DMRs10+0.6652\times DMRs11.$

And the risk score was calculated for each patient. In the training dataset, patients were categorized into high-risk group when a risk score was> −2.69 and into low-risk group when a risk score was ≤ −2.69, using the optimal cutoff value determined by X-lite software. A worse overall survival was observed for the high-risk group (HR = 2.58, 95% CI = (1.64, 4.05)) with reference to the low-risk group after multiple adjustment (Table 2). The survival probability was depicted in the Kaplan-Meier survival curve plot (Figure 2a), In addition, we used the time-dependent ROC curve to evaluate the prediction accuracy of risk model in the training dataset, and the calculated area under the ROC curve was 0.791(Figure 2b).

Table 2. Survival analysis for risk module with OS as endpoints in the training dataset and validation dataset

	Training dataset (n = 319)						Validation dataset (n = 80)
	Univariate			Adjusted ^a			Univariate			Adjusted ^a
Variable	HR	95% CI	P	HR	95% CI	P	HR	95% CI	P	HR	95% CI	P
Age ^b	2.05	1.22–3.46	0.007	1.53	0.90–2.60	0.120	1.99	0.83–4.75	0.121	1.49	0.59–3.73	0.400
Stage ^c	2.35	1.53–3.62	<.001	2.02	1.30–3.12	0.002	2.81	0.99–7.91	0.050	2.97	1.04–8.46	0.041
Gender ^d	0.90	0.62–1.30	0.565	0.90	0.62–1.32	0.591	0.92	0.44–1.92	0.815	0.70	0.32–1.51	0.360
Risk^e	2.99	1.91–4.66	<.001	2.58	1.64–4.05	<.001	2.87	1.36–6.06	0.006	2.77	1.25–6.15	0.012

a. Adjusted for Age, Pathological Stage, Gender b. ≥60 vs <60 (Ref)

c. Stages III and IV vs Stage II (Ref).

d. Female (Ref) vs Male.

e. High Risk vs Low Risk (Ref).

Figure 2. Kaplan-Meier curves and time-dependent receiver operating characteristic curves of methylation signature predicting the overall survival of MIBC. (a), (b) in training dataset. (c), (d) in validation dataset

Risk model validation

In order to determine the stability of the risk model, we used the validation dataset to evaluate the constructed risk model. Verification dataset were divided into two groups, that is, high risk (n = 50) and low risk (n = 30) according to the same cutoff values as defined in the training dataset. An unfavorable overall survival were observed for high-risk group before (HR = 2.87, 95% CI = (1.36, 6.06)) and after (HR = 2.77, 95% CI = (1.25, 6.15)) multiple adjustment in comparison with low-risk group (Table 2, Figure 2c). The prediction accuracy of the risk model was assessed (AUC = 0.668, Figure 2d).

Subgroup analyses

In order to observe the impact of gender on the prognosis of MIBC, we evaluated the relative risk of high-risk group compared with low-risk group for total dataset (HR = 2.56, 95% CI = (1.73, 3.78)), male dataset (HR = 2.72, 95% CI = (1.71, 4.34)), and female dataset (HR = 2.18, 95% CI = (1.07, 4.44)), respectively (data not shown). The prediction accuracy (AUC value) was basically similar for different groups (0.762 for total dataset, 0.775 for male, and 0.718 for female, P = .441, Figure 3)

Figure 3. Kaplan-Meier curves (a) in Total dataset, (b) in Female dataset, (c) in Male dataset and time-dependent receiver operating characteristic curves (d) of methylation signature predicting the overall survival of MIBC

Nomogram construction based on risk model and clinical covariates

In the training dataset, we found that pathologic stage and risk models were associated with overall survival in both univariate and multivariate analyses. Therefore, we developed a five year survival probability nomogram based on risk model and pathological stage (Figure 4).

Figure 4. Methylation signature nomogram to predict the risk of survival in patients with MIBC. (a) Nomogram of methylation signature predict the 5-years Survival probability of MIBC. (b) Calibration curve for the OS nomogram model in training dataset. (c) Calibration curve for the OS nomogram model in validation dataset

Discussion

The biological mechanism of bladder cancer has been widely studied for decades. Initially, most researches focused on the expression and biological function of a single gene.²⁸ With the introduction of high-throughput technology, researchers began to shift their attention to the association of DNA methylation with carcinogenesis.²⁹ Methylation is a dynamic chemical modification of DNA cytosine, which regulates biological development by inhibiting the transcription of genes.³⁰ Studies had shown that abnormal methylation of tumor suppressor gene and proto-oncogene promoter region played an important role in cancer occurrence and development.^11,13,31 In this study, we used publicly available data warehouses to perform genome-wide DNA methylation analysis at the gene region level and correlated it with clinical information to assess the relationship between aberrant methylation regions and the prognosis of MIBC.

In our study, we applied bumphunter method to compare methylation differences between tumor and normal samples and we identified 490 DMRs from a total of 58,449 DMPs. Further, univariate Cox regression was used to evaluate the relationship between each DMR and the prognosis of MIBC. In order to avoid collinearity and overfitting for our model, we therefore screened the DMRs that were strongly correlated with the prognosis of MIBC by performing the adaptive lasso method, an improvement of the traditional Lasso model for having overcome the limitation that the Oracle property was not satisfied.

Previous studies had explored the potential mechanism underlying the correlation between the identified single CpGs site with the unfavorable prognosis.^14,15 We, here, identified a serial of DMRs which had been associated with carcinogenesis and cancer prognosis including ARHGDIB, SHANK2, PROM1, RASGRF2, TRIM27, GPR75-ASB3, GALR1, and NRN1. To date, direct evidences with carcinogenesis have not been reported for other DMRs, namely, SEC31B, PXDNL, and PP12613.

Among these markers identified, several had been related to the development of nervous system as well as nervous system tumorigenesis, recurrence, and prognosis. NRN1 (neuritin 1) encodes a protein of the neuritin family, which plays an important role in the development of the nervous system, the proliferation of stem cells, and the aging process. It had been shown that NRN1 played a key role in the occurrence and development of cancer by promoting the morphological changes of cells.^32,33 PROM1 (prominin 1) maintains the properties of stem cells by inhibiting differentiation through encoding a pentaspan transmembrane glycoprotein. It was proved that the high expression of PROM1 was associated with the recurrence, metastasis and poor prognosis of nervous system tumors.^34,35

GALR1 (galanin receptor 1), ARHGDIB (Rho GDP dissociation inhibitor beta) and SHANK2 (SH3 and multiple ankyrin repeat domains 2) has been reported to be associated with the prognosis of head and neck cancer, muscle-infiltrating bladder and prostate cancer, and oral squamous cell carcinoma, respectively. GALR1 exerts biological effects by interacting with specific G protein-coupled receptors. Studies have shown that aberrant methylation in the promoter region of GALR1 is associated with the occurrence, development, and poor prognosis of head and neck cancer. It may be that the hyper-methylation of GALR1 promoter region inhibits the expression of GALR1 mRNA, which leads to the loss of control of the expression of cyclin D1.³⁶ In addition, it has been reported that GALR1 affects the sensitivity of patients to chemotherapeutic drugs by participating in the mitogen-activated protein kinase (MAPK) signaling.³⁷ ARHGDIB is a member of Rho (or ARH) protein family, which is involved in cell signaling, proliferation, cytoskeletal organization, and secretion.³⁸ ARHGDIB was reported to prevent the metastasis of cancer by inhibiting the activation of Rho GTP enzyme.³⁹ Several studies from the United States have confirmed that the low expression of ARHGDIB is associated with metastasis and poor prognosis in muscle-infiltrating bladder and prostate cancer.^40–42 SHANK2 was reported to encode a synaptic protein in the Shank family that was involved in the development of the nervous system.⁴³ It has been also reported that SHANK2 promotes the invasion of cancer cells through synergistic action with cortactin, which is a protein involved in the cytoskeleton, and its high expression was related to the poor prognosis of oral cancer.^44,45

Several other markers were previously observed to be relevant to tumorigenesis in different mechanisms. RASGRF2 (Ras protein specific guanine nucleotide releasing factor 2) encodes a calcium-regulated nucleotide exchange factor, which plays a key role in invasion, colonization, metastasis, and diffusion of tumor cells through the negative regulation of CDC42 signaling pathway.⁴⁶ TRIM27 (tripartite motif containing 27), a transcriptional inhibitor, mediates the regulation of thioredoxin binding protein 2 (TBP-2) by histone deacetylases (HDAC), which affects the susceptibility of cancer cells to oxidative stress and cisplatin.⁴⁷ GPR75-ASB3 and ABS3 encode the same protein by downregulating tumor necrosis factor R2 (TNF-R2) ubiquitin and degradation and affect tumor necrosis factor α (TNF-α) cell response and JNK signaling pathway, which is involved in the process of tumor cell proliferation, apoptosis, and metastasis.⁴⁸ To our knowledge, their association with cancer prognosis has not reported so far.

To date, SEC31B, PXDNL, and PP12613 methylation abnormality of promoter regions has not been reported to be related to carcinogenesis. We, among the first, observed in this study that aberrant methylation of SEC31B, PXDNL, and PP12613 promoter regions was associated with the prognosis of MIBC. Until recently, we know that SEC31B (SEC31 homolog B) encodes a protein of unknown function that mediates vesicular traffic from the endoplasmic reticulum, PXDNL (peroxidasin like) encodes a peroxidase that participates in basement membrane synthesis,⁴⁹ and PP12613 is a noncoding gene. The unknown mechanism underlying the process of carcinogenesis and pathological differentiation warrants to be further studies.

Compared with the study of Robertson et al, our study has the advantage of using multiple CpG sites over a single CpG site in the promoter region of the gene to better measure the association of the methylation level of that gene with MIBC survival. Meanwhile, we developed a new clinical survival risk prediction tool.⁵⁰ We also noted that Xie and colleagues had proposed a statistical method named GDP (Group lass regularized Deep learning for cancer Prognosis) which seems to be more efficient in some ways.⁵¹ However, disadvantages existed for their method compared to adaptive lasso used in our study, such as, unsatisfied Oracle property, lack of clinical prior group knowledge and low sparsity within the group. We therefore constructed a new prognostic risk model based on gene methylation by used adaptive lasso method. Furthermore, some studies have shown that the prognosis of MIBC was different between genders and the worse prognosis was usually among women.⁵² We therefore performed a gender-based subgroup analysis, and our results suggested that the risk model based on 11 DMRs was independent of gender. However, this study has its limitations. First, our data were abstracted from the public database, where treatment information was not available. Furthermore, lack of information about cancer recurrence prevents us from exploring the relationship between risk models and cancer recurrence. In addition, limited by the data acquisition, we only could verify the model internally, which limits the extrapolation of results.

Conclusion

We analyzed methylation data and clinical data from patients with MIBC. 11 abnormal methylation regions were identified and a prognostic risk model based on methylation markers was constructed and internally validated. Simultaneously, nomograms including methylated markers and clinical risk factors were established to predict the survival probability, which is believed to be valuable for the development of individualized treatment of MIBC.

Author contributions

Conception and design: Yueyi Feng, and Xiaochen Shu

Financial support: Xiaochen Shu

Collection and assembly of data: Yueyi Feng, Yiqing Jiang, Qingting Feng, Lingkai Xu, Yun Jiang

Data analysis and interpretation: Yueyi Feng, Yiqing Jiang, Qingting Feng, Lingkai Xu, Yun Jiang, Fang Meng, and Xiaochen Shu

Manuscript writing: Yueyi Feng, Yiqing Jiang, Fang Meng, and Xiaochen Shu

Final approval of manuscript: Yueyi Feng, Yiqing Jiang, Qingting Feng, Lingkai Xu, Yun Jiang, Fang Meng, and Xiaochen Shu

Disclosure of potential conflicts of interest

No potential conflicts of interest were disclosed.

Acknowledgments

The corresponding author would be grateful for co-workers’ contributions including collecting, managing, maintaining, and analyzing data as well as writing and proofing the manuscript submitted.

Supplementary material

Supplemental data for this article can be accessed on the publisher’s website.

Additional information

Funding

The present work was conducted by the support of Scientific Research Foundation for Talented Scholars in Soochow University, China (Q413900215).