Genome-Scale Methylation Analysis Identifies Immune Profiles and Age Acceleration Associations with Bladder Cancer Outcomes

Bladder cancer is a malignant urogenital neoplasm and is classified into non–muscle-invasive bladder cancer (NMIBC) and muscle-invasive bladder cancer. In 2022, an estimated 81,000 new cases of bladder cancer and 17,000 deaths from the disease occurred in the United States (1). The common risk factors of bladder cancer are age, sex, and smoking. Bladder cancer is four times more common in men compared with women (2). About 90% of patients with bladder cancer are age 55 or older, and patients younger than 60 have a higher 10-year overall survival (OS) rate than patients older than 60 (3). Around 50% to 60% of new cases are attributed to smoking, and current smoking has a positive association with the risk of recurrence (4). The conventional treatment for bladder cancer is surgery or surgery in combination with chemotherapy drugs or intravesical immunotherapy [Bacillus Calmette-Guérin (BCG)] (5). Even though transurethral resection and immunotherapy generally control the disease (1, 6), the tumor recurrence rate is about 40% after treatment (7, 8). Predictive biomarkers that alert clinicians to recurrence would help to improve the clinical management of bladder cancer.

Circulating immune cell profiles have been associated with outcomes in patients with bladder cancer. For example, CD8⁺ cell proportions were associated with a decreased risk of tumor recurrence (9). Also, an elevated neutrophil-to-lymphocyte ratio (NLR) has been associated with worse OS and higher recurrence rate (10, 11). Previously (12), we measured peripheral blood DNA methylation profiles of patients with NMIBC and applied cell-type deconvolution to estimate the proportions of six immune cell types (13). CD4T and CD8T cell proportions were associated with decreased risk of death and recurrence. Yet, there have been limited studies investigating the relationship of circulating immune profiles in bladder cancer with disease outcomes. Furthermore, subtypes of each major cell type have been shown to affect cancer development distinctly. For instance, cytotoxic CD4⁺ T cells can kill the bladder tumor cells, and in contrast, regulatory CD4⁺ T can suppress the activity of cytotoxic CD4⁺ T cells and lead to tumor growth indirectly (14). To broaden the scope of the effects of circulating immune profiles on bladder cancer outcomes, it is necessary to investigate the association between leukocyte subtypes and bladder cancer outcomes.

Recently, our group developed an enhanced method to perform high-resolution cell mixture deconvolution to resolve 12 immune cell types in blood using DNA methylation measures [naïve and memory B, CD4T, and CD8T, as well as regulatory T, monocyte, natural killer (NK) cells, neutrophils, basophils, and eosinophils; ref. 15]. Because DNA methylation involves gene regulation for cell lineage specification (16), cell-specific differentially methylated regions can be utilized to distinguish cell types with reference-based deconvolution (13, 17, 18). Also, DNA methylation cytometry is more efficient for immune profiling with high accuracy than flow cytometry and can be applied to archival specimens.

Epigenetic clocks have been developed to estimate chronologic age or physiologic age in regard to aging outcomes, such as cancers and all-cause mortality (19–21). Age acceleration derived from these clocks has been associated with prospective risk in lung, kidney, and pancreatic cancer (22–24). Moreover, age acceleration has been associated with outcomes in other cancers (25, 26). Though only a few studies have reported the association of age acceleration with bladder cancer risk, they have not exhibited consistent results nor mentioned the subtype of bladder cancer they investigated. For example, one study showed that Pheno and Grim age acceleration were positively associated with the prospective risk in bladder cancers (22). Another study reported that Horvath and Hannum age acceleration was not associated with bladder cancer risk (27). Because aging is one of risks of bladder cancer (28), we investigated the association of age acceleration with bladder cancer outcomes.

Here, we hypothesized that DNA methylation-derived immune cell proportions and age acceleration are associated with bladder cancer outcomes. We applied our new methylation cytometry approach for extended immune cell resolution to DNA methylation profiles of archival blood samples from a population-based study containing NMIBC (N = 601) patients. We then tested the association of cancer outcomes with each leukocyte subtype proportion and age acceleration. We also used partDSA (29), a classification and regression trees method, and a semi-supervised recursively partitioned mixture model (SS-RPMM; ref. 30), to group/cluster our subjects based on cancer outcomes, patients’ demographics, tumor characteristics, and methylation profiles of specific CpG loci.

