Skip to main content

Immune profiles and DNA methylation alterations related with non-muscle-invasive bladder cancer outcomes



Non-muscle-invasive bladder cancer (NMIBC) patients receive frequent monitoring because ≥ 70% will have recurrent disease. However, screening is invasive, expensive, and associated with significant morbidity making bladder cancer the most expensive cancer to treat per capita. There is an urgent need to expand the understanding of markers related to recurrence and survival outcomes of NMIBC.

Methods and results

We used the Illumina HumanMethylationEPIC array to measure peripheral blood DNA methylation profiles of NMIBC patients (N = 603) enrolled in a population-based cohort study in New Hampshire and applied cell type deconvolution to estimate immune cell-type proportions. Using Cox proportional hazard models, we identified that increasing CD4T and CD8T cell proportions were associated with a statistically significant decreased hazard of tumor recurrence or death (CD4T: HR = 0.98, 95% CI = 0.97–1.00; CD8T: HR = 0.97, 95% CI = 0.95–1.00), whereas increasing monocyte proportion and methylation-derived neutrophil-to-lymphocyte ratio (mdNLR) were associated with the increased hazard of tumor recurrence or death (monocyte: HR = 1.04, 95% CI = 1.00–1.07; mdNLR: HR = 1.12, 95% CI = 1.04–1.20). Then, using an epigenome-wide association study (EWAS) approach adjusting for age, sex, smoking status, BCG treatment status, and immune cell profiles, we identified 2528 CpGs associated with the hazard of tumor recurrence or death (P < 0.005). Among these CpGs, the 1572 were associated with an increased hazard and were significantly enriched in open sea regions; the 956 remaining CpGs were associated with a decreased hazard and were significantly enriched in enhancer regions and DNase hypersensitive sites.


Our results expand on the knowledge of immune profiles and methylation alteration associated with NMIBC outcomes and represent a first step toward the development of DNA methylation-based biomarkers of tumor recurrence.


In 2021, the estimated number of bladder cancer deaths is projected to be 17,200, with an estimated number of new cases of 83,730 in the USA. Bladder cancer is the fourth most common cancer among men and twelfth most common among women, which may in part be due to smoking prevalence rates being higher in men than women as cigarette smoking accounts for half of all cases (47%) in the USA. [1]. Seventy-five percent of bladder cancers are diagnosed as low-grade non-muscle invasive tumors, NMIBC [2]. Cystoscopy is used for diagnosis (biopsy) with transurethral excision of localized tumors as the primary treatment [3,4,5]. Although transurethral excisions can successfully control the disease and mortality from bladder cancer among patients with localized tumors is low, 45% of NMIBC cases have recurrences within 12 months of surgery [6]. In addition, frequent invasive follow-up via cystoscopy without prognostic markers leads to significant patient morbidity and comes with a considerable cost burden to the health care system, estimated at approximately $4 billion dollars annually in the USA [7]. To control the patient and healthcare burdens associated with NMIBC, there is an imminent need for biomarkers to identify those at the highest risk of tumor recurrence.

Peripheral blood immune profiles have been associated with different outcomes in bladder cancer patients and may have clinical utility for NMIBC prognosis [8,9,10,11,12,13]. For instance, in NMIBC, patients with elevated neutrophil-to-lymphocyte ratio (NLR) had poorer cancer-specific survival than patients with lower NLR [9,10,11]. Thus, elevated NLR could be a potential predictor of overall survival and cancer-specific survival for this disease. Other studies have indicated that bladder cancer patients with increased lymphocyte-to-monocyte ratio had poorer overall survival and cancer-specific survival [12, 13]. In addition, Bacillus Calmette–Guérin (BCG), a commonly used intravesical immunotherapy for NMIBC administered post-surgery, has been reported to reduce the proportion of natural killer T cells, memory CD4T, CD8T, and regulatory T cells in the peripheral blood of NMIBC patients [8]. Together, circulating immune profiles might be promising markers for reducing adverse outcomes in NMIBC patients. Previous studies have relied on complete blood count differential (CBC) tests to determine immune profile variables [14, 15]. The CBC test relies on fresh blood samples and is incapable of providing proportions of specific lymphocyte subtypes [16].

Our prior work has established methods to infer the immune profiles in archival samples using immune cell-type-specific DNA methylation [17, 18]. DNA methylation plays an essential role in gene regulation for cell lineage specification [19, 20]. Differentially methylated regions (DMRs) have been used to distinguish cell types, including leukocyte subtypes, and form the basis of reference-based deconvolution methods for estimating specific immune cell-type proportions [18, 21, 22]. Compared with cytological methods for determining cell type abundances and proportions, such as flow cytometry, DNA methylation-based cell-type deconvolution does not require a fresh substrate, intact cells, or batch-sensitive reagents, is reproducible and cost-effective relative to time-sensitive blood processing [23, 24]. With cell-type deconvolution approaches, it is possible to identify immune cell-type profiles and test their relationship with cancer outcomes in archived samples [25]. This study used the archival blood samples of NMIBC patients to test the association between the immune profiles and outcomes in bladder cancer patients.

In the present study, we sought to identify immune profiles and epigenetic features associated with disease recurrence in the hope that such information might help improve the management of NMIBC. Here, we hypothesized that CpG-specific DNA methylation and DNA methylation-derived immune cell profiles are associated with recurrence-free survival in NMIBC patients. We used archival blood samples from a population-based case–control study to obtain genome-scale DNA methylation profiles. We then investigated the association between the methylation-derived immune profiles and outcomes in NMIBC patients. Preliminary work from our group observed an association between methylation-derived NLR (mdNLR) and survival in bladder cancer patients using a smaller sample size (223 cases) and an early genome-scale methylation array (HumanMethylation27K array) [26]. In this study, we increased the sample size (603 cases), used a new more comprehensive array (HumanMethylationEPIC array) with 30 times as many measured features, and a new cell type deconvolution library [18] to estimate immune cell-type proportions. An epigenome-wide association study (EWAS) and enrichment analyses were used to determine possible CpG sites and gene sets associated with recurrence-free survival.


Profiles of DNA methylation were obtained from 685 peripheral blood samples using the Human MethylationEPIC array. Eighty-two subjects were excluded due to low-quality CpG value or bisulfite intensity (n = 11), or without muscle-invasive status, histopathology re-review, tumor grade, smoking status, and pack-years (n = 71) (Fig. 1). The remaining subjects (N = 603) included in the study group were 75.8% men (n = 457), 82.9% ever-smokers (n = 500), and had a median age of 66 (Table 1). We estimated the cell-type proportions for each patient using methylation values by performing FlowSorted.Blood.EPIC (see Additional file 7: Figure S1 for the distribution). Neutrophil-to-lymphocyte ratio (NLR) was then calculated according to the ratio of neutrophil proportion to lymphocyte proportion (B cell + CD4T cell + CD8T cell + NK cell), and the median methylation-derived NLR (mdNLR) was 1.97. Further details of study population characteristics are described in Table 1.

Fig. 1
figure 1

Flow chart of study

