Introduction
============

According to Global Cancer Statistics 2018, there will have been approximately 18.1 million new cancer cases and 9.6 million death cases worldwide [@B1]. Researchers around the world are constantly striving to improve the medical technology for providing more sensitive method for tumor diagnosis and effective approach for tumor therapy. However, due to the complexity of tumor formation mechanism, it is far from enough to understand the nature of cancer from the genetic layer in the traditional sense [@B2]. It is recognized that the expression of oncogene depends not only on the gene itself, but also on epigenetic modification without changing gene sequence [@B2],[@B3].

Epigenetics is a research hotspot in recent years, which is defined as no change in DNA sequence, but with heritable changes in gene expression [@B3]. Previously, the epigenetic researchers are mainly focused on DNA and histone modifications. It was even believed that, for example, mRNA only played the role in information transmission. However, with the rapid development of high-throughput sequencing technology, it is found that mRNA also undergo various modifications such as N6-methyladenosine (m^6^A), N1-methyladenosine (m^1^A) and pseudouridine methylation during the process of exon splicing, 5\'-capping and 3\'-tailing [@B4]-[@B6]. These modifications will then affect mRNA splicing, nucleation, stabilization, translation and other mRNA metabolism processes, thereby regulating gene expression. Up to now, 171 RNA modifications have been discovered [@B7]. Researches on m^6^A epigenetic modification have been increasingly attended. M^6^A is a methylation modification that occurs on RNA adenine (A), and is one of the most abundant modifications in most eukaryotic mRNAs and long-chain non-coding RNAs (lncRNAs) [@B8]-[@B10]. In transfer RNA (tRNA), ribosome RNA (rRNA), microRNA, the m^6^A methylation is likewise detected [@B11]. Similar to DNA and protein modification, the m^6^A methylation is dynamically and reversibly regulated by methyl-transferases (\"writers\"), binding proteins (\"readers\"), and demethylases (\"erasers\") [@B12]. The RNAs undergo the process of methyl group modification under the action of \"writers\" that such as METTL3, METTL14, WTAP, RBM15, KIAA1429, and ZC3H13 [@B13]-[@B19]. Then the \"readers\" including YTHDF1, YTHDF2, YTHDC1, YTHDC2, HNRNPC recognize those m^6^A-modified RNAs and function in RNA processing, nuclear export, translation and decay [@B12],[@B17],[@B20]. Relying on the \"Erasers\" (FTO, ALKBH5), the m^6^A is restored to adenosine and thus to achieve demethylation modification [@B21]. Once the component involved in the regulation of m^6^A modification has been lost, the physiological functions such as cell differentiation and embryo development would be affected and expression of genes would be abnormally regulated [@B4],[@B22].

At present, some diseases, such as obesity, diabetes, neuronal disease, infertility, autoimmune diseases, have been reported associated with m^6^A modification [@B23]-[@B28]. In addition, m^6^A methylation is likewise closely related to the tumor development. M^6^A-related protein is an important regulator of tumorigenesis, high or low expression level often directly determines the tumor pathological process. For example, in cervical squamous cell carcinoma (CSCC), FTO is significantly high expressed [@B29]. It would lower the m^6^A methyl level and in turn activate the beta-catenin pathway and affect ERCC1 expression and then lead to poor prognosis. In colorectal cancer, YTHDF1 is highly expressed associated with the tumor diameter (P=0.009), lymph node metastasis (P=0.044), distant metastasis (P=0.036) and clinical stage (P=0.0226) [@B30]. As for lung cancer, report has indicated that FTO enhances the expression of myeloid zinc finger 1 (MZF1) by reducing m^6^A levels, thereby promoting the development of lung squamous cell carcinoma (LUSC) [@B31]. Another evidence points that miR-33a inhibits the proliferation of non-small cell lung cancer (NSCLC) cells by targeting METTL3 mRNA 3\'UTR binding sites [@B32]. Thus, the fully understanding of the pivotal m^6^A RNA methylation regulators is vital for lung cancer treatment. In all categories of lung cancers, lung adenocarcinoma (LUAD) is one of the most common histological types [@B33]. However, it lacks a comprehensive analysis of the expression of m^6^A RNA methylation regulators in lung cancer, especially in LUAD. The clinicopathological characteristics, diagnosis and prognostic value of such regulators remains to be explored.