A detailed description of subjects who participated in the current study is available in prior studies (31–33). Briefly, bladder cancer subjects were recruited from three phases of a New Hampshire population-based case–control study (34). The first phase collected blood samples from 331 individuals diagnosed between July 1994 and June 1998 (phase I). The second phase collected blood samples from 243 individuals diagnosed between July 1998 and December 2001 (phase II). The third study phase collected blood samples from 194 individuals diagnosed between July 2002 and December 2004 (phase III). Patients with bladder cancer were identified using the New Hampshire State Cancer Registry and hospital cancer registry (patients in phase III were identified using the hospital cancer registry only). Patients’ OS data were from the National Death Index, and tumor recurrence data were ascertained through chart review. We performed four comprehensive National Death Index (NDI) searches for the years 2008, 2010, 2014, and 2018 to identify cases of death. In addition, the New Hampshire State Cancer Registry had previously reported some deaths to us through their conducted searches. During each NDI search, we included all bladder cancer case–control study participants who had not been previously matched with an NDI death. To ensure accurate matching, we followed the NDI-recommended method and utilized the code and algorithms provided by the NDI to score the matches. Furthermore, we applied the NDI score interpretation using the recommended NPCR algorithm (35) to enhance the accuracy of our findings. Subjects without histopathology re-review, muscle-invasive status, tumor grade, or smoking status were excluded from the study. The remaining 601 patients with NMIBC were used in downstream statistical analyses. In addition, 40 patients received BCG in phase I, 30 received BCG in phase II, and 19 received BCG in phase III. All patients with BCG treatment had their blood drawn after BCG treatment, and all blood samples were taken after the initial diagnosis. This study was approved by the Dartmouth Human Research Protection Program (Institutional Review Board; approval number STUDY00010107). The written informed consent was obtained from the patients and the studies were conducted in accordance with Belmont Report.

After the blood draw, blood samples were kept at 4°C and frozen within 24 hours. DNA was extracted from blood samples using the QIAamp DNA Blood Kit (Qiagen) according to the manufacturer's protocol. Extracted DNA quantity and quality were assessed with Qubit 3.0 Fluorometer (Life Technologies) and Fragment Analyzer (Advanced Analytical). Then, extracted DNA underwent bisulfite conversion using EZ DNA Methylation Kit (Zymo Research) according to the manufacturer's protocol. Approximately 750 ng of bisulfite-modified DNA was used as input for the DNA methylation array.

After DNA extraction, quantification, and bisulfite modification, Infinium MethylationEPIC Bead Chips (Illumina, Inc.) were used to measure the DNA methylation status of bisulfite-modified DNA samples. Raw probe intensity data, iDAT files, from the methylation array were processed through preprocessNoob using minfi (RRID:SCR_012830; ref. 36), and quality control was performed using ENmix (37) R package. To distinguish from background noise, samples with more than 5% of probes with a detection P > 1.0 × 10⁻⁶ were not included. In addition, we dropped 32,713 probes that were not detected in more than 10% of the samples. Then, the bias of type-2 probe values was corrected using BMIQ (RRID:SCR_003446; ref. 38) from watermelon (RRID:SCR_001296; ref. 39) R package, and the ComBat (RRID:SCR_010974; ref. 40) was used to adjust for batch effects. Next, 106,522 probes previously reported to be cross-reactive, SNP-associated, non-CpG (CpH) methylation, and sex-specific were excluded (41). After these exclusions, 726,856 CpGs remained for downstream statistical analysis. The annotation for CpG sites was from IlluminaHumanMethylationEPICanno.ilm10b4.hg19.

Methylation age was estimated with the function methyAge from the ENmix (37) R package implementation of the methylation age estimation. Age acceleration was defined as the residual from a regression of methylation age on chronologic age. Cell-type proportions were estimated with the projectCellType_CP from the FlowSorted.Blood.EPIC (RRID:SCR_022540; ref. 13) R package. The NLR was calculated according to the ratio of neutrophil proportion to lymphocyte proportion.