Table 1 Characteristics of subjects after excluding subjects with missing values

To characterize 10-year recurrence-free survival (RFS), we generated Kaplan–Meier curves for each covariate and fit a Cox proportional hazard regression model for univariate analyses and multivariable analyses, respectively. In the Kaplan–Meier analysis, the NMIBC patients aged > 65 had a worse probability of RFS compared with the NMIBC patients with age ≤ 65 (P = 0.0006). Females had a greater probability of RFS compared with males (P = 0.002), and the NMIBC patients with high-grade tumors (grade 3 + 4) had a lesser probability of RFS than those with low-grade tumors (grade 1 + 2) (P = 8.0 × 10−5). Ever-smokers had a worse probability of RFS than the NMIBC patients who were never smokers (P = 3.0 × 10−4). NMIBC patients with low mdNLR had a greater probability of RFS than the patients with high mdNLR (P = 0.002) (Fig. 2). Consistent with the Kaplan–Meier results, in a multivariable Cox model, age > 65 (HR = 1.01, 95% CI = 1.00–1.03), high tumor grade (HR = 1.48, 95% CI = 1.17–1.87), ever-smoking (HR = 1.65, 95% CI = 1.22–2.25), and mdNLR (HR = 1.12, 95% CI = 1.04–1.20) were significantly associated with an increased hazard of RFS (Table 2).

Fig. 2
figure 2

Kaplan–Meier analysis of 10-year recurrence-free survival (RFS). 10-year RFS curves stratified by A age, B sex, C tumor grade, D smoking status, E BCG treatment status or F mdNLR level. P-values for log-rank tests are shown

Table 2 Cox proportional hazards 10-year recurrence-free survival models

We next investigated the association between immune cell-type proportions and RFS in multivariable models. CD4T (HR = 0.98, 95% CI = 0.97–1.00) and CD8T cell proportion (HR = 0.97, 95% CI = 0.95–1.00) were significantly associated with the decreased hazard of RFS. Monocyte cell proportion (HR = 1.04, 95% CI = 1.00–1.07) was significantly associated with the increased hazard of RFS. Both B cell and NK cell proportion hazard estimates were < 1 but not statistically significant (Additional file 1). We also examined the 10-year overall survival (OS) in univariate models and the multivariable models. Observed associations with RFS were consistent for OS for: age, male, high tumor grade, ever-smoking, and mdNLR (Table 3 and Additional file 7: Figure S2). In addition, CD4T, CD8T, B cell, and NK cell proportion were significantly associated with the decreased hazard of death; neutrophil cell proportion was significantly associated with the increased hazard of death (Additional file 2).

Table 3 Cox proportional hazards 10-year overall survival models

Next, we assessed the relation of NMIBC patient outcomes with DNA methylation. First, Cox proportional hazards models were fit for each CpG controlling for age, stratified sex, tumor grade, smoking status, and stratified BCG receiving status. Without adjustment for immune cell proportions, we identified 27,575 CpGs whose methylation was associated with a significant (P < 0.005) difference in hazard of RFS (Fig. 3A, Additional file 3). We then conducted similar analyses controlling for age, sex, tumor grade, smoking status, BCG status, and winsorized immune cell proportions. As expected and demonstrating the importance of adjusting for cell-type proportions, the fully adjusted models were attenuated and identified 2528 CpGs whose methylation was associated with a significant difference (P < 0.005) in the hazard of RFS. The 10 CpGs most strongly (with smallest P-value) associated with hazard of RFS corresponded to 10 genes: TMCO4 (cg04738197), LENG9 (cg12057190), CDC42EP5 (cg12057190), LNP1 (cg02540094), TOMM70A (cg02540094), RUNX2 (cg08012149), TBXAS1 (cg01584377), SSH1 (cg16237760), SFXN2 (cg08609163), and COG3 (cg06172950). 1572 CpGs were associated with the increased hazard of RFS, and 956 CpGs were associated with the decreased hazard of RFS (Fig. 3B). The complete list of RFS-associated CpGs is provided in Additional files 4 and 5.

Fig. 3
figure 3

Volcano plots of recurrence-free survival (RFS) associated CpGs from the epigenome-wide association study (EWAS) analyses. The Cox multivariable model that was fitted in EWAS was shown in each plot. CpGs are colored in red (A) 12,105 CpGs were associated with the increased hazard of NMIBC 10-year RFS, and 15,470 CpGs were associated with the decreased hazard of NMIBC 10-year RFS. (B) 1572 CpGs were associated with the increased hazard of NMIBC 10-year RFS, and 956 CpGs were associated with the decreased hazard of NMIBC 10-year RFS

To gain a better understanding of the regions where the hazard-associated CpGs were located, we tested for enrichment of CpG island region context among CpG loci associated with a significant change in the hazard of RFS. We found that 1572 CpGs associated with the increased hazard of RFS were significantly enriched in the open sea (OR = 1.14, 95% CI = 1.03–1.27) and were significantly depleted in CpG island S Shore regions (OR = 0.82, 95% CI = 0.67–0.99). The 956 CpGs associated with a decreased hazard of RFS were significantly enriched in CpG island N Shore regions (OR = 1.30, 95% CI = 1.06–1.57) and were significantly depleted in CpG island (OR = 0.60, 95% CI = 0.46–0.78) (Fig. 4A). We also tested for enrichment of other gene regulatory regions among CpG loci associated with a significant change in the hazard of RFS.

Fig. 4
figure 4

Genomic context enrichment analysis of CpG sites whose methylation state is significantly associated with recurrence-free survival. Enrichment analysis of (A) relation to CpG island and (B) genomic context of NMIBC recurrence-free survival associated CpGs. The 2528 CpGs from EWAS (P < 0.005) were tested for enrichment versus all modeled CpGs. The bar represents the 95% confidence interval. Mantel–Haenszel was used to test RFS-associated CpGs enrichment of CpG island-related gnome context. An odds ratio larger than 1 means enrichment, and an odds ratio smaller than 1 indicates depletion

Location of CpGs in regulatory regions also was investigated and the CpGs associated with the increased hazard of RFS were significantly enriched in enhancer regions (OR = 1.78, 95% CI = 1.43–2.18), DNase hypersensitive sites (DHS) (OR = 1.25, 95% CI = 1.13–1.40), 5’UTR regions (OR = 1.37, 95% CI = 1.17–1.60) and gene body (OR = 1.15, 95% CI = 1.04–1.27); however, these CpGs were significantly depleted in regions 200–1500 bps upstream of the transcription start site (TSS1500) (OR = 0.85, 95% CI = 0.73–0.99). In addition, while the CpGs associated with a decreased hazard were strongly enriched in enhancer regions (OR = 4.18, 95% CI = 3.44–5.05) and DHS (OR = 2.46, 95% CI = 2.13–2.86), they were depleted for TSS200 (OR = 0.70, 95% CI = 0.48–0.99), gene body (OR = 0.85, 95% CI = 0.74–0.97) and transcription factor binding sites (TFBS) (OR = 0.68, 95% CI = 0.55–0.83) (Fig. 4B). Using the CpGs associated with hazard of RFS from cell-type unadjusted models in tests for enrichment gave results that were largely consistent with those above based on CpGs from fully adjusted models (Additional file 7: Figure S3).

