A five-m6A regulatory gene signature is a prognostic biomarker in lung adenocarcinoma patients

Results

Seventeen m6A regulatory genes are differentially expressed in lung adenocarcinoma samples

We systematically analyzed the transcript levels of twenty m6A regulatory genes in 59 normal lung and 535 LADC tissue samples. The heat map showed that the expression levels of 17 m6A regulatory genes were significantly upregulated or downregulated in LADC tissues compared to normal lung tissues (Figure 1A). HNRNPA2B1, HNRNPC, HNRNPG, IGF2BP1, IGF2BP2, IGF2BP3, KIAA1429, METTL3, RBM15, YTHDF1, YTHDF2, and YTHDF3 were significantly upregulated, whereas, FTO, METTL14, METTL16, WTAP, and ZC3H13 were significantly downregulated in LADC tissues compared to the normal lung tissues (Figure 1A, 1B). However, expression levels of YTHDC1, YTHDC2, and ALKBH5 were similar in LADC and normal lung tissue samples (Figure 1A, 1B). We then analyzed single nucleotide variations (SNVs) in the seventeen differentially expressed m6A regulatory genes and found that the mutation frequency was ≤3% (Figure 1C). This suggested that differential expression of the 17 m6A regulatory genes was less likely to be caused by SNVs.

Figure 1. Expression patterns of twenty m6A regulatory genes in LDAC tissues and protein-protein interaction network analysis of m6A regulatory genes. (A) Heat map shows expression levels of 20 m6A regulatory genes in 59 normal lung and 535 LADC tissues. Note: *0.01 ≤ P < 0.05, **0.001 ≤ P < 0.01, ***P < 0.001. (B) Boxplot shows differential expression of the 20 m6A regulatory genes in LADC tissues compared to normal lung tissues. Note: *0.01 ≤ P < 0.05, **0.001 ≤ P < 0.01, ***P < 0.001. (C) Single nucleotide variations in seventeen m6A regulatory genes based on analysis of 561 LADC patient tissues. (D) PPI network analysis of seventeen m6A regulatory genes using STRING database. (E) The interacting numbers of each gene with the other sixteen genes. If the value of one gene is equal to 1, it means that only one of the other sixteen genes was interacted with this gene. (F) Spearman correlation analysis between 17 m6A regulatory genes. Blue indicates negative correlation in comparison with red indicating positive correlation. All the correlation coefficients are shown in the squares, and the areas of circles in the squares are positive correlated with the absolute value of corresponding correlation coefficients. Squares containing crosses denotes the P-values of correlation analyses are above 0.05.

PPI network and correlation analysis of 17 differentially expressed m6A regulatory genes

Protein-protein interaction (PPI) network analysis of the 17 m6A regulatory genes suggested that METTL3 and METTL14 were hub genes (Figure 1D, 1E), whereas, YTHDF3, IGF2BP2, and METTL16 were not linked to any other m6A regulatory genes (Figure 1D, 1E). In general, the ‘writer’ m6A regulatory genes showed significantly higher number of links and the ‘reader’ m6A regulatory genes showed significantly lower number of interacting partners (Figure 1D, 1E). Moreover, all ‘writers’ except METTL16 interacted with each other, whereas, METTL16 did not show any interaction with the other writers (Figure 1D). Correlation analysis showed weak to moderate association between the seventeen m6A regulatory genes (Figure 1F). The highest correlation was between KIAA1429 and YTHDF3 (r = 0.66), whereas, the correlation between FTO and HNRNPC was most negative (r = -0.29) (Figure 1F). Besides, KIAA1429 was the only m6A regulatory gene that significantly correlated with all the other sixteen m6A regulatory genes in the PPI network (Figure 1F).

Construction of a m6A-gene based prognostic signature

Univariate Cox regression analysis of the 398 training set LADC patients showed that six out of seventeen m6A regulatory genes, namely, IGF2BP1, IGF2BP2, IGF2BP3, HNRNPA2B1, METTL3, and HNRNPC were associated with OS (P < 0.05) (Figure 2A). Then, we analyzed these six genes using the LASSO Cox regression model and selected five genes, namely, IGF2BP1, IGF2BP2, HNRNPA2B1, METTL3, and HNRNPC to construct a five-gene prognostic signature (Figure 2B, 2C). We calculated their coefficients based on the lambda.min via 10-fold cross validation (Figure 2B–2D). The risk score was calculated as 0.037*IGF2BP1 + 0.004*IGF2BP2 + 0.007*HNRNPA2B1 + 0.022*HNRNPC + (-0.125)*METTL3. Then, we classified the training set LADC patients into high-risk and low-risk groups based on the median risk score value of 0.060. We then analyzed the distributions of risk scores and survival status of the training set LADC patients and found that survival times were significantly shorter for the high-risk patients compared to the low-risk group (Figure 2E, 2F).

Figure 2. Development of a five-gene prognostic signature. (A) The forest plots show association between expression levels of seventeen m6A regulatory genes and OS of LADC patients as determined by the univariate Cox regression model. (B) LASSO coefficients profiles of the six potentially prognostic m6A regulatory genes. (C) The coefficients of five prognostic m6A regulatory genes based on lambda.min using 10-fold cross validation. (D) The co-efficient distributions of the five prognostic m6A regulatory genes. (E) The distribution of risk scores in the training set LADC patients from the TCGA database. (F) The distribution of risk scores of the training set TCGA-LADC patients relative to their OS status.

Kaplan-Meier survival curve analysis confirmed that OS of the high-risk LADC patients in the training set was significantly shorter than the low-risk LADC patients (P < 0.0001) (Figure 3A). This suggested that the five-gene prognostic risk signature successfully distinguished high-risk and low-risk LADC patients. Furthermore, we performed ROC curve analysis and found that the AUC values of the five-gene prognostic signature for 3-year and 5-year survival of LADC patients were 0.684 and 0.646, respectively (Figure 3B). These AUC values were slightly lower than those considered as acceptable [35]. However, the AUC values for the five-gene prognostic signature were similar to the AUC values for the AJCC stages, which is a commonly used clinical prognostic indicator (Figure 3B). We then analyzed the association between the risk groups and their clinicopathological characteristics. The heat map showed significant differences in gender (P=0.041) and AJCC stages (P=0.028) between the low- and high-risk groups (Figure 3C). We then performed functional enrichment analysis of the differentially expressed genes between the low-risk and high-risk groups. GSEA results showed that the high-risk group was strongly associated with cancer-related pathways including mismatch repair (NES = 2.045, normalized P = 0.004, FDR q-val = 0.005), cell cycle (NES = 2.382, normalized P < 0.0001, FDR q-val < 0.0001), nucleotide excision repair (NES = 2.017, normalized P = 0.006, FDR q-val = 0.006), p53 signaling pathway (NES = 2.069, normalized P < 0.0001, FDR q-val = 0.004), and pathways in cancer (NES = 1.585, normalized P = 0.021, FDR q-val = 0.109) (Figure 3D).

Figure 3. The association of risk score with OS and clinicopathological characteristics. (A) Kaplan-Meier survival curve analysis of high-risk and low-risk LADC patients from the training set. (B) ROC curves show prediction accuracy of the prognostic signature and AJCC stages in the training set TCGA-LADC patients. (C) Heat map shows association between the expression levels of the five prognostic m6A regulatory genes and clinicopathological characteristics of high- and low-risk LADC patient subgroups. *0.01 ≤ P < 0.05, **0.001 ≤ P < 0.01, ***P < 0.001. (D) The positive correlation between cancer-related pathways and high-risk LADC patients.

The five m6A prognostic proteins are overexpressed in LADC samples

Western blot analysis of ten pairs of matched LADC and adjacent normal lung tissues showed that the expression levels of IGF2BP1, IGF2BP2, HNRNPA2B1, HNRNPC, and METTL3 proteins were significantly upregulated in LADC tissues compared to the corresponding normal lung tissues (Figure 4). This was consistent with the bioinformatics results for the TCGA-LADC dataset.

