DNA Methylation Modification Map to Predict Tumor Molecular Subtypes and Efficacy of Immunotherapy in Bladder Cancer

Background: Considering the heterogeneity and complexity of epigenetic regulation in bladder cancer, the underlying mechanisms of global DNA methylation modification in the immune microenvironment must be investigated to predict the prognosis outcomes and clinical response to immunotherapy.

Methods: We systematically assessed the DNA methylation modes of 985 integrated bladder cancer samples with the unsupervised clustering algorithm. Subsequently, these DNA methylation modes were analyzed for their correlations with features of the immune microenvironment. The principal analysis algorithm was performed to calculate the DMRscores of each samples for qualification analysis.

Findings: Three DNA methylation modes were revealed among 985 bladder cancer samples, and these modes are related to diverse clinical outcomes and several immune microenvironment phenotypes, e.g., immune-desert, immune-inflamed, and immune-excluded ones. Then patients were classified into high- and low-DMRscore subgroups according to the DMRscore, which was calculated based on the expression of DNA methylation related genes (DMRGs). Patients with the low-DMRscore subgroup presented a prominent survival advantage that was significantly correlated to the immune-inflamed phenotype. Further analysis revealed that patients with low DMRscores exhibited less TP53 wild mutation, lower cancer stage and molecular subtypes were mainly papillary subtypes. In addition, an independent immunotherapy cohort confirmed that DMRscore could serve as a signature to predict prognosis outcomes and immune responses.

Conclusion: Global DNA methylation modes can be used to predict the immunophenotypes, aggressiveness, and immune responses of bladder cancer. DNA methylation status assessments will strengthen our insights into the features of the immune microenvironment and promote the development of more effective treatment strategies.

Introduction

DNA methylation modification is one of the most representative epigenetic modifications, which is indispensable in vertebrate development and illnesses (Ortega-Recalde and Hore, 2019; da Rocha and Gendrel, 2019). Besides, it has been demonstrated to associate with multiple biological functions in cancer, e.g., the formation and evolution of tumor microenvironment, as well as impairment restoration in the immune cycle (Chen et al., 2020; Zhang et al., 2020a). On the other hand, abnormal DNA methylation is also significantly related to the occurrence of multiple cancer types, such as sarcoma (Koelsche et al., 2021), bladder cancer (Liu et al., 2021a), and vulvar intraepithelial neoplasia (Thuijs et al., 2021).

Bladder cancer is an extremely malignant urogenital neoplasm (Siegel et al., 2021). Heterogeneous distributions of genome clusters lead to molecular and cellular heterogeneity in tumors, which affect clinical outcomes and treatment responses (Qiu et al., 2019; Craig et al., 2020). Despite pronounced progress in the treatment of bladder cancer, more effective therapeutic strategies are still in demand. Studies have demonstrated that several genes involved in the occurrence and progression of bladder cancer are regulated by promoter methylation. For example, Chen X et al. built a diagnostic model based on 2 DNA methylation markers for early detection and recurrent monitoring of bladder cancer. Wilhelm CS et al. discovered that LINE1 hypomethylation may contribute to bladder cancer tumorigenesis, especially in women (Wilhelm et al., 2010). Kandimalla R et al. summarized the biomarkers of DNA methylation and identified that methylated genes, including SFRP1, SOX9, FHIT, CDH1, PMF1, RUNX3, LAMC2, and RASSF1A, are related to the poor clinical outcomes in bladder cancer patients (Kandimalla et al., 2013). In short, DNA methylation is involved in carcinogenesis or tumor inhibition across varying scenarios.

In recent years, immune checkpoint blockade therapy has emerged as a promising therapeutic strategy. It aims to enhance the immune activity of T lymphocytes to kill tumor cells by inhibiting immune checkpoints, such as PD-1 and its ligand PD-L1 (Yi et al., 2018). Studies have shown that immunotherapy could improve clinical outcomes of numerous tumors, such as ovarian cancer (Wan et al., 2021), bladder cancer (Han et al., 2021), and colorectal cancer (Liu et al., 2021b). However, patients respond to immunotherapy differently, and the effective rate of immune checkpoint blockade therapy has been less than 20%. The expression level of PD-L1, tumor microenvironment (TME), and tumor mutation burden (TMB) have been reported as signatures to evaluate the clinical responses to immunotherapy (Samstein et al., 2019). Previous studies demonstrated that DNA methylation may contribute to the alteration of TME. For instance, Sasidharan Nair V et al. discovered that DNA hypomethylation shall alter the expression of CTLA-4, TIGIT, and PD-1 genes (Sasidharan Nair et al., 2018). Elashi AA et al. also revealed that an abnormal promoter methylation profile is correlated to the peripheral upregulation of TIGIT and PD-1 in many cancers. They speculated that a combined administration of anti-PD-1 agents and demethylation inhibitors could be a more effective immunotherapeutic strategy than the current ones (Elashi et al., 2019). However, the regulatory mechanisms of global DNA methylation on tumor microenvironment and immune response in bladder cancer remain unclear.

In this study, genomic data and clinical information of 985 samples from six independent bladder cancer cohorts were included. DNA methylation modes were clarified by analyzing the expression of fifteen DNA methylation regulators in these samples We investigate the DNA methylation regulators rather than DNA methylation itself, because the biology function of DNA methylation would be altered according to the genomic environment. Specifically, three DNA methylation modes were identified to meet the criteria of immune-desert, immune-inflamed, and immune-excluded immunophenotypes, respectively. Moreover, an evaluation system was built to qualify the DNA methylation modes in individual patients, and the patients’ clinical responses to immunotherapy were assessed based on their DMRscore. Our study provides a new perspective to observe the global DNA methylation status and the immunophenotype of individual tumors in bladder cancer so that more specified precision medicine could be achieved.

Results

The Landscape of DNA Methylation Regulators in Bladder Cancer

We executed systematic research that included 15 DNA methylation regulators and summarizes the mutation rates of all these regulators in bladder cancer. Among 412 samples, 52 samples experienced alteration of DNA methylation regulators, with frequency 12.62%. According to the waterfall diagram, alterations of the MBD1 gene were the most frequent, and these alterations have been reported to participate in tumorigenesis. Besides, DNMT1, DNMT3A, and DNMT3B genes also exhibited an alteration frequency of 2% (Figure 1A). Furthermore, co-occurrence mutation was observed in several DNA methylation regulators despite their functional differences, including NTHL1, MBD3, MECP2, UHRF2, and ZBTB33 (Supplementary Figure S1C).

FIGURE 1