To further understand the biological function of the hazard-associated CpGs, Gene Set Enrichment Analysis (GSEA) Molecular Signature Database (MSigDB) was used to explore the potential gene sets which might associate with the tumor recurrence or death of NMIBC patients. The input was 2528 RFS associated CpGs from the Cox model EWAS adjusting for immune cell composition. In the top 10 hazard-associated gene sets in gene ontology (GO) terms, some gene sets were related to neurological system processing (Additional file 7: Figure S4A), and most related genes were associated with the increased hazard of RFS in NMIBC (Additional file 7: Figure S4B). Since mdNLR was associated with the hazard of recurrence or death of NMIBC, we were interested in immune-related gene sets. We checked the immunologic signature gene set, and the top 10 gene sets and their genes were related to immune cell regulation. However, only one gene set (BCELL VS MDC UP: up-regulated genes in B cells compared with myeloid dendritic cells after vaccination for influenza) consisting of 41 genes was significantly associated with RFS (FDR < 0.05) (Additional file 7: Figure S4C-D). Results from GO term analyses using 27,575 CpGs from the cell-type unadjusted models identified 5 pathways related to immunologic regulation among the top 10 pathways. Further, in the immunologic signature gene set specifically, we observed the top 10 pathways were associated with monocytes and lymphocytes (Additional file 7: Figure S5).

Locus overlap analysis (LOLA) was used to test the enrichment of CpGs in genomic regions. As our analysis was conducted on blood samples, results focus on the genomic regions within hematopoietic stem cells. When controlling for immune cell profiles, the 2528 CpGs associated with the hazard of RFS in NMIBC patients were most significantly enriched in Histone H3 acetylated at lysine 9 and 14 (H3K9K14ac) (Q-value = 1.9 × 10−20) (Additional file 7: Figure S6A). LOLA results for the 27,575 CpGs associated with the hazard of RFS from cell type unadjusted models were most significantly enriched in cistrome of the promyelocytic leukemia protein (PML) (Q-value < 0.05) (Additional file 7: Figure S6B).

We also ran a differentially methylated regions analysis using DMRcate. In this study, we found 11 CpGs in a specific genome region overlapping with the gene BLCAP. NMIBC patients without tumor recurrence or death within 10 years had a higher mean of methylation levels for this region compared with NMIBC patients with tumor recurrence or death within 10 years, and one CpG in this region was significant (FDR = 2.63 × 10−14) (Additional files 6, 7: Figure S7).


In this study, we tested whether immune profiles and epigenetic features are associated with NMIBC recurrence. Although previous work observed the association between NLR and overall survival in NMIBC patients, it was in a smaller study sample and used DNA methylation data from a dated (second generation) array platform with ~ 27,000 CpGs. In this study, our sample size was nearly three times larger and DNA methylation data was collected using the current genome-scale platform (fourth generation), measuring ~ 860,000 CpGs for which an optimized cell-type deconvolution library exists to determine highly accurate immune cell-type proportions. We extended tests of association with the patient outcomes beyond the methylation-derived neutrophil-to-lymphocyte ratio (mdNLR) to include leukocyte-specific cell-type proportions. Our findings suggest that elevated mdNLR increased the hazard of RFS in NMIBC patients. These findings are consistent with previous studies demonstrating that NLR was significantly higher in high-risk NMIBC patients [27], and increased NLR was positively associated with poor prognosis [28, 29]. While past studies have shown a significant association between NLR and outcomes in NMIBC patients using the conventional method of CBC tests [30,31,32], this study uses NLR derived from blood methylation profiles. Compared with flow cytometry, estimation using the differentially methylated regions (DMR) library has several advantages: it does not require long sample processing time, large volume of blood, or intact cells as DNA and 5-methylcytosine are both stable [33]. The advantage of our approach is the ability to use archived samples that many investigators may already have. Moreover, utilizing blood samples to monitor patient outcomes is less invasive compared with a cystoscopy, the routine screening method. For NMIBC patients with a high risk of tumor recurrence, BCG is the standard intravesical immunotherapy to induce immune system eliminating bladder cancer cells that might be left after surgery [8, 34, 35], and therefore, blood immune profile is a potential prognostic factor. Immune profiling with DNA methylation data is a promising avenue for assessing NMIBC prognosis.

As NLR is composed of lymphocyte and neutrophil proportions, we also considered the association between each methylation-derived immune cell type proportion and patient outcomes. Interestingly, increased CD4T and CD8T cell proportions were associated with decreased NMIBC recurrence-free survival and overall survival. Prior work has shown that NMIBC patients with a high CD4T cell count in BCG pretreatment microenvironment have a significantly prolonged recurrence-free survival compared to patients with a low CD4T cell count [35]. In addition to the immune cell types explored in this study, other cell types have been shown to affect NMIBC development, such as GATA3 + T cells, regulatory T cells, and tumor-associated macrophages [35]. Other work showed an increased CD8+ILT2+ T cell proportion was associated with a significantly increased hazard of NMIBC recurrence [36]. Further, peripheral neutrophil x platelet/lymphocyte was inversely correlated with high-risk NMIBC recurrence-free survival [37]. These results demonstrate the potential of immune cell profiles in evaluating the prognosis of NMIBC patients, and future work mapping DNA methylation profiles of additional immune cell types and states could add detail to investigations of immune profiles and bladder cancer outcomes. In addition, peripheral immune cell distribution is affected by potentially residual confounding factors such as, infection [38], inflammation [39], lifestyles, treatments [40, 41], obesity [42], chronic alcohol consumption [43], and type 2 diabetes [44]. Since our study includes hundreds of subjects and has a decade or more of follow-up time, in the future, we will re-explore the effect of these residual confounding factors on the association of circulating immune cell distribution and NMIBC outcomes. In the future, we will perform higher resolution methylation cell mixture deconvolution to resolve additional immune cell types and employ blood count on prospectively collected samples. Together, these data will allow us to understand if specific subsets of neutrophils or lymphocytes contribute to high NLR in patients with poor outcomes.

Bacillus Calmette–Guérin (BCG) is the most commonly used immunotherapy for high-risk NMIBC. After transurethral resection of bladder tumors, NMIBC patients may receive BCG intravesical therapy to induce an immune response in the bladder to attack cancer cells. Previous studies have shown that tumor immune environments may interact with BCG and interfere with the efficacy of this therapy. For instance, IL-12 secreted by BCG-induced monocytes was increased in NMIBC patients without tumor recurrence, a phenomenon that may involve the innate immune memory of circulating monocytes [45]. Nevertheless, patients receiving BCG in our study did not have peripheral immune cell profiles that differed from patients who were BCG naïve, and BCG was not significantly associated with NMIBC outcomes. As the number of patients with BCG treatment was relatively low (N = 89; 14.8%) and limited power to observe possible therapy induced changes in peripheral immune profiles, future work is needed to address the potential associations of BCG with peripheral immune profiles and patient outcomes. Most BCG-treated patients were already high-risk. In addition, blood samples were collected only after surgery or initial BCG treatment. To assess the changes in immune profiles over time, having multiple blood draws is necessary for future work.