Figure 4. Expression levels of the five prognostic proteins in LADC patient samples. (A) Western blot analysis shows expression levels of IGF2BP1, IGF2BP2, HNRNPA2B1, HNRNPC, and METTL3 proteins in 10 paired normal lung and LADC tissue samples. (B–F) Relative expression levels of the five potentially prognostic m6A regulatory proteins, namely, IGF2BP1, IGF2BP2, HNRNPA2B1, HNRNPC, and METTL3 in 10 LADC tissues. GAPDH was used as loading control. The values were normalized by log₂ fold change (ratio of tumor to normal tissue expression) of the target proteins.

Evaluation of the five-m6A regulatory gene signature risk score as an independent prognostic predictor and construction of a predictive nomogram

Next, we performed univariate and multivariate Cox regression analysis to determine if the five-gene signature-derived risk score was an independent prognostic marker using the training cohort (n = 398) LADC patient data. The clinicopathological characteristics of the training cohort are shown in Table 1. Univariate Cox regression model showed that AJCC stage (HR, 1.629; 95% confidence interval (CI), 1.400-1.894; P < 0.0001) and the five-m6A regulatory gene signature risk score (HR, 3.332; 95% CI, 2.242-4.952; P < 0.0001) were prognostic factors for LADC patients (Figure 5A). Furthermore, multivariate Cox regression model also showed that AJCC stage (HR, 1.574; 95% CI, 1.353-1.831; P < 0.0001) and the five m6A regulatory gene signature risk score (HR, 3.119; 95% CI, 2.059-4.724; P < 0.0001) were independent prognostic factors for LADC patients (Figure 5B).

Table 1. Clinicopathological characteristics of patients in the training cohort and validation cohort.

Clinical characteristics	Training cohort		Validation cohort
	TCGA (No. (% ))		GSE72094 (No. (% ))	GSE41271 (No. (% ))
Total	398		320	168
Age at diagnosis
Median (IQR), years	65 (58-72)		71 (64-77)	64 (57-72)
Gender
Female	215 (54.0)		170 (53.1)	81 (48.2)
Male	183 (46.0)		150 (46.9)	87 (51.8)
Race
Black	45 (11.3)		9 (2.8)	11 (6.5)
White	353 (88.7)		311 (97.2)	157 (93.5)
AJCC Stages
Stage I	219 (55.0)		212 (66.2)	93 (55.4)
Stage II	92 (23.1)		52 (16.2)	27 (16.1)
Stage III	64 (16.1)		46 (14.4)	44 (26.2)
Stage IV	23 (5.8)		10 (3.1)	4 (2.4)
Smoking History
No	60 (15.1)		30 (9.4)	16 (9.5)
Yes	338 (84.9)		290 (90.6)	152 (90.5)
Survival Status
Alive	257 (64.6)		233 (72.8)	103 (61.3)
Dead	141 (35.4)		87 (27.2)	65 (38.7)
Survival Time (months)	21.0 (13.9-35.2)		28.1 (19.3-34.4)	39.0 (20.0-64.5)
Abbreviations: IQR: interquartile range; AJCC: American Joint Committee on Cancer.

Figure 5. Construction of a prognostic risk signature-based nomogram and evaluation of its performance in the TCGA dataset. (A, B) The forest plot shows association of clinicopathological parameters including risk score and the OS status of TCGA-LADC dataset as assessed by univariate and multivariate Cox regression models. (C) The nomogram with risk score and AJCC stages to predict 1-year, 3-year and 5-year OS of individual LADC patients. (D) The calibration curves and c-index values for the predicted 1-year, 3-year, and 5-year OS based on the nomogram for the TCGA-LADC patients from the training set. (E) ROC curve analysis shows the variations in the AUC values for the nomogram and AJCC stages in the training and validation cohorts for the 1-year to 5-year follow-up period. (F) The DCA analysis shows 1-year, 3-year and 5-year OS for LADC patients based on the nomogram.

We then constructed a nomogram based on coefficients from multivariate Cox regression model for the AJCC stage and the five m6A regulatory gene signature risk score and predicted the 1-year, 3-year and 5-year survival probabilities of LADC patients by calculating the total score for each patients by adding the variable values for AJCC stage and risk score (Figure 5C).

Next, we estimated the discrimination efficacy and prediction accuracy of the nomogram in the training set by evaluating c-index and AUC values as well as calibration curves. Our results showed that the nomogram was well calibrated because the curves were close to the diagonal line (Figure 5D). The c-index value for the nomogram was 0.71 (CI, 0.66-0.76) (Figure 5D), and the AUC value was 0.75 in the follow-up period (Figure 5E). Thus, the nomogram showed a higher discrimination power than the conventionally used prognostic index—AJCC stage (Figure 5E). Furthermore, decision curve analysis (DCA) showed that the threshold probabilities for 1-year, 3-year, and 5-year OS based on the nomogram were 0.06-0.48, 0.18-0.81, and 0.42-0.80, respectively (Figure 5F). The net clinical benefit gained from the nomogram was higher than hypothetical treat-all-patients or treat-none scenarios. These data showed that, the nomogram including AJCC stage and risk score was clinically beneficial for predicting 1-year, 3-year, and 5-year OS of LADC patients.

Validation of the prognostic model in LADC patient datasets from the GEO database

We further verified the reliability of the prognostic model using two validation LADC patient cohorts, GSE72094 (n=320) and GSE41271 (n=168). The clinicopathologic characteristics of the validation cohorts are shown in Table 1. We calculated the risk scores of all patients and classified them into high-risk and low-risk groups based on the previously determined cutoff value of 0.060. Kaplan Meier survival curve analyses showed that the OS was significantly shorter for LADC patients in the high-risk group compared to the low-risk group in both GSE72094 (P = 0.0006) and GSE41271 (P < 0.0001) datasets (Figure 6A, 6B). Furthermore, multivariate Cox regression analysis showed that the risk score was an independent prognostic factor for both the GSE72094 (HR, 6.772; 95% CI, 3.331-13.767; P < 0.0001) and GSE41271 (HR, 2.809; 95% CI, 1.791-4.406; P < 0.0001) datasets (Figure 6C, 6D). ROC curve analysis of the prognostic signature showed that the AUC values were 0.695 and 0.656, respectively, for the 3-year and 5-year OS of LADC patients in the GSE72094 dataset and 0.704 and 0.684, respectively, for the 3-year and 5-year OS of LADC patients in the GSE41271 dataset (Figure 6E, 6F). Furthermore, calibration curves for the 1-year, 3-year, and 5-year OS based on the nomogram were close to the diagonal line with a c-index value of 0.72 (CI, 0.67-0.77) for the GSE72094 dataset and 0.70 (CI, 0.63-0.77) for the GSE41271 dataset (Figure 6G, 6H). This suggested satisfactory performance of the prognostic model. The AUC values for both the GSE72094 and GSE41271 datasets were approximately 0.75 for the prognostic model in the 1-year, 3-year and 5-year follow-up period (Figure 5E). These values were higher than the AUC values of approximately 0.65 for the AJCC stage only.

Figure 6. Validation of the risk signature and nomogram using GEO datasets. (A, B) Kaplan-Meier survival curves show OS of high-risk and low-risk LADC patients in the GSE72094 and GSE41271 datasets. (C, D) The forest plot shows association between clinicopathological characteristics including risk score and OS status as assessed by the multivariate Cox regression model for the GSE72094 and GSE41271 datasets. (E, F) ROC curve analysis shows the prediction accuracy of the prognostic signature in the GSE72094 and GSE41271 datasets. (G, H) The calibration curves and c-index values for the predicted 1-year, 3-year, and 5-year OS in the LADC patients from the GSE72094 and GSE41271 datasets based on the nomogram.

Identification of co-expression modules and hub genes related to the five prognostic m6A regulatory genes using WGCNA

