The workflow of this study is presented in Figure S1A . Gene expression data and complete clinical annotations of LUAD samples were retrospectively collected from publicly available datasets of the Gene-Expression Omnibus (GEO) and TCGA databases. A total of six GEO lung adenocarcinoma cohort somatic mutation, copy number variation (CNV), and clinical data, including tumor stage, age, sex, and overall survival times/states were obtained from TCGA databases (GSE116959, GSE58772, GSE99995, GSE68571, GSE68465, and GSE26939) and one TCGA-LUAD cohort were enrolled for further analysis. For the GEO microarray data, the normalized matrix files were directly downloaded. As for TCGA cohort, the RNA sequencing data (FPKM value) for gene expression was downloaded from the Genomic Data Commons and then transformed into transcripts per kilobase million (TPM) format. Batch effects among GEO datasets were corrected using the “ComBat” algorithm of the sva package. The baseline information of all eligible datasets is summarized in Table S1 .

We collected DNA methylation regulator-related studies, and a total of 21 regulator genes (including “writers’: DNMT1, DNMT3A, DNMT3B; “easers”: TET1, TET2, TET3; and “readers”: MBD1, MBD2, MBD3, MBD4, TDG, SMUG1, UHRF1, UHRF2, ZBTB4, ZBTB24, ZBTB33, ZBTB38, NSUN2, MGMT, DMAP1) were extracted for DNA methylation modification pattern identification ( Table S2 ). The ConsensusClusterPlus algorithm was employed to conduct unsupervised clustering of individual tumor samples with the gene expression profiles of 21 DNA methylation regulators (18). The cluster assignments were stable when k = 2.

We performed Gene-Set Variation Analysis (GSVA) enrichment analysis with “GSVA” R packages to study the differed biological pathway between different DNA methylation modification patterns in cancer samples (19). The Hallmarker gene-set of “c2.cp.kegg.v7.4.symbols” was retrieved from the MSigDB database (20). Adjusted p-values <0.05 were considered statistically significant for GSVA analysis. Gene ontology (GO) and pathway annotation for DNA methylation pattern-related genes were performed using the R package “clusterProfiler” with a cutoff value of p < 0.05.