FIGURE 1. Landscape of genetic alteration and transcriptome variation of DNA methylation regulators in bladder cancer. (A) The alteration frequency of 15 DNA methylation regulators in 412 bladder cancer samples (TCGA-BLCA). the annotation of each variant types was displayed by the bagplots right barplots. Each cohort represented an individual sample. The stacked barplot below displayed conversion ratio for each sample. (B) The location of CNV alteration of DNA methylation regulators on 23 chromosomes was displayed by circular plot. (C) Principal component analysis (PCA) for the transcriptome characteristics of 15 DNA methylation regulators to distinguish tumors from normal samples in GSE13507 cohort. Tumor samples were labeled with blue color and normal samples were labeled with yellow. (D) The frequency of copy number variation in TCGA-BLCA cohort. Deletion frequency: the green dot and amplification frequency: red dot. The number represented the variation frequency. (E) The transcriptome characteristics of 15 DNA methylation regulators between normal and bladder cancer tissues. Tumor: red box; Normal: blue box. The median value: black lines in boxes, the outliers: black dots out boxes. The asterisks represented the statistical p value (*p < 0.05; **p < 0.01; ***p < 0.001).

In addition, a prevalent CNV alteration was observed in the fifteen regulators (Figure 1D). Specifically, DNMT3B, UHRF2, and MECP2 demonstrated a widespread frequency of amplification in samples while MBD3, UHRF1, and NTHL1 were frequently detected. The locations and circle sequences of the DNA methylation regulators along the chromosomes are depicted in Figure 1B. Moreover, the principal component analysis revealed that the bladder cancer can be distinguished from normal samples by observing the expression levels of the 15 regulators (Figure 1C). In addition, the mRNA expressions of the 15 DNA methylation regulators are also significantly different between BLCA tumors and normal tissues (Figure 1E). In a word, the genomic imbalance of DNA methylation regulators is vital for bladder cancer tumorigenesis and development.

DNA Methylation Regulator Clusters

The complete clinical information and transcriptome data of six GEO datasets (GSE13507, GSE31684, GSE32548, GSE48075, GSE48476, GSE80691) and TCGA-BLCA were enrolled into one cohort for further exploration. A univariate Cox regression analysis was performed to find the prognostic value of the 15 DNA methylation regulators in bladder cancer patients (Supplementary Figure S1B). Following that, the comprehensive landscape of the regulators’ intercorrelation and their prognostic attributes for bladder cancer were calculated by network planning (Figure 2B). From these results, we speculated that DNA methylation regulators may be related to the heterogeneity of bladder cancer. Therefore, unsupervised clustering was performed to explore ultramodern DNA methylation regulator clusters (DMRclusters) based on the expression levels of the regulators in the meta-cohort. Three DMRclusters were classified, including 306, 348, and 331 sample patients in DMRcluster A, B, and C, respectively, and these distinct DMRclusters could be distinguished via the principal component analysis (Figure 2D). Specifically, DMRcluster A presented a particularly prominent survival advantage, but DMRCluster B exhibited the worst clinical outcome in the integrated cohort (Figure 2A).

FIGURE 2

FIGURE 2. Clusters of DNA methylation modes and biological profiles of each cluster. A Kaplan-Meier curve with p value 0.032 displayed a remarkable difference among three DNA methylation modes, the DMRcluster B presented a remarkable poor clinical outcome. DMRcluster (A): 306 samples, DMRcluster (B): 348 samples and DMRcluster (C): 331 samples. The meta cohort including 985 samples (GSE13507, GSE31684, GSE32548, GSE48075, GSE48476, GSE80691 and TCGA-BLCA). (B) The interplay among DNA methylation regulators in bladder cancer. Red and gary represented readers and writers respectively. The size of circles displayed the influence of each regulator on clinical outcomes. The lines connecting regulators represented their interactions, and thickness represented the correlation strength. Negative correlation was labeled with blue and positive correlation was labeled with red. Risk factor: purple, favorable factor: green. (C) The proportion of molecular subtypes in the three DNA methylation modes (TCGA-BLCA). Basal squamous subtype, green; Luminal subtype, blue; luminal infiltrated subtype, red; luminal papillary subtype, yellow; and Neuronal subtype, olivedrab. (D) Principal component analysis (PCA) for the transcriptome characteristics of three DMRclusters. DMRcluster A was labeled with blue color, DMRcluster (B) was labeled with red color and DMRcluster (C) was labeled with green. (E,F) Gene set variant analysis displayed the activation status of biological pathways in diverse DMRclusters. The heatmap help us to observe the difference of biology pathway activity among three DMRclusters. Blue represented inhibited pathways and red represented activated pathway. (E) DMRcluster (A) vs. DMRcluster (B); F DMRcluster (B) vs. DMRcluster (C).

The Immune Features of Distinct DMRclusters

The GSVA enrichment analysis was performed to identify the biological processes in the DMRclusters. Judging from the results, DMRcluster A was enriched in immune activation pathways, such as the T/B cell receptor signaling pathway, toll-like receptor signature, as well as complement and coagulation cascades. On the other hand, DMRcluster B was prominently associated with immune suppression, and DMRcluster C was even more prominent in carcinogenic pathways, including the P53 signature pathway and the ERBB signature pathway (Figures 2E,F). As expected, subsequent analyses revealed that DMRcluster A was significantly enriched in cells related to acquired immunity, including activated B cells, central memory CD4/CD8 T cells, and activated dendritic cells (Figure 3B). Such a finding could well explain the results of the survival analysis. Meanwhile, stromal activity (e.g., epithelial-mesenchymal transition &, EMT) was remarkably enriched in DMRcluster C (Figure 3A). Based on all the previous results, we speculated that these DMRclusters had remarkably diverse features in terms of immune cell infiltration into the tumor microenvironment. Specifically, DMRclusters A, B, C were featured by immune-inflamed, immune-desert, and immune-excluded phenotypes, respectively.

FIGURE 3

FIGURE 3. Tumor microenvironment characteristics and transcriptome profile in three DNA methylation modification modes. (A) stromal activation pathways among three different DNA methylation modification modes include EMT, angiogenesis and Pan-F-TBRS. (B) The content of each tumor microenvironment immune infiltrating cells in three DNA methylation modification modes. The median value: black lines in boxes, the outliers: black dots out boxes. (C) The heatmap help us to observe the expression level of DNA methylation regulators among different DMRclusters (TCGA-BLCA cohort). DMRcluster subtypes, Molecular subtypes, Histology, Grade, Stage, Gender, Age and Survival status were used as patient annotations. Red represented high expression level of DNA methylation regulators and blue represented low expression level. (D) Functional annotation for DNA methylation-related genes. Red represented more enriched genes and blue represented a small number of enriched genes.

The Transcriptome Data and Clinical Features of DMRclusters

