DNA methylation regulator-mediated modification pattern defines tumor microenvironment immune infiltration landscape in colon cancer

Abstract

Emerging evidence implies a non-negligible role of DNA methylation in tumor immunity, however, its comprehensive impact on tumor microenvironment (TME) formation and immune activation remains unclear. In this study, we integrated 24 DNA methylation regulators among 754 colon cancer patients to distinguish different modification patterns via an unsupervised clustering method, and explore their TME immune characteristics. Three DNA methylation modification patterns with distinct prognosis and biological behaviors were identified, consistent with three known phenotypes of immune-inflamed, immune-excluded, and immune-desert. We then determined a DNA methylation gene signature and constructed a DNA methylation score (DMS) to quantify modification patterns individually through principal component analysis algorithms. DMS-low group had characteristics of specific molecular subtypes, including microsatellite instability, CpG island methylator phenotype positive, and mutant BRAF, presented by increased mutation burden, activation of DNA damage repair and immune-related pathways, highly TME immune cells infiltration, and hence, a preferable prognosis. Further, low DMS was also demonstrated to be correlated to better response and prolonged survival of anti-PD-L1 antibody, indicating that DMS could be considered as an effective predictive tool for immunotherapy. In conclusion, our work presented a landscape of different DNA methylation modification patterns, and their vital role in the formation of TME diversity and complexity, which could help to enhance understanding of TME immune infiltration characteristics and more importantly, guide immunotherapy strategies more effectively and personalized.

Introduction

Colon cancer is common worldwide and remains one of the leading causes of cancer-related mortality (Sung et al., 2021). As a biologically heterogeneous disease, colon cancer derives from the accumulation of a series of genetic and epigenetic changes that transform normal glandular epithelium into malignant invasive adenocarcinoma. Now, it is appreciated that there are multiple molecular pathways involved in these genetic mutations and epigenetic alterations during colon cancer development, including microsatellite instability (MSI), the CpG island methylator phenotype (CIMP), chromosome instability (CIN), and somatic mutations of critical oncogenes like BRAF and KRAS (Nguyen et al., 2020).

Epigenetics, referring to heritable alterations in gene expression that are not dependent on changes in the DNA sequence, plays an important role in the pathogenesis of colon cancer (Lao and Grady. 2011). Therein, aberrant DNA methylation, one of the most widely studied epigenetic modifications, could lead to the dysregulation of gene expression in colon cancer. DNA methylation is a reversible modification process mediated by DNA methyltransferases (DNMTs) that facilitate the catalytic addition of methyl groups to the fifth position of the cytosine of CpG dinucleotides to generate 5-methylcytosine (5mC). On the contrary, 5mC could be reversed to unmodified cytosine through TET dioxygenase-mediated oxidation (Bestor 2000). Over the past two decades, emerging studies have elucidated the epigenetically regulatory mechanism of DNA methylation in colon cancer-specific gene expression patterns. Hypermethylation in the promoter region could silence the expression of tumor-suppressor genes, and contribute to loss of function. For instance, MSI, one of the hallmarks of molecular subtypes of colon cancer (Dienstmann et al., 2017), is the consequence of a deficiency in the DNA mismatch repair (MMR) system, which is not only due to the genetic mutation of MMR-related genes, but also the results of hypermethylation of MLH1 gene (Herman et al., 1998). On the contrary, global hypomethylation of the genome has been demonstrated to influence colon cancer development through inducing CIN and global loss of imprinting (Suter et al., 2004). These insights have improved our understanding of colon cancer pathophysiology and provide clues to discover novel biomarkers and therapeutic targets.

Recently, immunotherapy, especially the inhibitor targeting immune checkpoints like CTLA-4, PD-1, or PD-L1, has achieved durable anti-tumor activity in a range of cancer types. However, there are many patients, particularly in colon cancer with microsatellite stable (MSS), do not benefit from this advanced treatment (Le et al., 2015). The major reason is thought to lack lymphocytes infiltration in the MSS tumor microenvironment (TME), forming a “immune-desert” phenotype and resulting a weak immunoreactivity to the immunotherapy (Chen and Mellman, 2017). Increasing evidence demonstrated DNA methylation regulators mediated regulation is critical in anti-tumor immune response through involving in many processes of the cancer-immunity cycle (Cao and Yan, 2020). To be specific, DNA methylation-associated mutagenesis could generate tumor neoantigens (Alexandrov et al., 2013). In tumor cells, DNMTs could suppress the expression of MHC-I to dysregulate the antigen-presenting machinery, and suppress the expression of pro-inflammatory chemokines, such as CXCL9 and CXCL10, which are required by effector T-cells to permeate the TME and execute an immune attack (Chen et al., 2017; Luo et al., 2018). Hypermethylation of the PD-L1 promotor region inhibits its expression and leads to an inferior prognosis in various cancer types (Goltz et al., 2017; Heiland et al., 2017). Meanwhile, in immune cells, MBD2, the “reader” of methylated DNA, is necessary to induce dendritic cells phenotypic activation and then initiate the T cell response (Cook et al., 2015). By contrast, TET2 and HDAC2, the “eraser” of DNA methylation and histone acetylation respectively, coordinate to suppress IL-6 expression of dendritic cells, inhibiting the inflammatory response (Zhang et al., 2015).

However, to date, the majority of studies focus on the function of one or two DNA methylation regulators, which cannot reflect the whole landscape of DNA methylation in the formation of tumor-permissive immune environment. Therefore, comprehensive recognition of the TME immune characteristics mediated by multiple DNA methylation regulators, including TME infiltrating immune cells and activity of immune/inflammatory-related pathways, could enhance our understanding of TME immune regulation, and further provide novel perspectives for cancer immunotherapy. In this study, we integrated the transcriptomic and clinical information of 754 colon cancer samples to identify DNA methylation modification patterns with distinct TME immune characteristics, which were highly consistent with three known immune phenotypes, including immune-inflamed, immune-excluded, and immune-desert phenotype, respectively. In addition, we determined the DNA methylation gene signature and constructed a scoring system to quantify modification patterns for individual patients, which could be served as an effective biomarker for predicting the efficacy and prognosis of immunotherapy.

Materials and methods

Study design, colon cancer datasets collection, and data processing