The Single Sample Gene-Set Enrichment Analysis (ssGSEA) was implemented to determine variations in immune leukocyte subtype abundance between different DMP clusters using the R package “GSEAbase.” Subsequently, the abundances of 22 immune cell types for each tumor specimen were further identified by estimating relative subsets of RNA transcripts (CIBERSORT; https://cibersort.stanford.edu/) using the gene expression profile of LUAD cancer (21).

We used three R packages (“limma,” “edgeR,” and “Deseq2”) to identify the DNA methylation modification-related differentially expressed genes (DEGs) across distinct DMP phenotypes. Univariate Cox model analysis was performed to calculate their association with overall survival and to extract prognostic DEGs to construct a scoring system. We then conducted principal component analysis (PCA) using the identified DEG prognostic genes, and PCA 1–2 components were selected to act as signature scores to construct DNA methylation-relevant gene signature, which we defined as the DMS. This method mainly focused on evaluating the score for each patient in the dataset with the largest group of well-correlated (or anticorrelated) genes. We constructed the DMS using an algorithm similar to the GGI: DMS = Σ(PC1 _i + PC2 _i ), where i stood for the relevant gene expression of the selected set (22, 23).

We searched the immunotherapeutic characteristics of the TCGA-LUAD patients using the Cancer Immunome Database (TCIA), which is a publicly available dataset containing corresponding clinical pathology information (24). The immunotherapeutic GEO datasets were included in this study: GSE135222 (anti-PD-1/PD-L1 treatment) (25) and GSE91061 (anti-CTLA-4/PD-1 treatment) (26). We also retrieved IMvigor210 datasets containing data on atezolizumab treatment and extracted the relative gene expression profiles and clinical notes using the R “IMvigor210CoreBiologies” package (27). The raw count data were then transformed into TPM value format.

In this study, 21 DNA methylation regulators were collected from the published data. Figure 1A displays the dynamic reversible process of RNA methylation mediated by regulators and their corresponding potential biological functions. We further analyzed the genetic alterations of DNA methylation regulators in LUAD. Among the 561 cases, 123 samples (21.93%) harbored genetic alterations of DNA methylation regulators, primarily including nonsense or missense mutations. TET1 (4%) had the highest mutation frequency, followed by DNMT3A, TET3, and DNMT3B ( Figure 1B ). A mutation cooccurrence pattern across all DNA methylation regulators was examined, and significant gene mutation relationships were identified between DNMT1 and ZBTB4, MBD4 and NSUN2, and DNMT1 and UHRF1 ( Figure S1B ). Further analysis of CNV alteration frequency in the 21 regulators showed a prevalence of CNV mutations. Among these, NSUN2, DMAP1, SMUG1, DNMT3B, ZBTB33/38, and DNMT3A showed relatively higher amplification frequency in terms of CNV, while MBD1/2/3, ZBTB4/24, UHRF1, TDG, and TET2 had a widespread frequency of copy number deletion variants ( Figure 1C ). The location of CNV among the DNA methylation regulators on chromosomes is shown in Figure 1D . To investigate whether these genetic variants influenced the gene expression of regulators in cancer samples, we compared the expression of the 21 DNA methylation regulators between normal and LUAD samples and found that regulators with significant genetic mutation as well as CNV amplification were significantly higher expressed in cancer groups ( Figures 1B–E ), indicating that the regulator gene mutation and the alternation of CNV could be risk factors leading to the multiperturbation of the translation of DNA methylation regulators. Univariable and multivariable Cox proportional hazards regression analyses were performed to evaluate the association between regulators and patient survival. The Forest plot and Kaplan-Meier curve showed that DMAP1 and MGMT were upregulated in tumor specimens and could be considered protective factors associated with prolonged survival time, while ZBTB38 and MBD2/3 downregulation in the tumor group was recognized as a risk factor with worse overall survival ( Figures S1C–H ), indicating that the unbalance expression of DNA methylation regulators could contribute to tumorigenesis. These analyses demonstrated that the crosstalk between DNA methylation regulators in the genomic and transcriptomic landscape has a potential role in LUAD occurrence and complexity.

To determine the pattern of DNA methylation modification mediated by 21 regulators, the TCGA-LUAD cohort and their corresponding clinical data were prepared for analysis ( Table S3 ). A comprehensive landscape of the interactions and the prognostic significance of the 21 regulators was visualized using the conetwork plot ( Figure 2A ). We dissected the relationship among those regulators and found that most regulators showed a positive correlation with each other, which is consistent with the above analysis. We also demonstrated that the expression of UHRF1, SMUG1, TDG, and MBD1/2/3 presented the most significant positive correlation with other regulators and might be the risk factors for carcinoma formation. Notably, the readers, such as MGMT and ZBTB4, showed a remarkably negative relationship with the other regulators. Thus, these findings indicated that crosstalk among DNA methylation regulators may play a critical role in the onset of LUAD.

We continued to conduct unsupervised clustering based on the expression of 21 DNA methylation regulators to identify DNA methylation modification patterns in LUAD samples, and two distinct clusters were accordingly obtained via the ConsensusClusterPlus package in R software, including 246 cases in cluster A and 267 cases in cluster B ( Figure 2B ). We renamed those clusters as DNA methylation pattern (DMP-A and DMP-B) ( Table S4 ). Survival analysis for the two DNA methylation modification clusters showed that DMP-B had a relatively better prognosis in the TCGA cohort ( Figure 2C ). In addition, we conducted heatmap analysis of the relative expression of 21 DNA methylation regulators between distinct DNA methylation modification patterns and observed that most regulators were markedly elevated in the DMP-A cluster ( Figure 2B ). Principal component analysis (PCA) also indicated that two DMP clusters were completely segregated, indicating that they could be easily distinguished via the gene expression pattern of the 21 regulators ( Figure 2D ). The heatmap analysis of the tumor methylation status determined that the DMP-B cluster presented a low-methylation epigenotype with worser clinical outcomes ( Figure 2F ). These results demonstrated that altered expression of DNA methylation regulators could affect patient survival by regulating tumor methylation levels, which contributes to the high heterogeneity of LUAD.

To further explore differences in biological behaviors among the two DNA methylation modification patterns, we performed GSVA enrichment analysis. As shown in Figure 2E and Table S5 , DMP-A was mainly presented a CD4+ immune regulator and stromal activation phenotype, associated with many related pathways, such as the cell cycle, DNA replication, Toll-like receptor, and Nod-like receptor signaling pathways in cancer. Whereas, DMP-B represented enriched pathways associated with immune/inflammation hyperactivation including the activation of the PPAR signaling pathway and the autoimmune response, including asthma, autoimmune thyroid disease, allograft rejection, and graft-versus-host disease. Furthermore, we conducted a TME cell infiltration analysis using the ssGSVA algorithm and demonstrated that DMP-A showed a relatively lower proportion of immune cell infiltration and exhibited a shorter survival time ( Figures 2B , 3A ). However, DMP-B did present a matching survival advantage and significantly elevated innate immune cells including activated B cells, MDSCs, macrophages, mast cells, natural killer cells, and activated CD8+ T cells and so on. The CIBERSORT analysis also indicated that DMP-B was markedly enriched in immunoactive cells, such as M2 macrophages, memory B cells, plasma cells, and Tregs, which is coherently indicative of a better prognosis ( Figures 2B , 3B ). Previous studies have reported that the immune-excluded phenotype in tumors also presented immune cells limited to the stroma surrounding the tumor cell nests (28). Therefore, we suspected that the stromal activation and low immune cell density in DMP-A pattern caused the suppression of antitumor response, which was also validated in the GEO LUAD cohort (GSE116959, GSE58772, and GSE99995) ( Figure S2I ). Altogether, these results indicated that LUAD patients could be classified into two immune phenotype groups, among which the DMP-A pattern represented a noninflamed phenotype (immune-excluded or immune-deserted trait) characterized by stromal activation and weakened immune regulator status, whereas the DMP-B pattern represented an immune-inflamed phenotype characterized by immune/inflammation hyperactivation and multicomponent immune cell infiltration in the TME.

To further explore the characteristics of the DMP in the LUAD cohort, a total of 2,720 DNA methylation regulator pattern related to the DEGs were obtained and subjected to conduct enrichment analysis. The KEGG pathway related to the cell cycle and the p53 signaling pathway were also significantly enriched in DMP-related gene set ( Figure 3C ). The results of GO terms showed that those genes were mainly enriched in the cell cycle process and DNA replication ( Figure 3D ). We subsequently selected 792 prognostic genes (p < 0.05) via univariate Cox model analysis ( Table S6 ). Based on the expression of the 792 representative genes, unsupervised clustering analysis was performed to classify LUAD patients into three distinct clusters, namely, DNA methylation gene cluster I/II/III ( Figures 4A, B ; Figures S2A–E ). The heatmap and survival analysis showed that patients in gene cluster III were characterized by earlier clinical stage and better prognosis, while patients in gene cluster II were characterized by advanced staging and worser survival outcome ( Figures 4B, C ). In addition, multivariate Cox regression analysis indicated that the overall survival between the gene cluster gene signatures remained statistically significant after considering multiple factors, including age, tumor stage (T and N), and sex ( Figure 4D ), indicating that the identified gene clusters could represent independent factors of prognosis. We also noted that the gene expression of the 21 regulators exhibited significant differences among the three gene clusters, and in particularly patients in gene cluster III were related to the higher expression of protective factors (TET2/MGMT), which is in line with the findings regarding DNA methylation modification patterns ( Figure 4E ). The patients in gene cluster III also presented higher immune cell infiltration, indicating a phototype of immune activation. While the levels of almost all immune cells in gene cluster II were relatively lower than in other groups, these presented a phenotype of immune suppression ( Figure 4F ). We further evaluated the DNA methylation status of the gene clusters, and results indicated that the gene cluster II exhibited a hypermethylation status, whereas the gene cluster III presented a weaker methylation status ( Figure 4G ). All these results also confirmed that the perturbed DNA methylation regulator genetic expression was highly correlated with different immune responses, genetic methylation levels, and clinical outcomes and contributed to LUAD heterogeneity, which was consistent with the previous findings.

To accurately evaluate the DNA methylation status of individual patients with LUAD, we developed a scoring scheme termed the DMS based on the expression of 792 DNA methylation-related genes.

The DMP-B pattern had a higher level of DMS compared with patients in DMP-A ( Figure 5A ). Notably, gene cluster III showed the highest level of DMS, followed by gene clusters I and II ( Figure 5B ). With an optimum cutoff value of 2.17 determined by the survminer package, we divided LUAD patients into high- and low-DMS groups. The alluvial diagram summarized the attribute changes of patients according to the DNA methylation regulator pattern, gene cluster, and DMS groups ( Figure 5C ). We then examined the correlation between biological processes, immune cell infiltration, and the level of DMS signatures using Spearman’s analysis. The DMS signatures were significantly negatively correlated with cell cycle and DNA/RNA repair signatures but were positively correlated with immune activation and stromal-relevant pathways ( Figure 5D ; Figure S2F ). In addition, we calculated the distribution differences of tumor genome somatic mutations between the DMS groups using the R “maftools” package. As we expected, the low-DMS group with shorter survival time presented significantly higher tumor mutational profiles than the high-DMS group, consistent with recent findings that gene alterations correlated highly with tumor invasion and cell proliferation. Furthermore, survival analysis revealed that the DMS-low patients were significantly correlated with a worse prognosis in the TCGA cohort (p < 0.001, Figure 5E ), which was further validated by the ROC curves (AUC = 0.652, Figure 5F ). The multivariate Cox regression model confirmed the DMS factor could stand as an independent prognostic biomarker for evaluating patient outcomes in the TCGA-LUAD cohort (HR, 2.03 [95% confidence interval, 1.34–3.08], p < 0.001, Figure 5G ). To further verify this DMS model, we also performed multivariate Cox regression model and prognosis analysis using meta-GEO cohorts (GSE68465, GSE68571, and GSE26939), which provided further support that DMS was a significant prognostic factor for predicting patient outcomes ( Figures S2G, H ). The analysis of tumor mutation burden (TMB) confirmed that the low-DMS cluster was significantly correlated with a higher TMB ( Figures 5H, I ). The mutational landscape of significantly mutated genes showed that most genes (including TP53 and TTN) had higher somatic alteration rates in the lower-DMS group, whereas KRAS had a higher somatic alteration rate in the higher-DMS group. We further observed that there was a markedly negative correlation between the TMB score and DMS (R = −0.48, p < 0.001), demonstrating the crosstalk between the DMS and genetic mutation evaluation ( Figure 5J ). Meanwhile, the methylation heatmap analysis revealed that the higher-DMS group exhibited hypomethylation compared with the lower-DMS group ( Figure S2J ). These results demonstrated that the DMS could represent the DNA methylation modification patterns and comprehensively could reflect genomic variation modifications, as well as effectively predict the prognosis of LUAD patients.

Significant progress has been made to identify the effective immune-related signatures that correlated with response to antitumor therapy, particularly represented by TMB and ICIs. Our analysis found that the gene expression of the CTLA-4/PD-1/PD-L1 in TCGA-LUAD cohort was significantly increased in the DMS-low group ( Figures 6A–C ). Considering the high correlation across DMS levels with the immune response, we explored whether the DMS system could predict the patients’ response to ICI treatment in two independent immunotherapy cohorts treated with CTLA-4/PD-1/PD-L1 antibody inhibitors. In the TCIA-LUAD cohort (anti-CTLA-4/PD-1 treatment), patients classified as the DMS-high group exhibited a significantly clinical response ( Figures 6D–G ), indicating the immunotherapeutic benefits of CTLA-4/PD-1 antibody treatment in DMS-high patients. In the anti-CTLA-4/PD-1/PD-L1 GEO cohort (GSE91061 and GSE135222), higher DMS in tumor patients was associated with a stronger immune response and higher therapeutic benefits ( Figures 6H, I ). The progression-free survival analysis showed that higher-DMS patients exhibited prolonged survival ( Figure 6I ). A similar result was also identified in the IMvigor210 cohort (anti-PD-L1 treatment) ( Figures 6J–M ). Those results demonstrated that the DMS signature system was strongly associated with the tumor immune response and might contribute to predicting the efficacy value of the anti-CTLA-4/PD-1/PD-L1 immunotherapy.