To further investigate these DMRclusters in diverse biological processes and clinical features, we focused on the TCGA-BLCA cohort containing 407 bladder cancer patients and their exhaustive clinical information. Similarly, the patients were classified into three clusters with unsupervised clustering (Supplementary Figures S2A–D). Judging from the results DNA methylation regulators’ transcriptional profiles among the three DMRclusters demonstrated significant difference, which was validated by one-way ANOVA analysis (Supplementary Figure S2E). Specifically, DMRcluster A revealed high expression of DNMT3B and DNMT3A, DMRcluster B was characterized by higher DNMT1 and UHRF1 expressions, and DMRcluster C exhibited lower contents of DNMT1, DNMT3A, DNMT3B, and UHRF1 at various extents (Figure 3C). Patients with the luminal infiltrated subtype were characterized by DMRcluster A, while the basal squamous subtype was featured by DMRcluster B (Figure 2C); besides, both DMRclusters B and C were enriched in the neuronal subtype (Figure 4D). In bladder cancer treatments, the neuronal subtype is particularly difficult because of its poor clinical outcome, while the luminal papillary subtype is prone to better survival. Thus, we performed the K-M analysis, and the results also validated our conjecture that patients characterized by DMRclusters B and C exhibited significantly more rapid disease progression and poorer clinical outcomes, while DMRcluster A presented a remarkable survival advantage (Supplementary Figure S2F). In addition, the luminal infiltration subtype in bladder cancer is characterized by low tumor cell purity and high lymphocytic infiltration. Most patients with the luminal infiltration subtype were categorized into DMRcluster A, and only a small amount of luminal infiltration was observed in DMRcluster B (Figure 4D), suggesting that DMRcluster A is related to immune activation and DMRcluster B is associated with the immune-desert phenotype.

FIGURE 4

FIGURE 4. Construction of DNA methylation signatures for individual sample. (A) The heatmap help us to observe the transcriptome landscape among different Gene.clusters (TCGA-BLCA cohort). Gene.cluster subtypes, Molecular subtypes, Histology, Grade, Stage, Gender, Age and Survival status were used as patient annotations. Red represented high expression level and blue represented low expression level. B-C Differences in DMRscore among three DMRclusters (B) or Gene.clusters (C) in TCGA-BLCA cohort (Kruskal-Wallis test, p < 0.001). D-E Kaplan-Meier curve displayed a remarkable difference among three Gene. clusters ((D), p < 0.001) or DMRscore subgroups ((E), p < 0.001) in TCGA-BLCA cohort. (F) Sankey diagram displayed the alteration of DMRclusters, molecular subtypes, Gene.cluster and DMRscore. (G) The transcriptome characteristics of 15 DNA methylation regulators among Gene.clusters. Gene.cluster (A): blue box; Gene.cluster (B): red box. Gene.cluster C: green box. The median value: black lines in boxes, the outliers: black dots out boxes. (H) Differences in the known gene signatures between high DMRscore and low DMRscore subgroups. APC: antigen-presenting cells. The asterisks represented the statistical p value (*p < 0.05; **p < 0.01; ***p < 0.001).

Functional Annotations of DMRGs

To further explore the potential biological processes in each DMRcluster, the R package named “limma” was performed to find DMRGs, and a total of 832 genes were selected (Supplementary Figure S2G). GO analysis was executed on the DMRGs using the R package “clusterProfiler.” DMRGs were prominently enriched in immunity activation pathways, DNA methylation, and cell proliferation, which verified that DNA methylation is vital in the immune regulation of tumor progression (Figure 3D).

To further investigate the mechanisms of DNA regulation, the patients were classified into three genomic subtypes based on the expression of the 832 DMRGs. Similarly, the genomic subtypes were identified via the unsupervised clustering algorithm. They were termed Gene cluster A, B, and C, respectively (Supplementary Figures S3A–D), and they were all related to DNA methylation in bladder cancer. A heat map also demonstrated that the three Gene clusters can be distinguished by their signature transcriptomes (Figure 4A). According to the K-M survival method, Gene cluster A presented a remarkable survival advantage, while Gene cluster B was proved to be associated with a poorer prognosis (Figure 4D). Moreover, the three Gene clusters revealed significant differences in the expression of DNA methylation regulators (Figure 4G).

TME Characteristics in the Three Gene.clusters

To identify the role of Gene clusters in the immune regulation of TME, we investigated the expression of cytokines and chemokines in Gene clusters. The targets of identification were chosen from the literature, among which ZEB1, TGFB2, PDGFRA, VIM, COL4A1, TGFBR2, TWIST1, ACTA2, and SMAD9 are related to transcripts of the transforming growth factor (TGF) b/EMT pathway. Besides, HAVCR2, CD80, LAG3, CD86, TIGIT, PDCD1, TNFRSF9, PD-L1, IDO1, CTLA4, and PD-L2 are associated with the transcripts of immune checkpoints, and CXCL10, PRF1, CD8A, CXCL9, GZMB, GZMA, TNF, IFNG, and TBX2 are associated with immune-activated transcripts (Sotiriou et al., 2006; Barbie et al., 2009; Ritchie et al., 2015; Zeng et al., 2019).

We found that the transcripts related to immune activation pathways were significantly up-regulated in Gene cluster B, but the patients in this cluster did not show an expected survival advantage. Previous studies revealed that high stromal activation was associated with limited immune activation (MacGregor et al., 2019). Therefore, we investigated the transcripts related to the (TGF)b/EMT pathway in this cluster and demonstrated stromal activation within. Based on these findings, we assumed that anti-tumor effects of immune cells in Gene cluster B are limited by stromal activation, indicating that Gene cluster B is the immune-excluded subtype. Besides, the transcripts of immune checkpoints were examined as highly expressed in Gene cluster B, suggesting that immunotherapy may bring unexpected outcomes (Supplementary Figures S3F–H).

To further investigate the functions of DMRGs, we examined the identified pathways in bladder cancer patients. Gene cluster A was found to enrich in CD8 T effector, DNA replication, mismatch repair, and antigen processing machinery pathways (Supplementary Figure S3E). Previous studies demonstrated that bladder cancer can be classified into five subtypes according to the molecular phenotype. Among them, the luminal-papillary subtype exhibits the best prognosis with a five-year survival rate of 60%. On the other hand, the five-year survival rate of neuronal bladder cancer is only 17% (Robertson et al., 2017). Our findings suggested that Gene cluster A was almost fully composed of the luminal-papillary subtype, which was relevant to survival advantage (Figure 4F).

Individual Modification Patterns of DNA Methylation

By now, the experimental results have confirmed that DNA methylation is irreplaceable in the formation of distinct TME landscapes. However, investigations above were not helpful to predict the DNA methylation status of an individual sample as they were conducted on a population. Since tumors are heterogeneous and complex, we built the DMRscore model to qualify the DNA methylation status based on the expression of DMRGs.

In this section, we attempted to assess whether DMRscore is effective in predicting clinical outcomes. Patients were classified into groups of high and low DMRscores according to the best cutoff value. The correlation results between DMRscores and clinical outcomes showed that patients in the low DMRscore group exhibited a remarkable clinical advantage, while those in the other group demonstrated less satisfactory clinical outcomes (TCGA-BLCA cohort) (Figure 4E, p = 7.158e-05). The Figure 4H shown that patients with high DMRscore were enriched in APC_co_inhibition, T_cell_coinhibition, which revealed that this subgroup presented immunosuppression. Subsequently, we examined whether the DMRscore can serve as an independent index to evaluate the clinical outcomes of bladder cancer. Multivariate Cox regression analysis was used to take the independent indices, including age and DMRscore, into the calculation, and the results confirmed DMRscore as an independent and robust prognostic index (HR = 1.05; Supplementary Figure S4A). Variations of individual patients are displayed by the Sankey diagram (Figure 4F).