Supplementary Figure S1A depicted the workflow of our present study. Generally, we searched publicly available transcriptomic data and clinical annotation of colon cancer samples from the Gene-Expression Omnibus (GEO) and the Cancer Genome Atlas (TCGA) database. Seven eligible colon cancer datasets with comprehensive survival information [GSE39582 (n = 562), GSE38832 (n = 122), GSE39084 (n = 70), GSE72970 (n = 124), GSE103479 (n = 155), GSE87211 (n = 196), and TCGA-colon adenocarcinoma cohort (COAD, n = 430)], including 1,659 patients, were collected for our further analysis (Supplementary Table S1). For the RNA sequencing data of the TCGA-COAD cohort, the fragments per kilobase of transcript per million mapped reads (FPKM) value of each sample was downloaded from the Genomic Data Commons (GDC, https://portal.gdc.cancer.gov/) by the “TCGAbiolinks” R package directly (Colaprico et al., 2016). And for the microarray data of the GSE cohort, we used the normalized matrix files downloaded from GEO. The somatic mutation data, copy number variation (CNV) data, and TCGA pan-cancer RNA sequencing data were acquired from the University of California Santa Cruz (UCSC) Xena browser (https://xenabrowser.net).

Moreover, we also included an immunotherapeutic cohort (IMvigor210, n = 348), urothelial carcinoma treated with anti-PD-L1 antibody atezolizumab, to evaluate the effect of DNA methylation modification in immunotherapy. The expression data and clinical information were available from the “IMvigor210” R package (Mariathasan et al., 2018).

Unsupervised clustering for 24 DNA methylation regulators

After a systematic review of published articles, we identified a total of 24 DNA methylation regulators, including 3 writers (DNMT1, DNMT3A, DNMT3B), 3 erasers (TET1, TET2, TET3), and 18 readers (MBD1, MBD2, MBD3, MBD4, ZBTB33, ZBTB38, ZBTB4, ZBTB24, UHRF1, UHRF2, MECP2, UNG, TDG, NTHL1, SMUG1, NSUN2, MGMT, DMAP1). The protein-protein interactions (PPI) network among 24 regulators were analyzed by the STRING interaction database (https://string-db.org/) (Szklarczyk et al., 2019) and visualized by the Cytoscape software (Shannon et al., 2003).

Three GEO datasets (GSE39582, GSE38832, and GSE39084) with the same microarray platform and no prognostic differences were integrated as meta-cohort (n = 754) to identify different DNA methylation modification patterns mediated by 24 regulators (Supplementary Figure S2A). The non-biological technical biases caused batch effect among each cohort was eliminated by the “ComBat” algorithm of the “sva” R package. Based on the expression of 24 regulators, the Nonnegative Matrix Factorization (NMF) method was used to determine different DNA methylation modification patterns through the “NMF” R package (Gaujoux and Seoighe. 2010).

Identification of differentially expressed genes and generation of DNA methylation gene signature

First, to identify DNA methylation phenotype-related genes, we used the empirical Bayesian approach of the “limma” R package to determine differentially expressed genes (DEGs) among three DNA methylation modification patterns (Ritchie et al., 2015). The adjusted p-value was set as < 0.01 to select significant DEGs. Then, we performed a univariate Cox regression analysis to confirm the DEGs that were significantly related to the prognosis (p-value < 0.01). Subsequently, based on the expression of these prognostic DEGs, we conducted the second NMF clustering algorithm to obtain DNA methylation gene clusters as well as validate their stability. Furthermore, through the principal component analysis (PCA) method, we used these prognostic DEGs to construct the DNA methylation gene signature, termed DNA methylation score (DMS), which could quantify the DNA methylation modification pattern for each patient. The procedure of establishing the DMS was similar to a previous study, and we added the principal component 1 and 2 to acquire the signature scores (Zhang et al., 2020; Gao et al., 2021). The formula of DMS was expressed as follow:where i is the expression of DNA methylation signature genes.

Gene set variation analysis, gene set enrichment analysis, and functional annotation

To investigate the difference in the biological processes among different DNA methylation modification patterns, gene clusters, and DMS groups, we performed gene set variation analysis (GSVA) and gene set enrichment analysis (GSEA) analyses through the “GSVA” and “clusterProfiler” R packages, respectively (Yu et al., 2012; Hanzelmann et al., 2013). The hallmark gene sets (h.all.v7.5.1.symbols) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) gene sets (c2.cp.kegg.v7.5.1.symbols) were downloaded from the MSigDB database for running enrichment analysis (http://www.gsea-msigdb.org/gsea/msigdb/). All enrichment p-values were adjusted by the Benjamini-Hochberg methods and less than 0.05 were considered statistically significant (Thissen et al., 2002). Moreover, we performed Gene Ontology (GO) functional annotation for DEGs via the “clusterProfiler” R package, with the cutoff value of false discovery rate (FDR) < 0.05.

We additionally collected 18 classical biological processes constructed by Mariathasan et al. (2018), including: 1) CD8 T effector; 2) DNA damage repair (DDR); 3) antigen-processing machinery (APM); 4) immune checkpoint; 5) cell cycle regulators; 6) Fanconi anemia; 7) pan-fibroblast TGFβ response signature (Pan-F-TBRS); 8) epithelial-mesenchymal transition (EMT) markers including EMT1, EMT2 and EMT3; 9) WNT targets; 10) fibroblast growth factor receptor 3 (FGFR3) related signature; 11) cell cycle signature; 12) mismatch repair 13) homologous recombination; 14) nucleotide excision repair; 15) DNA replication; 16) base excision repair. The detailed gene set of the corresponding pathway was provided in Supplementary Table S2.

Estimation of the tumor microenvironment infiltrating cell abundance and the immune infiltration score

We used three different algorithms, including the single-sample gene-set enrichment analysis (ssGSEA) algorithm of the “GSVA” R package, the CIBERSORT method, and the Tumor Immune Estimation Resource (TIMER) database, to evaluate the infiltrating abundance of various TME immune cells, such as B cell, CD8+ T cell, dendritic cell, and macrophage et al. The gene sets of each type TME infiltrating cell were extracted from the study of Charoentong and listed in Supplementary Table S3 (Charoentong et al., 2017).

Besides, through applying “xCell” and “ESTIMATE” methods, we calculated the TME stromal score, immune score, estimate score, and microenvironment score based on signature gene expression to infer the fraction of stromal and immune cells in colon cancer samples. The calculation was performed by the “xCell” and “estimate” R packages, respectively (Yoshihara et al., 2013; Aran et al., 2017).

Cell culture

The colorectal cancer cell line HCT116 and normal colonic epithelial cell line NCM460 were purchased from Type Culture Collection of the Chinese Academy of Science (Shanghai, China). HCT116 cells were cultured in McCoy’s 5A Medium (16600-082, Gibco) with 10% fetal bovine serum (10099141C, Gibco) in 5% CO₂ at 37°C. NCM460 cells were cultured in RPMI 1640 Medium (118575-093, Invitrogen) with same condition.

RNA extraction and real-time quantitative reverse transcription-polymerase chain reaction)

Total RNA was isolated using TRIzol™ Reagent (15596026, Invitrogen) and quantified with a NanoDrop 2000™ (Thermo Fisher Scientific, United States). 1 μg RNA was used for the reverse transcription reaction to generate cDNA through the PrimeScript™ RT Reagent Kit with gDNA Eraser (RR047A, TaKaRa) according to the manufacturer’s protocols. The mRNA expression was determined by rt-qPCR, which was performed using Ultra SYBR Mixture (CW0957M, CWBIO) and a LightCycler^® 480 II system (Roche, Shanghai, China). The mRNA expression of ACTB was used as a reference. The primers used in this study were listed in Supplementary Table S4.

Statistical analysis

The normality of data was tested by the Shapiro-Wilk test. For the comparison of two groups, we used the t-test to detect the significant difference between normally distributed data, and the Wilcoxon test for skewed distributed data. For the comparison of three or more groups, one-way ANOVA and Kruskal-Wallis tests were conducted to detect the significant difference between normal distributed and skewed distributed data, respectively. A chi-squared test was used to compare the frequency differences between the two groups. Correlation coefficients were calculated by the Spearman and distance correlation analyses.