We compared 199 high-risk and 199 low-risk LADC patient groups in the training set and identified 741 differentially expressed genes (DEGs; 335 up-regulated and 406 down-regulated) (Figure 7A). We then developed a co-expression network of the DEGs using WGCNA. We achieved a higher scale-free topology fit index (> 0.8) and higher mean connectivity by constructing a scale-free network with the soft threshold power set as 8 (Figure 7B, 7C). These genes were classified by average linkage hierarchical clustering into three modules, namely, blue (209 genes), brown (173 genes), and turquoise (337 genes) (Figure 7D). We then investigated the association between clinicopathological characteristics and modules by calculating the correlation coefficients between module eigengenes (MEs) and clinicopathological characteristics. The turquoise module showed highest correlation with the prognostic risk score (Figure 7E, 7F).

Figure 7. Identification of gene modules using WGCNA. (A) Volcano plot shows 741 differentially expressed genes or DEGs (335 up-regulated and 406 down-regulated) in 199 high-risk and 199 low-risk LADC patient groups. (B) Analysis of scale-free topology fit index for various soft threshold powers (β). (C) Analysis of mean connectivity for various soft threshold powers. (D) Dendrogram of all differentially expressed genes clustered into various gene modules on the basis of a topological overlap-derived dissimilarity measure. The modules are displayed with different colors in the horizontal bar below the dendrogram. (E) Distribution of gene significance values in different modules. (F) The module-trait relationship heat map shows the correlation between module eigengenes (MEs) and clinicopathological characteristics. The correlation coefficients and p values are shown in the column for each ME and the corresponding clinicopathological characteristic.

GO and KEGG enrichment analysis of the turquoise module genes showed enrichment of genes involved in biological processes such as cell division, cell proliferation, and cell cycle, as well as, signaling pathways related to cell cycle, and p53 signaling pathway (Figure 8A, 8B).

Figure 8. Gene enrichment analysis in the turquoise module and determination of hub genes in the PPI network. (A) GO enrichment analysis of genes in the turquoise module. (B) KEGG pathway enrichment analysis of genes in the turquoise module. (C) Correlation analysis of genes in the turquoise module and the risk score. The cutoff values were set as gene significance (GS) > 0.55 and module membership (MM) > 0.85. (D) PPI network of all genes in the turquoise module including top 15 hub genes with the highest degree. PPI network was constructed using the STRING database and the Cytoscape software. (E) Venn diagram shows six common genes in the three lists of candidate hub genes. (F) The forest plot shows association between expression levels of the six hub genes and the OS status of LADC patients as assessed by the univariate Cox regression model.

We identified 22 genes as potential hub genes because they showed highly significant correlation with the risk score (GS > 0.55) and turquoise module (MM > 0.85) (Figure 8C). Based on connectivity of the weighted network, the top fifteen genes (top 10%) of the turquoise module were designated as the hub gene candidates. Furthermore, we performed PPI analysis of the 337 genes in the turquoise module using STRING database and cytoscape software, and identified the top 15 genes (top 10%) with the highest degree as the hub gene candidates (Figure 8D). By intersecting the three lists of candidate hub genes, we identified six common genes (CCNA2, CCNB1, BUB1B, BUB1, KIF2C, and KIF11) (Figure 8E). Moreover, the expression levels of these six genes was associated with the overall survival of LADC patients according to the univariate Cox regression model, thereby suggesting that these six hub genes were potential prognostic factors for LADC (Figure 8F).

Association of m6A regulatory gene prognostic signature with tumor immunity

We used the ESTIMATE algorithm [36] to evaluate tumor immune microenvironment in all LADC patient samples, including immune scores, stromal scores, and tumor purity. In general, we did not observe any significant differences in the immune scores, stromal scores, and tumor purity between the low-risk and high-risk groups (all P > 0.05, Figure 9A). Moreover, the correlation coefficients between the prognostic signature-based risk scores and tumor immune microenvironment scores were all below 0.2 (Figure 9B), which are regarded as negligible [37, 38]. We then investigated the association between the m6A prognostic signature and immune infiltration by evaluating the proportions of 21 immune cell types. The results showed significantly lower proportions of M1 macrophages, activated memory CD4⁺ activated T cells, CD8⁺ T cells, and neutrophils (Figure 9C). We also explored the association between immune checkpoint proteins and the m6A prognostic signature. The results showed significantly higher levels of PD-L1 and PD-L2 in high-risk patients compared to the low-risk patients (Figure 9D). Tumor mutational burden (TMB) is an important biomarker that is useful to predict response to PD-1/PD-L1 targeted immunotherapy in multiple cancer types [39–42]. We observed that the correlation between risk score and TMB was low (r = 0.28) (Figure 9E). Moreover, high-risk patients demonstrated significantly higher TMB (Figure 9F). This suggested that ICB therapy was potentially beneficial for high-risk LADC patients. Finally, in order to assess potential responses of different risk groups to ICB therapy, including CTLA-4 and PD-1 targeted therapy, we performed TIDE algorithm analysis and found that ICB response was significantly higher in the high-risk patients compared to the low-risk patients (Figure 9G). In summary, these findings demonstrate that ICB therapy is potentially beneficial for the high-risk LADC patients.

Figure 9. The association between m6A regulatory gene risk signature and tumor immunity of LADC patients. (A) Boxplot shows differences in immune scores, stromal scores and tumor purity between high-risk and low-risk LADC patient groups. (B) Spearman analysis shows the correlation between risk scores and immune microenvironment parameters such as the immune scores, stromal scores and tumor purity. (C) Relative proportions of infiltrating immune cell types in the high-risk and low-risk LADC patient groups. (D) Boxplot shows differences in expression levels of the immune checkpoint proteins between high-risk and low-risk LADC patient groups. (E) Spearman analysis shows correlation between TMB and risk scores. (F) Violin plot shows differences in TMB between high-risk and low-risk LADC patient groups. (G) Box plot shows the response to immune checkpoint therapy in high-risk and low-risk LADC patient groups.

Materials and Methods

Human lung adenocarcinoma tissue samples

We obtained paired LADC and adjacent noncancerous lung tissues from ten patients that were diagnosed through pathological examination at the First Affiliated Hospital of Zhejiang University. We obtained written informed consent from all patients before the tissues were collected. This study was approved by the Ethics Committee of the First Affiliated Hospital of Zhejiang University. The patient samples were frozen in liquid nitrogen and stored at −80° C before further use.

Human LADC patient datasets

