Comprehensive Analysis Revealed the Potential Implications of m6A Regulators in Lung Adenocarcinoma

Background: The biological significance of RNA N6-methyladenosine (m6A) decoration in tumorigenicity and progression has been highlighted in recent studies, but whether m6A modification plays a potential role in tumor microenvironment (TME) formation and immune regulation in lung adenocarcinoma (LUAD) remains unclear.

Methods: m6A modification features were evaluated by analyzing the multi-omics features of 17 m6A regulators in over 1900 LUAD samples, and at the same time, the correlation between these modification patterns and TME characteristics was analyzed. An m6A score signature–based principal component analysis (PCA) algorithm was constructed to assess the prognosis and responses of individual patients to immunotherapeutic and targeted therapies.

Results: Three different m6A modification patterns were determined in 1901 LUAD samples, which were found to be related to diverse clinical outcomes via different biological pathways. Based on the m6A score extracted from the m6A-associated signature genes, LUAD patients were separated into high- and low-m6A score groups. It was discovered that patients with high m6A scores had longer survival, lower tumor mutation loads, and low PD-L1/PDCD1/CTLA4/TAG3 expression level. In addition, LUAD patients with high m6A scores displayed lower IC₅₀ to some targeted drugs, including nilotinib, erlotinib, imatinib, and lapatinib.

Conclusion: m6A modification was significantly associated with the TME and clinical outcomes. These findings may help gain more insights into the role of m6A decoration in the molecular mechanism of LUAD, thus facilitating the development of more effective personalized treatment strategies.

Introduction

Lung cancer is the most prevalent malignancy of all cancers, with top incidence and mortality rate in the world (Siegel et al., 2020). Based on the pathologic type, lung cancer can be separated into non–small cell lung cancer (NSCLC) and small cell lung cancer (SCLC). In addition, lung adenocarcinoma (LUAD) remains the most common subtype of NSCLC (Chen et al., 2019a). The development of LUAD is a long-term complicated process involving multiple steps, including interactions of multiple genes with external factors. Despite remarkable progress in the treatment of LUAD with respect to surgery, radiotherapy, chemotherapy, and immunotherapy, the prognosis of LUAD patients is still unsatisfactory due to undesirable response to treatment in addition to invasion and migration of tumor cells. Therefore, it is essential to further comprehend the underlying molecular mechanisms of the occurrence and progression of LUAD and find some clinically effective diagnostic and prognostic biomarkers for the sake of developing more promising individualized treatment strategies for LUAD.

N6 adenosine (m6A) methylation is extensively present in mRNAs, lncRNAs, and miRNAs. It is the most prevalent RNA modification type and plays an important role in various physiopathological processes (Zhao et al., 2017; Zaccara et al., 2019; Gu et al., 2020). m6A decoration is loaded with a reversible and dynamic process, which is regulated by various types of regulatory factors: the methyltransferases (writer), binding proteins (reader), and degrees (eraser) (Wang et al., 2020a). It is necessary to study these regulatory proteins and gain a better understanding about the mechanisms of m6A modification in gene regulation since these proteins have a great impact on m6A modification (Li et al., 2019a; Chen et al., 2020). Ample evidence has shown that abnormal expression and gene variation of m6A regulatory factors are related to the progression of malignant tumors and abnormal immune regulation (Chen et al., 2019b; Wang et al., 2020b; Shulman and Stern-Ginossar, 2020; Zhao et al., 2020). A comprehensive analysis of genetic variations and expression interference behind the heterogeneity of LUAD would promote the discovery of novel biomarkers and therapeutic targets based on m6A modification (Chen et al., 2019b; Li et al., 2019b; Wang et al., 2020b; Shulman and Stern-Ginossar, 2020; Zhao et al., 2020). Previous studies have reported abnormal expression patterns of m6A regulatory factors in LUAD (Liu et al., 2020a; Zhang et al., 2020a; Li et al., 2021a). Zhou et al. analyzed 21 potential m6A regulators involved in the tumor immune microenvironment of LUAD, based on which they constructed a risk signature and used it to define the tendency of immune cell infiltration in LUAD and predict the prognosis of LUAD patients. They found that m6A regulators played a critical role in the tumor immune microenvironment (TME) of LUAD (Zhou et al., 2021). Another study revealed that LUAD patients with high risk scores had poorer survival rates. The analysis of univariate and multivariate Cox regression showed that m6A-RPS was an independent prognostic risk element, and a significant association was observed between the expression of m6A-RPS and m6A modulators (Sun et al., 2021). More studies have demonstrated that m6A-related genes are efficacious biomarkers for early diagnosis, prognostic prediction, and efficacy evaluation of LUAD (Zhang et al., 2020a; Zhuang et al., 2020; Li et al., 2021a; Zhang et al., 2021). The m6A reader YTHDC2 suppressed LUAD tumorigenesis in vivo and in vitro by directly targeting the solute carrier 7A11 (SLC7A11) (Ma et al., 2021). In addition, LUAD tumor growth was effectively inhibited via targeting the METTL3-dependent m6A decoration of FBXW7 (Wu et al., 2021a).

Previous research studies have concentrated on single or a few m6A-related genes and included only a small sample size. However, the occurrence and progression of LUAD is a long-term gradually evolutionary process of numerous genes and multiple steps and methylation decoration patterns. The clinical significance of the tumor microenvironment (TME) landscape and m6A regulatory factors in LUAD remains largely unknown. The purpose of our present study was to systematically evaluate the multi-omics features of 17 m6A regulators and the m6A decoration patterns by integrating the genomic and transcriptomic data of more than 1900 LUAD samples to see whether the m6A decorative pattern was associated with the TME characteristics. In addition, based on m6A regulators and related genes, we built a score system for quantifying m6A decoration patterns in individual tumors and forecasting the clinical response of LUAD patients to immune-targeted therapies. Moreover, we developed a nomogram based on the scoring system (m6Ascore signature) and some important clinical trials to predict the prognosis of LUAD patients. The results obtained in our study suggested that m6A modification played a pivotal role in forming various TME profiles which could serve as a guide for planning therapeutic intervention schemes for LUAD.

Methods

Collection and Pretreatment of Publicly Available Expression Data Sets