For the survival analysis, we focused on the overall survival (OS) and recurrence-free survival (RFS), and we obtained the best cut-off value through the “survcutpoint” function of the “survminer” R package. Kaplan-Meier method was used to depict the survival curves, and the log-rank test was utilized to identify significant survival differences between groups. Univariate Cox proportional hazards regression model was applied to calculate the hazard ratio (HR) and 95% confidence interval (95% CI) for each DNA methylation regulator and DNA methylation related gene. And the multivariate Cox model was performed to determine the independent prognostic factors when adjusted by clinical characteristics. The results of Cox regression analyses were visualized by the “forestplot” R package. The prediction performance of DMS to evaluate the OS probability at distinct times was assessed by the receiver operating characteristic (ROC) curves and quantified by the area under the curve (AUC), which were conducted via the “timeROC” R package.

All statistical analyses were accomplished in R 3.6.1 software, and all reported p-values were two-sided, with p-value < 0.05 as statistically significant.

Results

Multi-omics landscape analysis of DNA methylation regulators in colon cancer

In this study, we identified 24 DNA methylation regulators, including 3 writers, 3 erasers, and 18 readers. Figure 1A summarized the dynamic reversible DNA methylation modification process of these regulators. Firstly, we analyzed the incidence of somatic mutation and CNV of 24 regulators in colon cancer from the genomic perspective. We visualized the mutation landscape of 24 regulators through the “waterfall” function of the “maftools” R package. Among 399 available samples in the TCGA-COAD cohort, a total of 113 (28.3%) mutations occurred in 24 regulators (Figure 1B). Therein, 2 erasers, TET3 and TET1, had the highest mutation frequency (6%) with missense mutation as a major mutation type, while MBD3 did not present any mutations. Besides, we found significant co-occurrence mutation relationships between several regulators, such as TET3 and ZBTB38, DNMT1 and MBD1, along with UHRF2 and MECP2 (Supplementary Figure S1B). In addition, the CNV analysis revealed an amplification of copy number in 6 regulators and deletion in other 6 regulators. The location of CNV alterations of 12 regulators on the chromosomes was shown in Figure 1C, which was accomplished by the “RCircos” R package.

FIGURE 1

Secondly, to ascertain whether the above genomic variations affected the transcription of 24 regulators, we further compared the mRNA expression level of these regulators between colon cancer and normal samples from the transcriptomic perspective. We observed most of the regulators were significantly differentially expressed between tumor and normal tissues. Interestingly, compared to normal tissues, we found that many regulators with amplified CNV had a markedly higher expression level in colon cancer tissues (e.g., DNMT1 and DNMT3B), and vice versa (e.g., MBD2 and ZBTB4), suggesting that the alteration of CNV might be the prominent factor resulting in the abnormal expression of DNA methylation regulators (Figures 1C,D). To further confirm the expression level of these regulators in colorectal cancer cell line and normal colonic epithelial cell line, we selected three most significantly differentially expressed regulators between tumor and normal tissues, DNMT3B (writer), TET2 (eraser), and UHRF1 (reader), to compare their mRNA expression in HCT116 and NCM460 cells. The results of quantitative reverse transcription-polymerase chain reaction (qRT-PCR) showed that the mRNA level of DNMT3B and UHRF1 were significantly higher in HCT116 cells while the mRNA level of TET2 was significantly higher in NCM460 cells (Supplementary Figures S1C–E).

Thirdly, we explored the interaction relationship between 24 regulators from the proteomics perspective. The PPI network depicted the extensive protein interactions among these regulators (Figure 1E). Comprehensively, the above analyses presented a highly genetic heterogeneity and expressional anomalism of DNA methylation regulators in colon cancer from a multi-omics landscape perspective, hinting at a critical role of these regulators in tumorigenesis of colon cancer.

Prognosis and immune characteristics of 24 DNA methylation regulators in colon cancer

To clarify the role of 24 DNA methylation regulators in colon cancer clinical prognosis and TME cell infiltration characterization, we gathered three GEO datasets (GSE39582, GSE38832, and GSE39084) without prognostic differences as meta-cohort for further analyses (Supplementary Figure S2A). The Cox regression analysis identified high expression of nearly half of the regulators were associated with a favorable prognosis (Supplementary Figure S2B). Moreover, the interaction network of regulators visualized the correlation of regulators’ expression and their prognostic significance (Figure 2A; Supplementary Table S5). We found that the expression of most regulators was positively related to each other, not only in regulators with the same functional category, but also among writers, erasers, and readers.

FIGURE 2

TME infiltrating immune cells had been widely reported to display an epigenomic reprogramming, especially in aberrant DNA methylation (Loo Yau et al., 2019). Therefore, we investigated the correlation between 24 regulators and TME infiltrating cells. The expression of MBD1, MBD2, and ZBTB4 was positively correlated with infiltrating abundance of majority immune cells, which could explain their favorable prognostic value to some extent (Figure 2B). A previous study demonstrated that colon cancer patients with a MSI status presented an active intra-tumoral immune environment, and hence, had a robust response to immunotherapy and superior prognosis (Andre et al., 2020). Similarly, we found the expression of MBD1 and MBD2 was upregulated in MSI colon cancer patients, implying both regulators might involve in immune activation (Supplementary Figure S2C). In consideration of the MBD2 expression were highly correlated to immune cells quantity, especially for activated CD4+ and CD8+ T-cells, as well as its preferable prognostic value, we next thoroughly analyzed the role of MBD2 in tumor immune environment formation. First, the MBD2 high expression group exhibited a higher infiltrating abundance of most immune cells and TME immune score (Figures 2C,D). Second, GSEA analysis revealed several activated immune and inflammation-related pathways, such as the TNFα, interferon- and interleukin-mediated signaling pathways were significantly enriched in the MBD2 high expression group (Figure 2E). Third, owing to the vital role of MBD2 in tumor immunity, we additionally investigated whether its expression could predict the efficacy and prognosis of immunotherapy. In the IMvigor210 cohort, we observed a marked survival benefit in patients with MBD2 high expression, although the therapeutic response was similar between two groups (Figure 2F; Supplementary Figure S2D).

Collectively, above results exhibited crosstalk among 24 DNA methylation regulators and their significant impact on tumor immunity. The MBD2 expression was positively correlated to TME immune cells infiltration and might be a potential prognostic biomarker in immunotherapy.

DNA methylation modification patterns mediated by 24 regulators in colon cancer

As the markedly different clinical outcomes and TME infiltrating cells characterizations among 24 DNA methylation regulators in colon cancer, we speculated that these regulators may contribute to forming different DNA methylation modification patterns in individual colon cancer. Based on the expression of 24 regulators, we applied the NMF unsupervised clustering to obtain three distinct clusters, including DNA methylation modification pattern-A (n = 306), pattern-B (n = 227), and pattern-C (n = 221), in meta-cohort (Supplementary Figures S3A,B). Accordingly, we found different regulators were differentially expressed in different patterns. Pattern-A was characterized by high expression of TET1, TET2, ZBTB4, ZBTB33, ZBTB38, MBD4, UHRF2, and MECP2; pattern-B was characterized by high expression of DNMT1, DNMT3B, ZBTB24, UHRF1, and UNG; and pattern-C was characterized by high expression of DMAP1, SMUG1 and MGMT (Figure 3A; Supplementary Figure S3C). In addition, the graph of principal component analysis (PCA) showed that three patterns were distinctly segregated, indicating a complete distinguishment of three patterns (Figure 3B). Moreover, the survival analysis revealed a notably favorable OS and RFS in pattern-B (Figures 3C,D).