In this study, the expression of 13 widely reported m^6^A RNA regulators in LUAD and normal samples were systematically analyzed. The different clinicopathological characters which related to each m^6^A modification regulator were then provided. It was found that the expression of m^6^A RNA methylation regulators played a responsible role in the progression of LUAD, which even were identified as the effective diagnostic and prognostic factors.

Methods
=======

Datasets
--------

Genotype-Tissue Expression (GTEx) dataset for normal lung tissues was downloaded from <https://www.gtexportal.org/> [@B34]. Data of lung adenocarcinoma and squamous cell lung carcinoma of TCGA [@B35] was downloaded from UCSC Xena (<https://xena.ucsc.edu/>) [@B36]. For GTEx and TCGA dataset, RNA-sequencing data (FPKM values) were normalized into log~2~ (FPKM+1). Clinical data of TCGA was downloaded from cBioPortal (<http://www.cbioportal.org/>) [@B37],[@B38]. Data of GSE75037 [@B39], GSE63459 [@B40], GSE29013 [@B41], GSE30219 [@B42], GSE37745 [@B43], and GSE50081 [@B44] were downloaded from the Gene Expression Omnibus (<https://www.ncbi.nlm.nih.gov/geo/>).

After normalization by \"LiMMA-normalizeBetweenArrays\" (<http://www.bioconductor.org/>), GSE75037 and GSE63459 were served as validation cohorts of diagnostic model. GSE29013, GSE30219, GSE37745, and GSE50081 were merged by Batch normalization to serve as validation cohorts of prognostic risk model. The details of expression data were shown in Table [1](#T1){ref-type="table"}. GTEx and TCGA dataset were acted as training cohort, and GSE75037, GSE63459, GSE29013, GSE30219, GSE37745, and GSE50081 as validation cohort. For GTEx and TCGA dataset, we defined expression value as log~2~ (FPKM+1), while for gene chips, expression value was defined as log~2~(expression).

Selection of investigative genes
--------------------------------

There were 13 m^6^A RNA methylation genes brought into the study, namely, METTL3, METTL14, WTAP, KIAA1429, RBM15, ZC3H13, YTHDC1, YTHDC2, YTHDF1, YTHDF2, HNRNPC, FTO, and ALKBH5 [@B45],[@B46]. We investigated their role in diagnosis, progression, and prognosis of lung adenocarcinoma of both training and validation cohort.

Diagnostic role of 13 m^6^A RNA methylation genes in lung adenocarcinoma
------------------------------------------------------------------------

After normalization by \"LiMMA-normalizeBetweenArrays\", 288 samples of normal lung tissues from GTEx were added into TCGA dataset of lung cancer to increase the number of the normal group. Then wilcox test was applied to differential analysis of 13 m^6^A RNA methylation genes between normal and tumorous samples. Package \"corrplot\" was used to analyze the correlation among each gene. App ClueGO (<http://www.ici.upmc.fr/cluego/version2.3.3>) of Cytoscape (<http://cytoscape.org>, version3.5.1) was to illuminate the relationship between m^6^A RNA methylation genes and related pathways [@B47]-[@B48]. Then the least absolute shrinkage and selection operator (lasso) regression model was made to diagnose lung adenocarcinoma in the training cohort by the "lars" package [@B49]-[@B50]. Eventually, 11 m^6^A RNA methylation related genes and coefficients were identified with smallest 10-fold cross validated mean square error in the training cohort.

Diagnostic score =

i\*Xi

(Coefi is the coefficient of each selected gene, Xi is the expression value.)

The diagnostic score model was further validated in the independent cohorts-GSE75037 and GSE63459. Besides, the distributions of diagnostic score were shown in different tumor stage in training and validation cohort to demonstrate the early-cancer discriminability of the diagnostic score model.

Construction of the clusters stratified by m^6^A RNA methylation genes and correlation analysis
-----------------------------------------------------------------------------------------------

We used \"ConsensusClusterPlus\" to cluster tumor samples into two to five groups [@B51]. Then Kaplan-Meier analysis and correlation analysis were applied to explain the correlation between clusters.

Prognostic evaluation of 13 m^6^A RNA methylation genes in lung adenocarcinoma
------------------------------------------------------------------------------

To assess the prognostic evaluation of each m^6^A RNA methylation gene, univariate Cox regression analysis was performed. Besides, lasso regression model was also made with 13 m6A RNA methylation genes to optimize the prognostic meaning.

Risk score =

i\*Xi

(Coefi is the coefficient of each selected gene, Xi is the expression value.)

