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”.

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.

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”.

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”.

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.

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:

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.

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.

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).

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.

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.

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.

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 ).

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 ).

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 ).

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 ).

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 ).

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.

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 ).

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.