FIGURE 3

Then, we performed a GSEA analysis to explore the different hallmark signaling pathways among three patterns (Supplementary Table S6). We found pattern-A presented the activated enrichment pathways associated with stromal and carcinogenic activation, such as the KRAS, TGFβ, EMT, and angiogenesis pathways, as well as the highest stromal score (Figure 3E; Supplementary Figure S3D). Immune regulation and tumor suppressor-related pathways were significantly enriched in pattern-B, including interferon-α, and -γ, TNFα, P53 signaling, G2M checkpoint, and DNA damage repair related pathways (Figure 3F). In contrast, these pathways were all inhibited in pattern-C (Figure 3G). We subsequently compared the abundance of TME infiltrating cells among three patterns. Not surprisingly, high infiltration of activated CD4+ and CD8+ T-cells were observed in pattern-B (Supplementary Figure S3E; Supplementary Table S7), which was consistent with the results of enrichment into activated immune-related pathways, and may explain the preferable prognosis (Figures 3C,D). Meanwhile, several innate immune cells were remarkably enriched in pattern-A, including eosinophil, mast cell, and plasmacytoid dendritic cell, while the content of the majority of immune cells was relatively lower in pattern-C (Supplementary Figure S3E).

From the above analyses, we were surprised to find three DNA methylation modification patterns exhibited significantly distinct biological behaviors and TME infiltrating cell characterizations, in accordance with the conception of three tumor immune phenotypes (Hegde and Chen, 2020). Pattern-A was classified as immune-excluded phenotype, characterized by innate immune cell infiltration and stromal activation; pattern-B was classified as immune-inflamed phenotype, characterized by a high population of adaptive immune cell infiltration and immune activation; pattern-C was classified as immune-desert phenotype, characterized by few infiltrations of immune cells and immune suppression.

Generation of DNA methylation gene signature

To further investigate the heterogeneity of three DNA methylation modification patterns, we identified 249 DEGs among three patterns (Supplementary Figures S4A,B). The GO functional annotation showed that these DEGs were significantly enriched in pathways related to DNA modification and damage repair events (Figure 4A; Supplementary Table S8), indicating the differences in clinical outcomes and TME characterizations among three patterns might result from these DEGs. Subsequently, we obtained 152 genes associated with prognosis through univariate Cox regression analysis (Figure 4B; Supplementary Table S9). Furthermore, based on the expression of 152 genes, we conducted unsupervised clustering to classify patients into different genomic subtypes. Surprisingly, consistent with three modification patterns, the NMF algorithm still clustered three distinct DNA methylation genomic phenotypes, named DNA methylation gene cluster-A (n = 272), cluster-B (n = 183), and cluster-C (n = 299), respectively (Figure 4C; Supplementary Figures S4C,D). The plot of PCA also presented a distinct separation of three gene clusters (Supplementary Figure S4E). As well, the prominent differences in the expression of 24 DNA methylation regulators were observed among three gene clusters, in line with the results of three modification patterns (Supplementary Figure S4F). In accordance with the feature of modification pattern-B, gene cluster-C had a superior prognosis and the highest CD8+ T-cells infiltrating abundance than other two clusters, demonstrating gene cluster-C was immune-inflamed phenotype (Figures 4D–F). Above all, these analyses verified the perspective that there were indeed three disparate immune phenotypes in colon cancer, representing different clinical and TME immune characterizations.

FIGURE 4

Characterization of DNA methylation score in different clinical and molecular traits

The above results demonstrated DNA methylation modification played a non-negligible role in TME formation and prognosis of colon cancer patients. Nevertheless, the remarkable heterogeneity between individual tumors limited the accurate application of DNA methylation modification patterns. Therefore, based on 152 prognostic DEGs, we constructed a DNA methylation gene signature, termed as DNA methylation score (DMS) to quantify the DNA methylation modification pattern of each patient.

We first evaluated the prognostic value of DMS in colon cancer, and determined −73.8 as the cut-off value to divide patients into low and high DMS groups. Prognostic analysis showed patients with low DMS had a prominently prolonged survival and recurrence-free time (Figures 5A,B). Besides, the multivariate Cox regression analysis demonstrated DMS was an independent risk prognostic biomarker for colon cancer (HR was 1.70 for OS and 1.66 for RFS; Supplementary Figure S5A). Additionally, to detect the stability of DMS, we enrolled other independent colon cancer cohorts to validate the prognostic value of DMS. Likewise, the DMS-low group had a preferable survival in the TCGA-COAD cohort, GSE72970, GSE103479, and GSE87211 cohorts (all p < 0.05, Supplementary Figures S5B–G). Furthermore, we continued to expand the application of DMS in other digestive system tumors, and we found rectum adenocarcinoma and stomach adenocarcinoma patients with low DMS presented a remarkably longer survival (both p < 0.05, Supplementary Figures S5H,I). The data above proved that the DMS we constructed could predict the survival of colon cancer effectively and stably.

FIGURE 5

Next, we explored the characteristics of DMS in different clinical and molecular subtype traits and fixed our attention on the GSE39582 cohort, which had comprehensive clinical information. We found DMS rose gradually by increasing the tumor TNM stage (p < 0.001, Figure 5C). Besides, patients with MSI, CIMP positive, CIN negative, and mutant BRAF showed a lower DMS, implying these subtypes may have a better prognosis (all p < 0.05, Figure 5C). Chemotherapy was a standard treatment strategy for colon cancer, and hence, we specifically examined the ability of DMS to predict the efficacy of chemotherapy. We revealed significantly longer recurrence-free time in the DMS-low group of colon cancer patients undergoing chemotherapy (Figure 5D), and the ROC curves to reflect the ability of DMS to predict the 1-, 3-, 5-, and 10-year survival probability were shown in Figure 5E.

Then, we investigated the distribution differences of somatic mutation between DMS-low and -high groups. As shown in Figures 5F,G, the DMS-low group exhibited a higher mutation frequency, and APC, TP53, and TTN were three common genes with the highest mutation frequency in both groups. As well, we quantified tumor mutation burden (TMB) value, and discovered a markedly negative correlation between DMS and TMB (Figure 5H). Same as the relation between MSI status and DMS, the TMB-high group presented a lower DMS, hinting low DMS may reflect genomic instability, massive tumor neoantigen generation, and induce immune activation (Figure 5I). Generally, these accumulated results manifested DMS was an effective predictive tool, and distinguished different molecular subtypes of colon cancer.

The role of DNA methylation score in the tumor microenvironment and biological processes