Ten-year OS was defined as the time interval from the date of initial diagnosis to death. Patients alive or lost to follow-up were censored at the last follow-up. Similarly, 10-year recurrence-free survival (RFS) was defined as the time interval from the date of initial diagnosis to the first tumor recurrence or death, whichever occurred first, and patients alive and without tumor recurrence or lost to follow-up were censored at the last follow-up. For OS and RFS, survival times were truncated at 10 years. In univariate and multivariable analyses, coxph from the survival (RRID:SCR_021137) R package was used to fit Cox proportional hazards models to evaluate the association between bladder cancer outcomes and each variable. Only immune cell profiles significantly associated with bladder cancer outcomes in the univariate Cox model were subjected to multivariable analyses. Cox.zph from the survival R package was employed to test the proportional hazards assumption. Predictors with assumption violations were included as strata in the Cox models. The linearity assumption was examined with the ggcoxfunctional from the survminer (RRID:SCR_021094) R package. We conducted 2% winsorization on immune cell profiles identified to violate the linearity assumption. FDR-corrected P value of < 0.05 was the significance threshold on multivariable analysis.

To explore interactions between clinical variables and immune cell proportions in survival analysis, we applied a partitioning deletion/substitution/addition algorithm [partDSA (29, 42); R package] for model building employing the inverse probability censoring weighted-L2 loss function. Those variables (age, Hannum or Pheno age acceleration, sex, tumor grade, smoking status, BCG treatment status, and immune cell type proportions) associated with 10-year OS were included in the multivariable model as input. The partDSA approach resulted in three groups of subjects that were based on neutrophil and CD8 naïve cell proportions. After the model was built, corresponding Kaplan–Meier curves were generated, and HRs and 95% confidence intervals (CI) were calculated using the Cox model.

To identify a novel set of blood DNA methylation profiles associated with cancer outcomes, we applied a SS-RPMM algorithm (30). This method uses the recursive partitioning mixture model (RPMM), demonstrating an effective and efficient unsupervised clustering procedure for methylation data (43–46). To avoid overfitting and provide for validation of the model, we randomly split the total population into a training and testing set at a 2:1 ratio, stratified by deceased status (whether subjects were deceased or censored within 10 years) to balance the distribution of outcome status between sets. We used the 10% most variable CpG loci in methylation beta values across all samples. After splitting subjects and subsetting CpGs, a series of Cox proportional hazards models were fit using the training set for each selected CpG loci adjusted for age, age acceleration, sex, tumor grade, smoking status, BCG treatment status, and immune cell-type proportions associated with 10-year OS of patients with NMIBC. Next, Cox-scores (|β|/se(β), where β = the proportional hazards estimate of the log-HR, and se = the standard error) were computed for each of the selected CpG loci, and Cpg loci were ranked based on the Cox scores. Subsequently, the top M (range: 5–50) loci with the largest absolute Cox scores were chosen using a x-fold cross-validation RPMM with the smallest median P value of the log-rank test for each potentially optimal number (M) of CpG loci in the training set. Then, RPMM was fit to the testing set for clustering subjects using the optimal M-selected CpG loci with the largest absolute Cox score, predicting the methylation class membership for the subjects. Then, all patients with NMIBC were clustered using RPMM based on the methylation levels of the optimal CpG sites. Finally, we evaluated the association of RPMM class membership with OS using Cox proportional hazards models.

All datasets generated and analyzed during this current study are available in the Gene Expression Omnibus repository at GSE183920.

DNA methylation profiles were obtained from 601 peripheral blood samples from patients with NMIBC using the Human MethylationEPIC array. The study group was 455 (75.7%) men, 306 (50.9%) former-smokers, 192 (32.0%) current-smokers, 89 (14.8%) with BCG treatment, and had a median age of 66 (Table 1). The distribution of chronologic age, methylation age, and age acceleration is shown in Supplementary Fig. S1A. Cell-type proportions were estimated for each patient using methylation cytometry (Supplementary Fig. S1B). Stratifying on median time from diagnosis to blood draw, patient and tumor characteristic summary statistics showed similar distributions. Furthermore, we observed no significant associations of immune profile variables with time from diagnosis to blood draw. To assess the potential modification of results by time to the blood draw, we performed an analysis testing the relation of immune profile variables with patient outcome, stratifying patients into two groups based on median time to the blood draw. We did not observe differences between associations of immune variables with outcomes between the groups based on time to blood draw.