This EWAS identified several CpGs associated with NMIBC recurrence-free survival when controlling for stratified sex, age, tumor grade, smoking status, stratified BCG receiving status, and immune cell profiles. With a P-value < 0.005, 2528 CpGs were found to be associated with the hazard of RFS.

Among the top CpGs whose methylation was associated with NMIBC recurrence or death in the fully adjusted models, some genes have been previously associated with bladder cancer. Slingshot homolog-1 (SSH1) had a positive association with tumor grade, tumor invasion, and tumor recurrence of bladder cancer patients [46]. Runt-related transcription factor 2 (RUNX2) is a key factor of osteoblast differentiation and has been reported to be associated with epithelial-mesenchymal transition in bladder tumors. Furthermore, RUNX2 could predict early recurrence in bladder cancer patients with high accuracy [47, 48]. Similar to past studies, SSH1 and RUNX2 were associated with the increased risk of tumor recurrence in our study. In addition, NMIBC patients with tumor recurrence within ten years had significantly higher methylation levels in the gene body of these two genes compared with patients without tumor recurrence.

When not adjusting for immune cell profiles, 27,575 CpGs were associated with the decreased hazard of RFS. The top 10 CpGs with the most significant P-value corresponded to 6 genes: BCL11A (cg24361098), TMCO4 (cg04738197), MICALCL (cg01518090), GRAP2 (cg21012238), TRAM2 (cg15085626), and KIRREL (cg10570484). While these genes have not been reported to be associated with bladder cancer outcomes, they have been reported to be involved in mechanisms promoting cancer development in tumors [49,50,51,52] or blood [53, 54]. Although we measured blood methylation in this study, we plan to measure tumor methylation and will explore the association of CpGs in these genes with NMIBC outcomes. What was more interesting is that the model, with or without controlling for immune cell profiles, led to different results. While immune cell profiles are usually not adjusted in the Cox model, we controlled for immune cell profiles since the immune system plays a key role in tumor development. This study presents a new perspective to demonstrate the difference between models with or without adjusting for immune profiles, indicating the need for further investigation on the involvement of immune profiles in outcome analyses.

Through DMR analysis, we found NMIBC patients with tumor recurrence or death within 10 years had a lower methylation level in the BLCAP region compared with patients without tumor recurrence or death. The bladder cancer-associated protein (BLCAP) gene encodes a protein that stimulates apoptosis. It has been reported that loss of protein expression is associated with bladder tumor progression, and the application of staining patterns for this protein could be a potential biomarker in bladder cancer [55]. In functional analysis, strong nuclear expression of BLCAP was associated with expression of p-STAT3 and overall poor disease outcome. Additionally, BLCAP was discovered to interact with STAT3 physically and may involve the STAT3-mediated progression of precancerous lesions to invasive bladder tumors [56]. These results were consistent with our finding that NMIBC patients with poor outcomes had lower methylation levels in the BLCAP region. Since the model we used for DMRs analysis was adjusted for immune cell profiles, we are curious whether immune cells may play roles in the interaction between BLCAP and STAT3 and will investigate this in the future. The CpG site with significantly lower methylation (cg10642330) is located in the 5’ UTR of BLCAP and the gene body of NNAT. The methylation beta value of this CpG site in NNAT had been reported significantly higher in prostate cancer tissue relative to adjacent normal tissues [57]. Though no CpG site in the BLCAP region was significantly associated with the hazard of tumor recurrence or death in the adjusted EWAS, five BLCAP CpGs (cg26083330, cg23757721, cg13790727, cg03061677, and cg04489586) were associated in the model not controlling for immune cell proportions (not including the DMR analysis site cg10642330). In the future, we will investigate the relation of BLCAP tumor methylation with survival outcomes.

Our results reveal several features of peripheral blood immune profiles that are associated with outcomes in NMIBC patients. We showed that higher CD4 or CD8 proportions were associated with decreased hazard of recurrence or death and further established that high NLR is associated with an increased hazard of RFS. The EWAS portion of our study also points to epigenetic reprogramming within the immune compartment being involved in tumor recurrence of NMIBC patients. Future study in a prospective setting will assess the clinical utility of incorporating methylation in predicting hazard of recurrence and shaping recommendations for disease surveillance. In addition to immune cells, future work examining cell-type proportions in tumor microenvironments of NMIBC patients is needed to understand the relationship between peripheral immune profiles with tumor-infiltrating immune profiles and patient outcomes. This work contributes to our understanding of associations of methylation-derived immune profiles and NMIBC patient outcomes and could further contribute to developments in epigenetic biomarkers of cancer.


Here we demonstrate the associations of non-muscle-invasive bladder cancer outcomes with immune profiles. In addition, we identify preliminary evidence of discrete and regional CpG methylation associations with bladder cancer outcomes. The findings could contribute to developing epigenetic biomarkers for recurrence-free survival in non-muscle-invasive bladder cancer.


Study subjects and samples

The subjects and data used in this study are described in more detail in prior publications [58,59,60]. Briefly, subjects were recruited from all three phases of a New Hampshire population-based bladder cancer case–control study [61]. The first wave of this study (phase 1) collected blood samples from 331 individuals diagnosed with incident bladder cancer between July 1994 and June 1998. The second study wave (phase 2) collected blood samples from 243 individuals diagnosed between July 1998 and December 2001. Finally, the third study wave (phase 3) obtained blood samples from 194 individuals recruited and diagnosed between July 2002 and December 2004. All the subjects were identified using the New Hampshire State Cancer Registry, hospital pathology departments, and hospital cancer registries, and all blood samples were collected after the time of diagnosis (time range: 20–1790 days). Among patients, 40 patients received BCG treatment in phase 1, 29 patients received BCG in phase 2, and 19 patients received BCG in phase 3. All patients with BCG treatment had blood drawn after receiving BCG (time range: 7–1542 days). An outline of data filtering and inclusion/exclusion criteria applied to these data are shown in Fig. 1. Briefly, subjects without muscle-invasive status, histopathology re-review, tumor grade, smoking status, or pack-years were removed from the study. Subjects that withstood the aforementioned exclusion criteria were retained and used in downstream statistical analyses.

DNA extraction, quantification, and bisulfite modification

Each blood sample was maintained at 4C and frozen within 24 h of blood draw. A hundred μl buffy coat was used to extract DNA. The QIAMP DNA blood & Tissue kit was used to extract DNA from blood samples according to the manufacturer’s protocol. Extracted DNA was quantified by using the Qubit 3.0 Fluorometer. After bisulfite modification, an established 5 mC microarray protocol optimized for Illumina methylation arrays was used to determine the genome-wide 5mC profile. DNA samples were subjected to bisulfite conversion (according to the manufacturer’s protocol of the Zymo EZ DNA methylation Kit) with an input of 750 ng per sample and whole-genome amplified prior to array hybridization. Recovered substrate ssDNA were submitted for DNA methylation array processing.