To better illustrate the underlying relevance between DMS and different molecular traits of colon cancer, we first compared DMS among different DNA methylation modification patterns and gene clusters. Pattern-B and gene cluster-C, representing the immune-inflamed phenotype, both had the lowest median DMS (Figures 6A,B). A similar distribution trend of DMS among three patterns and gene clusters was shown in Supplementary Figures S6A,B when we split the meta-cohort into three separate GSE cohorts. These results verified low DMS was related to better survival again, and the Sankey diagram visualized the distribution of individual colon cancer patients in different DNA methylation modification patterns, gene clusters, DMS, and clinical outcomes (Figure 6C).

FIGURE 6

Subsequently, we analyzed the correlation between DMS and several known biological processes signatures constructed by Mariathasan et al. (2018). The results revealed DMS was mightily negatively correlated with DNA damage repair-related signatures (Supplementary Figure S6C). Further comparative analysis demonstrated a higher activity of immune and DNA damage repair pathways in the DMS-low group, including CD8+ T effector, immune checkpoint, and mismatch repair (Figure 6D). Additional GSVA analysis revealed activated DNA damage repair and tumor suppresser-related pathways were significantly enriched in the DMS-low group, while stromal and carcinogenic activation pathways were markedly enriched in the DMS-high group (Figure 6E). Above analyses confirmed DNA methylation patterns participated in these biological processes and the DMS could reflect the activity of immune and stromal-related pathways.

We then examined the relationship between TME infiltrating cells and DMS using different immunocytes associated algorithms, including ssGSEA, CIBERSORT, and TIMER database, and we found the majority of immune cells, especially CD8+ T cells, presented a high infiltrating abundance in the DMS-low group (Figure 6F; Supplementary Figures S6G,H). Besides, the tumor microenvironment score, estimate score, and stromal score showed a positive correlation with DMS (Figure 6G; Supplementary Figures S6D,E). Furthermore, we compared the expression level of 24 DNA methylation regulators between DMS-low and -high groups and revealed that some regulators positively associated with immune cells were also highly expressed in the DMS-low group, such as TET2, MBD1, and MBD2 (Supplementary Figure S6F). At last, we investigated several molecules, including chemokines and cytokines associated with immune activation or suppression, which were referenced from published literature (Zeng et al., 2019). Of these molecules, CD8A, CXCL10, CXCL9, GZMA, GZMB, IFNG, PRF1, TBX2, and TNF were considered to be correlated to the transcripts of immune activation; and CD80, CD86, CTLA-4, HAVCR2, IDO1, LAG3, PD-1, PD-L1, PD-L2, TIGIT, and TNFRSF9 were related to the transcripts of immune checkpoints. We still found the expression of mRNAs relevant to immune activation was significantly upregulated in the DMS-low group (Supplementary Figure S6I). While two immune checkpoints, PD-1 and PD-L2, were highly expressed in the DMS-high group (Supplementary Figure S6J). Collectively, these results demonstrated DNA methylation involved in TME immunity response. Then combining the above conclusions that low DMS was related to higher TMB, we speculated that DNA methylation modification might influence the genetic mutation and generate a great number of tumoral neoantigens, inducing immune activation. These results provided a novel perspective to explore the mechanism of DNA methylation in tumor somatic mutation, and indirectly hinting DMS may play an important role in tumor immunotherapy.

The predictive value of DNA methylation score in anti-PD-1/PD-L1 immunotherapy

Immunotherapy represented by anti-PD-1/PD-L1 or CTLA-4 antibody had broadened the field of cancer treatment and brought huge clinical benefits in recent years. Here, we explored whether DMS could predict the therapeutic response and prognosis of patients treated with immunotherapy. In the IMvigor210 cohort, the DMS-low group presented a remarkably prolonged survival, and the multivariate Cox regression analysis also determined that higher DMS was an independent risk factor for prognosis (Figure 7A; Supplementary Figure S7A). The ROC curves for DMS in predicting 12-, 18-, 24-months survival probability were presented in Supplementary Figure S7B. The objective response rate (percentage of complete response and partial response) of patients with low DMS was significantly higher (Figure 7B). As well, patients in the response group showed a lower DMS compared to the non-response group, indicating a clinical benefit and treatment advantage of immunotherapy in patients with low DMS (Figure 7C; Supplementary Figure S7C).

FIGURE 7

Although we did not find a significant difference in PD-L1 expression level, a potential biomarker for immunotherapy, between DMS-low and -high groups in IMvigor210 cohort (Supplementary Figure S7D), the pan-cancer analysis demonstrated DMS was negatively related to PD-L1 and PD-L2 expression in the majority of cancer types including colon cancer, indirectly proving the reliability of DMS to predict the efficacy of immunotherapy (Figure 7D; Supplementary Figure S7E). In this immunotherapeutic cohort, the TMB value was also higher in the DMS-low group, which confirmed our above conclusions again (Figure 7E). In addition, we evaluated the prognostic value of TMB and DMS combination in immunotherapy, and observed a marked prognosis benefit in patients with low DMS and high TMB (Figure 7F).

More importantly, the biological processes signatures analyses showed that DNA damage repair-related pathways were significantly activated in the DMS-low group, while EMT and Pan-F-TBRS pathways were highly activated in the DMS-high group, indicating DMS was closely related to the DNA damage repair and stromal signatures in the setting of patients receiving immunotherapy (Figure 7G). Lastly, we investigated the difference of DMS among different immune phenotypes identified by the IMvigor210 study, and found the immune-excluded and -desert phenotypes had higher DMS, implying the response to immunotherapy of these subtypes was limited (Figure 7H). In summary, our work demonstrated a significant correlation between DNA methylation modification patterns and tumor immune phenotypes, and DMS could help to predict the response of immunotherapy.

Discussion

In this study, based on 24 DNA methylation regulators, we first identified three DNA methylation modification patterns with distinct TME infiltrating characteristics and biological behaviors in colon cancer. Moreover, we obtained prognostic DEGs among three modification patterns and established the DNA methylation gene signature, termed DMS, to quantify the DNA methylation modification profile of individual colon cancer, and more importantly, predict the efficacy and clinical outcome of immunotherapy.

