Comprehensive analysis of m7G modification patterns based on potential m7G regulators and tumor microenvironment infiltration characterization in lung adenocarcinoma

Background: The non-negligible role of epigenetic modifications in cancer development and tumor microenvironment (TME) has been demonstrated in recent studies. Nonetheless, the potential regulatory role of N7-methylguanosine (m⁷G) modification in shaping and impacting the TME remains unclear.

Methods: A comprehensive analysis was performed to explore the m⁷G modification patterns based on 24 potential m⁷G regulators in 817 lung adenocarcinoma (LUAD) patients, and the TME landscape in distinct m⁷G modification patterns were evaluated. The m⁷G score was established based on principal component analysis (PCA) to quantify m⁷G modification patterns and evaluate the TME cell infiltrating characteristics of individual tumors. Further, correlation analyses of m7Gscore with response to chemotherapy and immunotherapy were performed.

Results: We identified three distinct m⁷G modification patterns with the biological pathway enrichment and TME cell infiltrating characteristics corresponded to immune-desert, immune-inflamed and immune-excluded phenotype, respectively. We further demonstrated the m⁷Gscore could predict the TME infiltrating characteristics, tumor mutation burden (TMB), response to immunotherapy and chemotherapy, as well as prognosis of individual tumors. High m⁷Gscore was associated with increased component of immune cell infiltration, low TMB and survival advantage, while low m⁷Gscore was linked to decreased immune cell infiltration and increased TMB. Additionally, patients with lower m⁷Gscore demonstrated significant therapeutic advantages.

Conclusion: This study demonstrated the regulatory mechanisms of m⁷G modification on TME formation and regulation of lung adenocarcinoma. Identification of individual tumor m⁷G modification patterns will contribute to the understanding of TME characterization and guiding more effective immunotherapy strategies.

Introduction

Lung cancer is currently the second most frequently diagnosed cancer, and the leading cause of cancer death in the world, accounting for 11.4% of all new cancer diagnoses and causing 18.7% of cancer-related deaths (Sung et al., 2021). Lung adenocarcinoma (LUAD) is the most common subtype, accounting for approximately 40% of all lung cancers. Current clinical treatments of LUAD include surgery, chemotherapy, radiotherapy, immunotherapy and molecularly targeted agents. Unfortunately, the poor prognosis of patients highlights the urgent need for the development of novel and more specific therapeutic targets for the treatment of LUAD. Accumulated evidence suggests that aberrant epigenetic modifications, especially RNA methylation, play critical roles in cancer development and progression (Berdasco and Esteller, 2010).

N7-methylguanosine (m⁷G) is one of the most prevalent modifications occurring in transfer RNA (tRNA) (Gauss et al., 1979), ribosomal RNA (rRNA) (Motorin and Helm, 2011) and messenger RNA (mRNA) 5′cap (Capping, 1976), that plays an essential role in regulating RNA processing, exporting, metabolism and function. Meanwhile, as a universally conserved modified nucleosides, m⁷G was found widely among eubacteria, eukaryotes (Juhling et al., 2009), and a few archaea (Edmonds et al., 1991). Notably, recent researches have begun to demonstrate the existence of m⁷G modification within internal mRNAs in higher eukaryotes (Chu et al., 2018; Malbec et al., 2019), and identified the distribution features of the internal mRNA m⁷G using both methylated RNA immunoprecipitation sequencing and chemical modification-assisted BS-seq methods (Zhang et al., 2019). Accumulative evidences have unraveled part of the regulatory mechanisms of m⁷G modification within mRNA, for example, METTL1-WDR4 complex was demonstrated to act as the m⁷G methyltransferases for mRNAs and mediate their formation (Malbec et al., 2019). Thus, m⁷G has become a novel biological marker with critical regulatory roles with the rapid advancement of sequencing technology. To further investigate the m⁷G modification patterns in LUAD and elucidate their impact on tumor progression, we retrieved 24 potential m⁷G modification-related regulators by considering the previous research (Tomikawa, 2018) and exploring the Molecular Signatures Database (www.broadinstitute.org/gsea/msigdb/annotate.jsp).

Tumor progression depends not just on the genetic and epigenetic heterogeneity of tumor cells, but also on the tumor microenvironment (TME), a complex environment containing tumor cells, interstitial cells [e.g., fibroblasts, endothelial cells, tumor-associated macrophages (TAMs)], distant recruited cells [e.g., infiltrating immune cells and bone marrow-derived cells (BMDCs)], and non-cellular elements (e.g., extracellular matrix, cytokines, chemokines and new blood vessels) (Witz and Levy-Nissenbaum, 2006). Complex interactions between tumor and TME are critically involved in multiple malignant biological behaviors such as stimulating cells proliferation and angiogenesis, suppressing apoptosis, as well as inducing immune tolerance (Mantovani et al., 2008). Growing evidence reveals the pivotal role of TME in tumorigenesis, tumor progression and immune evasion (Quail and Joyce, 2013). In particular, TME significantly correlate with response to immune checkpoint blockade (ICB) therapy, and the evaluation of TME cell infiltrating characterization is crucial for the development of novel immunotherapeutic strategies (Ali et al., 2016). Thus, comprehensive analyses of the TME landscape facilitate the identification of distinct tumor immunophenotypes, and contribute to developing biomarkers of the response to immunotherapy and discovering novel targets for immunotherapy.

Compelling evidence has revealed the pivotal role of RNA methylation in shaping and impacting the TME (especially immune cells infiltrating) (Cao and Yan, 2020; Zhang et al., 2020). The m6A RNA methylation has been reported mediating the biological behavior of tumor cells and tumor-infiltrating immune cells by regulating RNA splicing, translation, initiation degradation and nuclear export (Pinello et al., 2018; Han et al., 2019; Wang et al., 2019). And m1A methylation modification has been confirmed to participate in the regulation of TME complexity and diversity based on immune cell infiltration (Gao et al., 2021; Liu et al., 2021). Chen et al. (2022) reported that METTL1-mediated m⁷G modification altered the immune characterization and dynamic interplay between tumor cells and surrounding stromal compartment. Galloway et al. (2021) reported that m⁷G cap methyltransferase RNMT increased translational capacity during T cell activation by coordinating mRNA processing.

However, due to methodological limitations, these studies have necessarily focused on one or several m⁷G regulators and cell types, while the antitumor effect of RNA modification is a highly coordinated process that regulated by numerous tumor suppressor factors. Therefore, a comprehensive analysis of the correlation between TME cell infiltration characterizations and multiple m⁷G regulators will further elucidated the mechanisms of m⁷G modification regulating the TME characterization and provide novel support for more effective immunotherapy.

In this study, we extracted and integrated the genomic data of 817 LUAD samples from the public databases to comprehensively analyze the m⁷G modification patterns, as well as explored the TME cell infiltrating characteristics under different patterns. We identified three m⁷G modification patterns which corresponded to immune-desert, immune-inflamed and immune-excluded phenotype, respectively, revealing that m⁷G modification played an indispensable role in shaping the TME characterization. Furthermore, a set of scoring system was established to quantify the individual tumor m⁷G modification patterns in LUAD patients. Improving the m⁷G modification patterns by targeting m⁷G regulators or m⁷G-related genes may alter TME cell-infiltrating characteristics, that may contribute to the development of novel immunotherapy target or optimization of combination therapy strategies.