First, we examined associations of three major methylation age clocks, Horvath age (19), Hannum age (21), and DNAmPhenoAge (20), with bladder cancer outcomes. Our findings showed that Horvath age acceleration was not associated with 10-year OS (Supplementary Table S1). Because the Horvath clock was developed using data from a subset of CpG loci on the Illumina HumanMethylation27 (27K) CpG BeadChip (∼27,000 features) compared with the Hannum and PhenoAge clocks (which were developed using 20 times more features with data from CpGs on the Illumina HumanMethylation450 (450K) and EPIC (850K) BeadChips), we focused on Hannum and Pheno age acceleration in subsequent analyses. Then, we fit multivariable Cox proportional hazards models for demographic and tumor characteristic variables to investigate associations with 10-year OS and RFS. For patients with NMIBC, age, age acceleration, smoking, and high tumor grade were associated with worse RFS and OS. Women had a better survival outcome compared with men (Table 2; Supplementary Table S2).

Next, the association between immune cell-type proportions and bladder cancer patient outcomes was investigated. We fit multivariable Cox models for immune cell-type proportions associated with bladder cancer outcomes in Cox univariate models (Table 3; Supplementary Table S3). CD4T memory, CD8T memory, and NK cell proportions were associated with a decreased risk of death. Whereas NLR, CD8T naïve, and neutrophil cell proportions were each associated with an increased risk of death. In analyses of RFS, CD4T memory was associated with a decreased risk of tumor recurrence, and NLR, monocyte, and regulatory T cell proportions were associated with an increased risk of tumor recurrence (Table 3; Supplementary Table S3). Hazard estimates for associations of demographic and tumor characteristic variables with outcomes were similar after adjusting for cell composition. Because there was a stronger (smaller P value) association of immune cell profiles with 10-year OS compared with 10-year RFS, subsequent analyses are focused on 10-year OS of patients with NMIBC. Because BCG treatment may affect immune cell composition, we performed a sensitivity analysis limiting our analysis to patients not receiving BCG treatment and observed results consistent with those obtained for all patients with NMIBC (Table 3; Supplementary Table S4).

To partition the covariate space and explore interactions of clinical variables with immune cell proportions in 10-year OS, we applied a partitioning deletion/substitution/addition [partDSA (29)] algorithm for model building (Fig. 1A). partDSA is an analytic algorithm employing recursive partitioning. It uses loss functions to build clinically interpretable predictors of risk for a given event based on the covariate information. In brief, this method divides the covariate space into mutually exclusive regions and stratifies patients into distinct risk groups with respect to an outcome. As a result, it can make guidelines for estimating a patient's prognosis from clinical and biological information (42). Age, Hannum or Pheno age acceleration, sex, tumor grade, smoking status, BCG treatment status, and immune cell type proportions associated with 10-year OS were included in the model as input. The outcome was 10-year OS. The partDSA analysis divided subjects into three groups based on neutrophil and CD8 naïve cell proportions. The patients with NMIBC with neutrophil proportion >76.46 had the worst 10-year OS (Group 2, n = 17, HR = 4.93, 95% CI = 2.78–8.71) compared with Group 1 patients (patients with neutrophil cell proportions ≤76.46 and CD8 naïve cell proportions ≤1.76; n = 454; Fig. 1B). Although the 10- and 5-year OS rates (Fig. 1C) for Groups 1 and 3 were not statistically significantly different, we observed that their corresponding Kaplan–Meier curves separated after 5 years. To further investigate the difference between 10-year OS in Group 1 and Group 3, we selected the patients with NMIBC in Groups 1 and 3 whose time intervals from the date of initial diagnosis to death or the last follow-up were greater than 5 years (sample size: Group 1 = 401; Group 3 = 117), using the Kaplan–Meier method to observe 10-year OS. As we expected, for patients with time intervals from the date of initial diagnosis to death or last follow-up of more than 5 years, Group 3 patients (neutrophil proportion ≤76.46, CD8T naïve proportion >1.76) had better 10-year OS (HR = 0.37, 95% CI = 0.19–0.71) compared with Group 1 patients (Fig. 1D). Next, we examined the distribution of immune cell-type proportions, methylation age, and age acceleration in the three groups. Among the groups, patients with NMIBC in Group 2 had lower B naïve, eosinophil, monocyte, NK, CD4T memory, CD4T naïve, CD8T memory, and CD8T naïve proportions compared with other groups. Also, patients in Group 2 had higher neutrophil proportion, NLR, chronologic age, methylation age, and age acceleration compared with other groups (Supplementary Fig. S2).