Increasing evidence revealed tumors commonly hijacked various epigenetic mechanisms to escape the supervision of the immune system. Particularly, certain regulator mediated DNA methylation and demethylation played an indispensable role in adaptive immune response, including generation of tumoral neoantigen, dysregulation of antigen-presenting machinery, and suppression of anti-tumor cytokine production (Cao and Yan, 2020). However, the analysis of the whole landscape of DNA methylation modification in colon cancer was limited, and its impact on TME immune response remained unclear. Here, we integrated the transcriptomic information of 24 DNA methylation regulators and revealed three distinct DNA methylation modification patterns, and their features were consistent with three classical tumor immunophenotypes: immune-inflamed, -excluded, and -desert (Hegde and Chen, 2020). Modification pattern-B was immune-inflamed phenotype and presented an activated immune status, characterized by activation of immune, tumor suppressor and DNA damage repair related pathways, as well as a high infiltrating abundance of activated CD8+ T cell. Accordingly, patients in pattern-B had a better prognosis than the other two patterns. In detail, abundant immune cells were positioned in proximity to the tumor cells, accompanied by many proinflammatory and effector chemokines and cytokines, suggesting the presence of pre-existing anti-tumor immunity in the tumor parenchyma (Harlin et al., 2009; Gajewski et al., 2013; Herbst et al., 2014). Mechanically, an inflammatory TME was the basis of immune-inflamed phenotype, also known as a hot tumor, containing pro-inflammatory cytokines which provided a more favorable condition for T cell activation and expansion, including type I and type II interferons, TNF-α, IL-2, and IL-12. By contrast, pattern-C was immune-desert phenotype and exhibited an immune suppression status, characterized by inhibition of multiple immune and inflammatory-related pathways and a low proportion of TME infiltrating immune cells. This phenotype might be the result of immunological ignorance, the induction of immune tolerance, and a lack of appropriate environment for T cell activation, and hence, reflected the absence of pre-existing anti-tumor immune response (Gajewski et al., 2013; Herbst et al., 2014; Hegde et al., 2016; Kim and Chen. 2016). Notably, we found several innate immune cells, such as eosinophil, mast cell, natural killer cell, and plasmacytoid dendritic cell, were enriched in the TME of pattern-A, while the patients in this pattern did not present a matching preferable prognosis as pattern B. Conceptually, the immune-excluded phenotype was characterized by stromal activation and the presence of vast immune cells, while they were retained in the stroma surrounding the tumor cells nests rather than penetrated the tumor parenchyma (Salmon et al., 2012; Herbst et al., 2014; Joyce and Fearon, 2015; Hegde et al., 2016). Previous study also provided evidence that T-cells proliferation and activation were observed after anti-PD1/PD-L1 agents but no infiltration (Herbst et al., 2014). Therefore, pattern-A was classified as immune-excluded phenotype, presenting a pre-existing anti-tumor immune response while being rendered ineffective by the retention of immune cells in the surrounding stroma. Generally, the immune-desert and -excluded phenotypes were both considered as the cold tumor, containing numerous immune-suppressive cytokines that contributed to impairing the anti-tumor response. These could explain the alike poorer prognosis of pattern-A and -C.

We next screened DEGs among three patterns, and GO functional annotation revealed they were significantly associated with DNA modification and damage repair related pathways, suggesting the different clinical and biological characteristics among three patterns might be the results of differentially expressed of these genes. We further identified prognostic DEGs, termed as DNA methylation signature genes, to perform unsupervised clustering. Likewise, we found three genomic subtypes, named DNA methylation gene clusters, whose clinical outcomes and immune cell infiltrating traits were similar to three modification patterns. Our comprehensive analyses strongly revealed three immune phenotypes in colon cancer with distinct clinical and TME immune characteristics, which enhanced our understanding of the non-negligible impact of DNA methylation in shaping different TME landscapes.

Whereas, above analyses were performed based on the patient population, which could not accurately predict the specific modification pattern in individual patients. As a consequence, based on the above signature genes, we developed the DNA methylation score, DMS, to quantify the certain modification pattern for each colon cancer patient. We found DMS could precisely discriminate three immune phenotypes, with the lowest median DMS in immune-inflamed type and the highest median DMS in immune-desert type. Additionally, we revealed a markedly positive relation between DMS and tumor stage, with DMS increasing gradually from the stage I to IV. Moreover, DMS was an independent prognostic risk factor, and patients with high DMS presented an inferior survival, which was validated in multiple colon cancer cohorts. Integrally, above results indicated DMS was a reliable tool to reflect the individual DNA methylation modification pattern, and predict clinical outcomes of colon cancer.

Colon cancer was a highly heterogeneous disease, resulting from a series of distinct genetic and epigenetic changes, and a subset of molecular alterations was considered to drive the cellular and clinical behavior of cancer, including MSI, CIMP, CIN, BRAF, and KRAS mutations (Marisa et al., 2013; Phipps et al., 2015). The CIMP subtype was referred to a distinct epigenome with a high frequency of methylated genes, and approximately 20% of colon cancer were CIMP+ tumors (Toyota et al., 1999). Previous study demonstrated that CIMP+ colon tumors had a unique association with BRAF^V600E oncogene mutation, and CIMP-associated methylation of MLH1 induced mismatch repair deficiency and resulted in a genomic instability status, also known as MSI, to generate more mutation burden and neoantigen (Weisenberger et al., 2006). In addition, aberrant epigenetic alterations as well contributed to the dysregulation of antigen-presenting machinery in tumor cells, leading to acquiring the adaptive immune response (Alexandrov et al., 2013). Here, we revealed DMS-low subtype was characterized as CIMP positive, mutant BRAF, MSI, and higher TMB, indicating this epigenotype of colon cancer identified by DMS had specific molecular alterations. Further analyses manifested the activation of immune and DNA damage repair related pathways, as well as abundant immune cells infiltration in the DMS-low group. Integrally, our findings substantiated DNA methylation modification was involved in the genomic instability, resulting in the accumulation of tumoral mutation burden and generation of neoantigen, which on the one hand activated the DNA damage repair pathway, and on the other hand enhanced immunogenicity to further activate the immune response.

Emerging evidence demonstrated that MSI and elevated TMB could heighten the anti-tumor activity of immunotherapy (Rizvi et al., 2015; Le et al., 2017; Mouw et al., 2017), and hence, we further investigated the ability of DMS to predict efficacy and prognosis of patients received anti-PD-1/PD-L1 antibody. In an immunotherapy cohort of advanced urothelial carcinoma patients treated with atezolizumab, we confirmed DMS-low group had a higher objective response rate and prolonged survival time than the DMS-high group. Higher PD-L1 expression and TMB value were considered to imply a favorable efficacy of immunotherapy (Gibney et al., 2016), and in this study, we discovered DMS was negatively correlated to TMB and PD-L1 expression in colon cancer, which validated the predictive value of DMS indirectly. Deeply, the pathway and immune phenotype analyses revealed activation of DNA damage repair and immune-inflamed environment in DMS-low patients, providing clues to illustrate the regulatory mechanism of DNA methylation in shaping TME immune landscape, and also confirmed patients with low DMS could benefit from immunotherapy. Collectively, we considered DMS had the potential to be an excellent predictive biomarker for immune checkpoint inhibitor, and might promote personalized colon cancer immunotherapy in the future.

A major limitation of this work was the public survival and transcriptomic data of colon cancer immunotherapy was not accessible yet. Therefore, the predictive performance of DMS needed to be further certified in the colon cancer immunotherapeutic cohort.

In conclusion, for the first time, we uncovered three distinct DNA methylation modification patterns in colon cancer, and illustrated their extensive regulatory mechanism in tumor immune environment formation, which was a non-negligible factor to cause individual TME heterogeneity and different clinical outcomes. Our integrated analyses of DNA methylation modification would contribute to enhancing the understanding of tumor immune characteristics, and providing novel insights to guide immunotherapy more effectively.

References

• 1

  AlexandrovL. B.Nik-ZainalS.WedgeD. C.AparicioS. A.BehjatiS.BiankinA. V.et al (2013). Signatures of mutational processes in human cancer. Nature500 (7463), 415–421. 10.1038/nature12477

  • CrossRef
  • Google Scholar

• 2

  AndreT.ShiuK. K.KimT. W.JensenB. V.JensenL. H.PuntC.et al (2020). Pembrolizumab in microsatellite-instability-high advanced colorectal cancer. N. Engl. J. Med.383 (23), 2207–2218. 10.1056/NEJMoa2017699

  • CrossRef
  • Google Scholar