Materials and methods

Data extraction and preprocessing

The flowchart of this study was shown in Figure 1. Gene expression profiles and matching clinical annotation were retrieved from The Cancer Genome Atlas (TCGA, https://www.cancer.gov/) and Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database. Patients with incomplete survival information were excluded from further analyses. Finally, 817 LUAD patients were included from four datasets (GSE50081, GSE37745, GSE30219, and TCGA-LUAD) for further analyses in this study. For consistency, Fragments Per Kilobase of transcript per Million reads sequenced (FPKM) values (TCGA-RNA sequencing data) were converted into transcripts per kilobase million (TPM) values. Batch effects in this cohort were removed using the ComBat algorithm (sva R package). The detailed information for each dataset was summarized in Supplementary Table S1. The somatic mutation profiles and Copy Number Variation (CNV) data of LUAD samples were obtained from TCGA database. These data were analyzed using the R (version 4.1.1) and R Bioconductor packages.

FIGURE 1

FIGURE 1. Overview of the study design and analytical flow. LUAD, lung adenocarcinoma; m⁷G, N7-methylguanosine; TME, tumor microenvironment; DEG, differentially expressed gene; PCA, principal component analysis.

Identification and unsupervised clustering for 7-methylguanosine regulators

Twenty four regulators related to m⁷G modification were identified from the published literature (Tomikawa, 2018) and three Gene Set (M26066, M26714, M18244) from Gene Set Enrichment Analysis (GSEA, http://www.gsea-msigdb.org/gsea/index.jsp), including METTL1, WDR4, NSUN2, DCP2, DCPS, NUDT3, NUDT4, NUDT10, NUDT11, NUDT16, AGO2, CYFIP1, EIF3D, EIF4E, EIF4E2, EIF4E3, EIF4G3, GEMIN5, LARP1, NCBP1, NCBP2, IFIT5, LSM1, and SNUPN. To identify distinct m⁷G modification patterns and categorize patients into subgroups, we performed unsupervised consensus clustering analysis (K-Means algorithm, Euclidean distance measure) using “ConsensusClusterPlus” R package (Hartigan and Wong, 1979; Wilkerson and Hayes, 2010), and conduct 1,000 repetitions to ensure the stability of classification.

Biological pathway enrichment analysis

To explore the difference in biological pathway between distinct m⁷G modification patterns, we used “GSVA” R package to conduct Gene set variation analysis (GSVA), which is a non-parametric, unsupervised method for estimating variation of gene set enrichment through the samples of an expression dataset (Hänzelmann et al., 2013). The gene sets of “c2. cp.kegg.v7.5.1,” “c5. go.bp.v7.5.1,” “c5. go.cc.v7.5.1,” and “c5. go.mf.v7.5.1” were downloaded from GSEA database for GSVA analysis, and adjusted p value less than 0.05 was considered statistically significant.

Comprehensive analysis of the tumor microenvironment characterization

The Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data (ESTIMATE) algorithm (Yoshihara et al., 2013) was performed to infer the fraction of stromal (defined as stromal score) and immune cells (defined as immune score) in each LUAD samples using “ESTIMATE” R package (v1.1.0, https://bioinformatics.mdanderson.org/estimate/rpackage.html). The ESTIMATE score is the sum of the immune score and the stromal score, and represents the comprehensive proportion of both components in the TME. Tumor purity was defined as the percentage of malignant cells in a solid tumor sample. The single-sample gene-set enrichment analysis (ssGSEA) algorithm (Barbie et al., 2009) was performed to quantify the immune cell infiltration using “GSVA” R package. The immune cell population were determined with reference to the study of Zhang (Zhang et al., 2020). Specific marker gene sets for each immune cell type (Supplementary Table S2) were derived from the published literatures (Barbie et al., 2009; Charoentong et al., 2017), which contained both innate immune cells (eg, eosinophils, neutrophils, macrophages) and adaptive immune cells (e.g., CD4⁺ T cell, CD8⁺ T cell, regulatory T cell). The gene set of immune-checkpoints was referred to Mariathasan et al. (2018).

N7-methylguanosine-related genes identification and N7-methylguanosine gene signature construction

Differentially expressed genes (DEGs) among distinct m⁷G modification patterns were determined by empirical Bayesian method using the “limma” R package (Smyth, 2004), and the selection criteria was set as adjusted p value <0.001. The intersections of distinct DEGs were defined as m⁷G-related DEGs. Based on m⁷G-related DEGs, all patients were classified into several subgroups for further analysis by performing unsupervised consensus clustering analysis (K-Means algorithm, Euclidean distance measure). This procedure was repeated 1,000 times to ensure the stability of classification. The “clusterProfiler” R package was used to perform GO enrichment analysis for the m⁷G-related DEGs, and significant enrichment pathways (adjust p value <0.05 and Q value <0.05) were displayed in the barplot.

Furthermore, to quantify the m⁷G modification patterns of individual tumors, a m⁷G gene signature (named as m⁷Gscore) was conducted according to the following steps. Firstly, univariate Cox regression analysis was performed to identified significant (p < 0.001) prognosis m⁷G-related DEGs. Following principal component analysis (PCA), both principal components 1 (PC1) and principal components 2 (PC2) were extracted to act as the gene signature score. Finally, we applied a method similar to gene expression grade index (GGI) (Sotiriou et al., 2006; Zhang et al., 2020) to define the m⁷Gscore of each patient: $m^{\mathit{7}}\mathrm{G}score=\sum \left({PC\mathit{1}}_i+{PC\mathit{2}}_i\right)$, where i is the expression of m⁷G-related DEGs.

Prediction of immunotherapy and chemotherapy response

We investigated the predictive capacity of m⁷Gscore in responding immunotherapy and four common first-line chemotherapy drugs (cisplatin, paclitaxel, docetaxel, gemcitabine) (Ettinger et al., 2021) for LUAD. The clinical response to immunotherapy was inferred by the Tumor Immune Dysfunction and Exclusion (TIDE, http://tide.dfci.harvard.edu/), an algorithm to simulate two primary mechanisms of tumor immune evasion: the induction of T cell dysfunction in tumors with high cytotoxic T lymphocytes (CTLs) infiltration and the prevention of T cell infiltration in tumors with low CTL (Jiang et al., 2018). Generally, a lower T cell dysfunction signature score predicts a better response to immunotherapy. The 50% inhibiting concentration (IC50) values of the four chemotherapy drugs were predicted using the pRRophetic algorithm (Geeleher et al., 2014) and the value was normally transformed.

Statistical analysis

The Student’s t test was used to compare the differences between two groups, and one-way ANOVA and Kruskal-Wallis tests were used to compare the differences among multiple groups. Spearman correlation coefficient was used for correlation analysis. Survival curves were constructed using the Kaplan-Meier method, and the log-rank test was used to identify the significance of differences. Univariate Cox regression analysis was performed to estimate the hazard ratios (HR) and 95% confidence intervals (CI). Multivariate Cox regression analysis was employed to identify independent prognostic factors, and only patients with complete clinical information were included in final multivariate analysis. The waterfall plots of a mutational landscape in TCGA-LUAD cohort were generated using “maftools” R package (Mayakonda et al., 2018). The copy number variation (CNV) landscape of m⁷G regulators in 23 pairs of chromosomes was visualized using “RCircos” R package (Zhang et al., 2013). All statistical analyses were performed using R software (version 4.1.1), and a p value < 0.05 was considered statistically significant.

Results

Analysis of genetic variation, expression and prognostic value of m⁷G regulators in lung adenocarcinoma

Genetic alteration is a critical factor influencing gene expression and function, we firstly explored the incidence of somatic mutations and copy number variations (CNV) of 24 m⁷G regulators in LUAD. As shown in Figure 2A, 80 of 561, (14.26%) samples experienced mutations of m⁷G regulators. Among them, the mutation frequency of EIF4G3 was the highest (3%), followed by LARP1 (2%), while nearly half of the regulators did not show any mutations. The summary of CNV showed that AGO2, NSUN2, METTL1, NCBP2 and NUDT3 were the top five regulators with highest CNV frequency, and amplification variations obviously higher than deletion (Figure 2B). The location of m⁷G regulators CNV on chromosomes was displayed in Figure 2C. To explore whether the genetic variation affect the expression level of m⁷G regulators in LUAD patients, the mRNA expression levels of these regulators were further analyzed between normal and LUAD samples, which indicated that most of these regulators were dysregulated in LUAD samples (Figure 2D). Moreover, univariate (Supplementary Figure S1A) and multivariate (Supplementary Figure S1B) Cox regression analyses further identified four independent poor prognostic factors (NUDT11, NUDT4, LARP1 and NCBP2) and two protective factors (EIF4E3, NUDT10) for LUAD patients. A complex regulatory network depicted the regulatory relationship (Supplementary Table S3) among m⁷G regulators and their prognostic significance for LUAD patients (Figure 2E). We found that most of regulators displayed a remarkably positive correlation in expression, whereas a few negative correlations among EIF4G3 and AGO2/NSUN2/WDR4/METTL1/LSM1, NUDT16 and LSM1/METTL1. In brief, the genetic alteration and expression level of m⁷G regulators are highly heterogeneous and significantly correlated with prognosis, indicating that the expressional alteration of m⁷G regulators played a crucial role in the LUAD occurrence and development.

FIGURE 2

FIGURE 2. Landscape of genetic and expression variation of m⁷G regulators in lung adenocarcinoma (LUAD). (A) The mutation frequency of 24 m⁷G regulators in 561 patients with LUAD from TCGA-LUAD cohort, and each column represents individual patients. The top barplot depicts tumor mutation burden and mutation frequency in each regulator is given on the right. The right barplot depicts the proportion of each variant type. The stacked barplot below depicts fraction of conversions in each sample. (B) The copy number variation (CNV) alteration frequency of m⁷G regulators in TCGA-LUAD cohort. The height of each column represents the alteration frequency (Red dot: the amplification frequency; green dot: the deletion frequency). (C) The location of CNV alteration of m⁷G regulators on 23 chromosomes using TCGA-LUAD cohort. (D) The expression of m⁷G regulators between normal tissues and LUAD tissues (Red: tumor; blue: normal). The horizontal line indicates the median, the lower and upper boundaries of the boxes the interquartile range, and the dots the outliers. Asterisks indicate statistical significance, *p < 0.05, **p < 0.01, ***p < 0.001. (E) The interaction between m⁷G regulators in LUAD. Each regulator is represented as a circle, where the size of the circle represents the effect of each regulator on the prognosis, and p values were calculated by Log-rank test (*p < 0.05, **p < 0.01, ***p < 0.001).

Identification of m⁷G methylation modification patterns mediated by 24 regulators

Based on the expression of 24 m⁷G regulators, three distinct m⁷G modification patterns were identified using unsupervised consensus clustering with optimal clustering stability, including 230 patients in pattern A, 262 patients in pattern B and 325 patients in pattern C (Supplementary Figures S2A–F). We named these patterns “m⁷Gcluster” A-C. Principal component analysis (PCA) confirmed the significant distinction existed on the m⁷G regulators expression among these m⁷G modification patterns (Supplementary Figures S2G). The expression patterns of m⁷G regulators and comparison of baseline clinicopathological characteristics in the three m⁷Gclusters were shown by the heatmap (Supplementary Figures S2H). m⁷Gcluster A exhibited high expression of METTL1, WDR4, NSUN2, DCPS, NUDT3, AGO2, EIF4E, LARP1, NCBP1, NCBP2, EIF4G3, LSM1; m⁷Gcluster B was characterized by decreased expression in almost regulators, except for METTL1, LSM1 and SNUPN; m⁷Gcluster C exhibited high expression of DCP2, NUDT16, CYFIP1, EIF4E3, and IFIT5.

Biological behaviors and the tumor microenvironment characterization in distinct m⁷G modification patterns

To identify the biological significance of distinct m⁷G modification patterns, we performed GSVA enrichment analysis (Supplementary Table S4). Specifically, m⁷Gcluster A was enriched in common oncogenic signaling pathways (e.g., mTOR, Notch and NSCLC signaling pathway), while lacked immune activation process (e.g., cytokine-cytokine receptor interaction, antigen processing and presentation), leading to the activation of abnormal biological characteristics including cell cycle, basal transcription factors, spliceosome, etc. (Figure 3A). On the contrary, in m⁷Gcluster B, various oncogenic signaling pathways such as mTOR, Notch and cancer associated pathways were strikingly suppressed, and immune activation process were significantly activated (Figures 3A,B). As expected, patients in m⁷Gcluster A showed the worst prognosis, whereas patients in m⁷Gcluster B had the best prognosis (Figure 3C). Significantly, m⁷Gcluster C was significantly enriched in innate immune activation process (such as Fc epsilon RI, Nod like receptor, Toll like receptor, and Fc gamma R-mediated phagocytosis signaling pathways) (Figure 3B), while patients with this m⁷G modification pattern did not show a matching survival advantage. We speculated that this result may be related to high-level enrichment of stromal interactions pathways (such as ERBB, Wnt, TGF beta, Adherens Junction, Focal Adhesion and Regulation of Actin Cytoskeleton pathways), as well as B cell receptor (BCR) and T cell receptor (TCR) which could upregulate B and T cells activation threshold.

FIGURE 3

FIGURE 3. Biological characteristics and tumor microenvironment characterization of each m⁷G modification pattern. (A,B) The heatmap visualizes the enrichment of biological processes using GSVA analysis in distinct m⁷G modification patterns; (A) m⁷Gcluster A vs. m⁷Gcluster B; (B) m⁷Gcluster B vs. m⁷Gcluster C. Activated pathways are colored red and inhibited pathways are colored blue. The LUAD cohorts are used as sample annotations. (C) Survival analyses for the three m⁷G clusters based on 817 patients with LUAD from four cohorts (TCGA-LUAD, GSE30219, GSE50081, GSE37745). (D) The abundance of 23 TME infiltrating cells in the three m⁷G clusters. The horizontal line indicates the median, the lower and upper boundaries of the boxes the interquartile range, and the dots the outliers. Asterisks indicate statistical significance, *p < 0.05, **p < 0.01, ***p < 0.001. (E–H) Violin plots show differences in the (E) immune score, (F) stromal score, (G) ESTIMATE score and (H) tumor purity between distinct m⁷G clusters.

To further explore the TME characterization, we compare the difference in immune cell infiltration among distinct m⁷G modification patterns (Supplementary Table S5; Figure 3D). Increased adaptive (e.g., B cell, CD8⁺ T cell, T helper cell) and innate (e.g., macrophage, NK cell, dendritic cell) immune cell infiltration was exhibited in m⁷Gcluster B and C as compared to m⁷Gcluster A. To our surprise, m⁷Gcluster C was remarkably rich in innate immune cell infiltration including dendritic cell, eosinophil, gamma delta T cell, macrophage, mast cell, natural killer cell, and cell types associated with immune suppression such as MDSC and regulatory T cell. Additionally, the ESTIMATE algorithm revealed the lowest level of stromal and immune cell infiltration in m⁷Gcluster A, while revealed the highest level in m⁷Gcluster C (Figures 3E–H). Based on the above results and previous research (Chen and Mellman, 2017), We could summarize that the three m⁷G modification patterns corresponded to three different immunophenotypes. Strikingly, m⁷Gcluster A was identified as immune-desert phenotype, characterized by suppressed immune-related pathways and deficient immune cell infiltration. m⁷Gcluster B was identified as immune-inflamed phenotype, characterized by immune activation and high level of adaptive immune cells infiltration. More importantly, we found that the TME characterization of m⁷Gcluster C was consistent with the immune-excluded phenotype described by Chen (Chen and Mellman, 2017), which was characterized by enhanced tumor stroma activity and abundant innate immune cells trapped in surrounding tumor cell nests.

Identification of m⁷G modification-related genomic subtypes and transcriptome characterization

In order to further explore the potential biological significance of distinct m⁷G modification patterns, we distinguished 2071 m⁷G-related differentially expressed genes (DEGs) (Supplementary Figure S3A, Supplementary Table S6). Subsequently GO enrichment analysis of m⁷G-related DEGs revealed significant enrichment of biological processes related to m⁷G modification and immune system (Supplementary Table S7, Supplementary Figure S3B), which confirmed that m⁷G modification played a critical role in immune modulation of the TME. To clarify the potential mechanisms, we further performed unsupervised consensus clustering based on the 401 prognosis m⁷G-related DEGs (Supplementary Table S8) and classify the entire LUAD cohort into three main m⁷G-related genomic subtypes (named as m⁷G gene cluster A-C), including 181 patients in cluster A, 268 patients in cluster B and 368 patients in cluster C (Supplementary Figures S4A–F). Signature genes expression level and baseline clinicopathological characteristics for the different clusters were displayed in Figure 4A. We found that m⁷G gene cluster A and B showed opposite gene expression patterns, and patients with alive status or clinical stage I-II were mainly concentrated in the cluster B. Kaplan-Meier analysis indicated that gene cluster A exhibited poorer prognosis, while gene cluster B exhibited better prognosis, and gene cluster C exhibited intermediate prognosis (Figure 4B). Moreover, prominent differences were observed in the expression of m⁷G regulators among the three gene clusters (Figure 4C), which demonstrated again that m⁷G ​​modification modulate the genomic phenotype of LUAD patients.

FIGURE 4

FIGURE 4. Clinicopathologic characteristics, genomic profiling and tumor microenvironment characteristics among distinct m⁷G modification-related genomic subtypes. (A) The heatmap visualizes the gene expression levels across the whole genome and comparison of baseline clinicopathological characteristics in each sample. Blue represents low expression and red represent high expression. (B) Survival analyses for the three gene clusters based on 817 patients with LUAD from four cohorts (TCGA-LUAD, GSE30219, GSE50081, GSE37745). (C) Difference in the m⁷G regulators expression among three gene clusters. (D) Difference in the abundance of 23 TME infiltrating cells among three m⁷G gene clusters. (E) Difference in the immune-checkpoint related gene expression among three gene clusters. The horizontal line indicates the median, the lower and upper boundaries of the boxes the interquartile range, and the dots the outliers. Asterisks indicate statistical significance, *p < 0.05, **p < 0.01, ***p < 0.001.

To explore the relationship between m⁷G-related genomic features and the tumor immune microenvironment, we examined the immune cell infiltrating characteristics and immune-checkpoint related gene (ICG) expression in three gene clusters. As shown in Figures 4D,E, gene cluster A was characterized by low levels of immune cell infiltration and upregulated ICG expression. Conversely, gene cluster B was characterized by high levels of immune cell infiltration and downregulated ICB expression. Remarkably, gene cluster C exhibited high immune cell infiltration level, while high immune checkpoint related mRNAs expression, which might be relevant to the poor prognosis.

Construction and evaluation of m⁷Gscore, one m⁷G modification quantification system

Regrettably, the above conclusions were generated based on group-level analyses, the characteristics of m⁷G modification in individual patients were limited. Considering the heterogeneity of m⁷G modification, we constructed a scoring system (named as m⁷Gscore) based on prognosis m⁷G-related DEGs to accurately quantify and predict the individual tumors m⁷G modification pattern (Supplementary Table S9). Patients were classified into high m⁷Gscore group (n = 348) and low m⁷Gscore group (n = 469) using the optimal cut-off value 2.984 identified by the “surv_cutpoint” function from the “survminer” R package. An alluvial diagram was generated to depict the distribution transitions of individual patients among the m⁷Gclusters, m⁷G gene clusters and m⁷Gscore groups (Figure 5A). Nonparametric Kruskal-Wallis (K-W) test was performed to reveal the difference in m⁷Gscore among distinct m⁷Gclusters and m⁷G gene clusters. We noticed that m⁷Gcluster A presented the lowest median score, whereas m⁷Gcluster C presented the highest (Figure 5B), which suggested that low m⁷Gscore might be closely associated with deficient immune cell infiltration while high m⁷Gscore might be linked to stromal activation. Additionally, m⁷G gene cluster B had significantly higher m⁷Gscore than the other two clusters, while m⁷G gene cluster A showed the lowest median score (Figure 5C). Moreover, K-M analysis indicated that high m⁷Gscore conferred a significant survival benefit (p < 0.001, Figure 5D). Subsequent subgroup analyses (Supplementary Figures S5E–J) showed that the prognostic value of m⁷Gscore remained statistically significant for each subgroup based on gender (male, female), age (≤65, >65), and clinical stage (I-II, III-IV). m⁷Gscore also showed the prognostic value in different datasets (Supplementary Figures S5A–D). Additionally, As shown in Supplementary Figures S5H–L, patients with high m⁷Gscore exhibited a significantly higher percentage of alive status (68%), and patients who died had remarkably lower m⁷Gscore (p < 0.001). Univariate (Figure 5E) and multivariate (Figure 5F) Cox analyses confirmed m⁷Gscore as a robust and independent prognostic biomarker for evaluating patient outcomes [HR: 0.495 (0.393–0.623), p < 0.001].

FIGURE 5

FIGURE 5. Construction of the m⁷G signature and correlation of m⁷Gscore with clinicopathological features. (A) Alluvial diagram shows the changes of m⁷G clusters, gene clusters, m⁷Gscore and survival state. (B) Differences in m⁷Gscore among three m⁷G clusters (p < 0.001, Kruskal-Wallis test). (C) Differences in m6Ascore among three gene clusters (p < 0.001, Kruskal-Wallis test). (D) Survival analyses for low (469 cases) and high (348 cases) m⁷Gscore groups in the four cohorts using Kaplan-Meier curves (p < 0.0001, Log-rank test). (E,F) Univariate (E) and multivariate (F) analyses for m⁷Gscore using the Cox regression model.

Furthermore, we assessed the relationship between m⁷Gscore and the TME cell infiltration to explore whether m⁷G modification quantification system can reflect the TME heterogeneity. The ESTIMATE algorithm indicated that high m⁷Gscores were significantly associated with enhanced levels of immune and stromal cell infiltration as well as low tumor purity (p < 0.001; Figures 6A–D). Correlation analysis of immune cell infiltration and m⁷Gscore (Figure 6E) indicated that m⁷Gscore was significantly positively correlated with most types of innate immune cells (e.g., eosinophil, dendritic cell, mast cell) and B cell, whereas negatively correlated with active CD4 T cell and CD56dim NK cell. In addition, m⁷Gscore was negatively correlated with immunosuppression-related ICGs (including IDO1, PDCD1, LAG3, TNFRSF9), whereas positively correlated with CD28⁻CD80/86 (Figure 6F) which provides co-stimulatory signals for T-cell activation (Esensten et al., 2016; Ma et al., 2021a). Considering the above results, high m⁷Gscore indicated increased immune and stromal cell infiltration, while low m⁷Gscore was correlated with decreased immune cell infiltration and great immunosuppression. m⁷Gscore could reflect the m⁷G modification pattern of individual LUAD patients to further evaluate the TME characterization.

FIGURE 6

FIGURE 6. Correlations between m⁷G signature and the tumor microenvironment characteristics. (A–D) Violin plots show differences in the (A) immune score, (B) stromal score, (C) estimate score and (D) tumor purity between low and high m⁷Gscore groups. (E) Correlations between m⁷Gscore and the abundance of 23 TME infiltrating cells. (F) Correlations between m⁷Gscore and immune-checkpoint related gene expression. Positive correlation was marked with red and negative correlation with blue. The circle color represents Spearman coefficient value, the size of circle is inversely proportional to the p-value, and the asterisk stands for p < 0.05.

Analysis of tumor somatic mutation in patients with different m⁷Gscores

Accumulated evidence demonstrated that tumor mutation burden (TMB) is an emerging biomarker of response to immune checkpoint blockade (ICB) therapy (Yarchoan et al., 2017; Ready et al., 2019), we analyzed the landscape of tumor somatic mutation among patients with different m⁷Gscore to indirectly reflect the immunotherapeutic outcomes in TCGA-LUAD cohort. The waterfall plots suggested that low m⁷Gscore group exhibited more extensive tumor somatic mutation than the high m⁷Gscore group, with the rate of the 20th most significant mutated gene 22% versus 9% (Figures 7A,B). As shown in Figures 7C,D, the m⁷Gscore exhibited a significant negative correlation to TMB. Subsequent Kaplan-Meier survival analysis indicated that patients with a high TMB level had better OS than those with a low TMB level in low m⁷Gscore group (Figure 7E). The above results demonstrated that m⁷Gscore could effectively reflect the TMB level of LUAD, which indirectly indicated the values of distinct m⁷G modification patterns in predicting ICB therapy outcomes.

FIGURE 7

FIGURE 7. The correlation of m⁷Gscore and the tumor somatic mutation. (A,B) The waterfall plot of tumor somatic mutation in LUAD patients with low m⁷Gscore (A) and high m⁷Gscore (B), each column represents individual patients. The top barplot depicts tumor mutation burden (TMB) and mutation frequency in each gene is given on the right. The right barplot depicts the proportion of each variant type. The stacked barplot below depicts fraction of conversions in each sample. (C,D) The relationship between the m⁷Gscore and TMB. (E) The Kaplan–Meier curves of the OS of subgroup patients stratified by the m⁷Gscore and TMB.

Prediction of immunotherapy and chemotherapy response

Based on the TIDE (Tumor Immune Dysfunction and Exclusion) algorithm, patients in TCGA-LUAD cohort were classified into insensitive and sensitive groups. As shown in Figure 8A, m⁷Gscores were significantly higher in the insensitive group than in the sensitive group (p < 0.001). Similarly, m⁷Gscores were significantly positively correlated with T cell dysfunction scores in GSE30219 (R = 0.39, p < 0.001; Figure 8B) and GSE37745 (R = 0.40, p < 0.001; Figure 8C) cohorts, which suggested that patients with higher m⁷Gscores had poor immunotherapy response rates. Moreover, there were marked increases in the IC50 to cisplatin (Figure 8D, p < 0.001), paclitaxel (Figure 8E, p < 0.001), docetaxel (Figure 8F, p < 0.001) and gemcitabine (Figure 8G, p < 0.01) in high m⁷Gscore group, which indicated the poor efficacy to these drugs in patients with high m⁷Gscores compared to patients with low m⁷Gscores. Together, m⁷Gscore could effectively predict the response to chemotherapy and immunotherapy for LUAD patients.

FIGURE 8

FIGURE 8. The correlation of m⁷Gscore and the therapeutic response. (A) The difference in m⁷Gscore between sensitive and insensitive groups divided by TIDE algorithm in TCGA-LUAD cohort. (B,C) The correlation of m⁷Gscore and the T cell dysfunction signature in GSE30219 (B) and GSE37745 (C) cohorts. (D–G) Differences in the IC50 of cisplatin (D), paclitaxel (E), docetaxel (F), gemcitabine (G) between low and high m⁷Gscore groups. IC50, 50% inhibitory concentration, which negatively correlated with drug responsiveness.

Discussion

With the rapid advancement of deep sequencing and large-scale profiling (Enroth et al., 2019; El Allali et al., 2021), accumulating evidence has demonstrated that m⁷G modification is critical for maintaining the physiological conditions of cells and organisms (Pei and Shuman, 2002; Lindstrom et al., 2003; Haag et al., 2015), while its aberrant distribution is closely related to tumor development and progression (Ma et al., 2021b). Moreover, recent studies have also confirmed that m⁷G may affect the distribution and function of immune cells (Zhang et al., 2021; Gao et al., 2022), such as T cells (Galloway et al., 2021). As most studies have focused on single regulator or single TME cell type, a comprehensive recognition of TME infiltration characterizations mediated by multiple m⁷G regulators is still lacking. Exploring the role of different m⁷G modification patterns in the TME cell infiltration will help to enhance our understanding of the TME antitumor immune response and guide novel immunotherapy strategies.

In this study, three m⁷G methylation modification patterns were identified based on 24 potential m⁷G regulators, which had significantly distinct TME cell infiltration characterizations. The m⁷G cluster A was characterized by suppressed immune-related functions and deficient immune cells infiltration, consistent with immune-desert phenotype; cluster B was characterized by immune activation and high level of adaptive immune cells infiltration, consistent with immune-inflamed phenotype; cluster C was characterized by enhanced tumor stroma activity and abundant innate immune cells, consistent with the immune-excluded phenotype. As mentioned in previous literatures (Gajewski et al., 2013; Joyce and Fearon, 2015; Chen and Mellman, 2017), the immune-inflamed tumors can demonstrate infiltration of large number of immune cells, especially T, B and monocytic cells, in the tumor parenchyma; the immune-excluded tumors are also characterized by the presence of abundant immune cells, but the immune cells are retained in the stroma surrounding nests of tumor cells rather than penetrate the parenchyma; the immune-desert tumors are associated with the immunological ignorance and paucity of immune cells in either the tumor parenchyma or the stroma. Significant enrichment of stromal interactions pathways in m⁷G cluster C and the characteristics of TME cell infiltration in each cluster corroborate the accuracy of our immunophenotypic classifications for distinct m⁷G modification patterns. Not surprisingly, m⁷G cluster C exhibits activated innate immunity but no matching survival advantage.

Furthermore, differentially expressed genes among the three modification patterns (named as m⁷G-related DEGs) were identified and demonstrated to be significantly associated with immune-related biological pathways. Three m⁷G-related genomic subtypes were identified based on 401 prognostic m⁷G-related DEGs, which were also significantly related to immune cell infiltration and activation. This further demonstrated the crucial role of m⁷G modification in modulating the TME landscape. Given the individual heterogeneity of m⁷G modification, it was critical to quantify m⁷G modification patterns in individual tumors. Therefore, we developed a novel scoring system (m⁷Gscore) to assess the m⁷G modification pattern of individual LUAD patients. The patients with immune-excluded and immune-inflamed tumor presented a higher m⁷Gscore, while the patients with immune-desert tumor presented a lower m⁷Gscore. Also, high m⁷Gscores were significantly associated with enhanced levels of immune and stromal cells infiltration. These results suggested m⁷Gscore was a reliable tool for assessing individual tumor m⁷G modification patterns, which could further indicate the immune phenotype in tumor environment. Additionally, m⁷Gscore was proved to be an independent prognostic factor, with lower m⁷Gscores indicating poorer prognosis.

Our results demonstrated the significantly negative correlations of m⁷Gscore with the expression of IDO1, PDCD1, and LAG3, which have been considered as important targets for cancer immunotherapy (Cyriac and Gandhi, 2018; Ruffo et al., 2019; Tang et al., 2021). There was also a markedly negative correlation between m⁷Gscore and tumor mutation burden (TMB). Growing evidence demonstrated that patients with high TMB had a greater clinical response to anti-PD-1/PD-L1 immunotherapy (Yarchoan et al., 2017; Ready et al., 2019). The Tumor Immune Dysfunction and Exclusion algorithm further showed that higher m⁷Gscore was associated with T cells dysfunction and exclusion, which directly reflect the efficacy to T cell-based immunotherapy (Jiang et al., 2018). Thus, the above results fully demonstrated that individual m⁷G modification pattern could be an effective indicator that estimate the responsiveness to Immune checkpoint blockade (ICB) therapy. Many reports have revealed the interaction between chemotherapy and immunotherapy, and the differences in immune and stromal cell infiltration in TME jointly affect resistance to chemotherapy (Wang et al., 2016; Zhu et al., 2021). Fu et al. (2018) have reported the immunotypes could predict the efficacy of patients to adjuvant chemotherapy. In this study, patients with higher m⁷Gscore exhibited poor efficacy to several first-line chemotherapy drugs (including cisplatin, paclitaxel, docetaxel and gemcitabine) for LUAD, and this might due to lower T cell infiltration and higher stromal cell infiltration.

In a word, m⁷Gscore can act as an effective tool to evaluate the individual m⁷G modification pattern and the corresponding immune phenotypes for LUAD patients, further to guide treatment decisions in clinical practice. m⁷Gscore can also be a potential prognostic biomarker for predicting survival. More importantly, m⁷Gscore may guide the clinicians in predicting the clinical response to ICB therapy and the efficacy of adjuvant chemotherapy. Modifying the m⁷G modification patterns by targeting m⁷G regulators or m⁷G-related genes may improve unfavorable TME cell infiltrating characterization, that may contribute to the development of novel immunotherapy target or optimization of combination therapy strategies. These findings offered new insights for identifying distinct immune phenotypes and developing individualized cancer immunotherapy.

Despite the important strengths of this study, several limitations should be noted. First, due to the limited clinicopathological parameters in public datasets, there would be potential bias when the m⁷Gscore acted as a prognosis biomarker. Second, we did not evaluate the location of immune and stromal cell infiltration in the TME. Thirdly, we could not directly analysis the correlation between m⁷Gscore and LUAD patients’ response to therapy due to the lack of treatment-related information. Finally, our findings were carried out based on a bioinformatics analysis, and further experimental validation is warranted.

Conclusion

In conclusion, this study revealed the non-negligible role of m⁷G modification in TME heterogeneity and complexity. Identifying m⁷G modification patterns helps predict clinical response to ICB therapy and efficacy of adjuvant chemotherapy. We believe that the assessment of individual tumor m⁷G modification pattern will contribute to a more comprehensive understanding of TME characterization and facilitate the development of novel immunotherapy strategies.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

Ethics statement

TCGA and GEO belong to public databases. The patients involved in the database have obtained ethical approval.

Author contributions

TJ and LL conceived and directed the project. SM and JuZ designed the study. SM, JuZ, and MW performed the statistical analysis. JiZ, WW, YX, and RJ provided data collection and visualization. SM drafted the manuscript. TJ, LL, JuZ, and MW helped to write the manuscript. All authors contributed to the article and approved the submitted version.

Acknowledgments

The authors would like to acknowledge the TCGA and GEO database for providing their platforms and those contributors for uploading their valuable datasets.

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.

Supplementary material

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

References

Ali, H. R., Chlon, L., Pharoah, P. D., Markowetz, F., and Caldas, C. (2016). Patterns of immune infiltration in breast cancer and their clinical implications: A gene-expression-based retrospective study. PLoS Med. 13 (12), e1002194. doi:10.1371/journal.pmed.1002194

PubMed Abstract | CrossRef Full Text | Google Scholar

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

Berdasco, M., and Esteller, M. (2010). Aberrant epigenetic landscape in cancer: How cellular identity goes awry. Dev. Cell 19 (5), 698–711. doi:10.1016/j.devcel.2010.10.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Cao, J., and Yan, Q. (2020). Cancer epigenetics, tumor immunity, and immunotherapy. Trends Cancer 6 (7), 580–592. doi:10.1016/j.trecan.2020.02.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Capping, A. J. S. (1976). Capping of eucaryotic mRNAs. Cell 9 (42), 645–653. doi:10.1016/0092-8674(76)90128-8

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. Cell Rep. 18 (1), 248–262. doi:10.1016/j.celrep.2016.12.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, D. S., and Mellman, I. (2017). Elements of cancer immunity and the cancer-immune set point. Nature 541 (7637), 321–330. doi:10.1038/nature21349

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, J., Li, K., Chen, J., Wang, X., Ling, R., Cheng, M., et al. (2022). Aberrant translation regulated by mettl1/wdr4-mediated trna N7-methylguanosine modification drives head and neck squamous cell carcinoma progression. Cancer Commun. 42, 223–244. doi:10.1002/cac2.12273

PubMed Abstract | CrossRef Full Text | Google Scholar

Chu, J. M., Ye, T. T., Ma, C. J., Lan, M. D., Liu, T., Yuan, B. F., et al. (2018). Existence of internal N7-methylguanosine modification in mrna determined by differential enzyme treatment coupled with mass spectrometry analysis. ACS Chem. Biol. 13 (12), 3243–3250. doi:10.1021/acschembio.7b00906

PubMed Abstract | CrossRef Full Text | Google Scholar

Cyriac, G., and Gandhi, L. (2018). Emerging biomarkers for immune checkpoint inhibition in lung cancer. Semin. Cancer Biol. 52 (2), 269–277. doi:10.1016/j.semcancer.2018.05.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Edmonds, C. G., Crain, P. F., Gupta, R., Hashizume, T., Hocart, C. H., Kowalak, J. A., et al. (1991). Posttranscriptional modification of trna in thermophilic archaea (archaebacteria). J. Bacteriol. 173 (10), 3138–3148. doi:10.1128/jb.173.10.3138-3148.1991

PubMed Abstract | CrossRef Full Text | Google Scholar

El Allali, A., Elhamraoui, Z., and Daoud, R. (2021). Machine learning applications in rna modification sites prediction. Comput. Struct. Biotechnol. J. 19, 5510–5524. doi:10.1016/j.csbj.2021.09.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Enroth, C., Poulsen, L. D., Iversen, S., Kirpekar, F., Albrechtsen, A., and Vinther, J. (2019). Detection of internal N7-methylguanosine (M7g) rna modifications by mutational profiling sequencing. Nucleic Acids Res. 47 (20), e126. doi:10.1093/nar/gkz736

PubMed Abstract | CrossRef Full Text | Google Scholar

Esensten, J. H., Helou, Y. A., Chopra, G., Weiss, A., and Bluestone, J. A. (2016). Cd28 costimulation: From mechanism to therapy. Immunity 44 (5), 973–988. doi:10.1016/j.immuni.2016.04.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Ettinger, D. S., Wood, D. E., Aisner, D. L., Akerley, W., Bauman, J. R., Bharat, A., et al. (2021). Nccn guidelines insights: Non-small cell lung cancer, version 2.2021. J. Natl. Compr. Canc. Netw. 19 (3), 254–266. doi:10.6004/jnccn.2021.0013

PubMed Abstract | CrossRef Full Text | Google Scholar

Fu, H., Zhu, Y., Wang, Y., Liu, Z., Zhang, J., Xie, H., et al. (2018). Identification and validation of stromal immunotype predict survival and benefit from adjuvant chemotherapy in patients with muscle-invasive bladder cancer. Clin. Cancer Res. 24 (13), 3069–3078. doi:10.1158/1078-0432.CCR-17-2687

PubMed Abstract | CrossRef Full Text | Google Scholar

Gajewski, T. F., Woo, S. R., Zha, Y., Spaapen, R., Zheng, Y., Corrales, L., et al. (2013). Cancer immunotherapy strategies based on overcoming barriers within the tumor microenvironment. Curr. Opin. Immunol. 25 (2), 268–276. doi:10.1016/j.coi.2013.02.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Galloway, A., Kaskar, A., Ditsova, D., Atrih, A., Yoshikawa, H., Gomez-Moreira, C., et al. (2021). Upregulation of rna cap methyltransferase rnmt drives ribosome biogenesis during T cell activation. Nucleic Acids Res. 49 (12), 6722–6738. doi:10.1093/nar/gkab465

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, Y., Wang, H., Li, H., Ye, X., Xia, Y., Yuan, S., et al. (2021). Integrated analyses of M(1)a regulator-mediated modification patterns in tumor microenvironment-infiltrating immune cells in colon cancer. Oncoimmunology 10 (1), 1936758. doi:10.1080/2162402X.2021.1936758

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, Z., Xu, J., Zhang, Z., Fan, Y., Xue, H., Guo, X., et al. (2022). A comprehensive analysis of Mettl1 to immunity and stemness in pan-cancer. Front. Immunol. 13, 795240. doi:10.3389/fimmu.2022.795240

PubMed Abstract | CrossRef Full Text | Google Scholar

Gauss, D. H., Grüter, F., and Sprinzl, M. S. (1979). Compilation of trna sequences. Nucleic Acids Res. 6 (1), 419–r19. doi:10.1093/nar/6.1.419-a

CrossRef Full Text | Google Scholar

Geeleher, P., Cox, N., and Huang, R. S. (2014). Prrophetic: An R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One 9 (9), e107468. doi:10.1371/journal.pone.0107468

PubMed Abstract | CrossRef Full Text | Google Scholar

Hartigan, J. A., and Wong, M. A. (1979). Algorithm as 136: A K-means clustering algorithm. Appl. Stat. 28 (1), 100–108. doi:10.2307/2346830

CrossRef Full Text | Google Scholar

Haag, S., Kretschmer, J., and Bohnsack, M. T. (2015). Wbscr22/Merm1 is required for late nuclear pre-ribosomal rna processing and mediates N7-methylation of G1639 in human 18s rrna. RNA 21 (2), 180–187. doi:10.1261/rna.047910.114

PubMed Abstract | CrossRef Full Text | Google Scholar

Han, D., Liu, J., Chen, C., Dong, L., Liu, Y., Chang, R., et al. (2019). Anti-tumour immunity controlled through mRNA m6A methylation and YTHDF1 in dendritic cells. Nature 566 (7743), 270–274. doi:10.1038/s41586-019-0916-x

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 Bioinforma. 14 (7), 7. doi:10.1186/1471-2105-14-7

CrossRef Full Text | Google Scholar

Jiang, P., Gu, S., Pan, D., Fu, J., Sahu, A., Hu, X., et al. (2018). Signatures of T Cell dysfunction and exclusion predict cancer immunotherapy response. Nat. Med. 24 (10), 1550–1558. doi:10.1038/s41591-018-0136-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Joyce, J. A., and Fearon, D. T. (2015). Cell exclusion, immune privilege, and the tumor microenvironment. Sci. Sci. 348 (6230), 74–80. doi:10.1126/science.aaa6204

PubMed Abstract | CrossRef Full Text | Google Scholar

Juhling, F., Morl, M., Hartmann, R. K., Sprinzl, M., Stadler, P. F., and Putz, J. (2009). Trnadb 2009: Compilation of trna sequences and trna genes. Nucleic Acids Res. 37, D159–D162. doi:10.1093/nar/gkn772

PubMed Abstract | CrossRef Full Text | Google Scholar

Lindstrom, D. L., Squazzo, S. L., Muster, N., Burckin, T. A., Wachter, K. C., Emigh, C. A., et al. (2003). Dual roles for Spt5 in pre-mrna processing and transcription elongation revealed by identification of spt5-associated proteins. Mol. Cell. Biol. 23 (4), 1368–1378. doi:10.1128/MCB.23.4.1368-1378.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J., Chen, C., Wang, Y., Qian, C., Wei, J., Xing, Y., et al. (2021). Comprehensive of N1-methyladenosine modifications patterns and immunological characteristics in ovarian cancer. Front. Immunol. 12, 746647. doi:10.3389/fimmu.2021.746647

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, J., Han, H., Huang, Y., Yang, C., Zheng, S., Cai, T., et al. (2021). Mettl1/Wdr4-Mediated M(7)G trna modifications and M(7)G codon usage promote mrna translation and lung cancer progression. Mol. Ther. 29 (12), 3422–3435. doi:10.1016/j.ymthe.2021.08.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, Y., Qiu, M., Guo, H., Chen, H., Li, J., Li, X., et al. (2021). Comprehensive analysis of the immune and prognostic implication of Col6a6 in lung adenocarcinoma. Front. Oncol. 11, 633420. doi:10.3389/fonc.2021.633420

PubMed Abstract | CrossRef Full Text | Google Scholar

Malbec, L., Zhang, T., Chen, Y. S., Zhang, Y., Sun, B. F., Shi, B. Y., et al. (2019). Dynamic methylome of internal mrna N(7)-methylguanosine and its regulatory role in translation. Cell Res. 29 (11), 927–941. doi:10.1038/s41422-019-0230-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Mantovani, A., Allavena, P., Sica, A., and Balkwill, F. (2008). Cancer-related inflammation. Nature 454 (7203), 436–444. doi:10.1038/nature07205

PubMed Abstract | 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

Mayakonda, A., Lin, D. C., Assenov, Y., Plass, C., and Koeffler, H. P. (2018). Maftools: Efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 28 (11), 1747–1756. doi:10.1101/gr.239244.118

PubMed Abstract | CrossRef Full Text | Google Scholar

Motorin, Y., and Helm, M. (2011). Rna nucleotide methylation. Wiley Interdiscip. Rev. RNA 2 (5), 611–631. doi:10.1002/wrna.79

PubMed Abstract | CrossRef Full Text | Google Scholar

Pei, Y., and Shuman, S. (2002). Interactions between fission yeast mrna capping enzymes and elongation factor Spt5. J. Biol. Chem. 277 (22), 19639–19648. doi:10.1074/jbc.M200015200

PubMed Abstract | CrossRef Full Text | Google Scholar

Pinello, N., Sun, S., and Wong, J. J. (2018). Aberrant expression of enzymes regulating M(6)a mrna methylation: Implication in cancer. Cancer Biol. Med. 15 (4), 323–334. doi:10.20892/j.issn.2095-3941.2018.0365

PubMed Abstract | CrossRef Full Text | Google Scholar

Quail, D. F., and Joyce, J. A. (2013). Microenvironmental regulation of tumor progression and metastasis. Nat. Med. 19 (11), 1423–1437. doi:10.1038/nm.3394

PubMed Abstract | CrossRef Full Text | Google Scholar

Ready, N., Hellmann, M. D., Awad, M. M., Otterson, G. A., Gutierrez, M., Gainor, J. F., et al. (2019). First-line nivolumab plus ipilimumab in advanced non-small-cell lung cancer (checkmate 568): Outcomes by programmed death ligand 1 and tumor mutational burden as biomarkers. J. Clin. Oncol. 37 (12), 992–1000. doi:10.1200/JCO.18.01042

PubMed Abstract | CrossRef Full Text | Google Scholar

Ruffo, E., Wu, R. C., Bruno, T. C., Workman, C. J., and Vignail, D. A. A. (2019). Lymphocyte-activation gene 3 (Lag3): The next immune checkpoint receptor. Semin. Immunol. 42, 101305. doi:10.1016/j.smim.2019.101305

PubMed Abstract | CrossRef Full Text | Google Scholar

Smyth, G. K. (2004). Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat. Appl. Genet. Mol. Biol. 3, 1027. doi:10.2202/1544-6115.1027

PubMed Abstract | 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

PubMed Abstract | CrossRef Full Text | Google Scholar

Sung, H., Ferlay, J., Siegel, R. L., Laversanne, M., Soerjomataram, I., Jemal, A., 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. doi:10.3322/caac.21660

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, K., Wu, Y. H., Song, Y., and Yu, B. (2021). Indoleamine 2, 3-dioxygenase 1 (Ido1) inhibitors in clinical trials for cancer immunotherapy. J. Hematol. Oncol. 14 (1), 68. doi:10.1186/s13045-021-01080-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Tomikawa, C. (2018). 7-Methylguanosine modifications in transfer rna (trna). Int. J. Mol. Sci. 19 (12), E4080. doi:10.3390/ijms19124080

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, H., Hu, X., Huang, M., Liu, J., Gu, Y., Ma, L., et al. (2019). Mettl3-Mediated mrna M(6)a methylation promotes dendritic cell activation. Nat. Commun. 10 (1), 1898. doi:10.1038/s41467-019-09903-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, W., Kryczek, I., Dostal, L., Lin, H., Tan, L., Zhao, L., et al. (2016). Effector T cells abrogate stroma-mediated chemoresistance in ovarian cancer. Cell 165 (5), 1092–1105. doi:10.1016/j.cell.2016.04.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilkerson, M. D., and Hayes, D. N. (2010). Consensusclusterplus: A class discovery tool with confidence assessments and item tracking. Bioinformatics 26 (12), 1572–1573. doi:10.1093/bioinformatics/btq170

PubMed Abstract | CrossRef Full Text | Google Scholar

Witz, I. P., and Levy-Nissenbaum, O. (2006). The tumor microenvironment in the post-paget era. Cancer Lett. 242 (1), 1–10. doi:10.1016/j.canlet.2005.12.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Yarchoan, M., Hopkins, A., and Jaffee, E. M. (2017). Tumor mutational burden and response rate to Pd-1 inhibition. N. Engl. J. Med. 377 (25), 2500–2501. doi:10.1056/NEJMc1713444

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoshihara, K., Shahmoradgoli, M., Martinez, E., Vegesna, R., Kim, H., Torres-Garcia, W., et al. (2013). Inferring tumour purity and stromal and immune cell admixture from expression data. Nat. Commun. 4, 2612. doi:10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, B., Wu, Q., Li, B., Wang, D., Wang, L., and Zhou, Y. L. (2020). M(6)a 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, H., Meltzer, P., and Davis, S. D. (2013). Rcircos: An R package for circos 2d track plots. BMC Bioinforma. 14, 244. doi:10.1186/1471-2105-14-244

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, L. S., Liu, C., Ma, H., Dai, Q., Sun, H. L., Luo, G., et al. (2019). Transcriptome-wide mapping of internal N(7)-methylguanosine methylome in mammalian mrna. Mol. Cell 74(6), 1304–1316. doi:10.1016/j.molcel.2019.03.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, M., Song, J., Yuan, W., Zhang, W., and Sun, Z. (2021). Roles of rna methylation on tumor immunity and clinical implications. Front. Immunol. 12, 641507. doi:10.3389/fimmu.2021.641507

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, S., Zhang, T., Zheng, L., Liu, H., Song, W., Liu, D., et al. (2021). Combination strategies to maximize the benefits of cancer immunotherapy. J. Hematol. Oncol. 14 (1), 156. doi:10.1186/s13045-021-01164-5

PubMed Abstract | CrossRef Full Text | Google Scholar