From this risk score model, each patient could get a risk score. We defined patient into high risk group with the risk score≥ median value of all patients, and patient into low risk group with the risk score\< median value. Then Kaplan-Meier analysis was performed to display the prognostic performance of Risk score model in both training and validation cohorts. And correlation analysis was applied to explain the correlation between subgroups stratified by the risk score model and clinicopathological characters, including smoking history, TNM stage, metastasis, lymph node stage, gender, diagnosed age, and status in both training and validation cohorts.

Univariate and multivariate Cox regression analyses of clinicopathological characters and risk score
----------------------------------------------------------------------------------------------------

Univariate and multivariate Cox regression analyses of risk score and clinicopathological characters, including smoking history, TNM stage, metastasis, lymph node stage, gender, diagnosed age, and tumor stage, were performed to identify the prognostic performance of these characters. \"Rms\" package of R was used to build a nomogram with 1-year, 2-year, and 3-year overall survival (OS) as endpoints.

Statistical analysis
--------------------

All figures and data were accomplished by R (Version 3.5.3) and Graphad Prism 5.0 (La Jolla, CA). To estimate the diagnostic role of diagnostic model, operating characteristic (ROC) curves were plotted. One-way ANOVA or t-test was conducted to compare different clinicopathological characters in subclusters or risk groups. Kaplan-Meier plots, displayed as hazard ratio (HR) with 95% confidence intervals (CI), were performed to compare the OS of patients of different subclusters and risk groups. All results with P\<0.05 were considered significant.

Results
=======

Differential expression of m^6^A RNA methylation genes in lung adenocarcinoma
-----------------------------------------------------------------------------