• 3

  AranD.HuZ.ButteA. J. (2017). xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol.18 (1), 220. 10.1186/s13059-017-1349-1

  • CrossRef
  • Google Scholar

• 4

  BestorT. H. (2000). The DNA methyltransferases of mammals. Hum. Mol. Genet.9 (16), 2395–2402. 10.1093/hmg/9.16.2395

  • CrossRef
  • Google Scholar

• 5

  CaoJ.YanQ. (2020). Cancer epigenetics, tumor immunity, and immunotherapy. Trends Cancer6 (7), 580–592. 10.1016/j.trecan.2020.02.003

  • CrossRef
  • Google Scholar

• 6

  CharoentongP.FinotelloF.AngelovaM.MayerC.EfremovaM.RiederD.et al (2017). Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep.18 (1), 248–262. 10.1016/j.celrep.2016.12.019

  • CrossRef
  • Google Scholar

• 7

  ChenD. S.MellmanI. (2017). Elements of cancer immunity and the cancer-immune set point. Nature541 (7637), 321–330. 10.1038/nature21349

  • CrossRef
  • Google Scholar

• 8

  ChenK.LiuJ.LiuS.XiaM.ZhangX.HanD.et al (2017). Methyltransferase SETD2-mediated methylation of STAT1 is critical for interferon antiviral activity. Cell170 (3), 492–506. 10.1016/j.cell.2017.06.042

  • CrossRef
  • Google Scholar

• 9

  ColapricoA.SilvaT. C.OlsenC.GarofanoL.CavaC.GaroliniD.et al (2016). TCGAbiolinks: An R/bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res.44 (8), e71. 10.1093/nar/gkv1507

  • CrossRef
  • Google Scholar

• 10

  CookP. C.OwenH.DeatonA. M.BorgerJ. G.BrownS. L.ClouaireT.et al (2015). A dominant role for the methyl-CpG-binding protein Mbd2 in controlling Th2 induction by dendritic cells. Nat. Commun.6, 6920. 10.1038/ncomms7920

  • CrossRef
  • Google Scholar

• 11

  DienstmannR.VermeulenL.GuinneyJ.KopetzS.TejparS.TaberneroJ. (2017). Consensus molecular subtypes and the evolution of precision medicine in colorectal cancer. Nat. Rev. Cancer17 (2), 79–92. 10.1038/nrc.2016.126

  • CrossRef
  • Google Scholar

• 12

  GajewskiT. F.WooS. R.ZhaY.SpaapenR.ZhengY.CorralesL.et al (2013). Cancer immunotherapy strategies based on overcoming barriers within the tumor microenvironment. Curr. Opin. Immunol.25 (2), 268–276. 10.1016/j.coi.2013.02.009

  • CrossRef
  • Google Scholar

• 13

  GaoY.WangH.LiH.YeX.XiaY.YuanS.et al (2021). Integrated analyses of m(1)A regulator-mediated modification patterns in tumor microenvironment-infiltrating immune cells in colon cancer. Oncoimmunology10 (1), 1936758. 10.1080/2162402X.2021.1936758

  • CrossRef
  • Google Scholar

• 14

  GaujouxR.SeoigheC. (2010). A flexible R package for nonnegative matrix factorization. BMC Bioinforma.11, 367. 10.1186/1471-2105-11-367

  • CrossRef
  • Google Scholar

• 15

  GibneyG. T.WeinerL. M.AtkinsM. B. (2016). Predictive biomarkers for checkpoint inhibitor-based immunotherapy. Lancet. Oncol.17 (12), e542–e551. 10.1016/s1470-2045(16)30406-5

  • CrossRef
  • Google Scholar

• 16

  GoltzD.GevenslebenH.GrunenS.DietrichJ.KristiansenG.LandsbergJ.et al (2017). PD-L1 (CD274) promoter methylation predicts survival in patients with acute myeloid leukemia. Leukemia31 (3), 738–743. 10.1038/leu.2016.328

  • CrossRef
  • Google Scholar

• 17

  HanzelmannS.CasteloR.GuinneyJ. (2013). Gsva: Gene set variation analysis for microarray and RNA-seq data. BMC Bioinforma.14, 7. 10.1186/1471-2105-14-7

  • CrossRef
  • Google Scholar

• 18

  HarlinH.MengY.PetersonA. C.ZhaY.TretiakovaM.SlingluffC.et al (2009). Chemokine expression in melanoma metastases associated with CD8+ T-cell recruitment. Cancer Res.69 (7), 3077–3085. 10.1158/0008-5472.CAN-08-2281

  • CrossRef
  • Google Scholar

• 19

  HegdeP. S.ChenD. S. (2020). Top 10 challenges in cancer immunotherapy. Immunity52 (1), 17–35. 10.1016/j.immuni.2019.12.011

  • CrossRef
  • Google Scholar

• 20

  HegdeP. S.KaranikasV.EversS. (2016). The where, the when, and the how of immune monitoring for cancer immunotherapies in the era of checkpoint inhibition. Clin. Cancer Res.22 (8), 1865–1874. 10.1158/1078-0432.CCR-15-1507

  • CrossRef
  • Google Scholar

• 21

  HeilandD. H.HaakerG.DelevD.MercasB.MasalhaW.HeynckesS.et al (2017). Comprehensive analysis of PD-L1 expression in glioblastoma multiforme. Oncotarget8 (26), 42214–42225. 10.18632/oncotarget.15031

  • CrossRef
  • Google Scholar

• 22

  HerbstR. S.SoriaJ. C.KowanetzM.FineG. D.HamidO.GordonM. S.et al (2014). Predictive correlates of response to the anti-PD-L1 antibody MPDL3280A in cancer patients. Nature515 (7528), 563–567. 10.1038/nature14011

  • CrossRef
  • Google Scholar

• 23

  HermanJ. G.UmarA.PolyakK.GraffJ. R.AhujaN.IssaJ. P.et al (1998). Incidence and functional consequences of hMLH1 promoter hypermethylation in colorectal carcinoma. Proc. Natl. Acad. Sci. U. S. A.95 (12), 6870–6875. 10.1073/pnas.95.12.6870

  • CrossRef
  • Google Scholar

• 24

  JoyceJ. A.FearonD. T. (2015). T cell exclusion, immune privilege, and the tumor microenvironment. Science348 (6230), 74–80. 10.1126/science.aaa6204

  • CrossRef
  • Google Scholar

• 25

  KimJ. M.ChenD. S. (2016). Immune escape to PD-L1/PD-1 blockade: Seven steps to success (or failure). Ann. Oncol.27 (8), 1492–1504. 10.1093/annonc/mdw217

  • CrossRef
  • Google Scholar

• 26

  LaoV. V.GradyW. M. (2011). Epigenetics and colorectal cancer. Nat. Rev. Gastroenterol. Hepatol.8 (12), 686–700. 10.1038/nrgastro.2011.173

  • CrossRef
  • Google Scholar

