The flowchart of our study is depicted in Figure 1A .

The transcriptome profiling data and corresponding clinical data of LUAD were acquired from the TCGA database (https://portal.gdc.cancer.gov) and GEO database (https://www.ncbi.nlm.nih.gov/gds). In the TCGA cohort, the transcriptome data were downloaded as FPKM (fragments per kilobase per million mapped reads) and further converted to TPM (transcripts per million) using the “limma” R package for analysis. Then the “normalizeBetweenArrays” function of the R package “limma” was performed for data standardization. The TCGA-LUAD dataset including 535 tumor samples and 59 tumor-adjacent samples was used to compare the difference in the expression of m⁷G-related genes between tumor and normal tissues. For the GEO datasets, probe IDs were converted to gene symbols according to platform annotation files. Normalized expression values were log2-transformed and scaled before being used in model validations. Using the “combat” algorithm in the “sva” package of the R software, we correct the batch effect between the TCGA and GEO datasets. Genome mutation data of TCGA-LUAD [including somatic mutation and copy number variation (CNV)] were downloaded from the TCGA database and the UCSC Xena platform (https://gdc.xenahubs.net/). For the analyses involving clinical data, samples with unknown survival times were deleted.

We identified m⁷G-related genes from published literature (28) and the gene sets named “m7G(5′)pppN diphosphatase activity”, “RNA 7-methylguanosine cap binding”, and “RNA cap binding” from the Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb/search.jsp) with the keyword “7-methylguanosine.” To screen differentially expressed m⁷G-related genes (DEMGs) with the threshold of p-value <0.05 between tumor and normal adjacent tissues from the TCGA-LUAD dataset, the “limma” R package was used. After screening, the “heatmap” R package was applied for generating heatmaps. The immunohistochemistry (IHC) results from the Human Protein Atlas (HPA, https://www.proteinatlas.org) were used to validate the protein level of DEMGs in normal and tumor tissues. Meanwhile, we observed genetic mutation frequency and types of m⁷G-related genes in TCGA-LUAD samples by using the “maftools” R package.

Among proteins with co-expression coefficients >0.4, the STRING database (https://string-db.org/) was used to construct the protein–protein interaction (PPI) network. Cytoscape software (version 3.9.1) was used to visualize the network; moreover, the MCC algorithm of the cytoHubba plugin was used to screen the hub genes. The “reshape2” R package was used to identify the correlation between the expression of m⁷G-related genes.

Overall survival (OS) was assessed by the Kaplan–Meier method. The datasets of TCGA-LUAD and GSE68465 (29) were merged to explore the overall survival predictive value of each m⁷G-related gene. R packages “survival” and “survminer” were used.

To classify LUAD samples into different subgroups based on the m⁷G-related gene set, the R package “ConsensusClusterPlus” was used. We combined transcriptome profiling data and corresponding clinical data from the TCGA-LUAD and GSE68465 datasets for this analysis. The maximum number of clusters was nine. We selected 80% item resampling (pItem), 100% gene resampling (pFeature), a maximum evaluated k of 9, 50 resamplings (reps), kmeans (clusterAlg), Euclidean (distance), and a specific random seed (seed = 123456) in the R package “ConsensusClusterPlus” for this analysis. Based on the consensus matrices and the cumulative distribution function (CDF) curves of the consensus index, the optimum number of clusters was determined. The differences of survival and distribution of clinicopathologic characteristics were compared between different clusters and visualized using the R packages “survival”, “survminer”, “limma”, “ggplot2”, and “pheatmap”. “Gene set variation analysis (GSVA)” R packages were further used to explore the difference in biological processes between different clusters. We downloaded the gene set “c2.cp.kegg.v7.4.symbols” from the MSigDB database for this analysis. An adjusted p-value of less than 0.05 was considered statistically significant. A single-sample gene set enrichment analysis (ssGSEA) was used to quantify the enrichment scores and to represent the relative abundance of 23 tumor-infiltrating immune cell types between different clusters. The “maftools” R package was used to present mutational differences of each cluster. The CNV of different clusters was analyzed and visualized using the RCircos package in R.

Also, transcriptome profiling data and corresponding clinical data from the TCGA-LUAD and GSE68465 datasets were merged for this analysis. The combined data were randomly divided into a training cohort and an internal validation cohort in a 1:1 ratio. Meanwhile, the 442 LUAD samples from GSE72094 (30) were used as an external validation cohort. Briefly, by using the “limma” R package, the differentially expressed genes (DEGs) with adjusted p-value <0.05 and |log2FC| >0.585 between clusters were detected. Subsequently, univariate Cox analysis was applied to explore the prognosis-related DEGs. Then, we used the LASSO regression with 10-fold cross-validation to narrow down the prognosis-related DEGs applying the R package “glmnet,” and further performed the multivariate Cox regression analysis to establish a signature for evaluating the relationships between the DEGs and the survival of LUAD patients. The IHC results from the HPA were used to validate the protein level of this gene signature in normal and tumor tissues. Finally, we calculated the risk scores of each patient based on the following model formula: risk score = Σi Coefficient (mRNA_i ) × Expression(mRNA_i ). According to the median risk scores in the training group, patients were separated into high- and low-risk groups among both training and validation cohorts. The Kaplan–Meier method was performed to compare the OS between the high- and low-risk groups. The predictive value of the prognostic model was assessed through time-related receiver operating characteristic (ROC). The heatmaps were used to compare and visualize the distribution of clinicopathologic characteristics between the risk cohorts. Multivariate Cox regression analysis was applied to test the prognostic independence of risk score. The principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) which can get a low-dimensional cluster distribution from high-dimensional gene sets were utilized for validating the classification results. The risk scores of each patient were further combined with clinical characteristics to construct a nomogram through the “rms” R package.

To compare our prognostic model with other models previously reported in LUAD, four studies primarily focusing on m⁶A-related signatures were selected (22, 31–33). We extracted the genes included in these prognostic models and used ROC curves and concordance index (C-index) values to compare the predictive accuracy between different models. The R packages “survcomp,” “survival,” “ggplot2,” “ggpubr,” “limma,” “survminer,” and “timeROC” were applied.

To explore the potential mechanisms and pathways between the high- and low-risk groups, the Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analysis, and gene set enrichment analysis (GSEA) were conducted among DEGs between the high- and low-risk groups using the R packages “clusterProfiler,” “enrichplot,” “limma,” “ggplot2,” and “org.Hs.eg.db.”

Using the ssGSEA through the “GSVA” R package, we compared the enrichment scores, represented the relative abundance of 23 tumor-infiltrating immune cell types between the high- and low-risk groups, and then visualized the results through the R packages “limma,” “ggpubr,” and”reshape2.” The differential analyses of stromal score, immune score, ESTIMATE score, and immune cells were performed based on the results of CIBERSORT and ESTIMATE using the R software packages “CIBERSORT” and “estimate.” Moreover, the mRNA expression-based stemness index (mRNAsi) scores of LUAD were obtained from previous research (34). Spearman’s correlation analyses were used to establish the relationship between risk scores and immune cells as well as between risk scores and mRNAsi scores. The tumor immune dysfunction and exclusion (TIDE) score (35) was calculated online (http://tide.dfci.harvard.edu/) to assess the immune checkpoint inhibitor (ICI) response between the high- and low-risk groups. We also compared the expressions of 47 immune checkpoint-related genes ( Supplementary Table 1 ) between the high- and low-risk groups through the R software packages “limma,” “ggpubr,” and “ggplot2.”

The public dataset, Genomics of Drug Sensitivity in Cancer (GDSC, https://www.cancerrxgene.org) (36), was chosen to evaluate the response to chemotherapeutic and small molecule drugs between the high- and low-risk groups, which was done by computing the half-maximum inhibitory concentration (IC₅₀). The analyses were conducted through the “pRRophetic” R package.

All analyses were completed by using R language (Version 4.1.2). Student’s t-test, chi-squared test, or Wilcoxon test was applied to compare the differences between groups. Spearman’s correlation test was performed to evaluate the association between variables. A p-value of <0.05 was considered statistically significant.

After screening, we found 29 m⁷G-related genes, namely, AGO2, CYFIP1, DCP2, DCPS, EIF3D, EIF4A1, EIF4E, EIF4E1B, EIF4E2, EIF4E3, EIF4G3, GEMIN5, IFIT5, LARP1, LSM1, METTL1, NCBP1, NCBP2, NCBP2L, NCBP3, NSUN2, NUDT10, NUDT11, NUDT16, NUDT3, NUDT4, NUDT4B, SNUPN, and WDR4. Of these, 12 DEMGs with the threshold of |log2FC| >0.585 and p <0.05 were observed between 535 tumor and 59 normal TCGA-LUAD tissues, which consisted of 2 downregulated (NCBP2L, EIF4E3) and 10 upregulated (DCPS, EIF4E1B, EIF4G3, LARP1, LSM1, METTL1, NCBP1, NCBP2, NSUN2, and WDR4) DEMGs in tumor samples ( Figure 1B and Table 1 ). The protein expressions of DEMGs in normal and tumor tissues were further validated using the IHC results from the HPA platform. As shown in Figure 2 , the protein expressions of METTL1, NSUN2, EIF4G3, LARP1, NCBP1, and NCBP2 were higher in tumor tissues than in normal tissues; however, the protein expressions of WDR4, EIF4E1B, and NCBP2L were negative in both tumor and normal tissues. The PPI network between 24 DEMGs with the threshold of p <0.05 is shown in Figure 1C ; of these DEMGs, EIF4E, EIF4E1B, EIF4E2, NCBP1, NCBP2, EIF4E3, AGO2, NCBP2L, CYFIP1, and EIF4A1 were the top 10 hub genes. Spearman’s correlation analysis suggested that EIF4E3 was most frequently and negatively correlated with other DEMGs, as shown in Figure 1D . A total of 561 samples in the TCGA-LUAD cohort were used for somatic mutation of m⁷G-related genes; among them, 80 (14.26%) samples were observed to have experienced mutation events. A missense mutation was the most common type of variant classification, and EIF4G3, LARP1, NSUN2, AGO2, and CYFIP1 were the top 5 most frequently mutated genes ( Figure 1E ).

To further verify the predictive role of each m⁷G-related gene, we conducted survival analyses using the Kaplan–Meier method. As shown in Figure 3 , the high expression of most genes (including CYFIP1, DCPS, EIF3D, EIF4E, EIF4E2, EIF4G3, LARP1, METTL1, NCBP1, NCBP2, NUDT4, NUDT11, SNUPN, and WDR4) portended a significantly poor OS. Conversely, a high expression of NUDT3 predicted an improved OS ( Figure 3 ).

Based on the expression similarity of the 29 m⁷G-related genes, the consensus clustering method was applied to cluster the combined LUAD samples of the TCGA and GSE68465 cohorts. The cluster number k ranged from 2 to 9; when k = 2, a relatively clear-cut boundary was shown in the heatmap of the consensus matrix ( Figure 4A , Supplementary Figure S1A ), and a flat slope was seen in the CDF curve of the consensus index score ( Supplementary Figure S1A ). Hence, we selected k = 2 as the appropriate number of clusters and divided 936 LUAD samples into two clusters—namely, cluster 1 (C1, n = 400) and cluster 2 (C2, n = 536). Next, Kaplan–Meier survival analysis was applied to evaluate the prognostic value of this clustering. A significant difference in OS was observed between the two clusters (log-rank p < 0.001). C1 had a worse median OS ( Figure 4B ). Accordingly, the samples of C1 had higher levels of gene expression compared with C2, as shown in the heatmap ( Figure 4D ). Also, as seen in the heatmap, however, variations in other clinicopathological characteristics did not show statistical significance between the samples of each cluster. Similarly, C1 had a higher expression of most m⁷G-related genes compared with C2 ( Supplementary Figure S1B ). The results of the GSVA enrichment analysis based on the KEGG gene set showed that the process of cell cycle and DNA repair, such as non-homologous end-joining, spliceosome, nucleotide excision repair, mismatch repair, and cell cycle, was enriched in C1; C2 was prominently enriched in metabolism-associated pathways, such as tyrosine metabolism, arachidonic acid metabolism, drug metabolism cytochrome P450, metabolism of xenobiotics by cytochrome P450, sulfur metabolism, and alpha-linolenic acid metabolism ( Figure 4E ). From the results of the ssGSEA, we found that the scores of some immune cells, such as activated B cells, activated CD8 T cells, activated dendritic cells (aDCs), CD56 dim nature killer cells, myeloid-derived suppressor cells (MDSCs), macrophages, mast cells, monocytes, neutrophils, T follicular helper cells, type 1 T helper cells, and type 17 T helper cells, were significantly enriched in C2 ( Figure 4C ). In addition, the somatic mutations were more frequent in C1 than in C2 ( Figure 4F ). The frequencies and locations of the CNVs were also different in C1 and C2 ( Figures 4G, H ), and the copy number losses were more frequent in C2 than in C1.

Using the “limma” R package, the DEGs were first screened between the two clusters in a combined LUAD dataset of the TCGA and GSE68465 cohorts, and 130 DEGs were finally obtained. Then, we used the univariate Cox analysis to explore 112 prognosis-related DEGs. To prevent model overfitting, LASSO penalized Cox regression modeling was conducted to screen the key DEGs associated with survival. With this method, a novel prognostic gene model with four genes was constructed ( Figures 5A, B ). Subsequently, risk scores per sample were calculated using the following model formula: risk score = (0.1606739599952 × expression value of KIF20B) + (0.207218473949824 × expression value of HMMR) + (0.157719134596455 × expression value of ARNTL2) + (0.0802509860697548 × expression value of DKK1). The combined LUAD dataset was randomly divided into a training cohort and internal validation cohort in a 1:1 ratio, and the 442 LUAD samples from the GSE72094 dataset were used as an external validation cohort. The samples were divided into high-risk and low-risk groups according to the median threshold of risk scores in the training group ( Figures 5J–L ). The expressions of most m⁷G-related genes were significantly higher in the high-risk group than in the low-risk group ( Figure 5C ). As shown by the Kaplan–Meier analyses, patients in the high-risk group had significantly worse OS than those in the low-risk group (p < 0.001, Figures 5D–F ); analogously, as the risk score increased, more patients died ( Figures 5J–L ). In the training cohort, the AUC values of the present risk model were 0.734, 0.691, and 0.676 for the 1-, 2-, and 3-year prognoses, respectively ( Figure 5G ), with similar results observed in the internal and external validation groups ( Figures 5H, I ). The distribution patterns from t-SNE and PCA analyses showed that samples could completely be distinguished into high- and low-risk groups ( Figures 5M–O ). Taken together, these findings demonstrated the prognostic robustness of the novel prognostic model in patients with LUAD.

To compare our prognostic model with other models previously reported in LUAD, we searched four studies primarily focusing on m⁶A-related signatures (22, 31–33). The number of genes included in these models varied from 3 to 27. We found that the AUC and C-index values of our model were higher than those of other models except for the model constructed by Ouyang et al. (32) ( Figure 6A ).

Univariate and multivariate Cox regression analyses were performed by introducing age, gender, TNM stage, and risk scores to assess the independence of risk scores in the survival prediction of LUAD patients. Those variables with p <0.1 in the univariate analysis were selected for multivariate analysis. Among the samples in the training cohort, the results showed that age, T stage, N stage, and risk score were identified as independent negative prognostic factors for patients with LUAD ( Figure 6B ), and similar results were observed in the internal and external validation groups ( Supplementary Figures S2A, B ). Meanwhile, in the internal validation cohort, risk score was of greater assistance than the other clinical characteristics in predicting prognosis (HR = 1.689, 95% CI: 1.376–2.075, p < 0.001; Supplementary Figure S2A ). As shown in Figure 6C , the heatmap presented the distribution of clinicopathological features between the high- and low-risk cohorts. There were significant differences in N stage (p < 0.01), T stage (p < 0.001), and gender (p < 0.01) between the different risk groups. To facilitate the utilization of our risk model, we further constructed a nomogram using risk scores combined with clinical characteristics, as shown in Figure 6D , and calibration curves further verified this nomogram as reliable and accurate in predicting 1-, 3-, and 5-year OS ( Supplementary Figure S2G ).

By considering the DEGs between the high- and low-risk groups from combined samples of the TCGA and GSE68465 cohorts, we conducted GO enrichment analysis, KEGG pathway analysis, and GSEA to explore the potential biological functions of these genes. As shown in Figure 6E and Supplementary Figure S2D , “nuclear division” and “organelle fission” were the most enriched terms among the biological process categories, and “spindle” and “tubulin binding” were the most enriched terms among the cellular component and molecular function categories, respectively. Cell cycle was identified to be the most enriched among the KEGG pathways of the DEGs ( Figure 6F and Supplementary Figure S2C ). GSEA was further performed to identify the differential pathways enriched in GO and KEGG between the high- and low-risk subgroups, and the results revealed that cell cycle and DNA replication-related biological processes and pathways were enriched in the high-risk group; however, the biological processes and pathways enriched in the low-risk group were weakly associated with tumor initiation and progression ( Figure 6G and Supplementary Figure S2E ).

ssGSEA was performed to quantify the enrichment scores of 13 immune cell-related functions and 16 immune cells between the two risk groups. Intriguingly, the scores of some immune cells, such as aDCs, immature dendritic cells (iDCs), mast cells, and neutrophils, were significantly enriched in the low-risk group of the training cohort ( Figure 7A ). However, the scores of some immune functions, such as cytolytic activity, inflammation promoting, major histocompatibility complex (MHC) class I, parainflammation, type-I interferon (IFN) response, and T-cell co-inhibition, were significantly enriched in the high-risk group ( Figure 7B ). Similar results were observed in the internal testing dataset ( Supplementary Figure S2F ). We further analyzed the correlations between immune cell infiltration and the expression of genes involved in the construction of the prognostic model, as shown in Figure 7D . The high-risk group had a lower TIDE score ( Figure 7C ) and higher expressions of most immune checkpoint-related genes ( Figure 7E ), which suggested that patients in the high-risk group may benefit from immunotherapy. In addition, we observed that patients in the low-risk group had higher ESTIMATE, immune, and stromal scores ( Figure 7F ). Additionally, by using Spearman’s correlation analysis, a positive and significant correlation was observed between risk score and tumor mutation burden (TMB, R = 0.22, p < 0.001, Figures 7G, H ) as well as between risk score and mRNAsi score (RNAss, R = 0.35, p < 0.001, Figure 7I ), which demonstrated that LUAD patients with higher risk scores had higher RNAss values and TMBs.

We further evaluated the response to chemotherapeutic and small molecule drugs between the high- and low-risk groups, as described previously. In the high-risk group, a total of 57 drugs with obviously lower IC₅₀ values were observed; concomitantly, a total of 28 drugs were observed in the low-risk group ( Table 2 ). Based on the drugs commonly used to treat LUAD, we noticed that patients in the high-risk group were more sensitive to chemotherapeutic agents such as cisplatin, docetaxel, etoposide, gemcitabine, paclitaxel, and vinorelbine.