Clinical features of LUAD and gene expression data were collected retrospectively from TCGA (https://cancergenome.nih.gov/) and NCBI GEO databases (https://www.ncbi.nlm.nih.gov/geo/). Twenty-three m6A regulatory factors were obtained (Supplementary Table S1) (Chen et al., 2019b; Li et al., 2021b; Zhou et al., 2020; Liu et al., 2020b; Arguello et al., 2017; Warda et al., 2017; Barros-Silva et al., 2020; Wu et al., 2021b). TCGA somatic mutation data had been captured through the utilization of TCGAbiolinks R package (Colaprico et al., 2016) and visualized through application of the maftools R package (Mayakonda et al., 2018). The copy number variation (CNV) datasets were obtained from the Xena Public Data Center (https://xenabrowser.net). The copy number variation landscape of these 23 m6A regulators in human chromosomes was identified using R package “Rcircos”. TCGA–LUAD RNA-seq data (FPKM format) were downloaded from the Genomic Data Commons (GDC, https://portal.gdc.cancer.gov/) and converted into transcripts per kilobase million (TPM) format. The Gene Expression Omnibus (GEO) database was used to make a comprehensive query on all qualified LUAD data sets. A total of seven datasets [GSE30219 (Xie et al., 2011), GSE30219 (Rousseaux et al., 2013), GSE31210 (Okayama et al., 2012), GSE37745 (Botling et al., 2013), GSE50081 (Der et al., 2014), GSE68465 (Shedden et al., 2008), and GSE72094 (Schabath et al., 2016)] were obtained to represent different LUAD-independent studies. All datasets contained clinical data and survival information. The platform files and survival data of these datasets are shown in Supplementary Table S2. In accordance with the corresponding annotation file, the probe was transformed into a gene symbol. For genes with set signals of multiple probes, their values were averaged to produce a single gene expression value. At last, they were further combined to a meta-queue through the “ComBat” algorithm using the “sva” package (Leek et al., 2012) to decrease the batch effect from abiotic deviations.

Unsupervised Clustering for 17 m6A Regulators

By integrating the eight datasets, 17 out of the 23 m6A modulators were used to identify different m6A decorative patterns mediated by the m6A modulators. They included five writers (RBM15, METTL3, WTAP, RBM15B, and ZC3H13), one eraser (FTO), and 11 readers (YTHDC2, YTHDC1, YTHDF2, YTHDF1, YTHDF3, IGF2BP2, IGF2BP3, HNRNPC, HNRNPA2B1, LRPPRC, and FMR1). Based on the expressions of the 17 m6A regulatory factors, we performed an unsupervised cluster analysis to discover different m6A decoration patterns, and the patients were categorized for subsequent analysis. The abovementioned steps were performed by applying the “ConsensusClusterPlus” software package (Wilkerson and Hayes, 2010) and repeated 1000 times to ensure stability of the classification.

Gene Set Variation Analysis and Functional Annotation

To investigate the discrepancy in biological processes between the m6A decorative patterns, GSVA enrichment analysis was performed by using “GSVA” R package knowing that GSVA is an unsupervised and nonparametric method usually used to appraise changes in pathway and biological process activities in expression datasets (Hänzelmann et al., 2013). The gene sets of “c2.cp.kegg.v7.2.symbols” were downloaded from MSigDB (http://www.gsea-msigdb.org/gsea/msigdb) to be used for analysis of GSVA. After adjustment, p < 0.05 indicated statistical significance. Functional annotation of m6A-associated genes was analyzed by using the “clusterProfiler” R package (Wu et al., 2021c), and the cutoff value of the FDR was <0.05.

Estimation of the TME

The level of infiltration of diverse immune cells, including the TME and adaptation of innate immune cell types, was quantified by single sample GSEA (ssGSEA) (Barbie et al., 2009; Charoentong et al., 2017). Immune cell and stromal cell score in malignant tumors were estimated by the expression data (estimation) algorithm (Becht et al., 2016) using the distinct properties of the transcription profile to deduce the nature of the tumor cell and purity of the tumor. The infiltration level of stromal and immune cells was predicted by the score of stromal and immune cells calculated by the ESTIMATE algorithm and used as the basis for inferring tumor purity. LUAD tissue with enriched infiltration of immune cells indicated a better immune score and poorer tumor purity.

Identification of Differentially Expressed Genes Between Distinct m6A Phenotypes

To find m6A-associated genes, patients were categorized into three different m6A decoration patterns according to the expression of the 17 m6A regulatory factors. The deg between different decorative patterns was demonstrated by the empirical Bayesian approach of limma R package (Ritchie et al., 2015) p value <0.001 was defined as the significant standard to determine DEGs.

Construction of the m6A-Related Gene Signature (m6A score)

To determine the m6A decoration pattern of individual tumors, we constructed a scoring system to assess the m6A decoration pattern of individual LUAD patients. The genetic marker associated with m6A is called m6Ascore. Specifically, overlapping DEGs discovered from different m6A clusters were standardized and extracted to facilitate analyzing the prognosis of each gene using the univariate Cox regression model. Genes with vital prognosis were collected for principal component analysis (PCA) to construct m6A-associated gene markers. Principal components 1 and 2 were regarded as signature fractions. The merit of this approach is that the score is focused on the largest correlated (or anti-correlated) gene block in the set, while the contribution of genes that are not interacted with other members of the set is underweighted. Then, we use a similar formula from previous research to determine m6Ascore (Zhang et al., 2020b):

$\text{m}6\text{Ascore}={\displaystyle \sum \left(\text{PC}1i+\text{PC}2i\right)},$

where i is the expression of terminally determined genes associated with the m6A phenotype.

Gene Set Enrichment Analysis and Functional Annotation

Using the clusterProfiler R software package, GSEA was executed to deduce biological processes associated with m6A regulators, and we regarded p value <0.05 as statistically significant.

Estimation of Drug Sensitivity

Drug sensitivity was evaluated by using the pRRophetic R software package (Geeleher et al., 2014a) and defined by the concentration essential for 50% inhibition of cell growth (IC₅₀) (Geeleher et al., 2014b).

Statistical Analysis

The limma R package was used for DEG expression analysis, and Spearman’s method was used for correlation analysis. The statistical difference between two groups was calculated by using the Wilcoxon rank sum test. The Kruskal–Wallis test was used for comparing more than two groups. The univariate Cox regression model was used to estimate the risk ratio (HR) of m6A regulatory factors. The multivariate Cox regression model was used to determine the independent prognostic factors. The correlation between the m6A decoration pattern and prognosis with the“survminer” package in R was analyzed by the Cox proportional hazard model and Kaplan–Meier survival analysis. Using surv-cutpoint function from the “survival” software package, samples were stratified into high- and low-m6Ascore groups. The mutation landscape of m6Ascore patients in the TCGA–LUAD cohort with high or low subgroups was presented by the waterfall function of the maftools software package. Patients with complete clinical data were enrolled in ultimate multivariate Cox analysis. The results were visualized using the “forestplot” software package in R. The receiver operating characteristic (ROC) curve was used to assess the prognosis classified performance and other clinical factors of m6Ascore, and the area under the curve (AUC) was determined by using the “survivalROC” package. Using the clinical characteristics and m6Ascore signature as input, multivariable Cox regression analysis and variable selection were used to determine the strong combination of these predictors. Then, we built up a quantitative line chart with the “rms” package in R to predict the individual three- and Figure 1A five-year survival probabilities. To assess the prediction performance of the nomogram, we used the bootstraps method to calculate the C-index. At the same time, the calibration curves, with the Hosmer–Lemeshow test, were also used to judge the consistency between the model prediction values and actual results. All statistical p values were two-tailed, and p < 0.05 was considered as statistical significant. The data processing was accomplished in R 3.6.2 software.

FIGURE 1

FIGURE 1. Landscape of genetic and expression variation of m6A regulators in LUAD. (A) Metascape enrichment network and heat map of enrichment items in the input gene list colored with p value. (B) Determination of transcription factors related to 23 m6A regulatory factors by Pearson’s correlation analysis (|Cor| > 0.5 and p < 0.001). (C) Mutation frequency of 23 m6A regulatory factors in 561 LUAD patients in TCGA–LUAD cohort. Each column represents an individual patient. The abovementioned barplot is TMB, and the number on the right represents the mutation frequency of each regulatory factor. The barplot on the right shows the proportion of each variation type. The stacked barplot below shows the transformation section in each example. (D) CNV variation frequency of m6A regulatory factors in TCGA–LUAD cohort. The height of the column represents the frequency of change. Delete frequency, green dot; amplified frequency, red dot. (E) Location of CNV change of m6A regulatory genes on 23 chromosomes by TCGA–LUAD cohort. (F) Expression of 23 m6A regulatory factors in normal and LUAD tissues. Tumors, red; normal, blue. The upper and lower ends of the box represent the quartile range of values. The lines in the box represent intermediate values, and the real points represent outliers. The asterisk represents the statistical p value. (*p < 0.05; **p < 0.01; ***p < 0.001).

Results

Genetic Variation Landscape of m6A Regulators in LUAD

Twenty-three m6A regulators including two erasers, eight writers, and 13 readers (“erasers”: ALKBH5 and FTO; “writers”: METTL14, METTL3, METTL16, WTAP, VIRMA, ZC3H13, RBM15B, RBM15; and “readers”: YTHDC1, YTHDF1, YTHDC2, YTHDF2, FMR1, YTHDF3, HNRNPC, LRPPRC, RBMX, IGF2BP1, HNRNPA2B1, IGF2BP2, and IGF2BP3) were initially obtained, and by integrating multiple datasets, we finally included 17 m6A regulators including five writers (METTL3, RBM15, RBM15B, WTAP, and ZC3H13), one eraser (FTO), and 11 readers (YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, IGF2BP2, IGF2BP3, HNRNPA2B1, HNRNPC, FMR1, and LRPPRC) in this study. GO enrichment and Metascape analyses were performed on the 17 m6A regulatory factors, and the biological processes of valuable enrichment are concluded in Figure 1A. The transcription factor (TF) subset was from Cistrome (http://cistrome.org/). We, thus, obtained the corresponding TF expression values in TCGA database. Through Pearson’s correlation analysis (|Cor| > 0.5 and p < 0.001), we identified the TFs related to the 17 m6A regulators (Supplementary Table S3). As indicated in Figure 1B, most TFs were involved in the expression of YTHDC1 (n = 41), followed by ZC3H13 (n = 27). Moreover, we concluded the incidence of CNVs and somatic cell mutations of the 17 m6A regulators in LUAD. m6A regulator mutation was observed in 113 (20.14%) of the 561 samples, with the mutation frequency of ZC3H13 being the highest, followed by FMR1, RBM15, YTHDC2, LRPPRC, and YTHDC1, although METTL3 showed no mutations in LUAD samples (Figure 1C). Taking the relatively high mutation frequency of the “writer” gene ZC3H13 into account, further exploration suggested that mutations in ZC3H13 did not affect the expression of other m6A regulators (Supplementary Figure S1). The subsequent analysis of 17 m6A regulators demonstrated that mutations in CNV are widespread. YTHDF1, FMR1, IGF2BP2, METTL3, HNRNPC, IGF2BP3, YTHDF3, IGF2BP1, HNRNPA2B1, LRPPRC, YTHDC1, and FTO extensive CNV amplification was shown. In sharp contrast to this, YTHDF2, WTAP, YTHDC2, ZC3H13, RBM15, and RBM15B displayed general CNV deficiency (Figure 1D). Figure 1E showed the CNV change positions of the 17 m6A regulatory genes on chromosomes. To determine if the aforementioned genetic variants affected the m6A regulatory factor expression in LUAD patients, a difference in the mRNA expression of regulators between normal and LUAD samples was examined in TCGA–LUAD cohort (Figure 1F), and it was found that CNV changes were probably the primary factor interfering with the expression of m6A regulators. To identify the correlation between these m6A regulatory factors and prognosis of LUAD patients, we used Cox regression analysis. The results revealed that METTL3, YTHDC1, YTHDC2, and FTO were protective factors significantly correlated with longer overall survival rates, while WTAP, RBM15, HNRNPC, LRPPRC, IGF2BP2, and IGF2BP3 were risk factors (Supplementary Figure S2A). Of them, METTL3 and HNRNPC were independent prognostic genes for LUAD (Supplementary Figure S2B). The abovementioned analysis revealed a highly heterogenous pattern of altered genetics and expression in the m6A regulators among normal and LUAD samples, suggesting that imbalance in the expression of m6A regulators plays a key role in the development and progression of LUAD.

m6A Methylation Decoration Patterns Mediated by 17 Regulators

Eight data sets containing operating system data and clinical information (TCGA–LUAD, GSE29013, GSE30219, GSE31210, GSE37745, GSE50081, GSE68465, and GSE72094) were included into one meta-cohort. To learn the comprehensive landscape of m6A regulator interactions, the m6A regulatory network was used to define the prognostic significance of connections between the m6A regulators in LUAD patients (Figure 2A). The results showed that the crosstalk between the regulators, readers, and erasers played a key role in the formation of various m6A decorative patterns and were related to the occurrence and development of cancer. Based on these findings, we used the “ConsensusClusterPlus” package in R to categorize the patients qualitatively into different m6A decoration modes according to the expression of the 17 m6A regulators. Finally, three different decoration modes were determined based on unsupervised clustering. We named them as m6A cluster A (n = 909), B (n = 582), and C (n = 410) (Figure 2B). Among them, the m6A cluster-A showed an obvious survival advantage, while m6A cluster-C showed the worst prognosis (Figure 2B). In addition, we also noticed that the expression levels of these m6A regulators differed significantly between the different m6A decoration patterns. WTAP, RBM15, HNRNPC, LRPPRC, IGF2BP2, and IGF2BP3 were elevated markedly in the m6A cluster-C subtype, while YTHDC1, YTHDC2, and FTO were increased in the m6A cluster-A subtype (Figure 3A).

FIGURE 2

FIGURE 2. m6A methylation modification pattern and biological characteristics of each pattern. (A) Interaction of m6A regulators in LUAD. The circle size represents the impact of each regulatory factor on prognosis, and the range of calculated value by log-rank test was p < 0.0001, p < 0.001, p < 0.01, p < 0.05, and p < 1. Left half of the circle: purple represents prognostic risk factors and green represents prognostic favorable factors. Right half of the circle: the types of m6A regulators. The line connecting regulatory factors represents the interaction between regulatory factors, and the thickness represents the correlation intensity between regulatory factors, showing a negative association with blue and a positive association with pink. (B) Survival analysis based on three m6A patterns in 1901 patients from seven GEO cohorts (GSE30219, GSE31210, GSE29013, GSE37745, GSE68465, GSE50081, and GSE72094) and one TCGA–LUAD cohort, including 909 cases in m6Acluster A, 582 in m6Acluster B, and 410 in m6Acluster C. Kaplan–Meier curves with log-rank p value less than 0.001 show significant survival differences in survival rates between the three m6A modification modes. The overall survival of m6Acluster B was significantly better than that of the other two m6Aclusters. (C,D) Activation state of biological pathways in different m6A modification patterns as shown by GSVA enrichment analysis. The biological processes were visualized by using heat maps, with red representing activated pathways and blue representing inhibited pathways. LUAD queue was used as sample annotations. C: m6Acluster C vs. m6Acluster A; D: m6Acluster C vs. m6Acluster B.

FIGURE 3

FIGURE 3. Relationship between the three m6A modification patterns and expression of m6A regulators and TME landscape. (A) Expression difference of 17 m6A regulators in the three m6A modification patterns. m6Acluster-A, blue; m6Acluster-B, yellow; m6Acluster-C, red. The top and bottom of the box represented the range of quartiles of values. The lines in the boxes represented median values, and the solid dots represented outliers. The asterisk indicated the statistical p value (*p < 0.05; **p < 0.01; ***p < 0.001). (B) Expression difference of PD-L1 (CD274) in the three m6A modification patterns in the GSE72094 cohort. (C) Expression difference of PD-L1 (CD274) in the three m6A modification patterns in TCGA cohort. (D) Difference of TMB in the three m6A modification patterns in TCGA–LUAD cohort. (E) Comparison of immune scores across the three m6A modification patterns. (F) Comparison of tumor purity across the three m6A patterns. (G) TME landscape of the three m6A modification patterns. Tumor purity, immune score, ESTIMATE score, stroma score, and m6A clustering were noted above. Red indicated a high degree of infiltration of TME infiltrating cells, and blue indicated a low degree of infiltration of TME infiltrating cells. (H) Abundance of each TME infiltrating cell in the three m6A decoration modes. The upper and lower ends of the box represented the quartile range of values. The lines in the box represented intermediate values, and the real points represented outliers. The asterisk represented the statistical p value. (***p < 0.001; **p < 0.01; *p < 0.05). m6Acluster-A, blue; m6Acluster-B, yellow; m6Acluster-C, red. For comparisons of the three groups, the Kruskal–Wallis test was used.

TME Landscape in Distinct m6A Decoration Patterns

With the aim of investigating the biological behavior between these different m6A decoration patterns, we executed an analysis of GSVA enrichment and found that m6A cluster-A exhibited enrichment pathways related to complete immune activation, which included drug metabolism cytochrome P450, fatty acid metabolism, tyrosine metabolism, and histidine metabolism (Figures 2C,D). Knowing that tumor mutation burden (TMB) and PD-L1 are mature biomarkers for predicting anti–PD-1/L1 treatment response, we compared the TMB and expression level of PD-L1 in the three different m6A decoration clusters, finding that PD-L1 (Figures 3B,C) and TMB (Figure 3D) were obviously downregulated in the m6A cluster-A subtype. In addition, we quantified the overall infiltration (immune score) and tumor cell purity (tumor purity) of immune cells under the three decoration modes using the ESTIMATE algorithm and found that m6Acluster-A and -B had the higher immune score than m6Acluster-C (Figure 3E). Conversely, m6A cluster-C had higher tumor purity than m6A cluster-A and m6A cluster-B (Figure 3F). Subsequently, we performed a ssGSEA analysis to explore the relative abundances of 23 immune infiltration cells in the three subtypes (Figures 3G,H). Antitumor lymphocyte cell sub-populations, such as activated CD8+ T cells, NK cells, and NK T cells, were mainly concentrated in the m6A cluster-B subtype. Among them, CD56 dim natural killer cell, immature B cell, MDSC, macrophage, monocyte, regulatory T cell, type 1 T helper cell, and type 17 T-helper cells were also significantly infiltrated in m6A cluster-B tumors. These results suggested that m6A cluster-B tumors were embraced by more nontumor components, such as stromal and immune cells.

m6A Phenotype–Associated DEGs in LUAD

Despite the classification of LUAD patients into three m6A-modified phenotypes based on a consensus clustering algorithm for m6A regulator expression, the underlying genetic changes and expression disorders in these phenotypes remained unclear. Therefore, we further observed changes in the potential transcriptional expression of the three m6A decoration patterns in LUAD. 372 overlapping DEGs between the three m6A decoration patterns were determined by the empirical Bayesian method (Figure 4A). Next, univariate Cox analysis was executed for these 372 DEGs to determine the genes associated with prognosis. Finally, 303 m6A phenotype–related DEGs associated with prognosis were regarded as m6A-associated signature genes (Supplementary Table S4). The GO analysis of the enrichment for these characteristic genes showed significant changes in biological processes associated with cell cycle (Figure 4B and Supplementary Table S5). Based on the 303 most representative genes associated with the m6A phenotype, we executed unsupervised consensus cluster analysis and acquired three stable transcriptome phenotypes (m6A geneCluster-A, m6A geneCluster-B, and m6A geneCluster-C) (Supplementary Figure.S2E). The results showed that patients in geneCluster-B displayed the worst prognosis (Figure 4C) and the highest PD-L1 expression (Figures 4D,E) and TMB (Figure 4F). As shown in Figures 4D–F, the patients in geneCluster-C had the lowest PD-L1 expression and TMB but the longest overall survival. Moreover, surprisingly, the prognostic risk genes (RBM15, HNRNPC, LRPPRC, IGF2BP2, and IGF2BP3) were all poorly expressed in geneCluster-C, and the prognostic-friendly genes (YTHDC1, FTO, and YTHDC2) were highly expressed (Supplementary Figure S2G). Similarly, the prognostic risk gene (WTAP) was highly expressed in geneCluster-B, while the prognostic-friendly gene (METTL3) was less expressed (Supplementary Figure S2G). The prognostic differences between geneCluster-A, B, and C and the corresponding expressions of YTHDC1, YTHDC2, RBM15, HNRNPC, LRPPRC, IGF2BP2, IGF2BP3, METTL3, and WTAP were consistent with the results of Cox regression analysis (Supplementary Figures S2A,B). All these results indicated that the three stable transcriptomic phenotypes (geneCluster-A, B, and C) could well-distinguish between different prognostic populations.

FIGURE 4

FIGURE 4. m6A phenotype–related DEGs and functional annotation. (A) Venn diagram of 372 m6A-related DEGs between the three m6A decoration patterns. (B) Functional annotation for phenotype-related DEGs associated with the prognosis using GO enrichment analysis. (C) Estimation of the survival curve of m6A phenotype–related gene markers by Kaplan–Meier plotter (p < 0.001, log-rank test). (D) Expression difference of PD-L1 (CD274) in the three geneclusters in the GSE72094 cohort. (E) Expression difference of PD-L1 (CD274) in the three gene Clusters in TCGA cohort. (F) Difference in TMB between the three gene clusters in TCGA–LUAD cohort. (G) Overlapping m6A phenotype–related DEGs associated with prognosis as unsupervised clustering, and the patients were divided into different genomic subtypes defined as gene clustering A–C. Gene signature subtypes, m6A clusters, items, age, gender, smoking status, and life status were used as patient annotations.

Construction of the m6Ascore Signature and Exploration of its Clinical Relevance

While we discovered that m6A modifications could be critical in regulating prognosis and immune infiltration, correlation analysis was conducted only on the basis of patient groups and was unable to accurately forecast the m6A methylation decoration pattern in a single tumor. Considering the individual heterogeneity and complexity of m6A decoration, we constructed a score system to quantify the m6A modification patterns of individual LUAD patients according to these phenotypically associated genes, which we called m6Ascore. Next, we tried to further define the value of m6Ascore in forecasting the prognosis of LUAD patients. We determined the critical value −1.244678 by the measurement package in R and separated the patients into a low-m6Ascore group and a high-m6Ascore group. Patients in the high-m6Ascore group showed significant survival superiority. The median OS was twice as long as that of patients in the high-m6Ascore group (8.4 vs. 4.1) (Figure 5A). We expressed the attribute changes of individual patients in an alluvial map (Figure 5B). Our subgroup analysis on different groups by age, sex, stage, smoking status, and EGFR/Kras/TP53/STK11 mutation status showed that the prognosis of patients with high m6Ascore was significantly better than that of patients with low m6Ascore (age ≤ 65 subset, age > 65 subset, female subset, male subset, ever smoking subset, never smoking subset, stage I/II subset, stage III/IV stage, EGFR wt subset, STK11 wt subset, Kras wt/mut subset, and TP53 wt/mut subset), which further showed that the predictive performance of m6Ascore was reliable in various subgroups of LUAD patients (Supplementary Figure S3). Subsequently, a further comparison of m6Ascore differences between different subgroups revealed that these subsets (female subset, never smoking subset, EGFR-Mut subset, Kras-WT subset, TP53-WT subset, and alive subset) had higher m6Ascore (Supplementary Figure S4). Given that m6Ascore could well-predict patient prognosis, we further compared m6Ascores between distinct m6A modification patterns. m6Acluster-A was found to have the highest m6Ascore, followed by m6Acluster-B and m6Acluster-C with the lowest m6Ascore (Figure 5C). This simply demonstrated the accuracy of m6A modification patterns in distinguishing various prognostic patients (Figure 2B). Similarly, Figure 5D also demonstrated the accuracy of the geneCluster in the classification of patients with different prognoses (Figure 4C). Furthermore, Spearman’s analysis was used to detect the relationship between the relative abundance of 23 immune infiltrating cells and m6Ascore. The correlation matrix heat map showed a significantly negative correlation of m6Ascore with activated CD4 T cells and a significantly positive connection with mast cells (Figure 5E).

FIGURE 5

FIGURE 5. Construction of the m6Ascore signature and tumor somatic mutation. (A) Survival analysis of patients with low m6Ascore (813 cases) and high m6Ascore (1088 cases) by Kaplan–Meier curve (p < 0.001, log-rank test). (B) Alluvial map showing changes in m6Aclusters, geneCluster, m6Ascore, and status. (C) Comparison of m6Ascore across the three m6A modification patterns. (D) Comparison of m6Ascore across three geneClusters. (E) Spearman’s analysis of correlations between m6Ascore and TME infiltrating cells. Blue represented a negative association and red represented a positive association. (F,G) Waterfall diagram of tumor somatic mutation established by high m6Ascore (F) and low m6Ascore (G). Each column indicated an individual patient. The figure above showed TMB, and the number on the right represented the mutation frequency of each gene. The bar graph on the right showed the proportion of each variation type. (H) Correlations between TMB, m6Ascore, and geneCluster using Spearman’s analysis. (I) Comparison of TMB between high- and low-m6Ascore groups. (J) Kaplan–Meier curves for survival analysis of patients with low (389) and high (91) TMB in TCGA–LUAD cohort (p = 0.061, log-rank test). (K) Kaplan–Meier curve analysis of survival in “TMB-H + m6Ascore-H”, “TMB-H + m6Ascore-L”, “TMB-L + m6Ascore-H”, and “TMB-L + m6Ascore-L” groups in TCGA–LUAD cohort (p < 0.001, Log-rank test). TMB-H, high TMB; TMB-L, low TMB; m6Ascore-H, high m6Ascore; m6Ascore-L, low m6Ascore.

The Association Between m6Ascore and TME Cell Infiltration and Functional Annotation

ssGSEA analysis showed that low-m6Ascore tumors were significantly infiltrated by CD4 T, CD8 T, activated B, CD56 bright NK cell, gamma delta T cell, MDSC, macrophage, immature B cell, NK cell, regulatory T cell, and type 1 T helper cell, while high-m6Ascore tumors were significantly infiltrated by eosinophil, mast cell, monocyte, plasmacytoid dendritic cell, and T follicular helper cell (Supplementary Figure S5A). The ESTIMATE algorithm was used to quantify the tumor purity and immune score of each sample. The Wilcoxon rank sum test revealed significant differences in immune score and tumor purity between high- and low-m6Ascore groups, indicating that tumor purity was higher (Supplementary Figure S5C) and the immune score was lower (Supplementary Figure S5B) in the high-m6Ascore group than the low-m6Ascore group. The GSEA analysis of the two groups showed that DNA replication, p53 signaling pathway, cell cycle, IL-17 signaling pathway, and progesterone-mediated oocyte maturation were markedly enriched in low-m6Ascore tumors (Supplementary Figure S5D). GSVA enrichment analysis was then conducted to investigate the biological behavior between m6Ascore groups.

Tumor Somatic Mutation and the Role of m6Ascore in Predicting Immunotherapeutic and Targeted Therapy Benefits

Then, we analyzed the differences of somatic mutation distributions in TCGA–LUAD cohort between high- and low-m6Ascore groups by maftools R package, and we found that the tumor mutation burden was more extensive in the low-m6Ascore group than in the high-m6Ascore group. The mutation rate of the fifth most significant mutation gene was 28% in the high-m6Ascore group vs. 41% in the low-m6Ascore group (Figures 5F,G). Quantitative analysis of TMB demonstrated that tumors with low m6Ascore were significantly associated with high TMB (Figure 5I), showing a significant negative correlation between m6Ascore and TMB (Figure 5H). We further explored the correlation between TMB and m6Ascore and the prognosis of patients and found that low TMB seemed to have a worse prognosis (Figure 5J). Notably, when TMB and m6Ascore were both low, the prognosis of the patients was the worst (Figure 5K). These findings may help better understand the impact of m6Ascore classification on genomic variation and revealed the potential complex interactions between single somatic mutations, m6A decoration, and patient outcomes. The Wilcoxon rank sum test showed that the expression of immune checkpoints was obviously different (PD-L1, PDCD1, CTLA4, LAG3) between the high- and low-m6Ascore groups. The low-m6Ascore group showed higher immune checkpoint expression (Figures 6A–D). Ample evidence suggested that patients at a high TMB status and high immune checkpoint expression exhibited sustained clinical reactions to anti–PD-1/PD-L1 immunotherapy. Hence, the results mentioned above may not only indirectly prove that variations in the tumor m6A modification mode are key factors mediating the clinical response to anti–PD-1/PD-L1 immunotherapy but also confirm the value of m6Ascore in forecasting the outcome of immunotherapy. As targeted therapy has been widely used in LUAD therapy, it is important to identify subgroups of patients who may be more sensitive to some drugs. Here, we predicted the response of patients in the high- and low-m6Ascore groups to several commonly used drugs. As shown in Figure 6E, patients with high m6Ascore were more sensitive to nilotinib, erlotinib, imatinib, and lapatinib in the LUAD cohort, which may provide useful evidence for guiding personalized treatment strategies for LUAD patients.

FIGURE 6

FIGURE 6. Potential of m6Ascore in predicting therapeutic response. (A–D) Wilcoxon rank sum test shows significant differences in the expression of PD-L1 (A), PDCD1 (B), CTLA4 (C) and LAG3 (D) between low- and high-m6Ascore groups. (E) Comparison of drug sensitivity between high- and low-m6Ascore groups. Distribution of IC₅₀ of nilotinib, imatinib, erlotinib, and lapatinib between high- and low-m6Ascore groups.

Construction and Evaluation of the Nomogram

In view of the value of m6Ascore in judging the prognosis of patients, we further used the Cox proportional hazard model to explore the independent prognostic factors of LUAD patients. As shown in Supplementary Figures S6A,B, age, smoking status, stage, and m6Ascore were all independent prognostic factors. Based on these findings, we constructed a nomogram to verify the role of various factors in forecasting the prognosis of patients (Supplementary Figure S6C). Based on the results of multivariate Cox analysis, the score was assigned to each factor in the rosette, and the total rosette score was obtained from the sum of the individual scores of all predictors. Using the total score, the 3- and 5-year survival rates of patients could be estimated by predicting the total score downward. Compared with the ideal model, the 3- and 5-year calibration curves also revealed good consistency, which further indicated that the nomogram was stable in predicting outcomes in patients with LUAD (Supplementary Figure S6D).

Discussion

As a dynamic and reversible process regulated by m6A regulatory factors, m6A modification promotes or suppresses malignant behavior mainly by modulating the expression of targeted oncogenes or oncogenes (Li et al., 2017; Ma et al., 2017; Bai et al., 2019). Growing evidence has clarified that m6A modification plays critically important roles in innate immune, inflammatory, and antitumor effects through interplay with different m6A regulators (Shulman and Stern-Ginossar, 2020; Zeng et al., 2020). To expand our understanding about the role of m6A modification in tumorigenesis, provide valuable biomarkers for diagnostic and prognostic evaluation, and propose new therapeutic targets, it is necessary to gain a better understanding about the molecular and biological characteristics of m6A regulatory factors.

It was found in our study that the mutation rate of m6A regulators was relatively low in LUAD, ranging from 0 to 3%, among which the mutation rate of ZC3H13 was the highest (3%), and this mutation status did not affect the expression of other m6A regulators, indicating that the genetic distortion of m6A regulators may be more than the state change of a single gene. Through the study of the CNV characteristics of m6A regulators, we found that YTHDF1, FMR1, IGF2BP2, METTL3, HNRNPC, IGF2BP3, YTHDF3, IGF2BP1, HNRNPA2B1, LRPPRC, YTHDC1, and FTO had the advantage of increasing the CNV, while YTHDF2, WTAP, YTHDC2, RBM15, ZC3H13, and RBM15B had the advantage of reducing the CNV. Subsequent differential analysis revealed that these m6A regulators were differentially expressed among lung tumors and normal tissues. Relating gene CNV to RNA expression, we found that CNV-amplified m6A regulators (METTL3, YTHDF1, HNRNPC, LRPPRC, HNRNPA2B1, IGF2BP1, IGF2BP2, and IGF2BP3) were highly expressed in LUAD tissues as compared with those in normal lung tissues, and vice versa for FTO, suggesting that CNV change may be a significant factor leading to perturbations of m6A regulator expressions.

Based on the expression characteristics of m6A regulators, we identified three subtypes with different TME characteristics and prognosis. Our study found that patients with m6Acuster-A had a significant survival advantage, while m6Acluster-C subtype had the worst prognosis (Figure 2B). To further explore whether m6A regulators play a role in prognosis among different patients, we analyzed the gene expression patterns of m6A regulators among different m6Aclusters. The results revealed that the expression of these m6A regulators differed significantly in different m6A decoration patterns. The expression levels of prognostic risk genes (WTAP, RBM15, HNRNPC, LRPPRC, IGF2BP2, and IGF2BP3) were significantly increased in the m6Acluster-C subtype, while the prognostic-friendly genes (YTHDC1, YTHDC2, and FTO) were apparently increased in the m6Acluster-A subtype. This suggested that the perturbation of m6A regulator expression had a great impact on patient prognosis, thereby distinguishing m6A regulatory patterns with different characteristics. Geeleher et al. (2014a) found that TME structure played a key role in tumor progression and affected the immunotherapy effect. Baseline levels of tumor-infiltrating CD4+/CD8+ T cells, NK cells, M1 macrophage, and inflammatory cytokine secretion were demonstrated to be associated with immune response (Topalian et al., 2016; Galon and Bruni, 2019). In addition, PD-L1 and TMB are recognized as biomarkers for forecasting anti–PD-1/L1 treatment response. In this study, we found that antitumor lymphocyte subsets including activated CD4+/CD8+ T cells and NK T cells were mainly enriched in m6Acluster-B and m6Acluster-C subtypes, both of which exhibited higher PD-L1 expression, TMB, and poorer prognosis. This supported the potential predictive value of the three m6A-related subtypes in terms of immunotherapeutic response and prognosis. Furthermore, from DEGs identified in different m6A decoration patterns, we obtained three transcriptome subtypes according to the m6A signature gene and found that they were obviously correlated with various survival outcomes and the TME landscape. Based on the DEGs associated with prognosis, we established a scoring system called “m6Ascore signature”. The m6Ascore was correlated with TME characterization, immune cell infiltration, and also with predictors of immune response (expression of immune checkpoints) in LUAD, suggesting that m6A modification may affect the efficacy of immunotherapy. Tumors with infiltration of immune cells, especially CD8+ T cells, and high expression of PD-L1 in tumor cells and stroma were usually determined as thermo-tumors, and they were often more closely associated with immune checkpoint inhibitors (Chen and Mellman, 2017; Wu et al., 2018). It was found in our study that tumors with low m6Ascore had more aggregation of antitumor lymphocyte subsets, such as activated CD4+/CD8+ T cells and NK T cells. Moreover, this group of tumors had higher tumor mutation load and immune checkpoint expression. In addition, they also showed lower sensitivity to targeted therapies and poorer prognosis. This not only proved that m6Ascore signature had the potential to predict the curative effect of immune targeted therapy and prognosis of patients, but more importantly, it highlighted the importance of m6A modification in shaping tumor immunity. Interestingly, multivariate Cox analysis identified WTAP as an independent prognostic gene for LUAD patients, while this gene showed no statistical significance for OS in LUSC patients, suggesting that m6A regulators displayed different clinical value in various tumors (Gu et al., 2021).

Integrating different independent studies on common key characteristics of disease has become a preferred strategy. Since there may be deviations in a single experiment, it is necessary to seek the findings supported by several evidences to improve the reliability. Based on a larger LUAD cohort, this study disclosed the methylation modification pattern, TME landscape, and clinical significance of m6A regulatory factors in LUAD. Although we reviewed the literature and planned 23 recognized RNA methylation regulators, only 17 genes integrating multiple data sets were finally included. It is, therefore, necessary to incorporate more newly identified m6A regulators into the model to optimize the accuracy of m6A refitted decorative patterns. In addition, although a large number of retrospective data sets were used in this study to identify different m6A modification patterns and m6Ascore, there is still a lack of appropriate LUAD data sets based on immunotherapy regimens to verify the predictive robustness of m6Ascore to further strengthen our conclusions. Therefore, a prospective study involving a cohort of LUAD patients receiving immunotherapy is required to confirm our findings.

Conclusion

Based on m6A regulators, the present study conducted an integrated assessment of m6A modification patterns in more than 1,900 LUAD samples to characterize the multidimensional profile of m6A regulators in LUAD. It was found that m6A regulators served an indispensable role in shaping heterogeneity and sophistication of the TME. According to these findings, we postulated that the evaluation of the m6A modification pattern of LUAD would help better understand the infiltration characteristics of the TME and provided important insights into the effectiveness of immunotherapy. m6A modification could be a promising strategy for the management of LUAD.

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

LX and RD contributed to designing and drafting the manuscript. XW, GX, ZG and XX critically revised the manuscript. All the authors gave final approval of this article.

Conflict of Interest

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

Publisher’s Note

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

Acknowledgments

We are especially grateful for the help provided by the Affiliated Hospital of Nantong University in this study. In addition, the work benefited from TCGA and GEO databases and other online analytical tools. We are grateful for the access to the resources and the efforts of the staff to expand and improve the databases.

Supplementary Material

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

Supplementary Figure S1 | Effect of the mutation status of ZC3H13 on the expression of other m6A regulators in TCGA–LUAD cohort.

Supplementary Figure S2 | (A) Univariate and (B) multivariate Cox regression analyses of 17 m6A regulatory factors. (C) Unsupervised clustering of the m6A regulators and consensus matrices for k = 3. (D) Unsupervised clustering of 17 m6A regulators. The m6Acluster, project, age, sex, smoking status, and living condition were used as patient notes. Red indicated the high expression of regulatory factors, while green represented the low expression. (E) Unsupervised clustering of the m6A phenotype–related genes correlated with prognosis and consensus matrices for k = 3. (F) GSVA enrichment analysis of the activation state of biological pathways in different m6A modification patterns (m6Acluster B vs. m6Acluster A). Heat maps of these biological processes, with red representing activated pathways and blue representing inhibited pathways, using LUAD queue as sample annotation. (G) Expression of 17 m6A regulatory factors in the three gene clusters. The upper and lower ends of the box represented the quartile range of values. The lines in the box represented intermediate values, and the real points represented outliers. The asterisk represented the statistical p value. (***p < 0.001; **p < 0.01; *p < 0.05). The statistical differences between the three gene clusters were tested by using one-way ANOVA.

Supplementary Figure S3 | Predictive performance of m6Ascore in various subgroups of patients with LUAD. (A) Age ≤ 65 subset. (B) Age > 65 subset. (C) Female subset. (D) Male subset. (E) Ever smoking subset. (F) Never smoking subset. (G) Stage I/II subset (H) Stage III/IV stage. (I) EGFR wt subset. (J) STK11 wt subset. (K) Kras wt subset. (L) Kras mut subset. (M) TP53 wt subset. (N) TP53 mut subset.

Supplementary Figure S4 | Relationship between m6Ascore and clinical parameters. (A) Age. (B) Sex. (C) Smoking status. (D) Stage. (E) EGFR. (F) STK11. (G) Kras. (H) TP53. (I) Survival status. The right side of each figure showed the proportion of high- and low-m6Ascore groups in each clinical characteristic, and the left side showed the comparison of m6Ascores in each clinical subgroup.

Supplementary Figure S5 | (A) Comparison of the TME infiltration cells abundances in the high- and low-m6Ascore groups. (B) Comparison of immune scores in the high- and low-m6Ascore groups. (C) Comparison of tumor purity in the high- and low-m6Ascore groups. (D) GSEA analysis of the high- and low-m6Ascore groups. (E,F) GSVA analysis of the high- and low-m6Ascore groups. E, KEGG; F, GO.

Supplementary Figure S6 | Construction and evaluation of the nomogram. (A) Univariate Cox analysis of m6Ascore, m6Acluster, and some clinical parameters (age, sex, smoking status, and stage). (B) Multivariate Cox analysis of m6Ascore, m6Acluster, and some clinical parameters (age, sex, smoking status, and stage). (C) Nomogram for predicting the 3- and 5-year survival of patients. (D) Calibration curves of 3 and 5 years.

References

Arguello, A. E., DeLiberto, A. N., and Kleiner, R. E. (2017). RNA Chemical Proteomics Reveals the N6-Methyladenosine (m6A)-Regulated Protein-RNA Interactome. J. Am. Chem. Soc. 139, 17249–17252. doi:10.1021/jacs.7b09213

PubMed Abstract | CrossRef Full Text | Google Scholar

Bai, Y., Yang, C., Wu, R., Huang, L., Song, S., Li, W., et al. (2019). YTHDF1 Regulates Tumorigenicity and Cancer Stem Cell-like Activity in Human Colorectal Carcinoma. Front. Oncol. 9, 332. doi:10.3389/fonc.2019.00332

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, 108–112. doi:10.1038/nature08460

PubMed Abstract | CrossRef Full Text | Google Scholar

Barros-Silva, D., Lobo, J., Guimarães-Teixeira, C., Carneiro, I., Oliveira, J., Martens-Uzunova, E. S., et al. (2020). VIRMA-dependent N6-Methyladenosine Modifications Regulate the Expression of Long Non-coding RNAs CCAT1 and CCAT2 in Prostate Cancer. Cancers (Basel) 12, 12. doi:10.3390/cancers12040771

PubMed Abstract | CrossRef Full Text | Google Scholar

Becht, E., Giraldo, N. A., Lacroix, L., Buttard, B., Elarouci, N., Petitprez, F., et al. (2016). Estimating the Population Abundance of Tissue-Infiltrating Immune and Stromal Cell Populations Using Gene Expression. Genome Biol. 17, 218. doi:10.1186/s13059-016-1070-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Botling, J., Edlund, K., Lohr, M., Hellwig, B., Holmberg, L., Lambe, M., et al. (2013). Biomarker Discovery in Non-small Cell Lung Cancer: Integrating Gene Expression Profiling, Meta-Analysis, and Tissue Microarray Validation. Clin. Cancer Res. 19, 194–204. doi:10.1158/1078-0432.ccr-12-1139

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, 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, 321–330. doi:10.1038/nature21349

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, R., Xu, X., Qian, Z., Zhang, C., Niu, Y., Wang, Z., et al. (2019). The Biological Functions and Clinical Applications of Exosomes in Lung Cancer. Cell. Mol. Life Sci. 76, 4613–4633. doi:10.1007/s00018-019-03233-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, X.-Y., Zhang, J., and Zhu, J.-S. (2019). The Role of m6A RNA Methylation in Human Cancer. Mol. Cancer 18, 103. doi:10.1186/s12943-019-1033-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Z., Wu, L., Zhou, J., Lin, X., Peng, Y., Ge, L., et al. (2020). N6-methyladenosine-induced ERRγ Triggers Chemoresistance of Cancer Cells through Upregulation of ABCB1 and Metabolic Reprogramming. Theranostics 10, 3382–3396. doi:10.7150/thno.40144

PubMed Abstract | CrossRef Full Text | Google Scholar

Colaprico, A., Silva, T. C., Olsen, C., Garofano, L., Cava, C., Garolini, D., et al. (2016). TCGAbiolinks: an R/Bioconductor Package for Integrative Analysis of TCGA Data. Nucleic Acids Res. 44, e71. doi:10.1093/nar/gkv1507

PubMed Abstract | CrossRef Full Text | Google Scholar

Der, S. D., Sykes, J., Pintilie, M., Zhu, C.-Q., Strumpf, D., Liu, N., et al. (2014). Validation of a Histology-independent Prognostic Gene Signature for Early-Stage, Non-small-cell Lung Cancer Including Stage IA Patients. J. Thorac. Oncol. 9, 59–64. doi:10.1097/jto.0000000000000042

PubMed Abstract | CrossRef Full Text | Google Scholar

Galon, J., and Bruni, D. (2019). Approaches to Treat Immune Hot, Altered and Cold Tumours with Combination Immunotherapies. Nat. Rev. Drug Discov. 18, 197–218. doi:10.1038/s41573-018-0007-y

PubMed Abstract | 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, e107468. doi:10.1371/journal.pone.0107468

PubMed Abstract | CrossRef Full Text | Google Scholar

Geeleher, P., Cox, N. J., and Huang, R. (2014). Clinical Drug Response Can Be Predicted Using Baseline Gene Expression Levels and In Vitro Drug Sensitivity in Cell Lines. Genome Biol. 15, R47. doi:10.1186/gb-2014-15-3-r47

PubMed Abstract | CrossRef Full Text | Google Scholar

Gu, C., Shi, X., Dai, C., Shen, F., Rocco, G., Chen, J., et al. (2020). RNA m6A Modification in Cancers: Molecular Mechanisms and Potential Clinical Applications. The Innovation 1, 100066. doi:10.1016/j.xinn.2020.100066

PubMed Abstract | CrossRef Full Text | Google Scholar

Gu, C., Shi, X., Qiu, W., Huang, Z., Yu, Y., Shen, F., et al. (2021). Comprehensive Analysis of the Prognostic Role and Mutational Characteristics of m6A-Related Genes in Lung Squamous Cell Carcinoma. Front. Cell Dev. Biol. 9, 661792. doi:10.3389/fcell.2021.661792

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Leek, J. T., Johnson, W. E., Parker, H. S., Jaffe, A. E., and Storey, J. D. (2012). The Sva Package for Removing Batch Effects and Other Unwanted Variation in High-Throughput Experiments. Bioinformatics 28, 882–883. doi:10.1093/bioinformatics/bts034

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, H., Zhang, Y., Guo, Y., Liu, R., Yu, Q., Gong, L., et al. (2021). ALKBH1 Promotes Lung Cancer by Regulating m6A RNA Demethylation. Biochem. Pharmacol. 189, 114284. doi:10.1016/j.bcp.2020.114284

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X.-C., Jin, F., Wang, B.-Y., Yin, X.-J., Hong, W., and Tian, F.-J. (2019). The m6A Demethylase ALKBH5 Controls Trophoblast Invasion at the Maternal-Fetal Interface by Regulating the Stability of CYR61 mRNA. Theranostics 9, 3853–3865. doi:10.7150/thno.31868

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., Gu, J., Xu, F., Zhu, Q., Chen, Y., Ge, D., et al. (2021). Molecular Characterization, Biological Function, Tumor Microenvironment Association and Clinical Significance of m6A Regulators in Lung Adenocarcinoma. Brief Bioinform 22, 22. doi:10.1093/bib/bbaa225

CrossRef Full Text | Google Scholar

Li, Y., Xiao, J., Bai, J., Tian, Y., Qu, Y., Chen, X., et al. (2019). Molecular Characterization and Clinical Relevance of m6A Regulators across 33 Cancer Types. Mol. Cancer 18, 137. doi:10.1186/s12943-019-1066-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Z., Weng, H., Su, R., Weng, X., Zuo, Z., Li, C., et al. (2017). FTO Plays an Oncogenic Role in Acute Myeloid Leukemia as a N 6 -Methyladenosine RNA Demethylase. Cancer Cell 31, 127–141. doi:10.1016/j.ccell.2016.11.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, S., Li, Q., Chen, K., Zhang, Q., Li, G., Zhuo, L., et al. (2020). The Emerging Molecular Mechanism of m6A Modulators in Tumorigenesis and Cancer Progression. Biomed. Pharmacother. 127, 110098. doi:10.1016/j.biopha.2020.110098

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Y., Guo, X., Zhao, M., Ao, H., Leng, X., Liu, M., et al. (2020). Contributions and Prognostic Values of M 6 A RNA Methylation Regulators in Non‐small‐cell Lung Cancer. J. Cell Physiol. 235, 6043–6057. doi:10.1002/jcp.29531

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, J. z., Yang, F., Zhou, C. c., Liu, F., Yuan, J. h., Wang, F., et al. (2017). METTL14 Suppresses the Metastatic Potential of Hepatocellular Carcinoma by Modulating N 6 ‐methyladenosine‐dependent Primary MicroRNA Processing. Hepatology 65, 529–543. doi:10.1002/hep.28885

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, L., Chen, T., Zhang, X., Miao, Y., Tian, X., Yu, K., et al. (2021). The m6A Reader YTHDC2 Inhibits Lung Adenocarcinoma Tumorigenesis by Suppressing SLC7A11-dependent Antioxidant Function. Redox Biol. 38, 101801. doi:10.1016/j.redox.2020.101801

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, 1747–1756. doi:10.1101/gr.239244.118

PubMed Abstract | CrossRef Full Text | Google Scholar

Okayama, H., Kohno, T., Ishii, Y., Shimada, Y., Shiraishi, K., Iwakawa, R., et al. (2012). Identification of Genes Upregulated in ALK-Positive and EGFR/KRAS/ALK-negative Lung Adenocarcinomas. Cancer Res. 72, 100–111. doi:10.1158/0008-5472.can-11-1403

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Rousseaux, S., Debernardi, A., Jacquiau, B., Vitte, A. L., Vesin, A., Nagy-Mignotte, H., et al. (2013). Ectopic Activation of Germline and Placental Genes Identifies Aggressive Metastasis-Prone Lung Cancers. Sci. Transl Med. 5, 186ra66–186r. doi:10.1126/scitranslmed.3005723

PubMed Abstract | CrossRef Full Text | Google Scholar

Schabath, M. B., Welsh, E. A., Fulp, W. J., Chen, L., Teer, J. K., Thompson, Z. J., et al. (2016). Differential Association of STK11 and TP53 with KRAS Mutation-Associated Gene Expression, Proliferation and Immune Surveillance in Lung Adenocarcinoma. Oncogene 35, 3209–3216. doi:10.1038/onc.2015.375

PubMed Abstract | CrossRef Full Text | Google Scholar

Shedden, K., Shedden, K., Taylor, J. M., Enkemann, S. A., Tsao, M. S., Yeatman, T. J., et al. (2008). Gene Expression-Based Survival Prediction in Lung Adenocarcinoma: a Multi-Site, Blinded Validation Study. Nat. Med. 14, 822–827. doi:10.1038/nm.1790

PubMed Abstract | CrossRef Full Text | Google Scholar

Shulman, Z., and Stern-Ginossar, N. (2020). The RNA Modification N6-Methyladenosine as a Novel Regulator of the Immune System. Nat. Immunol. 21, 501–512. doi:10.1038/s41590-020-0650-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Siegel, R. L., Miller, K. D., and Jemal, A. (2020). Cancer Statistics, 2020. CA A. Cancer J. Clin. 70, 7–30. doi:10.3322/caac.21590

CrossRef Full Text | Google Scholar

Sun, J., Ping, Y., Huang, J., Zeng, B., Ji, P., and Li, D. (2021). N6-Methyladenosine-Regulated mRNAs: Potential Prognostic Biomarkers for Patients with Lung Adenocarcinoma. Front. Cell Dev. Biol. 9, 705962. doi:10.3389/fcell.2021.705962

PubMed Abstract | CrossRef Full Text | Google Scholar

Topalian, S. L., Taube, J. M., Anders, R. A., and Pardoll, D. M. (2016). Mechanism-driven Biomarkers to Guide Immune Checkpoint Blockade in Cancer Therapy. Nat. Rev. Cancer 16, 275–287. doi:10.1038/nrc.2016.36

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, T., Kong, S., Tao, M., and Ju, S. (2020). The Potential Role of RNA N6-Methyladenosine in Cancer Progression. Mol. Cancer 19, 88. doi:10.1186/s12943-020-01204-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Wang, Y., Luo, W., Song, X., Huang, L., Xiao, J., et al. (2020). Roles of Long Non-coding RNAs and Emerging RNA-Binding Proteins in Innate Antiviral Responses. Theranostics 10, 9407–9424. doi:10.7150/thno.48520

PubMed Abstract | CrossRef Full Text | Google Scholar

Warda, A. S., Kretschmer, J., Hackert, P., Lenz, C., Urlaub, H., Höbartner, C., et al. (2017). Human METTL16 Is a N 6 ‐methyladenosine (M 6 A) Methyltransferase that Targets pre‐mRNAs and Various Non‐coding RNAs. Embo Rep. 18, 2004–2014. doi:10.15252/embr.201744940

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, 1572–1573. doi:10.1093/bioinformatics/btq170

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, H., Dong, H., Fu, Y., Tang, Y., Dai, M., Chen, Y., et al. (2021). Expressions of m6A RNA Methylation Regulators and Their Clinical Predictive Value in Cervical Squamous Cell Carcinoma and Endometrial Adenocarcinoma. Clin. Exp. Pharmacol. Physiol. 48, 270–278. doi:10.1111/1440-1681.13412

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, S.-P., Liao, R.-Q., Tu, H.-Y., Wang, W.-J., Dong, Z.-Y., Huang, S.-M., et al. (2018). Stromal PD-L1-Positive Regulatory T Cells and PD-1-Positive CD8-Positive T Cells Define the Response of Different Subsets of Non-small Cell Lung Cancer to PD-1/PD-L1 Blockade Immunotherapy. J. Thorac. Oncol. 13, 521–532. doi:10.1016/j.jtho.2017.11.132

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, T., Hu, E., Xu, S., Chen, M., Guo, P., Dai, Z., et al. (2021). clusterProfiler 4.0: A Universal Enrichment Tool for Interpreting Omics Data. The Innovation 2, 100141. doi:10.1016/j.xinn.2021.100141

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, Y., Chang, N., Zhang, Y., Zhang, X., Xu, L., Che, Y., et al. (2021). METTL3-mediated m6A mRNA Modification of FBXW7 Suppresses Lung Adenocarcinoma. J. Exp. Clin. Cancer Res. 40, 90. doi:10.1186/s13046-021-01880-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, Y., Xiao, G., Coombes, K. R., Behrens, C., Solis, L. M., Raso, G., et al. (2011). Robust Gene Expression Signature from Formalin-Fixed Paraffin-Embedded Samples Predicts Prognosis of Non-small-cell Lung Cancer Patients. Clin. Cancer Res. 17, 5705–5714. doi:10.1158/1078-0432.ccr-11-0196

PubMed Abstract | CrossRef Full Text | Google Scholar

Zaccara, S., Ries, R. J., and Jaffrey, S. R. (2019). Reading, Writing and Erasing mRNA Methylation. Nat. Rev. Mol. Cell Biol 20, 608–624. doi:10.1038/s41580-019-0168-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeng, D., Ye, Z., Wu, J., Zhou, R., Fan, X., Wang, G., et al. (2020). Macrophage Correlates with Immunophenotype and Predicts Anti-PD-L1 Response of Urothelial Cancer. Theranostics 10, 7002–7014. doi:10.7150/thno.46176

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, D., Zhang, D., Wang, C., Yang, X., Zhang, R., Li, Q., et al. (2021). Gene and Prognostic Value of N6-Methyladenosine (m6A) Modification Regulatory Factors in Lung Adenocarcinoma. Eur. J. Cancer Prev. doi:10.1097/cej.0000000000000717

CrossRef Full Text | Google Scholar

Zhang, Y., Liu, X., Liu, L., Li, J., Hu, Q., and Sun, R. (2020). Expression and Prognostic Significance of m6A-Related Genes in Lung Adenocarcinoma. Med. Sci. Monit. 26, e919644. doi:10.12659/MSM.919644

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, B. S., Roundtree, I. A., and He, C. (2017). Post-transcriptional Gene Regulation by mRNA Modifications. Nat. Rev. Mol. Cell Biol 18, 31–42. doi:10.1038/nrm.2016.132

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, Q., Zhao, Y., Hu, W., Zhang, Y., Wu, X., Lu, J., et al. (2020). m6A RNA Modification Modulates PI3K/Akt/mTOR Signal Pathway in Gastrointestinal CancerA RNA Modification Modulates PI3K/Akt/mTOR Signal Pathway in Gastrointestinal Cancer. Theranostics 10, 9528–9543. doi:10.7150/thno.42971

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, H., Zheng, M., Shi, M., Wang, J., Huang, Z., Zhang, H., et al. (2021). Characteristic of Molecular Subtypes in Lung Adenocarcinoma Based on m6A RNA Methylation Modification and Immune Microenvironment. Bmc Cancer 21, 938. doi:10.1186/s12885-021-08655-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, Z., Lv, J., Yu, H., Han, J., Yang, X., Feng, D., et al. (2020). Mechanism of RNA Modification N6-Methyladenosine in Human Cancer. Mol. Cancer 19, 104. doi:10.1186/s12943-020-01216-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhuang, Z., Chen, L., Mao, Y., Zheng, Q., Li, H., Huang, Y., et al. (2020). Diagnostic, Progressive and Prognostic Performance of m6A Methylation RNA Regulators in Lung Adenocarcinoma. Int. J. Biol. Sci. 16, 1785–1797. doi:10.7150/ijbs.39046

PubMed Abstract | CrossRef Full Text | Google Scholar