To reassure the predictive effects of DMRscore, we examined its relationship with the identified clusters by Kruskal-Wallis tests. The test results suggested that DMRscore could be used to predict DNA methylation clusters. Specifically, both DMRcluster B and Gene cluster B showed the highest median DMRscore (Figures 4B,C). In addition, patients suffering from neuronal bladder cancer also exhibited the highest median DMRscore among five molecular subtypes (Figure 5A). In a word, DMRscore has been proved as an effective index to assess the DNA methylation status of individual samples and predict clinical outcomes. In order to develop the accuracy of predictive performance, the prognostic nomogram included a DMRscore, and other clinical variables was constructed to evaluate the 1-, 3-, and 5-year overall survival probabilities (Supplementary Figure S6).

FIGURE5

FIGURE5. Characteristics of DMRscore in TCGA molecular subtypes and tumor mutation burden. (A) Differences in DMRscore among diverse bladder cancer molecular subtypes. Basal squamous subtype, blue; Luminal subtype, red; luminal infiltrated subtype, green; luminal papillary subtype, sapphire; and Neuronal subtype, yellow. (B) Kaplan-Meier curve showed the clinical prognosis of patients with combination of DMRscore and adjuvant chemotherapy stratification. H, high. L, low. ADJC, adjuvant chemotherapy (p = 0.028). (C) Kaplan-Meier curve showed the clinical prognosis of patients with combination of DMRscore and TP53 stratification. H, high. L, low. MT, mutation type; WT, wild type (p < 0.001). (D) Difference in tumor mutation burden between high DMRscore and low DMRscore (p = 0.019). (E) Correlation between DMRscore and tumor mutation burden. (R = 0.25, p < 0.001) (F,G) The landscape of tumor somatic mutation in TCGA-BLCA established by high (F) and low DMRscore (G). Each column represented individual patients. The upper barplot displayed tumor mutation burden.

Particularly, the capability of DMRscore to assess the efficacy of adjuvant chemotherapy (ADJC) in bladder cancer patients was evaluated. DMRscore prediction results were not disturbed by ADJC: whether receiving ADJC or not, the low-DMRscore group always presented significant survival advantages. However, DMRscore cannot be utilized to judge whether ADJC can be applied on a bladder cancer patient, and patients with low DMRscores had shorter survival after ADJC. (Figure 5B). In addition, patients with high grade, TP53MT, and non-papillary subtypes of the cancer showed significantly higher DMRscores, with a poorer survival prognosis (Figure 6A). This also validated in E-MTAB-4321 cohort (Supplementary Figure S4C). Furthermore, the capability of DMRscore to assess the efficacy of TP53 mutation in bladder cancer patients was examined as well. We found that the L. DMRscore-TP53. WT group exhibited a remarkably advantageous survival, while the H. DMRscore-TP53. MT group demonstrated the worst clinical outcome (Figure 5C). K-M survival analysis and multivariate Cox regression analysis for the E-MTAB-4321 cohort also verified that DMRscore can serve as an independent prognostic index in bladder cancer (Figure 6C; Supplementary Figure S4B).

FIGURE 6

FIGURE 6. Role of DNA methylation modification in clinical prediction. (A) Differences in DMRscore among different clinical status. The median value: black lines in boxes, the outliers: black dots out boxes. MT, mutation type; WT, wild type. (B) Kaplan-Meier curve showed the clinical prognosis of patients with combination of DMRscore and NEO stratification. H, high. L, low. NEO, Newantigen burden (p < 0.001). (C) Survival analyses for high (135 samples) and low (341 samples) DMRscore subgroups in the E-MTAB-4321 cohort using Kaplan-Meier curves (p < 0.001). (D) Survival analyses for high (430 samples) and low (151 samples) DMRscore subgroups in the GEO-metacohort cohort using Kaplan-Meier curves (p < 0.001). (E) Survival analyses for high (147 samples) and low (151 samples) DMRscore subgroups in the anti-PD-L1 cohort (IMvigor210 cohort) using Kaplan-Meier curves (p = 0.004). (F,G) The proportion of patients with response to PD-L1 blockade immunotherapy in low or high DMRscore subgroups. SD, stable disease; PD, progressive disease; CR, complete response; PR, partial response. Responser/Nonresponer: 32%/68% in the low m6Ascore groups and 14%/86% in the high m6Ascore groups. H Differences in PD-L1 expression between low and high m6Ascore groups (p < 0.0001). L Differences in DMRscore among distinct tumor immune phenotypes in IMvigor210 cohort. (p = 0.029).

The Role of DNA Methylation Mode in anti-PD-1/PD-L1 Immunotherapy

Studies have verified that patients’ response to immunotherapy is related to the TMB frequency, and higher TMB statuses lead to a persistent response to anti-PD-1/PD-L1 immunotherapy. We investigated the somatic mutation frequencies between high- and low-DMRscore groups in the TCGA-BLCA cohort. However, TMB quantification analysis verified that DMRscore is significantly and positively correlated to TMB (Figures 5D,E). Besides, the waterfall diagram showed that the high-DMRscore group was more susceptible to somatic mutations than the other group, with somatic mutation frequencies of 146/148 (98.65%) and 234/251 (93.23%), respectively (Figures 5F,G). Thus, our experimental results indicated that patients with a high DMRscore exhibit good response to anti-PD-1/PD-L1 immunotherapy, which is contradictory to previous findings. Consequently, we speculated that TMB frequency cannot be utilized to predict the effect of immunotherapy in this model.

In order to further examine the prediction performance of the DMRscore model, we applied the established DMRscore signature to other independent bladder cancer cohorts. Almost all cohorts presented survival differences as revealed by the DMRscore model except for two GEO datasets with few samples (Supplementary Figures S5A–E). The prediction performance of the DMRscore model for tumor stages was assessed by the receiver operating characteristic (ROC) curves, and the area under the curve (AUC) was 0.699 and 0.721at 3 and 5 years, respectively (Supplementary Figures S5F,G). These data indicated that the DMRscore signature could serve as a new biomarker to predict clinical outcomes.

Immune checkpoint blockade therapy has undoubtedly produced significant therapeutic benefits for many cancer patients. Based on the collected immunotherapy cohorts, we explored whether DMRscore can serve as a signature to predict patients ' response to immunotherapy. Remarkable survival advantages are seen in patients with low DMRscores in the GEO-metacohort and anti-PD-L1 cohorts (IMvigor210, advantaged urothelial cancer) (Figures 6D,E). In addition, these patients also exhibited much better clinical outcomes when receiving anti-PD-L1 immune checkpoint blockade therapy (Figure 6F), and they were characterized by a significantly higher expression of PD-L1, a potential clinical response to immunotherapy (Figure 6H). Besides, we evaluated the tumor neoantigen burden in bladder cancer patients, and those with high neoantigen burden and low DMRscore signatures presented a significant survival advantage. (Figure 6B). Judging from Figure 6G, the DMRscore signature was a robust and potential biomarker to estimate patient response and clinical outcomes in immunotherapy. The immunophenotypes of metastatic urothelial cancer have been distinguished in the IMvigor210 cohort, so we studied the differences of DMRscore among them (Figure 6I). Most patients with low DMRscores exhibited the inflamed immunophenotype, to which individualized immunotherapy is crucial in treatment. In a word, DNA methylation modes are significantly related to tumor immunophenotypes and patients’ clinical responses to immunotherapy.