DNA methylation data

Bisulfite-modified DNA samples were measured for their DNA methylation status using the MethylationEPIC array, which interrogates > 860,000 CpG sites. Probe intensity data (iDAT files) from the EPIC methylation array were processed for quality control via the R package minfi [62] and ENmix [63] in R version 3.6. After quality control, 11 samples with low-quality CpG values or bisulfite intensity (threshold: 7000) were excluded from the study. The data were then normalized and conducted background correction through preprocessNoob procedure from minfi. The ComBat [64] was used to adjust for potential batch effect. Probes with a detection P > 1.0 × 10−6 in more than 10% of the samples were excluded (32,414). Also, 98,826 probes, which are cross-reactive, SNP-associated, and non-CpG (CpH) methylation [65], as well as 17,120 probes on sex chromosomes, were excluded. In total, 726,856 probes were used in downstream statistical analyses in downstream statistical analyses. IlluminaHumanMethylationEPICanno.ilm10b4.hg19 [66] was used to annotate CpG sites. Relation to CpG island was defined by the “Relation_to_Island” as in the Illumina annotation used in the genomic context analysis. ‘5’UTR,’ ‘Exon,’ ‘Gene Body,’ and ‘3’UTR’ contexts were defined by having ‘5UTR,’ ‘ExonBnd,’ ‘Body’ and ‘3UTR’ in UCSC_RefGene_Group. ‘Enhancer’ context was defined by having a record in the Phantom5_Enhancers. ‘DHS’ context was defined by finding a record in the DNase_Hypersensitivity_NAME. ‘TFBS’ context was defined by having a record in the TFBS_NAME.

Statistical analysis

The estimation of cell-type proportions was processed through the estimateCellCounts2 from the FlowSorted.Blood.EPIC package in Bioconductor (version 3.9; [18]). Methylation-derived neutrophil-to-lymphocyte ratio (mdNLR) was calculated by performing cell-mixture deconvolution to estimate the proportion of leukocyte subtypes and the ratio of neutrophil proportion to lymphocyte proportion was then computed. Individual leukocyte cell-type proportions and mdNLR were included in outcome analyses as continuous variables for Cox proportional hazard regression. In addition, mdNLR was dichotomized based on the median mdNLR for the Kaplan–Meier method.

Ten-year overall survival was defined as the time interval from the date of initial diagnosis to death within 10 years (all deaths were related with bladder cancer). Subjects who were alive or lost to follow-up were censored at the last follow-up. Ten-year recurrence-free survival was defined as the time interval from the date of initial diagnosis to the first tumor recurrence or death (all causes), whichever occurred first within 10 years. Patients alive and free of the disease or lost to follow-up were censored at the last follow-up. For both overall and recurrence-free survival, patient survival times over 10 years were truncated at 10 years, and patients were censored if the first tumor recurrence or death occurred after 10 years. The median survival times for the two survival outcomes were estimated using the Kaplan–Meier method. In multivariable analyses, Cox proportional hazard regression models were used to examine the association of each variable on bladder cancer outcomes and were fit via coxph in the survival R package. The proportional hazards assumption was tested by using cox.zph from the survival R package. The cox.zph function tests the proportionality of all the predictors in Cox models by creating interactions with time. As sex and BCG treatment status violated the proportional hazards assumption, stratification on both variables was included. The linearity assumption was examined via ggcoxfunctional from the R survminer package, and methylation-derived immune cell profiles were found to violate the linearity assumption. Hence, winsorization was used on methylation-derived immune cell profiles. The winsorization cut-point of each immune cell profile is shown in Additional file 7: Figure S1. A P-value of < 0.05 was the significance threshold on multivariable analysis. Cox model results were presented using the stargazer R package.

Epigenome-wide association study (EWAS), enrichment analysis, and differentially methylated regions (DMRs)

