m6A regulators featured by tumor immune microenvironment landscapes and correlated with immunotherapy in non-small cell lung cancer (NSCLC)

Introduction: Recent research has confirmed the critical role that epigenetic factors play in regulating the immune response. Nonetheless, what role m⁶A methylation modification might play in the immune response of non-small cell lung cancer (NSCLC) remains vague.

Methods: Herein, the gene expression, copy number variations (CNVs), and somatic mutations of 31 m⁶A regulators in NSCLC and adjacent control samples from the GEO and TCGA databases were comprehensively explored. Using consensus clustering, m⁶A modification patterns were identified. Correlations between m⁶A modification patterns and immune cell infiltration traits in the tumor immune microenvironment (TME) were systematically analyzed. Differentially expressed genes were verified and screened by random forest and cox regression analysis by comparing different m⁶A modification patterns. Based on the retained gene panel, a risk model was built, and m⁶Ascore for each sample was calculated. The function of m⁶Ascore in NSCLC prognosis, tumor somatic mutations, and chemotherapy/immunotherapy response prediction were evaluated.

Results: Consensus clustering classified all NSCLC samples into two m⁶A clusters (m⁶A_clusterA and m⁶A_clusterB) according to the expression levels of 25 m⁶A regulator genes. Hierarchical clustering further divides the NSCLC samples into two m⁶A gene clusters: m⁶AgeneclusterA and m⁶AgeneclusterB. A panel of 83 genes was screened from the 194 differentially expressed genes between m⁶A gene clusters. Based on this, a risk score model was established. m⁶A modification clusters, m⁶A gene clusters, and m⁶Ascore calculated from the risk model were able to predict tumor stages, immune cell infiltration, clinical prognosis, and tumor somatic mutations. NSCLC patients with high m⁶Ascore have poor drug resistance to chemotherapy drugs (Cisplatin and Gemcitabine) and exhibit considerable therapeutic benefits and favorable clinical responses to anti-PD1 or anti-CTLA4 immunotherapy.

Discussion: In conclusion, methylation modification patterns mediated by the m⁶A regulators in individuals play a non-negligible role in prognosis prediction and immunotherapy response, which will facilitate personalized treatment and immunotherapeutic strategies for NSCLC patients in the future.

Introduction

As the second most frequent malignant oncologic disease worldwide, lung cancer accounts for the greatest number of mortality (1). 85% of lung cancers are non-small cell lung cancers (NSCLC). Due to the majority of NSCLC patients being diagnosed at an advanced stage, its overall 5-year survival rate is only 8% (2). Therefore, new approaches are urgently required to explore novel mechanisms of NSCLC that are susceptible to therapeutic inventions. Recent research has shown that tumor immune microenvironment (TME) immune cell infiltrating characteristics are closely correlated with m⁶A modifications, which might provide an alternative choice (3–5). N6-methyladenosine (m⁶A) RNA modification is the most prevalent kind of RNA modification in eukaryotic cells. It is essential for the regulation of epigenetic processes, a variety of physiological functions, and the development of disease (4, 6, 7).

Methyltransferases, also known as “writers”, promote m⁶A methylation modification in RNA; demethylases, also known as “erasers”, remove m⁶A methyl groups from RNA; and binding proteins, also known as “readers”, bind to the m⁶A methylation site in RNA and perform specific biological functions. Three types of proteins regulate the dynamic and reversible process of m⁶A modification (7, 8). Comprehending the functions of m⁶A modification in post-transcriptional regulation would be more accessible by comprehensively exploring the expression and function of m⁶A regulatory proteins (9, 10). The development of malignant tumors and immunomodulatory disorders are correlated with dysregulated expression and genetic alterations of m⁶A regulator genes (10–12), demonstrating that m⁶A regulators may be crucial in regulating the immunological microenvironment of malignancies.

An increasing number of research have revealed the relationships between m⁶A modifications and the immune cell infiltrating characteristics of the TME (3–5). Wang et al. reported that dendritic cell activation and maturation were aided by METTL3-mediated m⁶A modification. Co-stimulatory molecules CD80 and CD40 are expressed less when METTL3 is knocked out (13). Studies suggest the vital role of TME in cancer progression and therapeutic responses with increasing evidence (14, 15). The immune response and the benefit of chemotherapy are reflected in the TME context that was established at diagnosis (15, 16). Clinical outcomes in various cancers are correlated with changes in the compositions of CD8 and CD4 positive T cells, macrophages, and cancer-associated fibroblast infiltration in the TME (17, 18).

The immune checkpoint blockade (ICB) therapy, which specifically targets the cytotoxic T lymphocyte antigen 4 (CTLA-4) and programmed cell death 1 (PD-1) or its ligands (PD-L1), has been used for cancer immunotherapy and has shown promising clinical results (19). Only a small portion of patients, nevertheless, might benefit from ICB treatment. Thus, exploring the TME and its associated mechanisms is urgently needed to improve immunotherapy’s efficacy. As was previously stated, the microenvironment of malignancies and immune cells are closely correlated with m⁶A modifications. Therefore, our comprehension of immunological regulation in the TME and immunotherapeutic tactics development will be enhanced by fully exploiting the effects of the regulatory network of RNA m⁶A modification enzymes on TME cells.

By thoroughly analyzing the gene expression profile of m⁶A regulators in 1558 NSCLC samples, the present study could distinguish different m⁶A modification patterns. We systematically correlated the characteristics of TME cell infiltration with genomic traits as well as the clinical and pathologic characteristics of NSCLC. We estimated the patterns of TME infiltration in 1,558 NSCLC samples. Further, a risk model was constructed based on a panel of 83 genes with differential expression, and m⁶Ascore was calculated for each sample. Consequently, we developed a method for quantifying the m⁶Ascore and found that it is a robust prognostic biomarker and a significant predictor of response to chemotherapy (Cisplatin and Gemcitabine) and immunotherapy (anti-PD1 or anti-CTLA4 immunotherapy).

Materials and methods

Data collection and pre-processing