Discussion

DNA methylation is closely related to tumorigenesis and tumor progression. The extent of DNA methylation varies among cancer types and different stages of cancer progression. For example, the progression of prostate cancer has been related to DNA hypomethylation (Fraser et al., 2017; Wu et al., 2020), while bladder cancer pathology was characterized by global DNA hypermethylation (Osei-Amponsa et al., 2020). Thus, this observation revealed that DNA methylation may occur in a cancer-specific manner and alter the tumor microenvironment. Liu P et al. demonstrated that DNMT1 regulated the tumor growth in bladder cancer via modulating the status of DNA methylation in the promoter of PTEN (Liu et al., 2020). Zhu Y et al. demonstrated that MBD2 was a protective signature against bladder carcinoma according to the RNA data from the peripheral blood lymphocytes of 98 bladder cancer patients and 135 frequency-matched control patients (Zhu et al., 2004). Ying L et al. confirmed that epigenetic repression of RGS2 by UHRF1 contributes to bladder cancer progression (Ying et al., 2015). However, most researches only focused on the effect of a single DNA methylation regulator on the alteration of TME and tumor progression. As a result, the landscape of immune cell infiltration characteristics, which is mediated by the synergistic effect of multiple DNA methylation regulators, remained less understood. By clarifying the roles of diverse DNA methylation modes in immune cell infiltration, our knowledge about TME and anti-tumor response could advance, and foundations of more efficient immunotherapy strategies could be established.

In this study, we identified three DNA methylation modes based on expression level 15 DNA methylation regulators, and each DMRcluster was found to correlate with significantly different TMEs. Specifically, DMRclusters A, B, and C are characterized by immune-inflamed, immune-desert, and immune-excluded phenotypes, respectively. The immune-inflamed phenotype, or “hot tumor,” is characterized by the existence of a large number of immune cells in the TME (Zhang et al., 2020b; Gruber et al., 2020; Yu et al., 2021). The other two phenotypes, or “cold tumor,” show non-inflammatory infiltration. Despite the immune-excluded phenotype exhibits considerable immune cell infiltration, the immune cells are constrained by the stromal component that can be present either in the tumor capsule or throughout the whole tumor tissue to prevent the immune cells from exerting anti-tumor effects (Lambrechts et al., 2018; Kaymak et al., 2021). Such an idea is verified by the strong stromal activation in DMRcluster C, where the EMT pathway inhibited the activity of immune cells. Thus, our classifications of different DNA methylation modes were confirmed feasible and effective.

In addition, we confirmed that the transcriptomes in distinct DNA methylation modes are different, and obtained differentially expressed genes among thses DNA methylation patterns. Their actual compositions are related to DNA methylation and immune-related biological pathways. Therefore, we termed these differentially expressed genes as DMRGs. Three genomic subtypes were divided from the samples based on the expression of DMRGs, and these subtypes were also significantly related to distinct immunophenotypes. Therefore, DNA methylation is indeed irreplaceable in shaping the TME, and a systematic assessment of DNA methylation modes will contribute to understanding the mechanisms of tumorigenesis and to the advancements of personal medicine.

Since tumors are heterogeneous, we built a DMRscore model to evaluate DNA methylation features in individual tumors. The patients with high DMRscores were characterized as the immune-desert phenotype, while the patients with low DMRscores were characterized as the immune-inflamed phenotype. These results were further verified in the IMvigor210 cohort whose immunophenotypes have been identified (Necchi et al., 2017). Comprehensive analyses suggest that DMRscore signature is a robust and potential biomarker to assess patients’ response to immunotherapy, Patients with low DMRscore displayed higher expression of PD-L1 compared to patients with higher DMRscore, and had a better response to Atezolizumab. In addition, patients with low DMRscores exhibited less TP53 wild mutation, lower cancer grade, low tumor mutation burden, and molecular subtypes were mainly papillary subtypes.

In summary, DMRscore can systematically assess the DNA methylation landscape and detect the TME characteristics, thereby identifying the immunophenotypes of individual patients for more efficient immunotherapeutic strategies. Besides, DMRscore can be used to evaluate other features of bladder cancer patients, including molecular subtypes, genetic mutation, tumor stage, and clinical histology. Moreover, DMRscore could serve as an independent prognostic indicator for effective prediction of clinical outcomes, as well as a factor that reflects the efficacy of and clinical responses to immunotherapy. Our research uncovers that DNA methylation can alter the immune microenvironment, resulting in the emergence of a “cold tumor.” Herein, we propose a new hypothesis: targeting DNA methylation regulators or DMR-related biological pathways could be effective to alter the DNA methylation status so that the unfavorable factors could be removed, and “cold tumors” could transform into “hot” ones. If proved correct, this hypothesis may promote the development of immunotherapeutic agents and drug combinations. Our study provided a new perspective to reveal the global DNA methylation status in bladder cancer patients, to predict the immunophenotype of individual tumors, and to promote individualized medicine.

Compared with existing investigations on prognostic signatures of bladder cancer, this study has some noteworthy advantages and shortcomings. Firstly, our investigation contributed to demonstrate the effect of DNA methylation modification in shaping of tumor microenvironment complexity and diversity, and explored the potential role of DNA methylation status to predict the clinical response to Atezolizumab therapy in urothelial carcinoma. The global DNA methylation landscape was constructed as the observation object to systematically investigate the effect of DNA methylation modification on tumor microenvironment, which has not been clarified before this study. Our study is mainly based on bioinformatics analysis and requires further clinical verification. Basic experiments are needed to verify the relationship between prognostic characteristics and immune infiltration; In the future we will conduct multicenter, large sample size studies to prospectively validate the model in order to further test the predictive potential and clinical ability of our model.

Methods

Data Acquisition and Processing

The workflow in our study is displayed in Supplementary Figure S1A. 7 sets of transcriptome data and their corresponding clinical annotations were obtained from The Cancer Genome Alta (TCGA) and Gene-expression omnibus (GEO) databases, in which patients without complete clinical annotation were excluded. The “ComBat” algorithm in the R package “sva” was used to correct the batch effect of non-biological technical deviations. The comprehensive information of all alternative bladder cancer datasets is summarized in Table 1. The transcriptome data were downloaded from UCSC Xena database, and the somatic mutation information was obtained from TCGA database. We investigated numerous DNA methylation regulators, including the DNA methyltransferase family (DNMT1, DNMT3A, DNMT3B), the methyl-CpG-binding domain proteins (MeCP2, MBD1, MBD2, MBD3, MBD4), the ubiquitin-like proteins containing PHD and RING finger domains (UHRF1, UHRF2), zinc-finger domain proteins (ZBTB33, ZBTB4), NTHL1, SMUG1, and UNG (Jones, 2012; Moore et al., 2013; Koch et al., 2018). All the data were processed with the R package “Bioconductor” in R software (version 4.0.3).