To investigate whether we could identify a blood DNA methylation profile associated with NMIBC survival, we applied a SS-RPMM method. The workflow is illustrated in Fig. 2A and Supplementary Fig. S3A. The subjects in the testing set were assigned to cluster membership using the methylation profiles of the optimal CpG sites [Supplementary Table S5A (Pheno) and Supplementary Table S5B (Hannum)]. Then, all patients with NMIBC were clustered using RPMM based on the methylation levels of the optimal CpG sites, resulting in two classes, rR and rL [“R” and “L” corresponded with branches in the dendrogram; r stands for root; Fig. 2B (Pheno); Supplementary Fig. S3B (Hannum)]. Methylation class membership was significantly associated with 10-year OS; patients in cluster rR had a more favorable 10-year OS compared with those in cluster rL in both the testing set [HR = 0.35, 95% CI = 0.20–0.60; Fig. 2C (Pheno); HR = 0.38, 95% CI = 0.21–0.68; Supplementary Fig. S3C (Hannum)] and when using all patients with NMIBC [HR = 0.35, 95% CI = 0.25–0.48; Fig. 2D (Pheno); HR = 0.37, 95% CI = 0.27–0.52; Supplementary Fig. S3D (Hannum)]. Then, we compared the distribution of immune cell-type proportions, chronologic age, methylation age, and age acceleration. Consistent with the models in Table 3 and Supplementary Table S3, we observed that patients in cluster rR (n = 288) had significantly higher B memory, B naïve, CD4T memory, CD4T naïve, CD8T memory, and NK cell proportions and had significantly lower basophil, eosinophil, monocyte, neutrophil cell proportions, NLR, chronologic age, methylation age, and age acceleration compared with patients in cluster rL (n = 313; Supplementary Fig. S4A and S4B). The model using Hannum age acceleration had similar results shown in Supplementary Fig. S4C and S4D.

The two methods above used immune profiles (partDSA) and DNA methylation profiles (SS-RPMM) to explore the association of these profiles with OS, respectively. We were curious whether we could combine immune and methylation information to produce a guideline for estimating a patient's prognosis. To gain a deeper understanding of interactions between clinical variables, immune cell proportions, and blood DNA methylation profiles in 10-year OS, we allocated patients with NMIBC based on clustering results from partDSA and SS-RPMM analyses in five groups (Supplementary Table S6). Groups 1 and 3, from the partDSA analysis, were divided into two subgroups based on the SS-RPMM analysis (Fig. 3A; Supplementary Fig. S5A). In Kaplan–Meier analysis, patients in Group 2 still had the worst 10-year OS rate among all groups. Within the Group 1 patients, the group G1_rR had better 10-year OS than the group G1_rL [HR = 0.40, 95% CI = 0.27–0.58; Fig. 3B (Pheno); HR = 0.45, 95% CI = 0.31–0.65; Supplementary Fig. S5B (Hannum)]. Consistent with Group 1 patients, in Group 3 patients, the group G3_rR had better 10-year OS compared with the group G3_rL [HR = 0.42, 95% CI = 0.18–0.95; Fig. 3B (Pheno); HR = 0.35, 95% CI = 0.15–0.79; Supplementary Fig. S5B (Hannum)].

This study aimed to test the relation of bladder cancer outcomes with high-resolution immune profiles using methylation cytometry and cell-independent methylation states in blood. In prior work, we observed that CD4T and CD8T proportions were associated with the decreased risk of death and tumor recurrence in patients with NMIBC (12); however, the results were not generalizable to naïve and memory subtypes of CD4T and CD8T cells. With recent advances in methylation cytometry for immune profiling (15), we were able to examine the association between the proportions of regulatory T cells, eosinophils, basophils, naïve and memory subtypes of CD4T, CD8T, and B cell and bladder cancer outcomes. Consistent with our previous study, NLR was associated with an increased risk of death and tumor recurrence in patients with NMIBC. For CD4T and CD8T cell subsets, CD4T memory and CD8T memory cell proportions were associated with a decreased risk of death. However, only CD4T memory cell proportion was associated with a reduced risk of death and tumor recurrence. One possible explanation for the different observations is that we introduced age acceleration into the models. This might confound the association between CD8T cell proportions subtypes and bladder cancer outcomes. Here, we identified associations between immune cell subtypes and age acceleration with bladder cancer outcomes. These factors could be potentially prognostic biomarkers of bladder cancer.