The study’s workflow is depicted in Figure S1. TCGA Gene expression data, genomic mutation data (including somatic mutation and copy number variation), and corresponding clinical data of NSCLC samples were downloaded from UCSC Xena (https://xenabrowser.net/datapages/). From the NCBI GEO database (https://www.ncbi.nlm.nih.gov/geo/), two additional datasets of NSCLC samples were downloaded. This project gathered 1705 samples in total, including TCGA-LUAD (N=568), TCGA-LUSC (N=545), GSE68465 (N=462), and GSE4573 (N=130) datasets (Table S1). Table S2 provides a summary of the clinical details of these samples. Table S3 lists the clinical details of each sample from TCGA dataset. Among the 1579 NSCLC samples, survival status and survival time are available for 1558 samples. Transcripts per kilobase million (TPM) values were generated from FPKM values of RNA sequencing data downloaded from TCGA. The raw “CEL” files for the Affymetrix-produced GEO microarray data were downloaded. R packages “affy” and “simpleaffy” were employed to adjust the background and perform quantile normalization. The ComBat function from the “SVA” R package was used to remove the batch effect between TCGA and GEO datasets and the integrated data after removing the batch effect is provided in Table S4. The genomic mutation status of NSCLC patients from the TCGA database was displayed in an oncoplot generated with the R package “maftools”.

Unsupervised consensus clustering of 25 m⁶A regulators

A total of 31 m⁶A regulators were gathered from the papers on m⁶A methylation modification. Due to the lacking of six m⁶A regulator genes (IGF2BP1, KIAA1429, METTL16, METTL14, ALKBH5, and RBMX) in GEO datasets, the remaining 25 m⁶A regulators were curated, including seven writers (METTL3, WTAP, RBM15, RBM15B, CBLL1, ZC3H13, and ZCCHC4), one eraser (FTO), 15 readers (YTHDF1, YTHDF2, YTHDF3, YTHDC1, YTHDC2, IGF2BP2, IGF2BP3, EIF3A, HNRNPA2B1, HNRNPC, FMR1, LRPPRC, ELAVL1, PRRC2A, and SND1), and two repellers (G3BP1 and G3BP2). Unsupervised consensus clustering was carried out using the expression levels of 25 m⁶A regulator genes to discriminate different m⁶A modification patterns with the R package “ConsensusClusterPlus”, which is based on a computational method called consensus clustering (20). Consensus Cumulative Distribution Function (CDF) and Delta area (relative change of area under the CDF curve) were used to select the proper clustering numbers within the high-throughput RNA-seq data. We use the parameters of a maximum evaluated k of 20, an 80% resampling rate, 1000 iterations, and Euclidean distance to determine the optimal number of clusters and to guarantee robustness.

Gene set variation analysis (GSVA) and functional annotation

Using the R package “GSVA”, gene set variation analysis (GSVA), an unsupervised and non-parametric method, was used to compute the pathway enrichment scores in order to explore the biological process variations among different m⁶A modification patterns (21). To conduct GSVA analysis, the well-defined KEGG gene sets of “c2.cp.kegg.v6.2.symbols” were downloaded from the MSigDB database (https://www.gsea-msigdb.org/gsea/index.jsp). Gene set enrichment analysis with a cutoff value of false discovery rate (FDR) < 0.01 was used to examine biological processes correlated with m⁶A regulators using the R package “clusterProfiler”.

Immune cell infiltration estimation

The single-sample gene-set enrichment analysis (ssGSEA) function from the R package “GSVA” was used to estimate the levels of immune cell infiltration. ssGSEA evaluates a specific gene set, including the gene expression data of 28 immune cells that represent different immune cell types, immune-related functions, and pathways in NSCLC (22). The enrichment scores representing the relative level of immune cell infiltration were compared between samples that belong to different m⁶A clusters by the Wilcox test. To illustrate their prognostic values, significantly different immune cells between m⁶A clusters were further analyzed by cox regression and visualized by the R package “forestplot”.

Identification of differentially expressed genes (DEGs) between distinct m⁶A clusters

Differential expression analysis was carried out using the R package “limma”, an empirical Bayesian approach, to identify DEGs associated with m⁶A (23). Genes with adjusted p < 0.05 (Benjamini-Hochberg adjustment) and |fold change| > 1.5 in expression were regarded as DEGs. Hierarchical cluster analysis was used to divide NSCLC patients into genomic clusters based on the DEGs. We used a bottoms-up approach called agglomerative clustering, in which the data points were initially isolated as separate groups and then merged iteratively based on similarity until one cluster had been formed. The similarity was measured with Ward’s linkage, namely the Euclidean distance between two clusters was defined by the increase in the sum of squared after the clusters were merged.

Dimension reduction and generation of m⁶A gene signatures

For all the identified differentially expressed genes between m⁶A clusters, the supervised machine learning algorithm random forest was applied for dimensionality reduction. After removing the redundant genes, the remaining genes more relevant to m⁶A modification went through survival analysis with the R package “Survminer”. Genes with significant survival results (p < 0.05) were added to a Cox regression model in further analysis. To explore the similarity between gene expression profile and prognosis efficiency, the m⁶A score was introduced. The m⁶A score was defined refer to the definition of gene expression grade index (GGI) (24), and the formula is as follows:

$$m^{6}Ascore=scale\left({\displaystyle \sum X-{\displaystyle \sum Y}}\right)$$

Where scale represents the transformation parameter of standardization and X and Y are the expression of gene sets with positive and negative Cox coefficients, respectively. The optimal cutoff value was computed using the surv-cutpoint function from the “survival” R package. All samples were subsequently stratified into m⁶Ascore-high and m⁶Ascore-low subgroups, and their relationships with prognosis were evaluated as well.

Correlation between m⁶A score and other pertinent biological processes

Mariathasan et al. have constructed a collection of genes to store genes related to a sort of biological processes, including Angiogenesis; CD8 T effector; Antigen processing machinery; Cell cycle; Cell cycle regulators; KEGG discovered histones; DNA damage repair; DNA replication; Fanconi anemia; FGFR3-related genes; Homologous recombination; Immune checkpoint; EMT1, EMT2, and EMT3 epithelial-mesenchymal transition (EMT) markers; Mismatch repair; Nucleotide excision repair; Pan-F-TBRS; WNT target (25). GSVA was used to quantify the above-mentioned biological processes in each sample with an enrichment score. Pearson correlation analysis was carried out between m⁶Ascore and enrichment score to reveal the relationship between m⁶Ascore and certain associated biological pathways.

Copy number variation (CNV) analysis

The Genomic Identification of Significant Targets in Cancer (GISTIC) method was used to identify the common CNV regions across all samples with TCGA Copy Number Segment data. The significance threshold of GISTIC was: False Discovery Rate (FDR), namely q-value ≤ 0. 05. The peak region for each significant region was identified with a confidence interval of 0.95. The GISTIC analysis made use of the MutSigCV module from GenePattern, an online analysis tool provided by the Broad Institute (https://cloud.genepattern.org/gp/pages/index.jsf).

Half maximal inhibitory concentration (IC₅₀) prediction and tumor immune dysfunction and exclusion (TIDE) analysis

In order to predict the clinical chemotherapeutic response from tumor gene expression profiles, the IC₅₀ values of clinical drugs (Cisplatin, Gemcitabine) were estimated using the R package “pRRophetic” (26). Then IC₅₀ values between high m⁶Ascore samples and low m⁶Ascore samples were compared. In addition, signatures of T cell dysfunction and exclusion were analyzed using the online algorithm TIDE (http://tide.dfci.harvard.edu/) to predict the cancer immunotherapy response to immune checkpoint blockade (ICB) (27). A higher TIDE prediction score indicates a poor prognosis and a poor response to ICB therapy.

NSCLC cell line m⁶A score calculation and chemotherapy drug IC₅₀ validation

The Cancer Cell Line Encyclopedia was used to download the gene expression data for NSCLC cell lines (CCLE: https://sites.broadinstitute.org/ccle/). The m⁶A score for each cell line was calculated using the m⁶Ascore formula. NSCLC cell lines were stratified into m⁶Ascore-high and m⁶Ascore-low groups based on the cutoff value. Genomics of Drug Sensitivity in Cancer (GDSC: https://www.cancerrxgene.org/) provided information on the drug sensitivity of the chemotherapeutic drugs Cisplatin and Gemcitabine. In general, a dose titration of Cisplatin and Gemcitabine (6 nM-6 μM for Cisplatin; 0.1 nM-0.1 μM for Gemcitabine) was administered to STRs-verified cell lines for 72 hours in culture media after they had been seeded in 96-well plates and grown for 24 hours at 37°C in 5% CO₂. Cell viability was determined using either a metabolic test (Resazurin or CellTiter-Glo) or DNA dye (Syto60). Every screening plate was put through rigorous quality control procedures. GraphPad Prism 9 software was used to calculate the IC₅₀ of Cisplatin and Gemcitabine in NSCLC cell lines and to generate dose-response curves.

Statistical analysis

The Wilcoxon test was used to determine whether scores between the two sample groups were statistically significant. The log-rank test from the R package “Survminer” was used to assess the statistical significance between the prognostic survival curves, which were generated using the Kaplan-Meier method. The prediction performance of immunotherapy by m⁶Ascore was assessed using the receiver operating character (ROC) curve, and the area under the ROC curve (AUC) was computed using the R package “pROC”. Patients with high and low m⁶Ascores had different mutational landscapes, which were visualized using the “maftools” R package.

Results

Genetic and transcriptional alteration landscapes of four types of RNA m⁶A methylation regulators in NSCLC

Table S5 lists the 31 RNA m⁶A modification regulators that were used in this study, which included ten methyltransferases “writers”, two demethylases “Erasers”, 17 RNA binding proteins “Readers”, and two “Repellers”. We first summarized the occurrence frequency of somatic mutations and copy number variations in 31 m⁶A regulator genes in TCGA NSCLC samples. Specifically, ZC3H13 and KIAA1429 had the greatest mutation frequency, reaching 4%. The most frequent mutation type was missense mutation (Figure 1A). Copy number variation (CNV) frequency analysis showed that copy numbers were generally changed among 31 regulatory factors. Copy number amplification commonly occurred in genes such as IGF2BP2, KIAA1429, and YTHDC1, while copy number deletion commonly occurred in genes such as RBM15, YTHDF2, and ZC3H13 (Figure 1B, Table S6). According to the expression of these 31 m⁶A regulator genes, principle component analysis could differentiate TCGA NSCLC samples from adjacent normal samples (Figure 1C). Gene expression analysis of 31 m⁶A regulators between TCGA NSCLC samples and adjacent control samples showed that most regulator genes were significantly overexpressed in NSCLC tissues, especially reader genes (IGF2BP1, IGF2BP2, and IGF2BP3) from the IGF2BPs family (Figure 1D).

FIGURE 1

Figure 1 Genetic and transcriptional variation landscapes of RNA m⁶A modification regulator genes in the TCGA dataset. (A) The frequency of mutations and the distribution of mutation types of m⁶A regulator genes in NSCLC tissue samples. (B) The frequency of m⁶A regulator gene copy number variation in NSCLC tissue samples. Blue triangles denote deletion and red dots denote amplification. (C) Principle component analysis results to distinguish NSCLC samples from adjacent normal tissues based on the expression of m⁶A regulator genes. (D) The expression profile of m⁶A regulator genes in cancer samples and adjacent normal samples. Asterisks in red indicate genes belonging to the IGF2BPs family. NS p > 0.05; *p < 0.05; **p < 0.01; ***p < 0.001.

Unsupervised clustering of m⁶A regulator genes

An m⁶A regulatory network was constructed for 1557 NSCLC samples with expression data and survival information available to describe the spearman correlations within m⁶A regulator genes and the correlations between m⁶A regulator genes and NSCLC prognosis (Figure 2A). Results suggested that different m⁶A modification patterns might be significantly influenced by interactions between different functional types of m⁶A regulators. Due to the absence of IGF2BP1, KIAA1429, METTL16, METTL14, ALKBH5, and RBMX expression data in GEO datasets, the remaining 25 m⁶A regulator genes were included for consensus clustering. Two subgroups were identified using unsupervised consensus clustering, and their respective names were m⁶A_clusterA and m⁶A_clusterB (Figure 2B). Biological pathway differences between two m⁶A clusters were identified using GSVA enrichment analysis. m⁶A_clusterA significantly enriched metabolic-related biological processes like fatty acid metabolism and tryptophan metabolism, while m⁶A_clusterB considerably enriched replication- and transcription-related biological pathways such as DNA repair and mismatch repair (Figure 2C, Table S7). The two m⁶A clusters also had significantly different prognoses, as demonstrated by the Kaplan-Meier curve of overall survival (p = 0.007) (Figure 2D).

FIGURE 2

Figure 2 m⁶A regulator genes unsupervised clustering in NSCLC samples. (A) m⁶A regulator gene interactions in NSCLC. The circle size reflects how each gene affects the ability to predict survival. The stronger the association between gene expression and prognosis, the larger the circle. The red circle indicates prognostic protective factors and the blue circle indicates prognostic risk factors. The spearman correlations between genes are linked by lines connecting them, with positive correlations denoted by solid grey lines and negative correlations by purple dotted lines. The interaction strength between m⁶A regulators determines line thickness. (B) Consensus clustering of m⁶A regulator genes, 1 and 2 represent two subgroups. (C) GSVA enrichment analysis revealed the activation state of biological pathways with distinct m⁶A clusters. The heatmap was used to depict these biological processes, with red indicating activation and blue indicating inhibition. (D) Two m⁶A clusters showed a significant survival differences, as depicted in the Kaplan-Meier curve.

The heatmap of 25 m⁶A regulator genes, which was classified by two m⁶A clusters, showed the relationship between the expression level and matching clinical information, such as cancer type, smoking indicator, stage, sex, and age. The distribution of clinical information between the two m⁶A clusters does not significantly differ. Notably, NSCLC patients in m⁶A_clusterB were more likely to express the IGF2BPs family, including IGF2BP2 and IGF2BP3 (Figure 3A). Among the most common gene mutations in people with NSCLC, we analyzed the mutation status of nine genes: EGFR, ALK, ROS1, BRAF, KRAS, CD274, MET, RET, and ERBB2. Results showed that only the mutation status of KRAS (p = 7.74e-05) and RET (p = 0.015) were significantly different between m⁶A_clusterA and m⁶A_clusterB, while the other seven genes were not significant (Chi-square test, p > 0.05). Excluding KRAS, NSCLC samples without the above-mentioned gene mutations are more in m⁶A_clusterA than in m⁶A_clusterB (Figure 3B). Furthermore, to illustrate the impact of m⁶A regulators on immune cell infiltration, ssGSEA was conducted based on the sample expression data and obtained the proportion distribution of 28 immune cell types in two different m⁶A clusters of NSCLC samples (Figure 3C). Results showed that 22 of 28 immune infiltration cells had differential expression between the two m⁶A clusters and most of them were highly expressed in m⁶A_clusterA, except for activated CD4 T cells and memory B cells, which were highly expressed in m⁶A_clusterB (Figure 3C, Table S8). This suggests that m⁶A_clusterA has an immune microenvironment that is hot and suppressive. For the 22 differentially expressed immune infiltration cells, univariate Cox regression analysis revealed that activated CD4 T cells (p = 0.0085), monocytes (p = 0.024), and activated B cells (p = 0.0212) were significantly associated with the prognosis of two m⁶A clusters (Figure 3D, Table S9).

FIGURE 3

Figure 3 Transcriptome characteristics and immune cell infiltration features in two m⁶A regulator gene clusters. (A) The expression pattern of the m⁶A regulator genes between the two m⁶A clusters that matched the clinical information, including cancer type, smoking indicator, stage, sex, and age. (B) The nine most common gene mutation statuses in NSCLC. (C) The boxplot shows the abundance of 28 infiltrating immune cells among the two m⁶A clusters. NS p > 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001. (D) Forest plots of hazard ratios (HRs) for 22 immune infiltration cells associated with OS and meaningful immune cells were marked in red (p < 0.05).

Identification of m⁶A signature genes

Using the R package “limma”, 194 m⁶A phenotype-related differentially expressed genes between m⁶A clusters were identified in order to explore the probable biological functions of each m⁶Acluster (Table S10). Unsupervised hierarchical cluster analysis has classified the NSCLC patients into two genomic clusters, termed m⁶A_gene_clusterA and m⁶A_gene_clusterB, which are roughly in accordance with the m⁶A modification pattern, based on the expression of the 194 differentially expressed genes (Figure 4A, Table S11). The log-rank test and Kaplan-Meier curve indicate that the genomic phenotypes of m⁶A modification were significantly related to OS in NSCLC patients and patients in m⁶A_gene_clusterB had better prognoses (Figure 4B, p = 0.0014). Among the 25 m⁶A regulator genes, 17 were significantly more abundantly expressed in m⁶A_gene_clusterA than in m⁶A_gene_clusterB, while three were significantly highly expressed in m⁶A_gene_clusterB (Figure 4C).

FIGURE 4

Figure 4 Genes that differ in expression are associated with m⁶A. (A) Unsupervised analysis and hierarchical clustering of differentially expressed genes associated with the m⁶A phenotype to divide NSCLC patients into two genomic groups.Patient annotations include m⁶A cluster, cancer type, smoking indicator, stage, sex, and age. (B) Kaplan-Meier curve depicts the association between m⁶A modification genomic phenotypes and OS. (C) The two m⁶A gene clusters’ expression of 25 m⁶A regulator genes. NS p > 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001.

Clinical and transcriptome features of m⁶A-related phenotypes

Among the 194 differentially expressed genes derived from the previous analysis, redundant genes were removed using the random forest algorithm, leaving 83 feature genes that were most closely related to the m⁶A relationship between these genes and NSCLC patient survival. These 83 genes were classified into two groups, a positive group with 30 genes and a negative group with 53 genes, based on the coefficient value of genes obtained from the Cox regression (Table S12). According to the m⁶Ascore formula, m⁶A scores for all samples were calculated and the best cutoff of m⁶Ascore (cutoff = -0. 7437558) was established by the “surv_cutpoint” function of the R package to classify NSCLC samples into m⁶Ascore high and low groups (Figure 5A, Table S13). Results of the survival analysis demonstrated that the m⁶Ascore might accurately characterize the prognosis of NSCLC patients (p < 0.0001), with the m⁶Ascore_low group has a good prognosis while the m⁶Ascore_high group has a bad prognosis (Figure 5B). To better depict the function of m⁶Ascore, GSVA was performed with known gene signatures. A significant positive correlation between m⁶Ascore and biological processes, including the cell cycle and DNA replication, has been found using pearson correlation analysis. In contrast, the correlation between m⁶Ascore and other biological processes, such as angiogenesis and EMT3, is significantly negative (Figure 5C, Table S14). m⁶Ascores between m⁶A clusters and between m⁶A gene clusters were compared using the Wilcox test. Results showed that m⁶A_clusterA had a substantially lower m⁶Ascore than m⁶A_clusterB, while m⁶A_gene_clusterA had a significantly higher m⁶Ascore than m⁶A_gene_clusterB (Figure 5D). Furthermore, as shown in Figure S2, there are significant differences in m⁶Ascores across several clinical categories, including EGFR mutation status, age, and sex (Figures S2A, B). Additionally, in both the TCGA and GEO datasets, NSCLC patients with high and low m⁶Ascores had significantly different overall survival probabilities (Figure S2C).

FIGURE 5

Figure 5 Construction of m⁶A scores. (A) An alluvial diagram of the m⁶A clusters showing different m⁶A gene clusters, m⁶Ascores, and survival outcomes. (B) The Kaplan-Meier curve demonstrated a significant correlation between m⁶A scores in the high and low groups and overall survival. (C) GSVA and Pearson correlation analysis show correlations between m⁶Ascore and biological pathways constructed with known gene signatures in NSCLC. A positive correlation is indicated by red, while a negative correlation is indicated by blue. The circle size is inversely proportional to the levels of significance, and X indicates no significant correlation. (D) Boxplots show the differences in m⁶A scores between m⁶A clusters and between m⁶A gene clusters. ****p < 0.0001.

Molecular characteristics of m⁶Ascore groups in TCGA datasets

The distribution of tumor somatic mutations in TCGA NSCLC datasets is visualized by the R package “maftools” (Figures 6A, B, left), and the GISTIC algorithm was used to evaluate and visualize the distribution of somatic copy number alterations in groups with high and low m⁶Ascores, respectively (Figures 6A, B, right). The mutational landscape revealed that the top 15 genes in the high m⁶Ascore group had higher tumor mutation frequencies than those in the low m⁶Ascore group. The most mutational gene was TP53 (80% vs. 37%). GISTIC plots revealed that the m⁶Ascore low group had significantly fewer copy number alterations than the m⁶Ascore high group, which is consistent with somatic mutations. Based on these data, we were able to more fully illustrate the impact that m⁶A score classification has on genomic variation and to demonstrate the potentially intricate interactions between individual somatic mutations/alterations and m⁶A modifications.

FIGURE 6

Figure 6 Molecular features of m⁶Ascore groups in TCGA datasets. (A) The distribution of tumor somatic mutations in NSCLC with a high m⁶Ascore is depicted in the waterfall plot. Each column represents a patient. TMB is depicted in the upper barplot, and each gene’s mutation frequency is indicated by the number on the left. The horizontal barplot on the right depicts the percentage of each variant type, and overall survival status is shown as patient annotations (left). GISTIC plot depicts the distribution of somatic copy number alterations in NSCLC with a high m⁶Ascore, regions of copy number amplification are highlighted in red, and regions of copy number deletion are highlighted in blue. Several significant gene names are marked (right). (B) The waterfall plot displays the somatic mutation distribution (left), and the GISTIC plot shows the distribution of somatic copy number alterations (right) in NSCLC with a low m⁶Ascore.

The value of m⁶Ascore in predicting chemotherapy and immunotherapy response

To extend the potential therapeutic use of m⁶Ascore, we explored whether the intrinsic m⁶Ascore in cancer cells could predict their response to various drugs, which was inspired by the cross-talk between m⁶Ascore and many key cancer-related pathways. Using the R package “pRRophetic” and the expression profile from TCGA and GEO datasets, IC₅₀ values of chemotherapeutic drugs Cisplatin and Gemcitabine were calculated. According to a comparison of the relative distribution of Cisplatin and Gemcitabine IC₅₀ values, the IC₅₀ value in the low m⁶A score group was significantly higher than that in the high m⁶A score group, indicating that the high m⁶A score group had poor drug resistance (Figure 7A). Moreover, based on the mRNA expression profile in the TCGA data set, the TIDE algorithm was employed to assess the clinical effect of ICB treatment in the m⁶Ascore high and low groups. Results revealed that patients with high m⁶Ascores had TIDE scores that were significantly lower than those with low m⁶Ascores. That was to say, compared to patients with low m⁶Ascores, NSCLC patients with high m⁶Ascores exhibited a better therapeutic benefit and clinical response to anti-PD1 or anti-CTLA4 immunotherapy (Figure 7B). According to the ROC curve, the m⁶Ascore may be a reliable biomarker for predicting outcomes and evaluating the therapeutic efficacy of anti-PD1 and anti-CTLA4 treatments (AUC = 0.88) (Figure 7C).

FIGURE 7

Figure 7 The response of the high and low m⁶Ascore groups to chemotherapy/immunotherapy. (A) Differences in IC₅₀ of Cisplatin and Gemcitabine treatment. (B) Differences in TIDE prediction score. (C) The receiver operator characteristic curve (ROC) shows the predictive performance of m⁶Ascore in NSCLC patients receiving anti-PD1 and anti-CTLA4 therapy (AUC, 0.88).

To validate the chemotherapy response in high and low m⁶Ascore groups, we first classified NSCLC cells into m⁶Ascore-high (NCI-H23, NCI-H2023, COR-L32, NCI-H1781, and NCI-H2170) and m⁶Ascore-low (VMRC-LCD, NCI-H1838, NCI-H2342, NCI-H661, and NCI-H2347) groups with distinct differences. The dose-response curves and IC₅₀ values indicated that Gemcitabine treatment is generally more effective in NSCLC cell lines. The m⁶Ascore-high group of NSCLC cell lines has higher anti-tumor activity when treated with Cisplatin and Gemcitabine (Figures 8A, B), which is in line with the predicted result that the group with a high m⁶A score has poor drug resistance.

FIGURE 8

Figure 8 Validation of chemotherapy drugs’ IC₅₀. (A, B) Measurement of IC₅₀ by DNA dye or metabolic assay in NSCLC cell panels (VMRC-LCD, NCI-H1838, NCI-H2342, NCI-H661, NCI-H2347, NCI-H23, NCI-H2023, COR-L32, NCI-H1781, NCI-H2170) treated with Cisplatin and Gemcitabine.

Discussion

The immunological state of TME in diverse malignancies is regulated by RNA m⁶A modification (28, 29). Recent studies have uncovered the relationships between m⁶A regulators and immune cell infiltration. The accumulation of myeloid-derived suppressor cells necessary to maintain the activation and proliferation of CD4+ and CD8+ T cells was reduced when METTL3 was selectively depleted in colorectal cancer cells (30). Natural killer cell homeostasis and maturation, as well as their anti-tumor and antiviral activity, are positively regulated by YTHDF2 (31). YTHDF1 deficiency in classical dendritic cells could enhance antigen presentation, initiate anti-tumor responses, and improve the therapeutic effectiveness of PD-L1 checkpoint blockade (32). Nevertheless, these research concentrated on the single m⁶A modification regulators, and the integrated roles of various m⁶A regulators in modifying immune characteristics in NSCLC need to be further investigated. More efficient immunotherapy approaches will result from a knowledge of the involvement of distinct m⁶A modification patterns in TME cell infiltration and the TME anti-tumor immune response.

The global transcriptional and genetic profiles of m⁶A modification regulator genes and their mutual correlation in NSCLC were the focus of this study. When compared to healthy controls, the expression of METTL3 and METTL14 in NSCLC samples exhibited a contrary trend. As well-studied m⁶A writers, METTL3 and METTL14 are reasonably thought to have similar functions. Although a similar phenomenon exists in prostate cancer (33), further studies are required in the future to figure this out. m⁶A readers from the IGF2BPs family showed high expression levels and mutation frequency, especially IGF2BP2. According to Li et al., IGF2BP2 and IGF2BP3 was essential for lung cancer progression, and they could identify and stabilize m⁶A sites and function as ‘readers’ in the post-transcriptional regulation manner (34). A study found that IGF2BP2 regulated macrophage activation in an m⁶A-dependent manner, which indicated a potential therapeutic target of macrophages in inflammatory diseases (35). RBM15 and its paralogue RBM15B contain RNA-binding motifs, which make it easier to recruit m⁶A methyltransferase to specific sites in RNA (36, 37). Our results showed that RBM15 and RBM15B were significantly up-regulated in NSCLC samples and had frequent CNV alterations, indicating their potential role in promoting cancer cell migration and invasion (38). Our comprehension of how epigenetic regulation affects diverse physiological processes and TME cell-infiltrating characterization will be enhanced by a thorough evaluation of the patterns of m⁶A modification, which will highlight the heterogeneity of m⁶A modification. These regulators’ effects on immune infiltration mechanisms via m⁶A modification require additional biology studies utilizing cell culture and even PDX mice models.

With consistent clustering, we have identified two independent m⁶A modification patterns with significantly different TME immune cell infiltration traits, among which m⁶A_clusterA showed immune inflamed phenotype and enriched abundant immune cells, known as hot tumor. Immune checkpoint inhibitors are frequently associated with greater benefits in such tumors (39). Therefore, m⁶A_clusterA had a better prognosis. As defined by us, m⁶A signature genes are the genes that are expressed differently in different m⁶A modification patterns. Two genomic clusters were identified using m⁶A signature genes, which is consistent with m⁶A modification clusters. m⁶A regulator genes were also differentially expressed between genomic clusters. This phenomenon demonstrates once more how the m⁶A modification significantly shapes various TME landscapes. Taking into account the individual variability of m⁶A modification, an m⁶Ascore algorithm based on the m⁶A signature genes was developed to compute the m⁶Ascore of each sample. m⁶Agenecluster and m⁶Ascore groups have apparent differences in prognosis, clinical features, or molecular characteristics. As a result, we demonstrated that m⁶Ascore might be used to evaluate the clinicopathological traits of patients, such as tumor inflammation, prognosis, genetic variation, and so forth. This demonstrated that m⁶Ascore was robust and reliable and could be utilized to identify the tumor immune phenotypes by comprehensively assessing individual tumor m⁶A modification patterns. Therefore, m⁶Ascore could be utilized to predict the effectiveness of chemotherapy and the patient’s clinical response to anti-PD1 and anti-CTLA4 immunotherapy for NSCLC, in addition to being used as an independent prognostic biomarker to predict patients’ survival. Immunotherapy response results predicted from m⁶Ascore classification demonstrated that the predictive performance of the m⁶Ascore in NSCLC patients treated with drugs anti-PD1 and anti-CTLA4 reaches a meaningful result (AUC = 0.88). Several novel ideas for cancer immunotherapy that alter m⁶A modification patterns by targeting m⁶A regulators or m⁶A signature genes were presented in our study. By harnessing the immune system, firing up the TME cell infiltration characterization to turn tumors from “cold” to “hot” is a promising strategy to explore novel drug combinations or novel immunotherapeutic agents.

This study is the first to comprehensively and methodically analyze the relationships between m⁶A regulator modification patterns, immune infiltration, and treatment resistance in NSCLC. Their relationships attracted emerging attention in recent years (40). The numerous results generated from this study will provide hints for other researchers in deep mechanism study and direct m⁶A immunotherapy in lung cancer. In addition, in the analysis process, our study is the first to use public drug sensitivity experimental data to validate the value of m⁶Ascore in predicting chemotherapy response, making our prediction results more reliable. However, limitations exist in this study. Firstly, when merging two large TCGA cohorts and two GEO datasets of NSCLC patients, intra-tumor heterogeneity in different databases was not considered. A study has proven that tumor heterogeneity affects cancer immunotherapy (41). Secondly, this study confirmed the strong impact of m⁶A modification on the immune characteristics of NSCLC, but only theoretically valid. Bench works are warranted in the future to explore the underlying mechanisms. Thirdly, this study lacks an external clinical cohort to verify the results. Therefore, large external NSCLC cohorts are needed for further validation.

Conclusions

In this study, a thorough examination of m⁶A regulators in NSCLC was conducted with bioinformatics analysis. We initially screened DEGs of m⁶Aclusters, then separated NSCLC patients into two categories, and computed the m⁶Ascore in order to build a risk model with a good predictive value for prognosis. The results of this study may help us learn more about how m⁶A signaling influences the progression and prognosis of NSCLC. This work emphasizes the important clinical implications of RNA modifications’ cross-talks and contributes to developing individualized immune therapeutic approaches for NSCLC patients.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author contributions

YW and BY conceived and designed this study. BY performed the analysis, interpreted the data, and wrote the manuscript. HQ, MZ, YY, XT, HY, and WH helped review the data and manuscript. YW supervised the study. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by grants from The National Natural Science Foundation of China (42125707, 41931291), the Special Research Fund for Central Universities, Peking Union Medical College (3332022033), the Non-profit Central Research Institute Fund of Chinese Academy of Medical Sciences (2019PT310027), the Chinese Academy of Medical Sciences Innovation Fund for Medical Sciences (CIFMS: 2021-RC310-006, 2021-RC310-018, 2021-I2M-1-018), and the State Key Laboratory of Molecular Oncology (SKLMO-2021-21).

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.

The reviewer YZ declared a shared affiliation, with no collaboration, with the authors BY, HQ, JZ, MZ, YY, and YW to the handling editor at the time of review.

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/fonc.2022.1087753/full#supplementary-material

Abbreviations

1. AUC, area under the ROC curves; CNV, copy number variation; CTLA4, cytotoxic T-lymphocyte associated protein 4; DEGs, differentially expressed genes; EMT, epithelial-mesenchymal transition; FDR, false discovery rate; GEO, gene expression omnibus; GISTIC, genomic identification of significant targets in cancer; GSVA, gene set variation analysis; IC₅₀, half maximal inhibitory concentration; ICB, immune checkpoint blockade; m⁶A, N6-methyladenosine; NSCLC, non-small cell lung cancer; PD-L1, programmed death-ligand 1; ROC, receiver operating characteristic; ssGSEA, single sample gene set enrichment analysis; TCGA, the cancer genome atlas; TIDE, tumor immune dysfunction and exclusion; TME, tumor immune microenvironment.

References

1. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin (2021) 71:209–49. doi: 10.3322/caac.21660

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Thai AA, Solomon BJ, Sequist LV, Gainor JF, Heist RS. Lung cancer. Lancet (2021) 398:535–54. doi: 10.1016/s0140-6736(21)00312-3

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Guo L, Yang H, Zhou C, Shi Y, Huang L, Zhang J. N6-methyladenosine RNA modification in the tumor immune microenvironment: Novel implications for immunotherapy. Front Immunol (2021) 12:773570. doi: 10.3389/fimmu.2021.773570

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Dong M, Shen W, Yang G, Yang Z, Li X. Analysis of m6A methylation modification patterns and tumor immune microenvironment in breast cancer. Front Cell Dev Biol (2022) 10:785058. doi: 10.3389/fcell.2022.785058

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Zhong J, Liu Z, Cai C, Duan X, Deng T, Zeng G. m(6)A modification patterns and tumor immune landscape in clear cell renal carcinoma. J Immunother Cancer (2021) 9:e001646. doi: 10.1136/jitc-2020-001646

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Roundtree IA, Evans ME, Pan T, He C. Dynamic RNA modifications in gene expression regulation. Cell (2017) 169:1187–200. doi: 10.1016/j.cell.2017.05.045

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Zaccara S, Ries RJ, Jaffrey SR. Reading, writing and erasing mRNA methylation. Nat Rev Mol Cell Biol (2019) 20:608–24. doi: 10.1038/s41580-019-0168-5

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Yang Y, Hsu PJ, Chen YS, Yang YG. Dynamic transcriptomic m(6)A decoration: writers, erasers, readers and functions in RNA metabolism. Cell Res (2018) 28:616–24. doi: 10.1038/s41422-018-0040-8

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Li XC, Jin F, Wang BY, Yin XJ, Hong W, Tian FJ. The m6A demethylase ALKBH5 controls trophoblast invasion at the maternal-fetal interface by regulating the stability of CYR61 mRNA. Theranostics (2019) 9:3853–65. doi: 10.7150/thno.31868

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Yang S, Wei J, Cui YH, Park G, Shah P, Deng Y, et al. m(6)A mRNA demethylase FTO regulates melanoma tumorigenicity and response to anti-PD-1 blockade. Nat Commun (2019) 10:2782. doi: 10.1038/s41467-019-10669-0

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Tong J, Cao G, Zhang T, Sefik E, Amezcua Vesely MC, Broughton JP, et al. m(6)A mRNA methylation sustains treg suppressive functions. Cell Res (2018) 28:253–6. doi: 10.1038/cr.2018.7

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Zhang Q, Xu K. The role of regulators of RNA m6A methylation in lung cancer. Genes Dis (2022). doi: 10.1016/j.gendis.2021.12.017

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Zeng D, Zhou R, Yu Y, Luo Y, Zhang J, Sun H, et al. Gene expression profiles for a prognostic immunoscore in gastric cancer. Br J Surg (2018) 105:1338–48. doi: 10.1002/bjs.10871

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Jiang Y, Zhang Q, Hu Y, Li T, Yu J, Zhao L, et al. ImmunoScore signature: A prognostic and predictive tool in gastric cancer. Ann Surg (2018) 267:504–13. doi: 10.1097/SLA.0000000000002116

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Rosenberg JE, Hoffman-Censits J, Powles T, van der Heijden MS, Balar AV, Necchi A, et al. Atezolizumab in patients with locally advanced and metastatic urothelial carcinoma who have progressed following treatment with platinum-based chemotherapy: a single-arm, multicentre, phase 2 trial. Lancet (2016) 387:1909–20. doi: 10.1016/S0140-6736(16)00561-4

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Turley SJ, Cremasco V, Astarita JL. Immunological hallmarks of stromal cells in the tumour microenvironment. Nat Rev Immunol (2015) 15:669–82. doi: 10.1038/nri3902

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Fridman WH, Zitvogel L, Sautes-Fridman C, Kroemer G. The immune contexture in cancer prognosis and treatment. Nat Rev Clin Oncol (2017) 14:717–34. doi: 10.1038/nrclinonc.2017.101

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Wieder T, Eigentler T, Brenner E, Rocken M. Immune checkpoint blockade therapy. J Allergy Clin Immunol (2018) 142:1403–14. doi: 10.1016/j.jaci.2018.02.042

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics (2010) 26:1572–3. doi: 10.1093/bioinformatics/btq170

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Hanzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinf (2013) 14:7. doi: 10.1186/1471-2105-14-7

CrossRef Full Text | Google Scholar

22. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U.S.A. (2005) 102:15545–50. doi: 10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res (2015) 43:e47. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Sotiriou C, Wirapati P, Loi S, Harris A, Fox S, Smeds J, et al. Gene expression profiling in breast cancer: understanding the molecular basis of histologic grade to improve prognosis. J Natl Cancer Inst (2006) 98:262–72. doi: 10.1093/jnci/djj052

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, et al. TGFbeta attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature (2018) 554:544–8. doi: 10.1038/nature25501

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Geeleher P, Cox N, Huang RS. pRRophetic: an r package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One (2014) 9:e107468. doi: 10.1371/journal.pone.0107468

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Shulman Z, Stern-Ginossar N. The RNA modification n-6-methyladenosine as a novel regulator of the immune system. Nat Immunol (2020) 21:501–12. doi: 10.1038/s41590-020-0650-4

PubMed Abstract | CrossRef Full Text | Google Scholar

29. He XF, Tan LQ, Ni J, Shen GP. Expression pattern of m(6)A regulators is significantly correlated with malignancy and antitumor immune response of breast cancer. Cancer Gene Ther (2021) 28:188–96. doi: 10.1038/s41417-020-00208-1

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Chen H, Pan Y, Zhou Q, Liang C, Wong CC, Zhou Y, et al. METTL3 inhibits antitumor immunity by targeting m(6)A-BHLHE41-CXCL1/CXCR2 axis to promote colorectal cancer. Gastroenterology (2022) 163:891-907. doi: 10.1053/j.gastro.2022.06.024

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Ma S, Yan J, Barr T, Zhang J, Chen Z, Wang LS, et al. The RNA m6A reader YTHDF2 controls NK cell antitumor and antiviral immunity. J Exp Med (2021) 218:e20210279. doi: 10.1084/jem.20210279

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Han DL, Liu J, Chen CY, Dong LH, Liu Y, Chang RB, et al. Anti-tumour immunity controlled through mRNA m(6)A methylation and YTHDF1 in dendritic cells. Nature (2019) 566:270. doi: 10.1038/s41586-019-0916-x

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Liu Z, Zhong J, Zeng J, Duan X, Lu J, Sun X, et al. Characterization of the m6A-associated tumor immune microenvironment in prostate cancer to aid immunotherapy. Front Immunol (2021) 12:735170. doi: 10.3389/fimmu.2021.735170

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Li B, Zhu L, Lu C, Wang C, Wang H, Jin H, et al. circNDUFB2 inhibits non-small cell lung cancer progression via destabilizing IGF2BPs and activating anti-tumor immunity. Nat Commun (2021) 12:295. doi: 10.1038/s41467-020-20527-z

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Wang X, Ji Y, Feng P, Liu R, Li G, Zheng J, et al. The m6A reader IGF2BP2 regulates macrophage phenotypic activation and inflammatory diseases by stabilizing TSC1 and PPARgamma. Adv Sci (Weinh) (2021) 8:2100209. doi: 10.1002/advs.202100209

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Patil DP, Chen CK, Pickering BF, Chow A, Jackson C, Guttman M, et al. m(6)A RNA methylation promotes XIST-mediated transcriptional repression. Nature (2016) 537:369–73. doi: 10.1038/nature19342

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Cai X, Chen Y, Man D, Yang B, Feng X, Zhang D, et al. RBM15 promotes hepatocellular carcinoma progression by regulating N6-methyladenosine modification of YES1 mRNA in an IGF2BP1-dependent manner. Cell Death Discov (2021) 7:315. doi: 10.1038/s41420-021-00703-w

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Wang T, Kong S, Tao M, Ju S. The potential role of RNA N6-methyladenosine in cancer progression. Mol Cancer (2020) 19:88. doi: 10.1186/s12943-020-01204-7

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Chen DS, Mellman I. Elements of cancer immunity and the cancer-immune set point. Nature (2017) 541:321–30. doi: 10.1038/nature21349

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Du A, Wu X, Gao Y, Jiang B, Wang J, Zhang P, et al. m6A regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in acute myeloid leukemia. Front Immunol (2021) 12:789914. doi: 10.3389/fimmu.2021.789914

PubMed Abstract | CrossRef Full Text | Google Scholar

41. McGranahan N, Furness AJ, Rosenthal R, Ramskov S, Lyngaa R, Saini SK, et al. Clonal neoantigens elicit T cell immunoreactivity and sensitivity to immune checkpoint blockade. Science (2016) 351:1463–9. doi: 10.1126/science.aaf1490

PubMed Abstract | CrossRef Full Text | Google Scholar