TABLE 1

TABLE 1. The gene expression profiles of bladder cancer included in this study.

Unsupervised Clustering of the Fifteen DNA Methylation Regulators

The unsupervised clustering algorithm was utilized to find out the distinct DNA methylation patterns based on the expression of the fifteen alternative DNA methylation regulators. The R package “ConsensuClsterPlus” was run 1,000 repetitive times to ensure the stability of classification, and the clustering number was assigned according to the K value. Subsequently, to verify the differences in biological functions among the three DNA methylation patterns, we ran the R package “GSVA.” GSVA (gene set variation analysis) is an unsupervised method to evaluate variations of biological pathways in a sample population (Hänzelmann et al., 2013). The gene sets identified from GSVA were named “c2.cp.kegg.v6.2.-symbols.”

Estimation of Tumor Microenvironment Cell Infiltration

The ssGSEA (single-sample gene set enrichment analysis) was utilized to calculate the relative amounts of gene components in each TME cell infiltration. A gene set that labels each immune cell type was adopted from the published studies (Charoentong et al., 2017). We investigated several immune cell types, including activated B cells, activated CD4 T cells, macrophages, eosinophils, CD56dim natural killer cells, and neutrophils. Each type of immune cell was counted for its enrichment score, and a box diagram was used to compare the scores in different DNA methylation patterns.

Identification of DNA Methylation Related Genes Among Distinct DNA Methylation Modes

To reveal which genes are DMRGs, we classified the samples into three DNA methylation clusters based on the expression levels of fifteen DNA methylation regulators. The differentially expressed genes (DEGs) among the clusters were picked out by an R package named “limma” with the adjusted p-value <0.001.

Construction of DNA Methylation Regulator Score Groups

To qualify the DNA methylation status of individual samples from bladder cancer patients, we designed a DNA methylation signature termed DMRscore for assessments. The DNA methylation score groups were designed as follows: After identifying the DMRGs, we extracted the overlapping genes in them. Patients were divided into groups with distinct immune subtypes for further analysis, which was performed by an unsupervised clustering approach on the overlapping DMRGs. The number of gene clusters and model stability were validated by the consensus clustering algorithm. Furthermore, the univariate Cox regression analysis was executed for each overlapping DMRG to find out the ones related to prognosis. Subsequently, the established DMRscore model was subject to the principal component analysis (PCA) that combines the linear high-dimension indicators into their linear independent low-dimension counterparts. Moreover, in PCA, both types of indicators retain their original information, and the speed of data processing is accelerated. Both principal components (i.e., the two types of indicators) were extracted to calculate the DMRscore as DMRscore = ∑(PC1i + PC2i), where i is the expression level of the prognostic-related DMRGs.

To comprehensively evaluate the clinical outcome of each patient, a prognostic nomogram that contained the T stage, M stage, N stage, Gender, Age, clinical Stage and DMRscore was constructed. Subsequently, the 1-, 3-, 5- year overall survival probabilities were assessed by the calibration curve. A calibration curve close to 45° indicated the prominent prediction ability of the constructed model.

The Relationship Between DNA Methylation Features and Other Relevant Biological Functions

We obtained several gene sets that are involved in certain biological processes, e.g., DNA damage repair, homologous recombination, cell cycle, mismatch repair, DNA replication, nucleotide excision, carcinogenesis, Pan-F-TBRS, EMT, angiogenesis, immune checkpoint, actions of CD8 T effector cells, and antigen processing (Rosenberg et al., 2016; Şenbabaoğlu et al., 2016; Mariathasan et al., 2018). The correlations between DNA methylation features and these processes were further identified via correlation analysis.

The Genomic Profiles of Immune Checkpoint Blockage Effects and Corresponding Clinical Information

In order to explore the predictive effect of DNA methylation statuses in immunotherapy, we included an immunotherapeutic cohort in this study, advanced urothelial cancer treated with atezolizumab (IMvigor210 cohort) (Rosenberg et al., 2016). Atezolizumab is anti-PD-1 monoclonal antibody. The transcriptome profiles of immune checkpoint blockage effects and their corresponding clinical information were obtained from the public dataset.

Statistical Analysis

Data processing was conducted solely on R software (version 4.0.3). The R package named “limma” was run to analyze differential gene expressions among distinct subtypes. The Spearman analysis and distance correlation analysis were performed to calculate correlation coefficients between the DNA methylation regulators and the infiltration of immune cells. The survival curves of bladder cancer patients were plotted via the Kaplan-Meier method, and the curves’ area under the curve (AUC) was calculated to evaluate the specificity and sensitivity of DMRscores obtained by the R package “pROC”. The location and circle sequence of the DNA methylation regulators along the chromosomes were depicted by the R package “RCircos”. Moreover, the R package “DEseq2” was run to normalize the raw data and convert the normalized cell count to TPM in the “IMvigor 210” cohort. Finally, the mutation landscape was drawn and presented via the R package “maftools.” All statistic p numbers were bilateral, and p < 0.05 was considered statistically significant.

Data Availability Statement