The study design and strategy is shown in Figure 10. RNA-seq, clinicopathological, and single nucleotide variation data for 535, 398, and 561 LADC patients were downloaded from the TCGA database (https://portal.gdc.cancer.gov), respectively. RNA-seq data was normalized by the Fragment Per Kilobase Million (FPKM) method to obtain expression values of all genes. We used five LADC datasets from the GEO database (GSE72094, GSE41271, GSE31908, GSE26939, and GSE29013) for validation. Among these, the largest two datasets, namely, GSE72094 and GSE41271 were chosen for independent external validation. The gene expression profiles and clinicopathological data for these datasets were obtained from the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) database. Since only a small number of patients were classified as “Asian” or “American Indian/ Alaskan Native”, we focused our analysis on white and black patients.

Figure 10. Schematic representation of study strategy. The flow diagram shows five main sub-sections of this study—(1) bioinformatics analysis of m6A regulatory gene expression; (2) identification of a five-m6A regulatory gene prognostic signature and construction of a prognostic prediction nomogram; (3) external validation using two GEO patient datasets (GSE72094 and GSE41271); (4) WGCNA, and (5) correlation analysis of the prognostic signature with tumor immunity and immunotherapy response.

Identification of differentially expressed m6A regulatory genes

We systematically compared the mRNA expression levels of twenty m6A regulatory genes, namely, FTO, ALKBH5, ZC3H13, RBM15, METTL3, METTL14, METTL16, KIAA1429, WTAP, YTHDF1, YTHDF2, YTHDF3, YTHDC1, YTHDC2, IGF2BP1, IGF2BP2, IGF2BP3, HNRNPC, HNRNPG, and HNRNPA2B1 [18, 19, 50, 90] in 535 LADC samples and 59 normal lung tissue samples in the TCGA dataset. The differentially expressed genes were determined using the raw p-value corrected for false discovery rate (FDR) and visualized using heat maps and box plots.

Western blotting

Western blot analysis was performed as described previously [91] using anti-METTL3 (Cat. no. 15073-1-AP; 1:1,000), anti-IGF2BP1 (Cat. no. 22803-1-AP; 1:1,000), anti-IGF2BP2 (Cat. no. 11601-1-AP; 1:1,000), anti-HNRNPC (Cat. no. 11760-1-AP; 1:5,000), and anti-HNRNPA2B1 (Cat. no. 14813-1-AP; 1:2,000) antibodies, all of which were purchased from Proteintech. Anti-GADPH antibody (Affinity; Cat. no. #T0004; 1:3,000) was used to determine GAPDH protein levels. The target protein bands were quantified using the ImageJ v1.53a software (NIH, Bethesda, Maryland, USA). Differential expression was determined by calculating the ratios of target protein values for tumor tissues relative to the corresponding paired adjacent noncancerous lung tissues followed by log₂FoldChange for each patient.

Bioinformatics analysis

Gene set enrichment analysis (GSEA) was performed using the GSEA version 3.0 software (Broad Institute, Uc San Diego) to functionally annotate gene sets that show significant differences between high-risk and low-risk LADC patient groups which were classified based on the m6A signature. The cancer-related pathways were extracted based on the following criteria: (1) |normalized enrichment score (NES)|>1; (2) normalized P value < 0.05; (3) FDR q-value < 0.25 [92, 93]. The protein-protein interaction network between m6A regulatory genes was constructed and analyzed using the STRING database (https://string-db.org/).

Statistical analysis

We compared the expression levels of 20 m6A regulatory genes in 535 LADC samples and 59 normal lung tissue samples using Wilcoxon-Mann-Whitney test. We also used Wilcoxon-Mann-Whitney test to compare expression levels of m6A regulatory genes based on age, gender, race, and smoking history of LADC patients. Kruskal-wallis test was to compare expression levels of m6A regulatory genes based on American Joint Committee on Cancer (AJCC) stages. Spearman correlation analysis was performed to investigate the relationships between different m6A regulatory genes.

Univariate Cox regression model was used to determine prognostic m6A regulatory genes. Least absolute shrinkage and selection operator (LASSO) Cox regression model was used to construct optimal prognostic risk signature. We calculated the lambda value corresponding to minimum mean error (lambda.min) for the five prognostic m6A regulatory genes and evaluated their coefficients using10-fold cross validation. The risk score of all patients in the training and validation cohorts was calculated as the sum of expression level of each gene multiplied by its corresponding coefficient. The patients were dichotomized into low-risk and high-risk groups according to the median value of risk scores. The prediction accuracies of the prognostic risk signature and AJCC stages were assessed by the Receiver operating characteristic (ROC) curve analysis.

Chi-square tests were used to compare the frequency distributions of grouped variables (including age, gender, race, smoking history, and AJCC stage) between the low- and high-risk groups. Kaplan-Meier survival curves and log-rank test were used to compare survival times of low- and high-risk groups. Univariate Cox regression model was used to assess the association between variables including risk score, clinicopathological characteristics and OS in the training set. Multivariate Cox regression model was used to identify variables independently associated with overall survival (P < 0.05). For continuous variables such as the prognostic risk score, the restricted cubic splines (RCS) with three knots located at the empirical quantiles (10%, 50%, and 90 %) were fitted to relax the linearity assumption of the model [94]. Then, forest plots were drawn to better visualize the association between each prognostic variable and OS. A nomogram was constructed to predict the 1-year, 3-year, and 5-year survival probability of LADC patients. Furthermore, prognostic model performance in the training and validation sets was evaluated by concordance index (c-index) values, area under the ROC curve (AUC) values, and calibration curves. A decision curve analysis (DCA) was performed to assess the clinical utility of the nomogram [95]. All statistical analyses were performed using the R statistical software version 3.5.2 (R Foundation, Vienna, Austria). A two-tailed P<0.05 was considered statistically significant.

Weighted gene co-expression network analysis

Differentially expressed genes (DEGs) between the high-risk and low-risk groups were screened using false discovery rate (FDR) < 0.05 and |logFoldChange| > 1 as cut-off parameters. A co-expression network of the DEGs was then constructed using WGCNA. A weighted adjacency matrix was constructed based on the power value β, and transformed into a topological overlap matrix (TOM). Subsequently, gene modules were identified by the dynamic shear approach. The most risk-related modules were determined according to module significance (MS) and correlation coefficients between module eigengenes (MEs) and risk scores. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using the DAVID database (https://david.ncifcrf.gov/). The hub genes were determined according to the following criteria: (1) module membership (MM) > 0.85 and gene significance (GS) > 0.55); (2) top ten percent genes based on connectivity of the weighted network; (3) top ten percent genes based on degree of protein-protein interactions (PPI); (4) genes significantly associated with overall survival (P < 0.05).

Tumor immunity analysis

ESTIMATE algorithm was used to calculate immune score, stromal score, and tumor purity of all LADC patients [36]. CIBERSORT was used to determine proportions of 22 immune cell subtypes in each sample [96] using P-value < 0.05 as the threshold value. Wilcoxon-Mann-Whitney test was used to evaluate the association between m6A regulatory gene prognostic signature and expression levels of ten immune checkpoint proteins, namely, PD-L1, PD-1, PD-L2, CTLA-4, IDO1, LAG3, TIM-3, TIGIT, CD27, and ICOS. TMB was defined as the total number of the mutations in the somatic coding region per million bases. In this study, the mutational frequency of each sample was computed based on the number of variants divided by 38. As the length of exons was 38 million, the number of variants divided by 38 was just equal to the number of variants per million bases (TMB). Tumor Immune Dysfunction and Exclusion (TIDE) algorithm was used to evaluate the potential response to CTLA-4 and PD-1 targeting ICBs in the high- and low-risk patient groups [97].

Abbreviations

ADC: adenocarcinoma; AJCC: American Joint Committee on Cancer; AUC: area under the ROC curve; C-index: concordance index; DCA: decision curve analysis; DEGs: differentially expressed genes; FDR: false discovery rate; GEO: Gene Expression Omnibus; GO: Gene Ontology; GS: Gene significance; GSEA: Gene set enrichment analysis; ICB: immune checkpoint blockade; KEGG: Kyoto Encyclopedia of Genes and Genomes; LASSO: least absolute shrinkage and selection operator; lncRNAs: long non-coding RNAs; m6A: N6-methyladenosine; MEs: Module eigengenes; miRNAs: microRNAs; MM: Module membership; MS: module significance; OS: overall survival; PPI: protein-protein interaction; RCS: restricted cubic splines; ROC: receiver operating characteristic; TCGA: The Cancer Genome Atlas; TIDE: Tumor Immune Dysfunction and Exclusion; TMB: tumor mutational burden; US: United States; WGCNA: Weighted gene co-expression network analysis.

Author Contributions

Wu X, Sheng HX, and Hu J designed the study and wrote the manuscript; Wu X, Xia PH, Wang YQ, Yu L, Lv W, and Hu J participated in literature search and data acquisition, analysis, and interpretation; Wu X, Sheng HX, and Hu J contributed to the revision of manuscript. All authors approved the final manuscript.

Acknowledgments

The authors are sincerely grateful to the TCGA and GEO databases for the access to data that was used in this study.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Funding

This work was supported by the General Project of Zhejiang Natural Science Foundation (Grant No. LY19H160039), National Key Research and Development Program of China (Grant No. 2017YFC0113500), Zhejiang Province Major Science and Technology Special Project (Grant No. 2014C03032), Zhejiang Provincial Key Discipline of Traditional Chinese Medicine (Grant No. 2017-XK-A33), and Zhejiang Lung Cancer Center (Grant No. JBZX-202007).

References

• 1. Siegel RL, Miller KD, Jemal A. Cancer statistics, 2017. CA Cancer J Clin. 2017; 67:7–30. https://doi.org/10.3322/caac.21387 [PubMed]
• 2. American Cancer Society. Cancer facts and figures 2020. http://www.cancer.org/content/dam/cancer-org/research/cancer-facts-and-statistics/annual-cancer-facts-and-figures/2020/cancer-facts-and-figures-2020.pdf.
• 3. National Cancer Institute and Surveillance Epidemiology and End Results (SEER) Program. Cancer of the Lung and Bronchus. https://seer.cancer.gov/csr/1975_2016/sections.html.
• 4. Lewis DR, Check DP, Caporaso NE, Travis WD, Devesa SS. US lung cancer trends by histologic type. Cancer. 2014; 120:2883–92. https://doi.org/10.1002/cncr.28749 [PubMed]
• 5. Meza R, Meernik C, Jeon J, Cote ML. Lung cancer incidence trends by gender, race and histology in the United States, 1973-2010. PLoS One. 2015; 10:e0121323. https://doi.org/10.1371/journal.pone.0121323 [PubMed]
• 6. Goldstraw P, Chansky K, Crowley J, Rami-Porta R, Asamura H, Eberhardt WE, Nicholson AG, Groome P, Mitchell A, Bolejack V, and International Association for the Study of Lung Cancer Staging and Prognostic Factors Committee, Advisory Boards, and Participating Institutions, and International Association for the Study of Lung Cancer Staging and Prognostic Factors Committee Advisory Boards and Participating Institutions. The IASLC lung cancer staging project: proposals for revision of the TNM stage groupings in the forthcoming (eighth) edition of the TNM classification for lung cancer. J Thorac Oncol. 2016; 11:39–51. https://doi.org/10.1016/j.jtho.2015.09.009 [PubMed]
• 7. Wang BY, Huang JY, Chen HC, Lin CH, Lin SH, Hung WH, Cheng YF. The comparison between adenocarcinoma and squamous cell carcinoma in lung cancer patients. J Cancer Res Clin Oncol. 2020; 146:43–52. https://doi.org/10.1007/s00432-019-03079-8 [PubMed]
• 8. Desrosiers R, Friderici K, Rottman F. Identification of methylated nucleosides in messenger RNA from novikoff hepatoma cells. Proc Natl Acad Sci USA. 1974; 71:3971–75. https://doi.org/10.1073/pnas.71.10.3971 [PubMed]
• 9. Alarcón CR, Lee H, Goodarzi H, Halberg N, Tavazoie SF. N6-methyladenosine marks primary microRNAs for processing. Nature. 2015; 519:482–85. https://doi.org/10.1038/nature14281 [PubMed]
• 10. Patil DP, Chen CK, Pickering BF, Chow A, Jackson C, Guttman M, Jaffrey SR. m(6)A RNA methylation promotes XIST-mediated transcriptional repression. Nature. 2016; 537:369–73. https://doi.org/10.1038/nature19342 [PubMed]
• 11. Chen XY, Zhang J, Zhu JS. The role of m⁶A RNA methylation in human cancer. Mol Cancer. 2019; 18:103. https://doi.org/10.1186/s12943-019-1033-z [PubMed]
• 12. Meyer KD, Saletore Y, Zumbo P, Elemento O, Mason CE, Jaffrey SR. Comprehensive analysis of mRNA methylation reveals enrichment in UTRs and near stop codons. Cell. 2012; 149:1635–46. https://doi.org/10.1016/j.cell.2012.05.003 [PubMed]
• 13. Ke S, Alemu EA, Mertens C, Gantman EC, Fak JJ, Mele A, Haripal B, Zucker-Scharff I, Moore MJ, Park CY, Vågbø CB, Kusśnierczyk A, Klungland A, et al. A majority of m6A residues are in the last exons, allowing the potential for 3' UTR regulation. Genes Dev. 2015; 29:2037–53. https://doi.org/10.1101/gad.269415.115 [PubMed]
• 14. Jia G, Fu Y, Zhao X, Dai Q, Zheng G, Yang Y, Yi C, Lindahl T, Pan T, Yang YG, He C. N6-methyladenosine in nuclear RNA is a major substrate of the obesity-associated FTO. Nat Chem Biol. 2011; 7:885–87. https://doi.org/10.1038/nchembio.687 [PubMed]
• 15. Yang Y, Hsu PJ, Chen YS, Yang YG. Dynamic transcriptomic m⁶A decoration: writers, erasers, readers and functions in RNA metabolism. Cell Res. 2018; 28:616–24. https://doi.org/10.1038/s41422-018-0040-8 [PubMed]
• 16. Haussmann IU, Bodi Z, Sanchez-Moran E, Mongan NP, Archer N, Fray RG, Soller M. m⁶A potentiates Sxl alternative pre-mRNA splicing for robust drosophila sex determination. Nature. 2016; 540:301–04. https://doi.org/10.1038/nature20577 [PubMed]
• 17. Liu N, Dai Q, Zheng G, He C, Parisien M, Pan T. N(6)-methyladenosine-dependent RNA structural switches regulate RNA-protein interactions. Nature. 2015; 518:560–64. https://doi.org/10.1038/nature14234 [PubMed]
• 18. Liu L, Wang Y, Wu J, Liu J, Qin Z, Fan H. N⁶-methyladenosine: A potential breakthrough for human cancer. Mol Ther Nucleic Acids. 2020; 19:804–13. https://doi.org/10.1016/j.omtn.2019.12.013 [PubMed]
• 19. Sun T, Wu R, Ming L. The role of m6A RNA methylation in cancer. Biomed Pharmacother. 2019; 112:108613. https://doi.org/10.1016/j.biopha.2019.108613 [PubMed]
• 20. Schöller E, Weichmann F, Treiber T, Ringle S, Treiber N, Flatley A, Feederle R, Bruckmann A, Meister G. Interactions, localization, and phosphorylation of the m⁶A generating METTL3-METTL14-WTAP complex. RNA. 2018; 24:499–512. https://doi.org/10.1261/rna.064063.117 [PubMed]
• 21. Tang C, Klukovich R, Peng H, Wang Z, Yu T, Zhang Y, Zheng H, Klungland A, Yan W. ALKBH5-dependent m6A demethylation controls splicing and stability of long -UTR mRNAs in male germ cells. Proc Natl Acad Sci USA. 2018; 115:E325–33. https://doi.org/10.1073/pnas.1717794115 [PubMed]
• 22. Liu N, Pan T. N6-methyladenosine–encoded epitranscriptomics. Nat Struct Mol Biol. 2016; 23:98–102. https://doi.org/10.1038/nsmb.3162 [PubMed]
• 23. Wang S, Sun C, Li J, Zhang E, Ma Z, Xu W, Li H, Qiu M, Xu Y, Xia W, Xu L, Yin R. Roles of RNA methylation by means of N⁶-methyladenosine (m⁶A) in human cancers. Cancer Lett. 2017; 408:112–20. https://doi.org/10.1016/j.canlet.2017.08.030 [PubMed]
• 24. Chen M, Wei L, Law CT, Tsang FH, Shen J, Cheng CL, Tsang LH, Ho DW, Chiu DK, Lee JM, Wong CC, Ng IO, Wong CM. RNA N6-methyladenosine methyltransferase-like 3 promotes liver cancer progression through YTHDF2-dependent posttranscriptional silencing of SOCS2. Hepatology. 2018; 67:2254–70. https://doi.org/10.1002/hep.29683 [PubMed]
• 25. Deng R, Cheng Y, Ye S, Zhang J, Huang R, Li P, Liu H, Deng Q, Wu X, Lan P, Deng Y. m⁶A methyltransferase METTL3 suppresses colorectal cancer proliferation and migration through p38/ERK pathways. Onco Targets Ther. 2019; 12:4391–402. https://doi.org/10.2147/OTT.S201052 [PubMed]
• 26. He L, Li H, Wu A, Peng Y, Shu G, Yin G. Functions of N6-methyladenosine and its role in cancer. Mol Cancer. 2019; 18:176. https://doi.org/10.1186/s12943-019-1109-9 [PubMed]
• 27. Huang Y, Su R, Sheng Y, Dong L, Dong Z, Xu H, Ni T, Zhang ZS, Zhang T, Li C, Han L, Zhu Z, Lian F, et al. Small-molecule targeting of oncogenic FTO demethylase in acute myeloid leukemia. Cancer Cell. 2019; 35:677–91.e10. https://doi.org/10.1016/j.ccell.2019.03.006 [PubMed]
• 28. Huang Y, Yan J, Li Q, Li J, Gong S, Zhou H, Gan J, Jiang H, Jia GF, Luo C, Yang CG. Meclofenamic acid selectively inhibits FTO demethylation of m6A over ALKBH5. Nucleic Acids Res. 2015; 43:373–84. https://doi.org/10.1093/nar/gku1276 [PubMed]
• 29. Sun L, Liu WK, Du XW, Liu XL, Li G, Yao Y, Han T, Li WY, Gu J. Large-scale transcriptome analysis identified RNA methylation regulators as novel prognostic signatures for lung adenocarcinoma. Ann Transl Med. 2020; 8:751. https://doi.org/10.21037/atm-20-3744 [PubMed]
• 30. Xu F, Chen JX, Yang XB, Hong XB, Li ZX, Lin L, Chen YS. Analysis of lung adenocarcinoma subtypes based on immune signatures identifies clinical implications for cancer therapy. Mol Ther Oncolytics. 2020; 17:241–49. https://doi.org/10.1016/j.omto.2020.03.021 [PubMed]
• 31. Li F, Wang H, Huang H, Zhang L, Wang D, Wan Y. m6A RNA methylation regulators participate in the Malignant progression and have clinical prognostic value in lung adenocarcinoma. Front Genet. 2020; 11:994. https://doi.org/10.3389/fgene.2020.00994 [PubMed]
• 32. Zhuang Z, Chen L, Mao Y, Zheng Q, Li H, Huang Y, Hu Z, Jin Y. Diagnostic, progressive and prognostic performance of m⁶A methylation RNA regulators in lung adenocarcinoma. Int J Biol Sci. 2020; 16:1785–97. https://doi.org/10.7150/ijbs.39046 [PubMed]
• 33. Zhang Y, Liu X, Liu L, Li J, Hu Q, Sun R. Expression and prognostic significance of m6A-related genes in lung adenocarcinoma. Med Sci Monit. 2020; 26:e919644. https://doi.org/10.12659/MSM.919644 [PubMed]
• 34. Liu Y, Guo X, Zhao M, Ao H, Leng X, Liu M, Wu C, Ma J, Zhu J. Contributions and prognostic values of m⁶A RNA methylation regulators in non-small-cell lung cancer. J Cell Physiol. 2020; 235:6043–57. https://doi.org/10.1002/jcp.29531 [PubMed]
• 35. Mandrekar JN. Receiver operating characteristic curve in diagnostic test assessment. J Thorac Oncol. 2010; 5:1315–16. https://doi.org/10.1097/JTO.0b013e3181ec173d [PubMed]
• 36. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, Treviño V, Shen H, Laird PW, Levine DA, Carter SL, Getz G, Stemke-Hale K, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013; 4:2612. https://doi.org/10.1038/ncomms3612 [PubMed]
• 37. Crane AK, Geller EJ, Myers EM, Fenderson JL, Wells E, Jannelli M, Connolly A, Matthews CA. Implementation of a standardized digital rectal exam to improve the accuracy of rectocele diagnosis. Int Urogynecol J. 2015; 26:107–11. https://doi.org/10.1007/s00192-014-2479-6 [PubMed]
• 38. Prasad K, Dash D, Kumar A. Validation of the hindi version of national institute of health stroke scale. Neurol India. 2012; 60:40–44. https://doi.org/10.4103/0028-3886.93587 [PubMed]
• 39. Goodman AM, Kato S, Bazhenova L, Patel SP, Frampton GM, Miller V, Stephens PJ, Daniels GA, Kurzrock R. Tumor mutational burden as an independent predictor of response to immunotherapy in diverse cancers. Mol Cancer Ther. 2017; 16:2598–608. https://doi.org/10.1158/1535-7163.MCT-17-0386 [PubMed]
• 40. Hellmann MD, Ciuleanu TE, Pluzanski A, Lee JS, Otterson GA, Audigier-Valette C, Minenza E, Linardou H, Burgers S, Salman P, Borghaei H, Ramalingam SS, Brahmer J, et al. Nivolumab plus ipilimumab in lung cancer with a high tumor mutational burden. N Engl J Med. 2018; 378:2093–104. https://doi.org/10.1056/NEJMoa1801946 [PubMed]
• 41. Wang Z, Duan J, Cai S, Han M, Dong H, Zhao J, Zhu B, Wang S, Zhuo M, Sun J, Wang Q, Bai H, Han J, et al. Assessment of blood tumor mutational burden as a potential biomarker for immunotherapy in patients with non-small cell lung cancer with use of a next-generation sequencing cancer gene panel. JAMA Oncol. 2019; 5:696–702. https://doi.org/10.1001/jamaoncol.2018.7098 [PubMed]
• 42. Mutation burden predicts anti-PD-1 response. Cancer Discov. 2018; 8:258. https://doi.org/10.1158/2159-8290.CD-NB2018-005 [PubMed]
• 43. Hirsch FR, Scagliotti GV, Mulshine JL, Kwon R, Curran WJ Jr, Wu YL, Paz-Ares L. Lung cancer: current therapies and new targeted treatments. Lancet. 2017; 389:299–311. https://doi.org/10.1016/S0140-6736(16)30958-8 [PubMed]
• 44. Shao Q, Liu S, Wang W, Zhang Y, Li F, Li J. Clinical investigation into the initial diagnosis and treatment of 539 patients with stage IV lung cancer. Onco Targets Ther. 2017; 10:535–41. https://doi.org/10.2147/OTT.S125108 [PubMed]
• 45. Chen Y, Peng C, Chen J, Chen D, Yang B, He B, Hu W, Zhang Y, Liu H, Dai L, Xie H, Zhou L, Wu J, Zheng S. WTAP facilitates progression of hepatocellular carcinoma via m6A-HuR-dependent epigenetic silencing of ETS1. Mol Cancer. 2019; 18:127. https://doi.org/10.1186/s12943-019-1053-8 [PubMed]
• 46. Bansal H, Yihua Q, Iyer SP, Ganapathy S, Proia DA, Penalva LO, Uren PJ, Suresh U, Carew JS, Karnad AB, Weitman S, Tomlinson GE, Rao MK, et al. WTAP is a novel oncogenic protein in acute myeloid leukemia. Leukemia. 2014; 28:1171–74. https://doi.org/10.1038/leu.2014.16 [PubMed]
• 47. Jin DI, Lee SW, Han ME, Kim HJ, Seo SA, Hur GY, Jung S, Kim BS, Oh SO. Expression and roles of Wilms’ tumor 1-associating protein in glioblastoma. Cancer Sci. 2012; 103:2102–09. https://doi.org/10.1111/cas.12022 [PubMed]
• 48. Zhu H, Gan X, Jiang X, Diao S, Wu H, Hu J. ALKBH5 inhibited autophagy of epithelial ovarian cancer through miR-7 and BCL-2. J Exp Clin Cancer Res. 2019; 38:163. https://doi.org/10.1186/s13046-019-1159-2 [PubMed]
• 49. Bai Y, Yang C, Wu R, Huang L, Song S, Li W, Yan P, Lin C, Li D, Zhang Y. YTHDF1 regulates tumorigenicity and cancer stem cell-like activity in human colorectal carcinoma. Front Oncol. 2019; 9:332. https://doi.org/10.3389/fonc.2019.00332 [PubMed]
• 50. Chai RC, Wu F, Wang QX, Zhang S, Zhang KN, Liu YQ, Zhao Z, Jiang T, Wang YZ, Kang CS. m⁶A RNA methylation regulators contribute to Malignant progression and have clinical prognostic impact in gliomas. Aging (Albany NY). 2019; 11:1204–25. https://doi.org/10.18632/aging.101829 [PubMed]
• 51. Wang J, Zhang C, He W, Gou X. Effect of m⁶A RNA methylation regulators on Malignant progression and prognosis in renal clear cell carcinoma. Front Oncol. 2020; 10:3. https://doi.org/10.3389/fonc.2020.00003 [PubMed]
• 52. Su Y, Huang J, Hu J. m⁶A RNA methylation regulators contribute to Malignant progression and have clinical prognostic impact in gastric cancer. Front Oncol. 2019; 9:1038. https://doi.org/10.3389/fonc.2019.01038 [PubMed]
• 53. Zhou Y, Yin Z, Hou B, Yu M, Chen R, Jin H, Jian Z. Expression profiles and prognostic significance of RNA N6-methyladenosine-related genes in patients with hepatocellular carcinoma: evidence from independent datasets. Cancer Manag Res. 2019; 11:3921–31. https://doi.org/10.2147/CMAR.S191565 [PubMed]
• 54. Liu Y, Wu L, Ao H, Zhao M, Leng X, Liu M, Ma J, Zhu J. Prognostic implications of autophagy-associated gene signatures in non-small cell lung cancer. Aging (Albany NY). 2019; 11:11440–62. https://doi.org/10.18632/aging.102544 [PubMed]
• 55. Xu Z, Wu Z, Zhang J, Zhou R, Ye L, Yang P, Yu B. Development and validation of an oxidative phosphorylation-related gene signature in lung adenocarcinoma. Epigenomics. 2020; 12:1333–48. https://doi.org/10.2217/epi-2020-0217 [PubMed]
• 56. Mo Z, Yu L, Cao Z, Hu H, Luo S, Zhang S. Identification of a hypoxia-associated signature for lung adenocarcinoma. Front Genet. 2020; 11:647. https://doi.org/10.3389/fgene.2020.00647 [PubMed]
• 57. Yu X, Zhang X, Zhang Y. Identification of a 5-gene metabolic signature for predicting prognosis based on an integrated analysis of tumor microenvironment in lung adenocarcinoma. J Oncol. 2020; 2020:5310793. https://doi.org/10.1155/2020/5310793 [PubMed]
• 58. Song C, Guo Z, Yu D, Wang Y, Wang Q, Dong Z, Hu W. A prognostic nomogram combining immune-related gene signature and clinical factors predicts survival in patients with lung adenocarcinoma. Front Oncol. 2020; 10:1300. https://doi.org/10.3389/fonc.2020.01300 [PubMed]
• 59. Yan M, Sun L, Li J, Yu H, Lin H, Yu T, Zhao F, Zhu M, Liu L, Geng Q, Kong H, Pan H, Yao M. RNA-binding protein KHSRP promotes tumor growth and metastasis in non-small cell lung cancer. J Exp Clin Cancer Res. 2019; 38:478. https://doi.org/10.1186/s13046-019-1479-2 [PubMed]
• 60. Gong F, Ren P, Zhang Y, Jiang J, Zhang H. MicroRNAs-491-5p suppresses cell proliferation and invasion by inhibiting IGF2BP1 in non-small cell lung cancer. Am J Transl Res. 2016; 8:485–95. [PubMed]
• 61. Gutschner T, Hämmerle M, Pazaitis N, Bley N, Fiskin E, Uckelmann H, Heim A, Groβ M, Hofmann N, Geffers R, Skawran B, Longerich T, Breuhahn K, et al. Insulin-like growth factor 2 mRNA-binding protein 1 (IGF2BP1) is an important protumorigenic factor in hepatocellular carcinoma. Hepatology. 2014; 59:1900–11. https://doi.org/10.1002/hep.26997 [PubMed]
• 62. Zhu J, Wang M, Hu D. Deciphering N⁶-methyladenosine-related genes signature to predict survival in lung adenocarcinoma. Biomed Res Int. 2020; 2020:2514230. https://doi.org/10.1155/2020/2514230 [PubMed]
• 63. Cheng M, Sheng L, Gao Q, Xiong Q, Zhang H, Wu M, Liang Y, Zhu F, Zhang Y, Zhang X, Yuan Q, Li Y. The m6A methyltransferase METTL3 promotes bladder cancer progression via AFF4/NF-κB/MYC signaling network. Oncogene. 2019; 38:3667–80. https://doi.org/10.1038/s41388-019-0683-z [PubMed]
• 64. Du M, Zhang Y, Mao Y, Mou J, Zhao J, Xue Q, Wang D, Huang J, Gao S, Gao Y. MiR-33a suppresses proliferation of NSCLC cells via targeting METTL3 mRNA. Biochem Biophys Res Commun. 2017; 482:582–89. https://doi.org/10.1016/j.bbrc.2016.11.077 [PubMed]
• 65. Zheng W, Dong X, Zhao Y, Wang S, Jiang H, Zhang M, Zheng X, Gu M. Multiple functions and mechanisms underlying the role of METTL3 in human cancers. Front Oncol. 2019; 9:1403. https://doi.org/10.3389/fonc.2019.01403 [PubMed]
• 66. Huang RS, Zheng YL, Li C, Ding C, Xu C, Zhao J. MicroRNA-485-5p suppresses growth and metastasis in non-small cell lung cancer cells by targeting IGF2BP2. Life Sci. 2018; 199:104–11. https://doi.org/10.1016/j.lfs.2018.03.005 [PubMed]
• 67. Xu X, Yu Y, Zong K, Lv P, Gu Y. Up-regulation of IGF2BP2 by multiple mechanisms in pancreatic cancer promotes cancer proliferation by activating the PI3K/Akt signaling pathway. J Exp Clin Cancer Res. 2019; 38:497. https://doi.org/10.1186/s13046-019-1470-y [PubMed]
• 68. Xuan Y, Wang J, Ban L, Lu JJ, Yi C, Li Z, Yu W, Li M, Xu T, Yang W, Tang Z, Tang R, Xiao X, et al. hnRNPA2/B1 activates cyclooxygenase-2 and promotes tumor growth in human lung cancers. Mol Oncol. 2016; 10:610–24. https://doi.org/10.1016/j.molonc.2015.11.010 [PubMed]
• 69. Zwicker J, Lucibello FC, Wolfraim LA, Gross C, Truss M, Engeland K, Müller R. Cell cycle regulation of the cyclin A, cdc25C and cdc2 genes is based on a common mechanism of transcriptional repression. EMBO J. 1995; 14:4514–22. https://doi.org/10.1002/j.1460-2075.1995.tb00130.x [PubMed]
• 70. Roberts JM. Evolving ideas about cyclins. Cell. 1999; 98:129–32. https://doi.org/10.1016/s0092-8674(00)81007-7 [PubMed]
• 71. Bharadwaj R, Yu H. The spindle checkpoint, aneuploidy, and cancer. Oncogene. 2004; 23:2016–27. https://doi.org/10.1038/sj.onc.1207374 [PubMed]
• 72. Logarinho E, Bousbaa H. Kinetochore-microtubule interactions “in check” by Bub1, Bub3 and BubR1: the dual task of attaching and signalling. Cell Cycle. 2008; 7:1763–68. https://doi.org/10.4161/cc.7.12.6180 [PubMed]
• 73. Sacristan C, Kops GJ. Joined at the hip: kinetochores, microtubules, and spindle assembly checkpoint signaling. Trends Cell Biol. 2015; 25:21–28. https://doi.org/10.1016/j.tcb.2014.08.006 [PubMed]
• 74. Ganguly A, Bhattacharya R, Cabral F. Cell cycle dependent degradation of MCAK: evidence against a role in anaphase chromosome movement. Cell Cycle. 2008; 7:3187–93. https://doi.org/10.4161/cc.7.20.6814 [PubMed]
• 75. Rath O, Kozielski F. Kinesins and cancer. Nat Rev Cancer. 2012; 12:527–39. https://doi.org/10.1038/nrc3310 [PubMed]
• 76. Huang T, Gao Q, Feng T, Zheng Y, Guo J, Zeng W. FTO knockout causes chromosome instability and G2/M arrest in mouse GC-1 cells. Front Genet. 2019; 9:732. https://doi.org/10.3389/fgene.2018.00732 [PubMed]
• 77. Wu R, Liu Y, Yao Y, Zhao Y, Bi Z, Jiang Q, Liu Q, Cai M, Wang F, Wang Y, Wang X. FTO regulates adipogenesis by controlling cell cycle progression via m⁶A-YTHDF2 dependent mechanism. Biochim Biophys Acta Mol Cell Biol Lipids. 2018; 1863:1323–30. https://doi.org/10.1016/j.bbalip.2018.08.008 [PubMed]
• 78. Song T, Zheng Y, Wang Y, Katz Z, Liu X, Chen S, Singer RH, Gu W. Specific interaction of KIF11 with ZBP1 regulates the transport of β-actin mRNA and cell motility. J Cell Sci. 2015; 128:1001–10. https://doi.org/10.1242/jcs.161679 [PubMed]
• 79. Wang SS, Lv Y, Xu XC, Zuo Y, Song Y, Wu GP, Lu PH, Zhang ZQ, Chen MB. Triptonide inhibits human nasopharyngeal carcinoma cell growth via disrupting lnc-RNA THOR-IGF2BP1 signaling. Cancer Lett. 2019; 443:13–24. https://doi.org/10.1016/j.canlet.2018.11.028 [PubMed]
• 80. Li L, Yang Y, Wu M, Yu Z, Wang C, Dou G, He H, Wang H, Yang N, Qi H, Xu X. β-Asarone Induces Apoptosis and Cell Cycle Arrest of Human Glioma U251 Cells via Suppression of HnRNP A2/B1-Mediated Pathway In Vitro and In Vivo. Molecules. 2018; 23:1072. https://doi.org/10.3390/molecules23051072 [PubMed]
• 81. Han J, Wang JZ, Yang X, Yu H, Zhou R, Lu HC, Yuan WB, Lu JC, Zhou ZJ, Lu Q, Wei JF, Yang H. METTL3 promote tumor proliferation of bladder cancer by accelerating pri-miR221/222 maturation in m6A-dependent manner. Mol Cancer. 2019; 18:110. https://doi.org/10.1186/s12943-019-1036-9 [PubMed]
• 82. Piñol-Roma S, Dreyfuss G. Cell cycle-regulated phosphorylation of the pre-mRNA-binding (heterogeneous nuclear ribonucleoprotein) C proteins. Mol Cell Biol. 1993; 13:5762–70. https://doi.org/10.1128/mcb.13.9.5762 [PubMed]
• 83. Hellmann MD, Paz-Ares L, Bernabe Caro R, Zurawski B, Kim SW, Carcereny Costa E, Park K, Alexandru A, Lupinacci L, de la Mora Jimenez E, Sakai H, Albert I, Vergnenegre A, et al. Nivolumab plus ipilimumab in advanced non-small-cell lung cancer. N Engl J Med. 2019; 381:2020–31. https://doi.org/10.1056/NEJMoa1910231 [PubMed]
• 84. Escudier B, Motzer RJ, Sharma P, Wagstaff J, Plimack ER, Hammers HJ, Donskov F, Gurney H, Sosman JA, Zalewski PG, Harmenberg U, McDermott DF, Choueiri TK, et al. Treatment beyond progression in patients with advanced renal cell carcinoma treated with nivolumab in CheckMate 025. Eur Urol. 2017; 72:368–76. https://doi.org/10.1016/j.eururo.2017.03.037 [PubMed]
• 85. Paramasivam A, Priyadharsini JV, Raghunandhakumar S. Implications of m6A modification in autoimmune disorders. Cell Mol Immunol. 2020; 17:550–51. https://doi.org/10.1038/s41423-019-0307-0 [PubMed]
• 86. Zhang C, Fu J, Zhou Y. A review in research progress concerning m6A methylation and immunoregulation. Front Immunol. 2019; 10:922. https://doi.org/10.3389/fimmu.2019.00922 [PubMed]
• 87. Garon EB, Rizvi NA, Hui R, Leighl N, Balmanoukian AS, Eder JP, Patnaik A, Aggarwal C, Gubens M, Horn L, Carcereny E, Ahn MJ, Felip E, et al, and KEYNOTE-001 Investigators. Pembrolizumab for the treatment of non-small-cell lung cancer. N Engl J Med. 2015; 372:2018–28. https://doi.org/10.1056/NEJMoa1501824 [PubMed]
• 88. Boussiotis VA. Molecular and biochemical aspects of the PD-1 checkpoint pathway. N Engl J Med. 2016; 375:1767–78. https://doi.org/10.1056/NEJMra1514296 [PubMed]
• 89. Leng C, Li Y, Qin J, Ma J, Liu X, Cui Y, Sun H, Wang Z, Hua X, Yu Y, Li H, Zhang J, Zheng Y, et al. Relationship between expression of PD-L1 and PD-L2 on esophageal squamous cell carcinoma and the antitumor effects of CD8⁺ T cells. Oncol Rep. 2016; 35:699–708. https://doi.org/10.3892/or.2015.4435 [PubMed]
• 90. Zhao X, Cui L. Development and validation of a m⁶A RNA methylation regulators-based signature for predicting the prognosis of head and neck squamous cell carcinoma. Am J Cancer Res. 2019; 9:2156–69. [PubMed]
• 91. Yuan XS, Wang ZT, Hu YJ, Bao FC, Yuan P, Zhang C, Cao JL, Lv W, Hu J. Downregulation of RUVBL1 inhibits proliferation of lung adenocarcinoma cells by G1/S phase cell cycle arrest via multiple mechanisms. Tumour Biol. 2016; 37:16015–27. https://doi.org/10.1007/s13277-016-5452-9 [PubMed]
• 92. Zeng L, Fan X, Wang X, Deng H, Zhang K, Zhang X, He S, Li N, Han Q, Liu Z. Bioinformatics analysis based on multiple databases identifies hub genes associated with hepatocellular carcinoma. Curr Genomics. 2019; 20:349–61. https://doi.org/10.2174/1389202920666191011092410 [PubMed]
• 93. Xu Z, Wang N, Xu Y, Hua L, Zhou D, Zheng M, Deng X. Effects of chronic PM_2.5 exposure on pulmonary epithelia: transcriptome analysis of mRNA-exosomal miRNA interactions. Toxicol Lett. 2019; 316:49–59. https://doi.org/10.1016/j.toxlet.2019.09.010 [PubMed]
• 94. Harrell F Jr. Regression modeling strategies. Springer International Publishing Switzerland. 2015. https://doi.org/10.1007/978-3-319-19425-7
• 95. Vickers AJ, Cronin AM, Elkin EB, Gonen M. Extensions to decision curve analysis, a novel method for evaluating diagnostic tests, prediction models and molecular markers. BMC Med Inform Decis Mak. 2008; 8:53. https://doi.org/10.1186/1472-6947-8-53 [PubMed]
• 96. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, Hoang CD, Diehn M, Alizadeh AA. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015; 12:453–57. https://doi.org/10.1038/nmeth.3337 [PubMed]
• 97. Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, Li Z, Traugh N, Bu X, Li B, Liu J, Freeman GJ, Brown MA, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med. 2018; 24:1550–58. https://doi.org/10.1038/s41591-018-0136-1 [PubMed]