• 27

  LeD. T.DurhamJ. N.SmithK. N.WangH.BartlettB. R.AulakhL. K.et al (2017). Mismatch repair deficiency predicts response of solid tumors to PD-1 blockade. Science357 (6349), 409–413. 10.1126/science.aan6733

  • CrossRef
  • Google Scholar

• 28

  LeD. T.UramJ. N.WangH.BartlettB. R.KemberlingH.EyringA. D.et al (2015). PD-1 blockade in tumors with mismatch-repair deficiency. N. Engl. J. Med.372 (26), 2509–2520. 10.1056/NEJMoa1500596

  • CrossRef
  • Google Scholar

• 29

  Loo YauH.EttayebiI.De CarvalhoD. D. (2019). The cancer epigenome: Exploiting its vulnerabilities for immunotherapy. Trends Cell Biol.29 (1), 31–43. 10.1016/j.tcb.2018.07.006

  • CrossRef
  • Google Scholar

• 30

  LuoN.NixonM. J.Gonzalez-EricssonP. I.SanchezV.OpalenikS. R.LiH.et al (2018). DNA methyltransferase inhibition upregulates MHC-I to potentiate cytotoxic T lymphocyte responses in breast cancer. Nat. Commun.9 (1), 248. 10.1038/s41467-017-02630-w

  • CrossRef
  • Google Scholar

• 31

  MariathasanS.TurleyS. J.NicklesD.CastiglioniA.YuenK.WangY.et al (2018). TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature554 (7693), 544–548. 10.1038/nature25501

  • CrossRef
  • Google Scholar

• 32

  MarisaL.de ReyniesA.DuvalA.SelvesJ.GaubM. P.VescovoL.et al (2013). Gene expression classification of colon cancer into molecular subtypes: Characterization, validation, and prognostic value. PLoS Med.10 (5), e1001453. 10.1371/journal.pmed.1001453

  • CrossRef
  • Google Scholar

• 33

  MouwK. W.GoldbergM. S.KonstantinopoulosP. A.D'AndreaA. D. (2017). DNA damage and repair biomarkers of immunotherapy response. Cancer Discov.7 (7), 675–693. 10.1158/2159-8290.CD-17-0226

  • CrossRef
  • Google Scholar

• 34

  NguyenL. H.GoelA.ChungD. C. (2020). Pathways of colorectal carcinogenesis. Gastroenterology158 (2), 291–302. 10.1053/j.gastro.2019.08.059

  • CrossRef
  • Google Scholar

• 35

  PhippsA. I.LimburgP. J.BaronJ. A.Burnett-HartmanA. N.WeisenbergerD. J.LairdP. W.et al (2015). Association between molecular subtypes of colorectal cancer and patient survival. Gastroenterology148 (1), 77–87. 10.1053/j.gastro.2014.09.038

  • CrossRef
  • Google Scholar

• 36

  RitchieM. E.PhipsonB.WuD.HuY.LawC. W.ShiW.et al (2015). Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res.43 (7), e47. 10.1093/nar/gkv007

  • CrossRef
  • Google Scholar

• 37

  RizviN. A.HellmannM. D.SnyderA.KvistborgP.MakarovV.HavelJ. J.et al (2015). Cancer immunology. Mutational landscape determines sensitivity to PD-1 blockade in non-small cell lung cancer. Science348 (6230), 124–128. 10.1126/science.aaa1348

  • CrossRef
  • Google Scholar

• 38

  SalmonH.FranciszkiewiczK.DamotteD.Dieu-NosjeanM. C.ValidireP.TrautmannA.et al (2012). Matrix architecture defines the preferential localization and migration of T cells into the stroma of human lung tumors. J. Clin. Invest.122 (3), 899–910. 10.1172/JCI45817

  • CrossRef
  • Google Scholar

• 39

  ShannonP.MarkielA.OzierO.BaligaN. S.WangJ. T.RamageD.et al (2003). Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Res.13 (11), 2498–2504. 10.1101/gr.1239303

  • CrossRef
  • Google Scholar

• 40

  SungH.FerlayJ.SiegelR. L.LaversanneM.SoerjomataramI.JemalA.et al (2021). Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. Ca. Cancer J. Clin.71 (3), 209–249. 10.3322/caac.21660

  • CrossRef
  • Google Scholar

• 41

  SuterC. M.MartinD. I.WardR. L. (2004). Hypomethylation of L1 retrotransposons in colorectal cancer and adjacent normal tissue. Int. J. Colorectal Dis.19 (2), 95–101. 10.1007/s00384-003-0539-3

  • CrossRef
  • Google Scholar

• 42

  SzklarczykD.GableA. L.LyonD.JungeA.WyderS.Huerta-CepasJ.et al (2019). STRING v11: Protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res.47 (D1), D607-D613–D613. 10.1093/nar/gky1131

  • CrossRef
  • Google Scholar

• 43

  ThissenD.SteinbergL.KuangD. (2002). Quick and easy implementation of the benjamini-hochberg procedure for controlling the false positive rate in multiple comparisons. J. Educ. Behav. Statistics27 (1), 77–83. 10.3102/10769986027001077

  • CrossRef
  • Google Scholar

• 44

  ToyotaM.AhujaN.Ohe-ToyotaM.HermanJ. G.BaylinS. B.IssaJ. P. (1999). CpG island methylator phenotype in colorectal cancer. Proc. Natl. Acad. Sci. U. S. A.96 (15), 8681–8686. 10.1073/pnas.96.15.8681

  • CrossRef
  • Google Scholar

• 45

  WeisenbergerD. J.SiegmundK. D.CampanM.YoungJ.LongT. I.FaasseM. A.et al (2006). CpG island methylator phenotype underlies sporadic microsatellite instability and is tightly associated with BRAF mutation in colorectal cancer. Nat. Genet.38 (7), 787–793. 10.1038/ng1834

  • CrossRef
  • Google Scholar

• 46

  YoshiharaK.ShahmoradgoliM.MartinezE.VegesnaR.KimH.Torres-GarciaW.et al (2013). Inferring tumour purity and stromal and immune cell admixture from expression data. Nat. Commun.4, 2612. 10.1038/ncomms3612

  • CrossRef
  • Google Scholar

• 47

  YuG.WangL. G.HanY.HeQ. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS16 (5), 284–287. 10.1089/omi.2011.0118

  • CrossRef
  • Google Scholar

• 48

  ZengD.LiM.ZhouR.ZhangJ.SunH.ShiM.et al (2019). Tumor microenvironment characterization in gastric cancer identifies prognostic and immunotherapeutically relevant gene signatures. Cancer Immunol. Res.7 (5), 737–750. 10.1158/2326-6066.CIR-18-0436

  • CrossRef
  • Google Scholar

• 49

  ZhangB.WuQ.LiB.WangD.WangL.ZhouY. L. (2020). m(6)A regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in gastric cancer. Mol. Cancer19 (1), 53. 10.1186/s12943-020-01170-0

  • CrossRef
  • Google Scholar

• 50

  ZhangQ.ZhaoK.ShenQ.HanY.GuY.LiX.et al (2015). Tet2 is required to resolve inflammation by recruiting Hdac2 to specifically repress IL-6. Nature525 (7569), 389–393. 10.1038/nature15252

  • CrossRef
  • Google Scholar