All datasets generated for this study are included in the article material, including the TCGA. database (https://portal.gdc.cancer.gov/), and GEO dataset (https://www.ncbi.nlm.nih.gov/gds/).

Author Contributions

Study design and concept, FY. Supervision, HJ. Adviser, CX. All authors contributed towards data collection and analysis.

Funding

This study was supported by the National Natural Science Foundation of China (Grant Numbers: 81872102 and 81803900).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We sincerely appreciate all members who participated in the data collection and analysis.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.760369/full#supplementary-material

References

Barbie, D. A., Tamayo, P., Boehm, J. S., Kim, S. Y., Moody, S. E., Dunn, I. F., et al. (2009). Systematic RNA Interference Reveals that Oncogenic KRAS-Driven Cancers Require TBK1. Nature 462 (7269), 108–112. doi:10.1038/nature08460

PubMed Abstract | CrossRef Full Text | Google Scholar

Charoentong, P., Finotello, F., Angelova, M., Mayer, C., Efremova, M., Rieder, D., et al. (2017). Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cel Rep. 18 (1), 248–262. doi:10.1016/j.celrep.2016.12.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, X., Pan, X., Zhang, W., Guo, H., Cheng, S., He, Q., et al. (2020). Epigenetic Strategies Synergize with PD-L1/pd-1 Targeted Cancer Immunotherapies to Enhance Antitumor Responses. Acta pharmaceutica Sinica B 10 (5), 723–733. doi:10.1016/j.apsb.2019.09.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Craig, A. J., von Felden, J., Garcia-Lezana, T., Sarcognato, S., and Villanueva, A. (2020). Tumour Evolution in Hepatocellular Carcinoma. Nat. Rev. Gastroenterol. Hepatol. 17 (3), 139–152. doi:10.1038/s41575-019-0229-4

PubMed Abstract | CrossRef Full Text | Google Scholar

da Rocha, S. T., and Gendrel, A.-V. (2019). The Influence of DNA Methylation on Monoallelic Expression. Essays Biochem. 63 (6), 663–676. doi:10.1042/ebc20190034

PubMed Abstract | CrossRef Full Text | Google Scholar

Elashi, A. A., Sasidharan Nair, V., Taha, R. Z., Shaath, H., and Elkord, E. (2019). DNA Methylation of Immune Checkpoints in the Peripheral Blood of Breast and Colorectal Cancer Patients. Oncoimmunology 8 (2), e1542918. doi:10.1080/2162402x.2018.1542918

PubMed Abstract | CrossRef Full Text | Google Scholar

Fraser, M., Sabelnykova, V. Y., Yamaguchi, T. N., Heisler, L. E., Livingstone, J., Huang, V., et al. (2017). Genomic Hallmarks of Localized, Non-indolent Prostate Cancer. Nature 541 (7637), 359–364. doi:10.1038/nature20788

PubMed Abstract | CrossRef Full Text | Google Scholar

Gruber, T., Kremenovic, M., Sadozai, H., Rombini, N., Baeriswyl, L., Maibach, F., et al. (2020). IL-32γ Potentiates Tumor Immunity in Melanoma. JCI insight 5 (18). doi:10.1172/jci.insight.138772

PubMed Abstract | CrossRef Full Text | Google Scholar

Han, H. S., Jeong, S., Kim, H., Kim, H.-D., Kim, A. R., Kwon, M., et al. (2021). TOX-expressing Terminally Exhausted Tumor-Infiltrating CD8+ T Cells Are Reinvigorated by Co-blockade of PD-1 and TIGIT in Bladder Cancer. Cancer Lett. 499, 137–147. doi:10.1016/j.canlet.2020.11.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Hänzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: Gene Set Variation Analysis for Microarray and RNA-Seq Data. BMC Bioinformatics 14, 7. doi:10.1186/1471-2105-14-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Jones, P. A. (2012). Functions of DNA Methylation: Islands, Start Sites, Gene Bodies and beyond. Nat. Rev. Genet. 13 (7), 484–492. doi:10.1038/nrg3230

PubMed Abstract | CrossRef Full Text | Google Scholar

Kandimalla, R., van Tilborg, A. A., and Zwarthoff, E. C. (2013). DNA Methylation-Based Biomarkers in Bladder Cancer. Nat. Rev. Urol. 10 (6), 327–335. doi:10.1038/nrurol.2013.89

PubMed Abstract | CrossRef Full Text | Google Scholar

Kaymak, I., Williams, K. S., Cantor, J. R., and Jones, R. G. (2021). Immunometabolic Interplay in the Tumor Microenvironment. Cancer Cell 39 (1), 28–37. doi:10.1016/j.ccell.2020.09.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Koch, A., Joosten, S. C., Feng, Z., de Ruijter, T. C., Draht, M. X., Melotte, V., et al. (2018). Analysis of DNA Methylation in Cancer: Location Revisited. Nat. Rev. Clin. Oncol. 15 (7), 459–466. doi:10.1038/s41571-018-0004-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Koelsche, C., Schrimpf, D., Stichel, D., Sill, M., Sahm, F., Reuss, D. E., et al. (2021). Sarcoma Classification by DNA Methylation Profiling. Nat. Commun. 12 (1), 498. doi:10.1038/s41467-020-20603-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Lambrechts, D., Wauters, E., Boeckx, B., Aibar, S., Nittner, D., Burton, O., et al. (2018). Phenotype Molding of Stromal Cells in the Lung Tumor Microenvironment. Nat. Med. 24 (8), 1277–1289. doi:10.1038/s41591-018-0096-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, C., Liu, R., Wang, B., Lian, J., Yao, Y., Sun, H., et al. (2021). Blocking IL-17A Enhances Tumor Response to Anti-PD-1 Immunotherapy in Microsatellite Stable Colorectal Cancer. J. Immunother. Cancer 9 (1). doi:10.1136/jitc-2020-001895

CrossRef Full Text | Google Scholar

Liu, P., Wu, L., Chand, H., Li, C., Hu, X., and Li, Y. (2020). Silencing of miR-152 Contributes to DNMT1-Mediated CpG Methylation of the PTEN Promoter in Bladder Cancer. Life Sci. 261, 118311. doi:10.1016/j.lfs.2020.118311

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Z., Sun, T., Zhang, Z., Bi, J., and Kong, C. (2021). An 18-gene Signature Based on Glucose Metabolism and DNA Methylation Improves Prognostic Prediction for Urinary Bladder Cancer. Genomics 113 (1 Pt 2), 896–907. doi:10.1016/j.ygeno.2020.10.022

PubMed Abstract | CrossRef Full Text | Google Scholar

MacGregor, H. L., Sayad, A., Elia, A., Wang, B. X., Katz, S. R., Shaw, P. A., et al. (2019). High Expression of B7-H3 on Stromal Cells Defines Tumor and Stromal Compartments in Epithelial Ovarian Cancer and Is Associated with Limited Immune Activation. J. Immunotherapy Cancer 7 (1), 357. doi:10.1186/s40425-019-0816-5

CrossRef Full Text | Google Scholar

Mariathasan, S., Turley, S. J., Nickles, D., Castiglioni, A., Yuen, K., Wang, Y., et al. (2018). TGFβ Attenuates Tumour Response to PD-L1 Blockade by Contributing to Exclusion of T Cells. Nature 554 (7693), 544–548. doi:10.1038/nature25501

PubMed Abstract | CrossRef Full Text | Google Scholar

Moore, L. D., Le, T., and Fan, G. (2013). DNA Methylation and its Basic Function. Neuropsychopharmacol 38 (1), 23–38. doi:10.1038/npp.2012.112

PubMed Abstract | CrossRef Full Text | Google Scholar

Necchi, A., Joseph, R. W., Loriot, Y., Hoffman-Censits, J., Perez-Gracia, J. L., Petrylak, D. P., et al. (2017). Atezolizumab in Platinum-Treated Locally Advanced or Metastatic Urothelial Carcinoma: post-progression Outcomes from the Phase II IMvigor210 Study. Ann. Oncol. 28 (12), 3044–3050. doi:10.1093/annonc/mdx518

PubMed Abstract | CrossRef Full Text | Google Scholar

Ortega-Recalde, O., and Hore, T. A. (2019). DNA Methylation in the Vertebrate Germline: Balancing Memory and Erasure. Essays Biochem. 63 (6), 649–661. doi:10.1042/ebc20190038

PubMed Abstract | CrossRef Full Text | Google Scholar

Osei-Amponsa, V., Buckwalter, J. M., Shuman, L., Zheng, Z., Yamashita, H., Walter, V., et al. (2020). Hypermethylation of FOXA1 and Allelic Loss of PTEN Drive Squamous Differentiation and Promote Heterogeneity in Bladder Cancer. Oncogene 39 (6), 1302–1317. doi:10.1038/s41388-019-1063-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Qiu, Z., Li, H., Zhang, Z., Zhu, Z., He, S., Wang, X., et al. (2019). A Pharmacogenomic Landscape in Human Liver Cancers. Cancer Cell 36 (2), 179–193. e11. doi:10.1016/j.ccell.2019.07.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). Limma powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies. Nucleic Acids Res. 43 (7), e47. doi:10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