Few studies have shown age acceleration from multiple age clocks to be associated with bladder cancer outcomes, and even then, they do not show consistent results (22, 27). In our study, the direction of hazard estimates among cell type proportions was similar. However, outcome associations with immune profile differed in Cox multivariable models controlling for Pheno and Hannum age acceleration. The observed difference was potentially due to the training methods of clocks. Unlike Hannum age using chronologic age as a surrogate, Pheno age mainly focuses on aging outcomes, such as cancers, diet, and all-cause mortality, hence, Pheno age can capture age-related outcomes and perform well in predicting survival compared with chronologic methylation clocks, such as Hannum and Horvath clocks. Horvath age acceleration was not associated with bladder cancer outcomes. One potential explanation is that the Horvath clock was built mainly using a subset of CpGs on the 27K methylation array platform, which had approximately 20 times fewer CpGs than the 450K and EPIC array platforms used for the developing of the Hannum and PhenoAge clocks. Furthermore, various cell and tissue types were used to develop the Horvath clock. However, Hannum and PhenoAge clocks were built based on data from blood DNA measures.

Our findings suggest that elevated NLR, neutrophil, basophil, regulatory T, and decreased CD4T memory cell proportions increased the risk of death and tumor recurrence in patients with bladder cancer. These findings met our expectations, consistent with previous studies demonstrating that peripheral blood NLR levels were associated with an increased risk of NMIBC recurrence after surgery (47–49). Moreover, basophil count was significantly associated with an increased risk of recurrence in patients with BCG-treated bladder cancer (50). One study indicated that peripheral (neutrophil × platelet)/(lymphocyte) was inversely correlated with a high risk of tumor recurrence in NMIBC (51). Even though some peripheral immune profiles we found had not been revealed to be associated with bladder cancer outcomes, these immune profiles in the tumor microenvironment had been reported to be associated with bladder cancer outcomes. For instance, higher regulatory T cell infiltration in the tumor microenvironment was associated with a shorter RFS (52), and regulatory T cell frequency within the tumor was inversely correlated with RFS in patients with NMIBC (53). Future work that integrates the assessment of cell-type proportions in the tumor microenvironment and periphery associated with bladder cancer outcomes would be of value.

The optimal 15 CpG loci selected by SS-RPMM using the model adjusting for Pheno age acceleration, and the optimal 50 CpG loci selected by SS-RPMM using the model adjusting for Hannum age acceleration track to several genes that have been reported to be involved in bladder cancer development. Sprouty-related EVH1 domain-containing protein 2 (SPRED2), is a negative regulator of the ERK-MAPK pathway, and has been reported to have increased mRNA and protein expression in NMIBC compared with carcinoma in situ and infiltrating urothelial carcinoma (54). In addition, patients with higher SPRED2 mRNA levels had better OS compared with low expression group (54). Peroxisome proliferator-activated receptor gamma (PPARG) high-activation has been reported to promote cell-cycle G₂ arrest and apoptosis, leading to suppression of tumor growth and better prognosis in patients with bladder cancer (55). Fibrous sheath interacting protein 1 (FSIP1) was overexpressed in protein and mRNA levels of bladder tumor tissues and cancer cell lines (56, 57). Besides, FSIP1 overexpression was associated with worse outcomes (56). Phosphorylated MAPK 14 (MAPK14) was overexpressed in bladder cancer cell lines and tissues (58). Furthermore, phosphorylated MAPK14 could combine and regulate RUNX2, which was identified in our previous study through epigenome‐wide association study (12), to promote the proliferation and migration of bladder cancer (58). Consistent with our findings, patients with a worse survival rate (rL group) had lower methylation levels in the CpG site, cg16145324, located in the MAPK14 gene region (Supplementary Table S5B).

To further investigate the interactions between clinical variables, immune profiles, and DNA methylation levels, we performed both partDSA and SS-RPMM. Though patients with NMIBC were divided into three groups using partDSA, the Kaplan–Meier curves for Group 1 and 3 patients were not significantly different. Interestingly, when we applied the optimal CpG loci selected by SS-RPMM, patients in Groups 1 and 3 were grouped into two groups, with a significant difference in 10-year OS. This finding illustrated the importance of methylation levels of specific CpG sites in evaluating bladder cancer prognosis.