An epigenome-wide association study (EWAS) was performed using ewaff R package ( to investigate the association of CpG-specific DNA methylation and bladder cancer recurrence. We fit Cox proportional hazards models independently to each CpG, controlling for age, sex, tumor grade, smoking status, Bacillus Calmette–Guérin (BCG) receiving status, and immune cell profiles. Since all P-values adjusted for false discovery rate from EWAS results are higher than 0.05, we relaxed the threshold for genomic context and enrichment analyses using a P-value of < 0.005.

For CpG sites associated (P < 0.005) with bladder tumor recurrence, we examined whether those CpGs were enriched in CpG island-related genomic context or regulatory regions via using Mantel–Haenszel tests, adjusted for Illumina probe type to eliminate the difference in distributions on genomic context. CpG island-related genomic context includes open sea, north shelves, north shores, islands, south shores, and south shelves. Regulatory regions include enhancers, DNase hypersensitivity sites, 5’UTR, TSS1500, TSS200, 1st Exon, Exon, gene body, 3’UTR, and transcription factor binding sites. Next, the Locus Overlap analysis (LOLA) R package in Bioconductor (version 1.20.0; [67]) was used to investigate enrichment of genomic regions limited to tissue equal to “hematopoietic stem cell.” Finally, gometh and gsameth from the missMethyl package in Bioconductor (version 1.6.2; [68]) were used to test the enrichment of gene sets for Gene Ontology (GO) terms and the C7: immunologic signature gene set in the Gene Set Enrichment Analysis (GSEA) Molecular Signature Database (MSigDB) cnetplot from the enrichplot package in Bioconductor (version 1.10.1; was used to plot gene-concept network plots.

Differentially methylated regions were identified and extracted through the dmrcate and extractRanges from the DMRcate R package [69]. The inputs were logit-transform of beta values (M-values). The phenotype of interest for comparison was “NMIBC patients with tumor recurrence within ten years or not” in our designed model. In addition, the designed model was adjusted for sex, age, tumor grade, smoking status, BCG receiving status, and immune cell profiles. We relaxed the threshold for DMRs analysis using an FDR-corrected P-value of < 0.1. Then, the visualizeGene from the sesame package [70] was used for observing the methylation levels of CpGs in regions identified by DMRcate.

Availability of data and materials

The datasets generated and analyzed during this current study are available in the gene expression omnibus repository at GSE183920.



Non-muscle-invasive bladder cancer


Neutrophil-to-lymphocyte ratio


Methylation-derived neutrophil-to-lymphocyte ratio


Epigenome-wide association study


Bacillus Calmette–Guérin


Complete blood count differential


Differentially methylated regions


Recurrence-free survival


Overall survival


Hazard ratio


Odds ratio


Confidence interval




DNase hypersensitive sites


Transcription factor binding sites


False discovery rate


Locus overlap analysis


Gene ontology


Gene set enrichment analysis


Molecular signature database


  1. Atlanta: American Cancer Society. Cancer Facts & Figures 2020. Am. Cancer Soc. 2020.

  2. Martinez Rodriguez RH, Buisan Rueda O, Ibarz L. Bladder cancer: present and future. Med Clín. 2017;149:449–55.

    Article  Google Scholar 

  3. Svatek RS, Lotan Y. Is there a rationale for bladder cancer screening? Curr Urol Rep. 2008;9:339–41.

    PubMed  PubMed Central  Google Scholar 

  4. Larré S, Catto JWF, Cookson MS, Messing EM, Shariat SF, Soloway MS, et al. Screening for bladder cancer: rationale, limitations, whom to target, and perspectives. Eur Urol. 2013;63:1049–58.

    PubMed  PubMed Central  Google Scholar 

  5. Degeorge KC, Holt HR, Hodges SC. Bladder cancer: diagnosis and treatment. 2017.

  6. Hall MC, Chang SS, Dalbagni G, Pruthi RS, Seigne JD, Skinner EC, et al. Guideline for the management of nonmuscle invasive bladder cancer (stages Ta, T1, and Tis): 2007 update. J Urol. 2007;178:2314–30.

    PubMed  PubMed Central  Google Scholar 

  7. Mossanen M, Gore JL. The burden of bladder cancer care: direct and indirect costs. Curr Opin Urol. 2014;24:487–91.

    PubMed  PubMed Central  Google Scholar 

  8. Lim CJ, Nguyen PHD, Wasser M, Kumar P, Lee YH, Nasir NJM, et al. Immunological hallmarks for clinical response to BCG in bladder cancer. Front Immunol. 2021;11:1–13.

    Google Scholar 

  9. Ogihara K, Kikuchi E, Yuge K, Yanai Y, Matsumoto K, Miyajima A, et al. The preoperative neutrophil-to-lymphocyte ratio is a novel biomarker for predicting worse clinical outcomes in non-muscle invasive bladder cancer patients with a previous history of smoking. Ann Surg Oncol. 2016;23:1039–47.

    PubMed  PubMed Central  Google Scholar 

  10. Tan YG, Eu E, Lau Kam On W, Huang HH. Pretreatment neutrophil-to-lymphocyte ratio predicts worse survival outcomes and advanced tumor staging in patients undergoing radical cystectomy for bladder cancer. Asian J Urol. 2017;4:239–46.

    Article  PubMed  PubMed Central  Google Scholar 

  11. Kang M, Jeong CW, Kwak C, Kim HH, Ku JH. Preoperative neutrophil-lymphocyte ratio can significantly predict mortality outcomes in patients with non-muscle invasive bladder cancer undergoing transurethral resection of bladder tumor. Oncotarget. 2017;8:12891–901.

    PubMed  PubMed Central  Google Scholar 

  12. Yoshida T, Kinoshita H, Yoshida K, Mishima T, Yanishi M, Inui H, et al. Prognostic impact of perioperative lymphocyte–monocyte ratio in patients with bladder cancer undergoing radical cystectomy. Tumor Biol. 2016;37:10067–74.

    CAS  Article  Google Scholar 

  13. Ma JY, Hu G, Liu Q. Prognostic significance of the lymphocyte-to-monocyte ratio in bladder cancer undergoing radical cystectomy: a meta-analysis of 5638 individuals. Dis Mark. 2019;2019.

  14. Bhindi B, Hermanns T, Wei Y, Yu J, Richard PO, Wettstein MS, et al. Identification of the best complete blood count-based predictors for bladder cancer outcomes in patients undergoing radical cystectomy. Br J Cancer. 2016;114:207–12.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  15. Ojerholm E, Smith A, Hwang W, Baumann BC, Tucker KN, Lerner SP, et al. Neutrophil-to-lymphocyte ratio as a bladder cancer biomarker: assessing prognostic and predictive value in SWOG 8710. Cancer. 2018;123:794–801.

    Google Scholar 

  16. Dixon LR. The complete blood count: physiologic basis and clinical usage. J Perinat Neonatal Nurs. 1997;11:1–18.

    CAS  PubMed  PubMed Central  Google Scholar 

  17. Accomando WP, Wiencke JK, Houseman EA, Nelson HH, Kelsey KT. Quantitative reconstruction of leukocyte subsets using DNA methylation. Genome Biol. 2014;15:1–12.

    Google Scholar 

  18. Salas LA, Koestler DC, Butler RA, Hansen HM, Wiencke JK, Kelsey KT, et al. An optimized library for reference-based deconvolution of whole-blood biospecimens assayed using the Illumina HumanMethylationEPIC BeadArray. Genome Biol. 2018;19:64.

    PubMed  PubMed Central  Google Scholar 

  19. Suelves M, Carrió E, Núñez-Álvarez Y, Peinado MA. DNA methylation dynamics in cellular commitment and differentiation. Brief Funct Genom. 2016;15:443–53.

    CAS  Google Scholar 

  20. Salas LA, Wiencke JK, Koestler DC, Zhang Z, Christensen BC, Kelsey KT. Tracing human stem cell lineage during development using DNA methylation. Genome Res. 2018;28:1285–95.

    CAS  PubMed  PubMed Central  Google Scholar 

  21. Baron U, Türbachova I, Hellwag A, Eckhardt F, Berlin K, Hoffmuller U, et al. DNA methylation analysis as a tool for cell typing. Epigenetics. 2006;1:55–60.

    PubMed  PubMed Central  Google Scholar 

  22. Houseman EA, Accomando WP, Koestler DC, Christensen BC, Marsit CJ, Nelson HH, et al. DNA methylation arrays as surrogate measures of cell mixture distribution. BMC Bioinform. 2012;13.

  23. Li Y, Pan X, Roberts ML, Liu P, Kotchen TA, Cowley AW, et al. Stability of global methylation profiles of whole blood and extracted DNA under different storage durations and conditions. Epigenomics. 2018;10:797–811.

    CAS  PubMed  PubMed Central  Google Scholar 

  24. Titus AJ, Gallimore RM, Salas LA, Christensen BC. Cell-type deconvolution from DNA methylation: a review of recent applications. Hum Mol Genet. 2017;26:R216–24.

    CAS  PubMed  PubMed Central  Google Scholar 

  25. Wiencke JK, Koestler DC, Salas LA, Wiemels JL, Roy RP, Hansen HM, et al. Immunomethylomic approach to explore the blood neutrophil lymphocyte ratio (NLR) in glioma survival. Clin Epigenet. 2017;9:1–11.

    Google Scholar 

  26. Koestler DC, Usset J, Christensen BC, Marsit CJ, Karagas MR, Kelsey KT, et al. DNA methylation-derived neutrophil-tolymphocyte ratio: an epigenetic tool to explore cancer inflammation and outcomes. Cancer Epidemiol Biomarkers Prev. 2017;26:328–38.

    CAS  PubMed  PubMed Central  Google Scholar 

  27. Aydın M, Bitkin A, Kadıhasanoğlu M, İrkılata L, Akgüneş E, Keleş M, et al. Correlation of neutrophil-lymphocyte ratio and risk scores in non-muscle invasive bladder cancer. Actas Urol Españolas. 2019;43:503–8.

    Google Scholar 

  28. Yuk HD, Jeong CW, Kwak C, Kim HH, Ku JH. Elevated neutrophil to lymphocyte ratio predicts poor prognosis in non-muscle invasive bladder cancer patients: initial intravesical bacillus Calmette-Guerin treatment after transurethral resection of bladder tumor setting. Front Oncol. 2019;9:1–8.

    Google Scholar 

  29. Getzler I, Bahouth Z, Nativ O, Rubinstein J, Halachmi S. Preoperative neutrophil to lymphocyte ratio improves recurrence prediction of non-muscle invasive bladder cancer. BMC Urol BMC Urology. 2018;18:1–10.

    Google Scholar 

  30. Vartolomei MD, Porav-Hodade D, Ferro M, Mathieu R, Abufaraj M, Foerster B, et al. Prognostic role of pretreatment neutrophil-to-lymphocyte ratio (NLR) in patients with non–muscle-invasive bladder cancer (NMIBC): a systematic review and meta-analysis. Urol Oncol Semin Orig Investig. 2018;36:389–99.

    Article  Google Scholar 

  31. Zhang Q, Lai Q, Wang S, Meng Q, Mo Z. Clinical value of postoperative neutrophil-tolymphocyte ratio change as a detection marker of bladder cancer recurrence. Cancer Manag Res. 2021;13:849–60.

    PubMed  PubMed Central  Google Scholar 

  32. Cantiello F, Russo GI, Vartolomei MD, Farhan ARA, Terracciano D, Musi G, et al. Systemic inflammatory markers and oncologic outcomes in patients with high-risk non–muscle-invasive urothelial bladder cancer. Eur Urol Oncol. 2018;1:403–10.

    PubMed  PubMed Central  Google Scholar 

  33. Teschendorff AE, Zheng SC. Cell-type deconvolution in epigenome-wide association studies: a review and recommendations. Epigenomics. 2017;9:757–68.

    CAS  PubMed  PubMed Central  Google Scholar 

  34. Witjes JA, Hendricksen K. Intravesical pharmacotherapy for non-muscle-invasive bladder cancer: a critical analysis of currently available drugs, treatment schedules, and long-term results. Eur Urol. 2008;53:45–52.

    PubMed  PubMed Central  Google Scholar 

  35. Pichler R, Fritz J, Zavadil C, Schäfer G, Culig Z, Brunner A. Tumor-infiltrating immune cell subpopulations influence the oncologic outcome after intravesical bacillus calmette-guérin therapy in bladder cancer. Oncotarget. 2016;7:39916–30.

    PubMed  PubMed Central  Google Scholar 

  36. Desgrandchamps F, LeMaoult J, Goujon A, Riviere A, Rivero-Juarez A, Djouadou M, et al. Prediction of non-muscle-invasive bladder cancer recurrence by measurement of checkpoint HLAG’s receptor ILT2 on peripheral CD8+ T cells. Oncotarget. 2018;9:33160–9.

    PubMed  PubMed Central  Google Scholar 

  37. Akan S, Ediz C, Sahin A, Tavukcu HH, Urkmez A, Horasan A, et al. Can the systemic immune inflammation index be a predictor of BCG response in patients with high-risk non-muscle invasive bladder cancer? Int J Clin Pract. 2021;75:1–9.

    Google Scholar 

  38. Douglas RG, Alford RH, Cate TR, Couch RB. The leukocyte response during viral respiratory illness in man. Ann Intern Med. 1966;64:521–30.

    PubMed  PubMed Central  Google Scholar 

  39. Woodell-May JE, Sommerfeld SD. Role of inflammation and the immune system in the progression of osteoarthritis. J Orthop Res. 2020;38:253–7.

    Article  PubMed  Google Scholar 

  40. Bedel C, Korkut M, Armag ̆an HH. NLR, d-NLR and PLR can be affected by many factors. Int Immunopharmacol. 2021;90:107154.

  41. Lin BD, Hottenga JJ, Abdellaoui A, Dolan CV, De Geus EJC, Kluft C, et al. Causes of variation in the neutrophil-lymphocyte and platelet-lymphocyte ratios: a twin-family study. Biomark Med. 2016;10:1061–72.

    CAS  PubMed  PubMed Central  Google Scholar 

  42. Elisia I, Lam V, Cho B, Hay M, Li MY, Kapeluto J, et al. Exploratory examination of inflammation state, immune response and blood cell composition in a human obese cohort to identify potential markers predicting cancer risk. PLoS ONE. 2020;15:1–21.

    Google Scholar 

  43. Laso FJ, Vaquero JM, Almeida J, Marcos M, Orfao A. Chronic alcohol consumption is associated with changes in the distribution, immunophenotype, and the inflammatory cytokine secretion profile of circulating dendritic cells. Alcohol Clin Exp Res. 2007;31:846–54.

    CAS  PubMed  PubMed Central  Google Scholar 

  44. Grossmannm V, Schmitt VH, Zeller T, Panova-Noeva M, Schulz A, Laubert-Reh D, et al. Profile of the immune and inflammatory response in individuals with prediabetes and type 2 diabetes. Diabetes Care. 2015;38:1356–64.

    Google Scholar 

  45. Graham CH, Paré JF, Cotechini T, Hopman W, Hindmarch CCT, Ghaffari A, et al. Innate immune memory is associated with increased disease-free survival in bladder cancer patients treated with bacillus Calmette-Guérin. Can Urol Assoc J. 2021;15:1–13.

    Google Scholar 

  46. Luo Q, Liu Y, Zhao H, Guo P, Wang Q, Li W, et al. Slingshot homolog-1 expression is a poor prognostic factor of pT1 bladder urothelial carcinoma after transurethral resection. World J Urol. 2020;38:2849–56.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  47. Abdelzaher E, Kotb AF. High coexpression of runt-related transcription factor. Appl Immunohistochem Mol Morphol. 2016;2:24.

    Google Scholar 

  48. Liu B, Pan S, Liu J, Kong C. Cancer-associated fibroblasts and the related Runt-related transcription factor 2 (RUNX2) promote bladder cancer progression. Gene. 2021;775:145451.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  49. Li L, Ugalde AP, Scheele CLGJ, Dieter SM, Nagel R, Ma J, et al. A comprehensive enhancer screen identifies TRAM2 as a key and novel mediator of YAP oncogenesis. Genome Biol. 2021;22:1–28.

    Google Scholar 

  50. Deng W, Wang Y, Gu L, Duan B, Cui J, Zhang Y, et al. MICAL1 controls cell invasive phenotype via regulating oxidative stress in breast cancer cells. BMC Cancer. 2016;16:1–11.

    CAS  Article  Google Scholar 

  51. Chen K, Zhao R, Yao G, Liu Z, Shi R, Geng J. Overexpression of kin of IRRE-Like protein 1 (KIRREL) as a prognostic biomarker for breast cancer. Pathol Res Pract. 2020;216:153000.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  52. Lee I, Yeom SY, Lee SJ, Kang WK, Park C. A novel senescence-evasion mechanism involving Grap2 and cyclin D interacting protein inactivation by Ras associated with diabetes in cancer cells under doxorubicin treatment. Cancer Res. 2010;70:4357–65.

    CAS  PubMed  PubMed Central  Google Scholar 

  53. Weniger MA, Pulford K, Gesk S, Ehrlich S, Banham AH, Lyne L, et al. Gains of the proto-oncogene BCL11A and nuclear accumulation of BCL11AXL protein are frequent in primary mediastinal B-cell lymphoma [1]. Leukemia. 2006;20:1880–2.

    CAS  PubMed  PubMed Central  Google Scholar 

  54. Satterwhite E, Sonoki T, Willis TG, Harder L, Nowak R, Arriola EL, et al. The BCL11 gene family: involvement of BCL11A in lymphoid malignancies. Blood. 2001;98:3413–20.

    CAS  PubMed  PubMed Central  Google Scholar 

  55. Moreira JMA, Ohlsson G, Gromov P, Simon R, Sauter G, Celis JE, et al. Bladder cancer-associated protein, a potential prognostic biomarker in human bladder cancer. Mol Cell Proteom. 2010;9:161–77.

    CAS  Google Scholar 

  56. Gromova I, Svensson S, Gromov P, Moreira JMA. Identification of BLCAP as a novel STAT3 interaction partner in bladder cancer. PLoS ONE. 2017;12:1–18.

    Google Scholar 

  57. Jacobs DI, Mao Y, Fu A, Kelly WK, Zhu Y. Dysregulated methylation at imprinted genes in prostate tumor tissue detected by methylation microarray. BMC Urol. 2013;13:1.

    Google Scholar 

  58. Baris D, Karagas MR, Verrill C, Johnson A, Andrew AS, Marsit CJ, et al. A case-control study of smoking and bladder cancer risk: emergent patterns over time. J Natl Cancer Inst. 2009;101:1553–61.

    PubMed  PubMed Central  Google Scholar 

  59. Schned AR, Andrew AS, Marsit CJ, Zens MS, Kelsey KT, Karagas MR. Survival following the diagnosis of noninvasive bladder cancer: WHO/International Society of Urological pathology versus WHO classification systems. J Urol. 2007;178:1196–200.

    PubMed  PubMed Central  Google Scholar 

  60. Kelsey KT, Hirao T, Schned A, Hirao S, Devi-Ashok T, Nelson HH, et al. A population-based study of immunohistochemical detection of p53 alteration in bladder cancer. Br J Cancer. 2004;90:1572–6.

    CAS  PubMed  PubMed Central  Google Scholar 

  61. Karagas MR, Tosteson TD, Blum J, Morris JS, Baron JA, Klaue B. Design of an epidemiologic study of drinking water arsenic exposure and skin and bladder cancer risk in a U.S. population. Environ Health Perspect. 1998;106:1047–50.

  62. Aryee MJ, Jaffe AE, Corrada-Bravo H, Ladd-Acosta C, Feinberg AP, Hansen KD, et al. Minfi: a flexible and comprehensive bioconductor package for the analysis of Infinium DNA methylation microarrays. Bioinformatics. 2014;30:1363–9.

    CAS  PubMed  PubMed Central  Google Scholar 

  63. Xu Z, Niu L, Li L, Taylor JA. ENmix: a novel background correction method for Illumina HumanMethylation450 BeadChip. 2016;44:1–6.

  64. Johnson WE, Li C, Rabinovic A. Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics. 2007;8:118–27.

    Google Scholar 

  65. Zhou W, Laird PW, Shen H. Comprehensive characterization, annotation and innovative use of Infinium DNA methylation BeadChip probes. Nucl Acids Res. 2017;45:e22.

    PubMed  PubMed Central  Google Scholar 

  66. KD H. IlluminaHumanMethylationEPICanno.ilm10b4.hg19: annotation for Illumina’s EPIC methylation arrays. R Packag version 060. 2017;

  67. Sheffield NC, Bock C. LOLA: enrichment analysis for genomic region sets and regulatory elements in R and bioconductor. Bioinformatics. 2016;32:587–9.

    CAS  PubMed  PubMed Central  Google Scholar 

  68. Phipson B, Maksimovic J, Oshlack A. MissMethyl: an R package for analyzing data from Illumina’s HumanMethylation450 platform. Bioinformatics. 2016;32:286–8.

    CAS  PubMed  PubMed Central  Google Scholar 

  69. Peters TJ, Buckley MJ, Statham AL, Pidsley R, Samaras K, Lord VR, et al. De novo identification of differentially methylated regions in the human genome. Epigenet Chromatin. 2015;8:1–16.

    Google Scholar 

  70. Zhou W, Triche TJ, Laird PW, Shen H. SeSAMe: reducing artifactual detection of DNA methylation by Infinium BeadChips in genomic deletions. Nucl Acids Res. 2018;46:1–15.

    CAS  Google Scholar 

Download references


Not applicable


Work at Dartmouth College was supported by the National Institutes of Health Grant Numbers R01CA216265, P20GM104416, R01CA057494, R01CA253976, and CDMRP/Department of Defense W81XWH-20–1-0778. Work at the University of California, San Francisco was supported by R01CA207360, P50CA097257, JKW is supported by the Robert Magnin Newman Endowed Chair in Neuro-oncology.

Author information




JC, BCC, and LAS contributed to the analysis. JC, BCC, LAS, DCK, KTK, JKW, AMM, and MRK contributed to the conception, design, and writing of the manuscript. ASA and JDS contributed to the data collection. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Brock C. Christensen.

Ethics declarations

Ethics approval and consent to participate

This case study was approved by the Dartmouth IRB.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests. JKW and KTK are co-founders of Cellintec which had no role in this work.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1

: Cox proportional hazards 10-year recurrence-free survival. Created via using stargazer() function

Additional file 2

: Cox proportional hazards 10-year overall survival. Created via using stargazer() function

Additional file 3

: The 27,575 CpGs associated with the hazard of 10-year RFS when the model without adjusting for immune cell profiles. The list of the 27,575 CpGs which were associated with the hazard of 10-year RFS from the EWAS with Cox model without adjusting for immune cell profiles.

Additional file 4

: The 1,572 CpGs associated with the increased hazard of 10-year RFS. The list of the 1,572 CpGs which were associated with the increased hazard of 10-year RFS from the EWAS with Cox model adjusting for immune cell profiles.

Additional file 5

: The 956 CpGs associated with the decreased hazard of 10-year RFS. The list of the 956 CpGs which were associated with the decreased hazard of 10-year RFS from the EWAS with Cox model adjusting for immune cell profiles.

Additional file 6

: The results of differentially methylated regions. no.cpgs: Number of CpG sites constituting the significant region.

Additional file 7

: Supplemental tables and figures.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Chen, JQ., Salas, L.A., Wiencke, J.K. et al. Immune profiles and DNA methylation alterations related with non-muscle-invasive bladder cancer outcomes. Clin Epigenet 14, 14 (2022).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:


  • DNA methylation
  • Non-muscle-invasive bladder cancer
  • Immune profile
  • Immunomethylomic
  • Recurrence
  • Survival