Robertson, A. G., Kim, J., Al-Ahmadie, H., Bellmunt, J., Guo, G., Cherniack, A. D., et al. (2017). Comprehensive Molecular Characterization of Muscle-Invasive Bladder Cancer. Cell 171 (3), 540–e25. doi:10.1016/j.cell.2017.09.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosenberg, J. E., Hoffman-Censits, J., Powles, T., van der Heijden, M. S., Balar, A. V., Necchi, A., et al. (2016). Atezolizumab in Patients with Locally Advanced and Metastatic Urothelial Carcinoma Who Have Progressed Following Treatment with Platinum-Based Chemotherapy: a Single-Arm, Multicentre, Phase 2 Trial. The Lancet 387 (10031), 1909–1920. doi:10.1016/s0140-6736(16)00561-4

CrossRef Full Text | Google Scholar

Samstein, R. M., Lee, C.-H., Shoushtari, A. N., Hellmann, M. D., Shen, R., Janjigian, Y. Y., et al. (2019). Tumor Mutational Load Predicts Survival after Immunotherapy across Multiple Cancer Types. Nat. Genet. 51 (2), 202–206. doi:10.1038/s41588-018-0312-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Sasidharan Nair, V., Toor, S. M., Taha, R. Z., Shaath, H., and Elkord, E. (2018). DNA Methylation and Repressive Histones in the Promoters of PD-1, CTLA-4, TIM-3, LAG-3, TIGIT, PD-L1, and Galectin-9 Genes in Human Colorectal Cancer. Clin. Epigenet 10 (1), 104. doi:10.1186/s13148-018-0539-3

CrossRef Full Text | Google Scholar

Şenbabaoğlu, Y., Gejman, R. S., Winer, A. G., Liu, M., Van Allen, E. M., de Velasco, G., et al. (2016). Tumor Immune Microenvironment Characterization in clear Cell Renal Cell Carcinoma Identifies Prognostic and Immunotherapeutically Relevant Messenger RNA Signatures. Genome Biol. 17 (1), 231. doi:10.1186/s13059-016-1092-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Siegel, R. L., Miller, K. D., Fuchs, H. E., and Jemal, A. (2021). Cancer Statistics, 2021. CA A. Cancer J. Clin. 71 (1), 7–33. doi:10.3322/caac.21654

CrossRef Full Text | Google Scholar

Sotiriou, C., Wirapati, P., Loi, S., Harris, A., Fox, S., Smeds, J., et al. (2006). Gene Expression Profiling in Breast Cancer: Understanding the Molecular Basis of Histologic Grade to Improve Prognosis. J. Natl. Cancer Inst. 98 (4), 262–272. doi:10.1093/jnci/djj052

CrossRef Full Text | Google Scholar

Thuijs, N. B., Berkhof, J., Özer, M., Duin, S., Splunter, A. P., Snoek, B. C., et al. (2021). DNA Methylation Markers for Cancer Risk Prediction of Vulvar Intraepithelial Neoplasia. Int. J. Cancer 148 (10), 2481–2488. doi:10.1002/ijc.33459

CrossRef Full Text | Google Scholar

Wan, C., Keany, M. P., Dong, H., Al-Alem, L. F., Pandya, U. M., Lazo, S., et al. (2021). Enhanced Efficacy of Simultaneous PD-1 and PD-L1 Immune Checkpoint Blockade in High-Grade Serous Ovarian Cancer. Cancer Res. 81 (1), 158–173. doi:10.1158/0008-5472.CAN-20-1674

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilhelm, C. S., Kelsey, K. T., Butler, R., Plaza, S., Gagne, L., Zens, M. S., et al. (2010). Implications of LINE1 Methylation for Bladder Cancer Risk in Women. Clin. Cancer Res. 16 (5), 1682–1689. doi:10.1158/1078-0432.ccr-09-2983

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, A., Cremaschi, P., Wetterskog, D., Conteduca, V., Franceschini, G. M., Kleftogiannis, D., et al. (2020). Genome-wide Plasma DNA Methylation Features of Metastatic Prostate Cancer. J. Clin. Invest. 130 (4), 1991–2000. doi:10.1172/jci130887

CrossRef Full Text | Google Scholar

Yi, M., Jiao, D., Xu, H., Liu, Q., Zhao, W., Han, X., et al. (2018). Biomarkers for Predicting Efficacy of PD-1/pd-L1 Inhibitors. Mol. Cancer 17 (1), 129. doi:10.1186/s12943-018-0864-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Ying, L., Lin, J., Qiu, F., Cao, M., Chen, H., Liu, Z., et al. (2015). Epigenetic Repression of Regulator of G-Protein Signaling 2 by Ubiquitin-like with PHD and Ring-finger Domain 1 Promotes Bladder Cancer Progression. Febs J. 282 (1), 174–182. doi:10.1111/febs.13116

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, X., Liu, W., Chen, S., Cheng, X., Paez, P. A., Sun, T., et al. (2021). Immunologically Programming the Tumor Microenvironment Induces the Pattern Recognition Receptor NLRC4-dependent Antitumor Immunity. J. Immunother. Cancer 9 (1). doi:10.1136/jitc-2020-001595

CrossRef Full Text | Google Scholar

Zeng, D., Li, M., Zhou, R., Zhang, J., Sun, H., Shi, M., et al. (2019). Tumor Microenvironment Characterization in Gastric Cancer Identifies Prognostic and Immunotherapeutically Relevant Gene Signatures. Cancer Immunol. Res. 7 (5), 737–750. doi:10.1158/2326-6066.cir-18-0436

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, B., Wu, Q., Li, B., Wang, D., Wang, L., and Zhou, Y. L. (2020). m6A Regulator-Mediated Methylation Modification Patterns and Tumor Microenvironment Infiltration Characterization in Gastric cancerA Regulator-Mediated Methylation Modification Patterns and Tumor Microenvironment Infiltration Characterization in Gastric Cancer. Mol. Cancer 19 (1), 53. doi:10.1186/s12943-020-01170-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Yang, M., Ng, D. M., Haleem, M., Yi, T., Hu, S., et al. (2020). Multi-omics Data Analyses Construct TME and Identify the Immune-Related Prognosis Signatures in Human LUAD. Mol. Ther. - Nucleic Acids 21, 860–873. doi:10.1016/j.omtn.2020.07.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, Y., Spitz, M. R., Zhang, H., Grossman, H. B., Frazier, M. L., and Wu, X. (2004). Methyl-CpG-binding Domain 2. Cancer 100 (9), 1853–1858. doi:10.1002/cncr.20199

PubMed Abstract | CrossRef Full Text | Google Scholar