While this study carefully evaluated the association of cancer outcomes with peripheral immune cell type proportions, there were potential study limitations. For instance, BCG treatment has been reported to affect immune cell composition (50, 51, 59) and methylation profiles (60). We do not have detailed information, such as cycle numbers or responsiveness for each patient, though only a few patients (14.8%) in our study received BCG treatment. In addition, our sensitivity analysis that limited to patients who did not receive BCG treatment demonstrated consistent results with the overall NMIBC group. In addition, case ascertainment for our population-based study resulted in a median time from diagnosis to study blood draw is 319 days in our dataset. Furthermore, some confounding factors, such as obesity (61), alcohol consumption (62), and type 2 diabetes (63), have been reported to affect peripheral immune cell distribution. However, we had incomplete information on these potential covariates.

Taken together, our analysis applied the latest differentially methylated region library and highlighted several peripheral immune profiles that were associated with bladder cancer outcomes. We assessed interactions between clinical variables, immune cell proportions, and blood DNA methylation profiles in 10-year OS using partDSA and SS-RPMM analyses, clustering patients with NMIBC into five groups according to the methylation levels of the optimal CpG loci, neutrophil, and CD8T naïve cell proportions. While few studies have investigated the association between bladder cancer outcomes and peripheral immune profiles, our findings provide insight into the potential of peripheral immune profiles to serve as prognostic biomarkers in bladder cancer. Future work examining immune profiles in tumor tissues as well as DNA methylation profiles of patients with bladder cancer is needed to integrate interactions between bladder cancer outcomes, methylation levels, peripheral immune environment, and tumor microenvironments to further validate the feasibility of methylation-derived immune profiles for epigenetic biomarkers of bladder cancer.

L.A. Salas reports grants from CDMRP/Department of Defense (W81XWH-20-1-0778) and NIGMS (P20 GM104416) during the conduct of the study; in addition, has a patent for U.S. Provisional Patent Application Enhanced DNA Methylation Library for Deconvoluting Peripheral Blood (63/148,695) pending. J.D. Seigne reports grants from NIH and CDMRP/Department of Defense during the conduct of the study; other support from Johnson & Johnson outside the submitted work. K.T. Kelsey reports other support from Cellintec outside the submitted work; in addition, has a patent 10,619,211 issued. B.C. Christensen reports grants from NIH during the conduct of the study; personal fees from University of California San Fransico outside the submitted work; in addition, has a patent for Enhanced DNA Methylation Library for Deconvoluting Peripheral Blood pending. No disclosures were reported by the other authors.

The funders had no role in study design, data collection, data analysis, or data interpretation.

J.-Q. Chen: Conceptualization, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. L.A. Salas: Conceptualization, formal analysis, supervision, investigation, methodology, writing–review and editing. J.K. Wiencke: Conceptualization, methodology, writing–review and editing. D.C. Koestler: Conceptualization, methodology, writing–review and editing. A.M. Molinaro: Conceptualization, methodology, writing–review and editing. A.S. Andrew: Data curation. J.D. Seigne: Data curation. M.R. Karagas: Conceptualization, methodology, writing–review and editing. K.T. Kelsey: Conceptualization, methodology, writing–review and editing. B.C. Christensen: Conceptualization, resources, formal analysis, supervision, funding acquisition, validation, methodology, project administration, writing–review and editing.

B.C. Christensen is supported by the NIH grant numbers R01CA216265. B.C. Christensen, M.R. Karagas, and L.A. Salas are supported by the NIH grant number P20GM104416. M.R. Karagas is supported by the NIH grant number R01CA057494. B.C. Christensen and K.T. Kelsey are supported by the NIH grant number R01CA253976. L.A. Salas is supported by the CDMRP/Department of Defense W81XWH-20-1-0778. J.K. Wiencke and A.M. Molinaro are supported by the NIH grant number P50CA097257. J.K. Wiencke is supported by the Robert Magnin Newman Endowed Chair in Neuro-oncology and the NIH grant numbers R01CA207360. D.C. Koestler is supported by the NCI Cancer Center Support Grant P30 CA168524 and the Kansas Institute for Precision Medicine COBRE (supported by the National Institute of General Medical Science award P20 GM130423).

The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.