Comparing 13 m^6^A RNA methylation genes between 347 normal samples and 526 tumorous samples of TCGA and GTEx, we found that 12 of these genes had statistic difference. Among these genes, KIAA1429, HNRNPC, YTHDC2, METTL3, WTAP, YTHDC1, and FTO were down expressed in lung adenocarcinoma, RBM15, ZC3H13, YTHDF1, YTHDF2, and ALKBH5 were up expressed (Figure [1](#F1){ref-type="fig"}A-B). These 13 m^6^A RNA methylation genes were mainly associated with RNA modification, dosage compensation by inactivation of X chromosome, and RNA destabilization (Figure [1](#F1){ref-type="fig"}C). In spearman correlation analysis, YTHDC1, YTHDC2, FTO, WTAP, HNRNPC, and METTL3 were positively related with each other. Besides, ZC3H13, ALKBH5, YTHDF1, and YTHDF2 were also positively related with each other in lung adenocarcinoma. Furthermore, YTHDF1, and YTHDF2 were negatively related with FTO, WTAP, HNRNPC, and METTL3 (Figure [1](#F1){ref-type="fig"}D).

Construction of diagnostic model
--------------------------------

The ROC curves of 13 m^6^A RNA methylation genes were shown in [Figure S1](#SM0){ref-type="supplementary-material"}A. To improve the diagnostic accuracy, 13 genes were served as candidate genes into lasso regression analysis. Figure [2](#F2){ref-type="fig"}A-B showed the process of model building in the training cohort. Diagnostic score = 0.033\*KIAA1429+0.116\*HNRNPC+0.115\*RBM15-0.067\*METTL3-0.048\*ZC3H13-0.221\*WTAP+0.213\*YTHDF1-0.132\*YTHDC1-0.135\*FTO+0.078\*YTHDF2+0.014\*ALKBH5. Each people could get a diagnostic score according to the diagnostic score model. The diagnostic scores of the tumor group were significantly higher than the normal group (Figure [2](#F2){ref-type="fig"}C). Area under the curve (AUC) was 0.996, P\<0.0001, indicating high sensitivity and specificity (Figure [2](#F2){ref-type="fig"}D). The diagnostic score model was further validated in two independent cohorts-GSE75037 and GSE63459. Figure [2](#F2){ref-type="fig"}E and [Figure S1](#SM0){ref-type="supplementary-material"}B showed that the variation trends of HNRNPC, RBM15, YTHDF1, YTHDC1, FTO, and YTHDF2 of GSE75037 were accord with these in the training cohort. Meanwhile, people of tumor group also had extraordinarily higher diagnostic scores than these of normal group (Figure [2](#F2){ref-type="fig"}F). The AUC of diagnostic score model was 0.971, which was higher than single gene (Figure [2](#F2){ref-type="fig"}G and [Figure S1](#SM0){ref-type="supplementary-material"}C). For GSE63459, the trends of RBM15, WTAP, YTHDF1, FTO, and YTHDF2 were consistent with these in the training cohort (Figure [2](#F2){ref-type="fig"}H and [Figure S1](#SM0){ref-type="supplementary-material"}D). Likewise, the diagnostic score model had its superiority in diagnosing lung adenocarcinoma (Figure [2](#F2){ref-type="fig"}I-J and [Figure S1](#SM0){ref-type="supplementary-material"}E). Besides, diagnostic score model had discrimination of each stage of lung adenocarcinoma, even the very early stage-stage IA in the training and validation cohorts (Figure [2](#F2){ref-type="fig"}K-L).

Correlation of m^6^A RNA methylation genes and clinicopathological characters
-----------------------------------------------------------------------------

In both training and validation cohorts, YTHDC2 was associated with tumor stage (P\<0.01, Figure [3](#F3){ref-type="fig"}A-C). Besides, HNRNPC was up expressed in progressed tumor group of training and merged validation cohorts (P\<0.05, Figure [3](#F3){ref-type="fig"}D-F). For mutation, there was no correlation between m^6^A RNA methylation genes and KRAS or EGFR mutation ([Figure S2](#SM0){ref-type="supplementary-material"}). For TP53 mutation, WTAP, RBM15, KIAA1429, YTHDF1, and YTHDF2 were all up expressed in mutation group (Figure [3](#F3){ref-type="fig"}G-H).

Construction of risk score model
--------------------------------

Before we build the risk score model, we clustered the patients into two to five subclusters (K=2-5) stratified by 13 m^6^A RNA methylation genes. Regretfully, subclusters stratified by m^6^A RNA methylation genes were not related with OS ([Figure S3](#SM0){ref-type="supplementary-material"}A-H). Due to subclusters that didn\'t have a good prognostic performance, we assessed the prognostic performance by a univariate Cox regression analysis in the training cohort. In these genes, HNRNPC was hazard of OS (HR: 1.8, 95%CI (1.22-2.656), P=0.003) (Figure [4](#F4){ref-type="fig"}A).

To better predict survival of lung adenocarcinoma with m^6^A RNA methylation genes, we conducted lasso regression analysis (Figure [4](#F4){ref-type="fig"}B-C). Then a ten-gene risk score model was built. Risk score=0.169\*ALKBH5-0.159\*FTO+0.581\*HNRNPC-0.348\*YTHDF2-0.265\*YTHDF1-0.123\*YTHDC2+0.434\*RBM15+0.143\*KIAA1429-0.200\*WTAP-0.310\*METTL3 (Figure [4](#F4){ref-type="fig"}D). We defined patient into high risk group with the risk score ≥ median value of all patients, and patient into low risk group with the risk score \< median value. The training cohort (n=500) and merged validation cohort (n=350) were divided into high- or low-risk group respectively. Figure [5](#F5){ref-type="fig"}A-B showed that risk score model had its superior performance in predicting outcome of lung adenocarcinoma in both training (HR: 1.828, 95%CI (1.347, 2.481)) and validation cohort (HR: 1.812, 95%CI (1.075, 2.563)). Figure [5](#F5){ref-type="fig"}C showed the expression of 10 m^6^A RNA methylation genes involved in risk score model. There existed correlation between risk group and TNM stage (P\<0.01), lymph node stage (P\<0.05), gender (P\<0.05), status (P\<0.001). The distributions of risk score, smoking history, tumor stage, metastasis, lymph node stage, tumor stage, gender, diagnostic age, and status of training cohort were shown in Figure [5](#F5){ref-type="fig"}D-L. And the correlations between risk score and gender, tumor stage, status were verified in merged validation cohort (Figure [5](#F5){ref-type="fig"}M-P).

Univariate and multivariate Cox regression analyses of clinicopathological characters and risk score
----------------------------------------------------------------------------------------------------

To investigate whether risk score was an independent risk factor, we implemented univariate and multivariate Cox regression analyses of relevant clinicopathological characters and risk score. In the univariate Cox regression analysis, lymph node stage, TNM stage, tumor stage, and risk score were all associated with bad OS (Figure [6](#F6){ref-type="fig"}A). While in the multivariate Cox regression analysis, only the risk score was related with bad OS (HR: 2.181, 95%CI (1.594-2.984), P\<0.001) (Figure [6](#F6){ref-type="fig"}B), indicating risk score was an independent risk factor of lung adenocarcinoma. A nomogram was built to facilitate clinicians to predict outcome. Diagnostic age, gender, TNM stage, and risk score were given points according to their impacts to the outcome (Figure [6](#F6){ref-type="fig"}C). By summing up all the points, we could get a total point for each patient. Then we could be aware of 1-year, 2-year, and 3-year survival of patients respectively.

Role of the diagnostic score model and prognostic risk model in squamous cell lung carcinoma
--------------------------------------------------------------------------------------------

We investigated the expression of 13 m^6^A RNA methylation genes in squamous cell lung carcinoma. From [Figure S4](#SM0){ref-type="supplementary-material"}A-B, the variation trends of 12 m^6^A RNA methylation genes were in line with these in lung adenocarcinoma. Besides, the result of correlation analysis of squamous cell lung carcinoma was roughly consistent with that in lung adenocarcinoma ([Figure S4](#SM0){ref-type="supplementary-material"}C). The diagnostic score could diagnose squamous cell lung carcinoma as well (AUC=0.992, P\<0.0001) ([Figure S4](#SM0){ref-type="supplementary-material"}D). However, the risk score that enabled clinicians to predict outcome of lung adenocarcinoma couldn\'t be applied to squamous cell lung carcinoma ([Figure S4](#SM0){ref-type="supplementary-material"}E).

Discussion
==========

Lung cancer is the most common malignant tumor with high morbidity and mortality worldwide, of which NSCLC accounts for 85-90%, including LUSC, LUAD and large cell lung cancer (LCLC) [@B52]. The underlying cause of the low lung cancer survival rate is its difficulty in early diagnosis. More than 70% patients are diagnosed at an advanced and usually fatal stage, lacking the effective therapeutic measures [@B53]. Thus, it becomes great significance to understand the development mechanism for its well diagnosis and treatment. Among these subtypes, LUAD is responsible for approximately 40% of all NSCLC cases [@B33]. The early diagnosis and treatment can greatly improve the survival rate of patients with LUAD. Recently, emerging evidence has indicated that the development of lung cancer is both affected by genetic variation and epigenetic variation [@B54]. Epigenetic regulates gene expression from multiple levels, including DNA methylation, RNA regulation, histone modification, and chromosome remodeling [@B55]. However, up to now, in lung cancers, only FTO and METTL3 have been reported as potential targets for its diagnosis and treatment [@B31],[@B32]. In this study, except FTO and METTL3, we also found that KIAA1429, HNRNPC, YTHDC2, WTAP, YTHDC1, RBM15, ZC3H13, YTHDF1, YTHDF2, and ALKBH5 had differential expression in LUAD.

In exploring the clinical value of these differential expression regulators, it was found that the joint detection of these regulators was more significant for the diagnosis of LUAD in comparison to the detection of mostly single regulator. Such combination provided higher predictive accuracy of 99.6% (AUC) for LUAD, while the most single regulator was less than 90%. The diagnostic score model we constructed involved 11genes. From the results of [Figure S1](#SM0){ref-type="supplementary-material"}, in the training and validation cohorts, YTHDF1 had a good diagnostic performance. So far, many studies in progress with regard toYTHDF1 in colorectal carcinoma and hepatocellular carcinoma, which was usually regarded as an independent prognostic factor [@B56]-[@B58]. However, there has been no related reports on LUAD, and no has been conducted as a diagnostic factor. Though the diagnostic score model had a unique superiority in the diagnosis of LUAD, we could use YTHDF1 alone to diagnose disease in consideration of operability and simplification.

In the absence of diagnostic characteristics, YTHDC2 was then found to be associated with tumor stage. Accurate pathological staging is the basis for setting on ideal treatment. Increasing researchers are also devoted in the research of molecular staging using gene or protein expression. Previously, in colon cancer, it has been found a significantly positive correlation between YTHDC2 expression and the tumor stage [@B59]. However, there is still lack of information for LUAD. The finding thus will provide new insight into our understanding of the function of YTHDC2 in LUAD. HNRNPC, except as a joint regulator in the diagnostic score model, it even exhibited differential expression in progressed tumor group and disease free group. Many studies have suggested that several types of solid tumor cells including breast cancer, gastric cancer, and esophageal squamous cell carcinoma acquired generation and development characteristics through HNRNPC disorder [@B60]-[@B62]. In lung epithelial cells, the stability of urokinase receptor (uPAR) mRNA was regulated by HNRNPC [@B63]. Increased uPAR expression as well as stabilization of uPAR mRNA would contribute to the pathogenesis of lung inflammation and neoplasia [@B63]. In this study, we preliminary confirmed that HNRNPC played a critical role in LUAD progression, which was expected to provide a potential therapeutic target.

The consensus clustering stratified by 13 m^6^A RNA methylation genes was not associated with the most clinicopathologic characters and survival. The inconsistent aberrated trends of these m^6^A RNA methylation genes might account for the frustrated utilization of the stratification. Then, we performed univariate Cox regression analysis of each selected m^6^A RNA methylation genes. We found that HNRNPC also was hazard of OS (HR: 1.8, 95%CI (1.22-2.656), P=0.003). As we have demonstrated, for breast cancer cells, the repression of HNRNPC could inhibit cell proliferation and tumor growth [@B60]. In gastric cancer, HNRNPC has been identified as a prognostic and therapeutic marker [@B61]. However, to the best of our knowledge, there has been rarely reported that HNRNPC was involved in the prognosis of LUAD.

As reported, a number of monogenic prognostic markers that were found significantly associated with prognosis, such as p53, HER2, and ERBB3 [@B64]. However, due to the heterogeneity of NSCLC, the genes involved might vary greatly among those different individuals [@B64]. Thus, it greatly prompted the researchers to construct gene expression profiles composed of multiple prognostic genes for using in the patients risk stratification. For instance, Kratz et al. ever have constructed a prognostic model using 14 genes which divided the patients with NSCLC into low-risk, intermediate-risk, and high-risk groups [@B65]. The results showed that the 5-year survival rate among three groups was significantly different, and the model exhibited satisfactory predictability in prognosis than traditional clinical pathology staging in multiple validation cohorts. Similarly, in this study, A ten-gene risk score model (Risk score=0.169\*ALKBH5-0.159\*FTO+0.581\*HNRNPC-0.348\*YTHDF2-0.265\*YTHDF1-0.123\*YTHDC2+0.434\*RBM15+0.143\*KIAA1429-0.200\*WTAP-0.310\*METTL3)was built. It stratified the OS of patients with LUAD into high and low risk categories, which exhibited entirely differences in survival. Following the discovery of these prognostic genes, the risk group distinguished through risk score was found existed correlation with TNM stage, lymph node stage, gender, status, and tumor stage in both training and validation cohort. Besides, these clinicopathological features were associated with bad OS, and risk score was an independent risk factor of LUAD. At present, TNM stage is commonly performed to assess the prognosis of patients. Due to the limited risk factors included in the TNM staging system, it is not possible to conduct a precise prediction for NSCLC. Validating and combining more clinical pathology and molecular biology risk factors to establish an effective prognosis model becomes an important direction for future research in the development of multidisciplinary comprehensive treatment. In this study, based on characteristics of diagnostic age, gender, TNM stage, and risk score, a nomogram was built to facilitate clinicians to predict outcome. To some extent, it was expected to address the issue of prognostic heterogeneity caused by single factor or insufficient risk factors analysis.

In spite of these indelible contributions, as the study went on, it was found that the above regulators exhibited no significant differential expression in KRAS and EGFR mutant LUAD. In TP53 mutant LUAD, the regulators including WTAP, RBM15, KIAA1429, YTHDF1, and YTHDF2 all were highly expressed. Thus, we inferred that those regulators were likely only effective in some types of LUAD, while the specific reasons required following elucidation. Besides, it was note-worthy that YTHDF2 was merely reported to be involved with TP53 in peripheral T-cell lymphoma and gastric cancer yet [@B66],[@B67]. YTHDF1 was functionally interactional with TP53 in gastric cancer [@B67]. To the best of our knowledge, it was not clear yet that TP53 mutant impacted the expression of WTAP, RBM15, and KIAA1429, in which more evidences were required to elucidate their mechanistic correlation.

In addition to LUAD, our work also indicated that m^6^A shared connections with LUSC, and the joint detection also provided higher predictive accuracy of 99.2% (AUC). Unfortunately, as we observed, these m^6^A RNA methylation regulators exhibited no prognostic value in LUSC. LUAD is more likely to occur in women and non-smokers [@B68]. LUSC is more common in older men and is closely related to smoking [@B69]. It is likely that different pathological features determine the different prognostic effects exhibited by the same regulators, which is worth investigating in future studies.

In this study, due to the limited sample size of normal samples, we normalized the lung tissues of GTEx and increased the credibility and accuracy of differential expression analysis. The robust analyses of correlation were performed among m^6^A RNA methylation regulators and diagnosis, pathological features or prognosis. The sample number of the study approximately reached 1900, which we thought to be enough to offset some confounding factors.

Conclusion
==========

In conclusion, our work systemically elucidated m^6^A\'s progressive role in LUAD, and provided insights into its diagnostic and prognostic function. However, m^6^A\'s functional details and its relationship with specific gene mutations and other types of NSCLC in controlling tumor progression merit further investigation.

Supplementary Material {#SM0}
======================

###### 

Supplementary figures and tables.

###### 

Click here for additional data file.

We are very grateful to the TCGA project team. Thanks for Chengyu Wu\'s help in the visualization of Graphical abstract.

Funding
=======

This work was supported by Wenzhou Science & Technology Bureau Project (grant number Y20170322).

m^6^A

:   N6-methyladenosine

LUAD

:   lung adenocarcinoma

ROC

:   receiver operating characteristic

AUC

:   area under the curve

lasso

:   least absolute shrinkage and selection operator

m^1^A

:   N1-methyladenosine

lncRNAs

:   long-chain non-coding RNAs

tRNA

:   transfer RNA

rRNA

:   ribosome RNA

CSCC

:   cervical squamous cell carcinoma

MZF1

:   myeloid zinc finger 1

LUSC

:   lung squamous cell carcinoma

NSCLC

:   non-small cell lung cancer

GTEx

:   Genotype-Tissue Expression

OS

:   overall survival

HR

:   hazard ratio

CI

:   confidence intervals

LCLC

:   large cell lung cancer

uPAR

:   urokinase receptor

![**Expression of m^6^A RNA methylation genes in lung adenocarcinoma.** (**A-B**) The differential expression of 13 m^6^A RNA methylation genes between 347 normal samples and 526 tumorous samples in lung adenocarcinoma. (**C**) Visualization of relation between m^6^A RNA methylation genes and related pathways. (**D**) Correlation analysis of 13 m^6^A RNA methylation genes. \*\*\*, P\<0.001.](ijbsv16p1785g001){#F1}

![**Diagnostic role of m^6^A RNA methylation genes related model.** (**A-B**) The solution paths of lasso regression model and the relationship between cross-validated mean square error (CV MSE) and model size. (**C**) The difference in the distribution of diagnostic score of training cohort. (**D**) The ROC curve of diagnostic score of training cohort. (**E**) The expression of 13 m6A RNA methylation genes in lung adenocarcinoma of validation cohort-GSE75037. (**F-G**) The distribution and ROC curve of diagnostic score of validation cohort-GSE75037. (**H**) The expression of 13 m6A RNA methylation genes in lung adenocarcinoma of validation cohort-GSE63459. (**I-J**) The expression of 13 m6A RNA methylation genes in lung adenocarcinoma of validation cohort-GSE63459. (**K**) The difference in the distribution of diagnostic score of tumor stage of training cohort. (**L**) The difference in the distribution of diagnostic score of tumor stage of validation cohort-GSE75037. \*, P\<0.05; \*\*, P\<0.01; \*\*\*, P\<0.001.](ijbsv16p1785g002){#F2}

![**Role of m^6^A RNA methylation genes in lung adenocarcinoma with different clinicopathological characters.** (**A-B**) The expression of m^6^A RNA methylation genes with different tumor stage of training cohort. (**C**) The expression of m^6^A RNA methylation genes with different tumor stage of merged validation cohort. (**D-E**) The expression of m^6^A RNA methylation genes with different tumor status of training cohort. (**F**) The expression of m^6^A RNA methylation genes with different tumor status of merged validation cohort. (G-H).The expression of m^6^A RNA methylation genes with TP53 mutation of training cohort.\*, P\<0.05; \*\*, P\<0.01; \*\*\*, P\<0.001.](ijbsv16p1785g003){#F3}

![**Construction of risk score model with m^6^A RNA methylation genes.** (**A**) Univariate Cox regression analysis of 13 m^6^A RNA methylation genes. (**B-C**)The process of building lasso model with size and coefficients by multivariate Cox regression. (**D**) The coefficients of 10 m^6^A RNA methylation genes involve in lasso risk model. Hazard ratio (95% confidence intervals).](ijbsv16p1785g004){#F4}

![**Relationship between clinicopathological characters and risk score model in lung adenocarcinoma of training cohort and merged validation cohort.** (**A**) Kaplan-Meier plot of high risk and low risk subgroup in training cohort. (**B**) Kaplan-Meier plot of high risk and low risk subgroup in merged validation cohort. (**C**).The heatmap of 10 m^6^A RNA methylation genes involve in lasso risk model and relationship between clinicopathological characters and risk subgroup. (**D-L**) The distribution of risk score in different risk (D), smoking history (E), stage (F), metastasis (G), lymph node (H), tumor stage (I), gender (J), age (K), status (L) subgroup of training cohort. (**M-P**) The validation of distribution of risk score in different risk (M), gender (N), tumor stage (O), status (P) subgroup. \*, P\<0.05; \*\*, P\<0.01; \*\*\*, P\<0.001, ns, no sense; HR, hazard ratio; 95%CI, 95% confidence intervals.](ijbsv16p1785g005){#F5}

![**Univariate and multivariate Cox regression analysis of clinicopathological characters and risk score.** (**A**) Univariate Cox regression analysis for evaluating the prognostic role of clinicopathological characters and risk score in lung adenocarcinoma. (**B**) Multivariate Cox regression analysis for evaluating the prognostic role of clinicopathological characters and risk score in lung adenocarcinoma. (**C**) The nomogram plot to predict 1-year, 2-year, 3-year overall survival. Summing up each point of clinicopathological characters and risk score to predict overall survival. Hazard ratio (95% confidence intervals).](ijbsv16p1785g006){#F6}

###### 

The basic information of series in the study

  --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
  Series accession numbers   Platform used                                   No. of normal samples   No. of tumorous samples   AJCC stage   Gender                Mean age, \[min, max\]   Region   Survival Outcome
  -------------------------- ----------------------------------------------- ----------------------- ------------------------- ------------ --------------------- ------------------------ -------- ------------------
  GSE75037                   Illumina Human WG-6 v3.0 expression beadchip    83                      83                        I:50;\       Female:118 Male:48    68,\[39,90\]             USA      NA
                                                                                                                               II:20;\                                                              
                                                                                                                               III:11;\                                                             
                                                                                                                               IV:2                                                                 

  GSE63459                   Illumina Human Ref-8 v3.0 expression beadchip   32                      33                        I:28;\       Female:34\            66,\[47,88\]             USA      NA
                                                                                                                               II:5         Male:31                                                 

  GSE29013                   Affymetrix Human Genome U133 Plus 2.0 Array     0                       30                        I:16;\       Female:10 Male:20     64,\[32,76\]             USA      OS
                                                                                                                               II:6;\                                                               
                                                                                                                               III:8                                                                

  GSE30219                   Affymetrix Human Genome U133 Plus 2.0 Array     0                       85                        NA           Female:19\            61,\[44,84\]             France   OS
                                                                                                                                            Male:66                                                 

  GSE37745                   Affymetrix Human Genome U133 Plus 2.0 Array     0                       106                       I:70;\       Female:60 Male:46     63,\[40,83\]             Sweden   OS
                                                                                                                               II:19;\                                                              
                                                                                                                               III:13;\                                                             
                                                                                                                               IV:4                                                                 

  GSE50081                   Affymetrix Human Genome U133 Plus 2.0 Array     0                       129                       I:93;\       Female:62\            69,\[40,86\]             Canada   OS
                                                                                                                               II:36        Male:67                                                 

  TCGA-LUAD                  Illumina RNAseq                                 59                      526                       I:172\       Female:274\           65,\[33,88\]             NA       OS
                                                                                                                               II:273\      Male:237\                                               
                                                                                                                               III:45\      NA:15                                                   
                                                                                                                               IV:18\                                                               
                                                                                                                               NA:18                                                                

  TCGA-LUSC                  Illumina RNAseq                                 49                      501                       I:114\       Female:131 Male:370   67,\[39,90\]             NA       OS
                                                                                                                               II:295\                                                              
                                                                                                                               III:71\                                                              
                                                                                                                               IV:21                                                                

  GTEX                       NA                                              288                     0                         NA           NA                    NA                       NA       NA
  --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

[^1]: \*These authors contributed equally to this work.

[^2]: ✉ Corresponding author: Yi Jin. Department of Rheumatology and Immunology, The Second Affiliated Hospital and Yuying Children\'s Hospital of Wenzhou Medical University, Wenzhou 325000, Zhejiang Province, China. Tel: +86-0577-88002611. E-mail: jinyi1990\@wzhealth.com.

[^3]: Competing Interests: The authors have declared that no competing interest exists.
