The flowchart of the present study was summarized in Fig. 1.
Identification of differentially expressed m6A regulators and lncRNA
To assess the m6A regulators that may be involved in the development of LUAD, we systematically investigated the expression patterns of m6A regulators between LUAD and normal lung tissues based on the available TCGA and GEO cohorts. As shown in Table 2, m6A related genes RBM15B, IGF2BP1, IGF2BP3, YTHDF1, HNRNPC, and HNRNPA2B1 were significantly differentially expressed and up-regulated in all four LUAD datasets (P < 0. 05). Especially, the expression of YTHDC1 was up-regulated in LUAD only in TCGA, while it was down-regulated in other datasets. So YTHDC1 was ruled out. Besides, 5907 DElncRNAs were identified. These results indicated that m6A RNA methylation regulators possessed essential biological roles in LUAD development.
Identification of mRlncRNA and prognostic lncRNA
We obtained 1035 mRlncRNAs by evaluating the correlation between 6 m6A regulators and 5907 DElncRNAs. Then 99 mRlncRNAs associated with prognosis were confirmed according to univariate Cox analysis.
Identification of the lncRNA prognostic signature in training cohort
We identified 13 lncRNAs which closely related to the OS from 5907 lncRNAs in the training cohort by LASSO analysis (Fig. 2a, Fig. 2b). Fig. 2c showed the differential expression of these 13 lncRNAs. Four of them were lowly expressed in LUAD (ADPGK-AS1, RP11-667K14.3, RP11-378A13.1, RP11-23J9.4) and nine were highly expressed in LUAD (RP11-563N4.1, CTD-3214H19.6, RP11-667K14.3, RP11- 303E16.2, MIR4435-1HG, CTD-2510F5.4, TMPO-AS1, Z83851.4, RP11-501C14.5, OGFRP1). The forest plot showed the hazard ratio (HR) with 95% confidence interval (CI) of 13 lncRNAs (Fig. 2d), whose detailed information were summarized in Table 4. Among these 13 lncRNAs, 7 were deleterious lncRNAs with LASSO coefficient > 0 (OGFRP1, RP11-501C14.5, Z83851.4, TMPO-AS1, CTD-2510F5.4, MIR4435-1HG, RP11-303E16.2), the other 6 were protective lncRNAs with LASSO coefficient < 0 (RP11-378A13.1, ADPGK-AS1, RP11-667K14.3, CTD-3214H19.6, RP11-23J9.4, RP11-563N4.1).
To assess the ability of the lncRNA signature in predicting the survival of LUAD, we established a risk score of the 13-lncRNA signature with the following formula:
Risk score = (-0.00181 × expression value of ADPGK-AS1) + (0.00047 × expression value of CTD-2510F5.4) + (-0.00349 × expression value of CTD-3214H19.6) + (0.00011 × expression value of MIR4435-1HG) + (0.00600 × expression value of OGFRP1) + (-0.00413 × expression value of RP11-23J9.4) + (0.00007 × expression value of RP11-303E16.2) + (-0.00002 × expression value of RP11-378A13.1) + (-0.00710 × expression value of RP11-563N4.1) + (-0.00256 × expression value of RP11-667K14.3) + (0.00052 × expression value of TMPO-AS1) + (0.00088 × expression value of Z83851.4) + (0.00282 × expression value of RP11-501C14.5).
Evaluating the predictive power of the 13-lncRNA prognostic signature in training cohort
The grouping was determined based on the best cut-off value of the 13-lncRNA signature, and the 356 LUAD patients of the training cohort were thus divided into high-risk group (n = 197) and low-risk group (n = 159). Fig. 3a and Fig. 3b illustrated the risk score, survival time and survival status of each patient in the training cohort. A heat map of 13-lncRNA expression was shown in Fig. 3c. Fig. 3d showed the Kaplan-Meier survival curves for the low- and high-risk groups in the training cohort. The survival time was significantly lower in the high-risk group than in the low-risk group (P < 0.001). In addition, time-dependent receptor operating characteristic (ROC) curves of the 13-lncRNA prognostic signature were shown in Fig. 3e. The area under the curve (AUC) of 1-,3-,5-year OS for the 13-lncRNA prognostic signature risk score was 0.66, 0.69, and 0.71, respectively.
Validation of the 13-lncRNA prognostic signature in testing cohort and entire TCGA LUAD cohort
To further confirm the predictive power and stability of the 13-lncRNA signature in predicting OS of LUAD patients, we validated the 13-lncRNA signature in the testing cohort (n = 90) and the entire TCGA LUAD cohort (n = 389) (Fig. 4 and Fig. 5). The risk score of 13-lncRNA signature was calculated according to the above formula. The KM survival curves of the high-risk group and low-risk group in the test cohort were shown in Figure 4D. Survival time was significantly higher in the low-risk group than that in the high-risk group (P = 0.009). Multivariate Cox regression analysis revealed that the risk score of 13-lncRNA signature (HR = 5.8,95% CI = 1.6-22) may be an independent prognostic indicator of LUAD. In addition, time-dependent ROC curves for 13-lncRNA prognostic signature were plotted. Consistent with the findings in the training cohort, 13-lncRNA prognostic signature had a good discrimination performance on the prognosis of patients with LUAD. The 1-year, 3-year, and 5-year AUC values of the prognostic signature risk scores were 0.61, 0.68, and 0.76, respectively (Fig. 4E). Similar validation results were obtained in the entire TCGA cohort of LUAD. Therefore, risk score of the 13-lncRNA signature might be considered as an independent prognostic indicator of LUAD OS.
Comparison of immune cell infiltration in high-risk and low-risk group
The modification of mRNA and lncRNA by m6A has been shown to contribute to the post-transcriptional regulation of the immune response [25]. Therefore, we explored the differences in the 22 tumor-infiltrating immune cells between the two groups.
In the high-risk group, the proportions of dendritic cells activated, macrophages M1, macrophages M2, neutrophils, T cells CD4 memory activated were increased, while the proportions of B cells memory, T cells CD4 memory resting were decreased (Fig. 6a, all P < 0. 05). These findings strongly suggested that this 13-lncRNA signature was associated with prognosis by interfering with immune cell infiltration in LUAD.
The sensitivity to immunotherapy and chemotherapy among the high-risk and low-risk group
Chemotherapy and immunotherapy are currently essential part of the treatment for advanced NSCLC [26]. The chemotherapy response for cisplatin, mitomycin and docetaxel of each LUAD patient in the TCGA cohort was calculated by using the R package of pRRophetic [28]. The estimated IC50 values showed that high risk group demonstrated had a better response to mitomycin and docetaxel, while low risk group had a better response to cisplatin (Fig. 6b, all P < 0. 05). THE expression of PD-L1 on tumor cells has been shown to be positively correlated with the efficacy of PD-1/PD-L1 in NSCLC [27].Fig. 6c showed that the expression levels of PD-L1, LAG-3, TIM-3 in the high-risk group were higher than that in the low-risk group (all P < 0.05). Above findings demonstrated that the mRlncRNA risk model had potential predictive value for chemotherapy and immunotherapy.
Pathway enrichment analysis
GSEA was further used to verify the potential molecular mechanism between high-risk and low-risk groups. The results showed that compared with the low-risk group, the high-risk group was more enriched in cell cycle, DNA replication, P53 signaling pathway, mismatch repair and nucleoside exception repair, all of which were highly related to tumorigenesis, tumor development and immunotherapy (Fig. 6d).
Identification of risk score as a key regulator to promote LUAD progression
The proximal-inflammatory (PI) and the proximal-proliferative (PP) subtypes were associated with higher risk scores, while the terminal-respiratory unit (TRU) subtype was associated with lower risk scores (Fig.7 a). The correlation between risk score and clinicalopathological characteristics of LUAD patients were shown in Fig.7 b-g, indicating that risk score was positively correlated with several categories, namely pathologic stage, tumor stage, lymph node status (all P <0.05).
Construction and validation of nomogram
To create a clinically applicable scoring tool to predict the OS of LUAD patients at 1-, 3-, and 5 years, we used risk score (based on 13-lncRNA signature), age, and pathological stage to create a nomogram (Fig. 8a). The 1-year, 3-year, and 5-year OS calibration curves showed that the nomogram's prediction of OS did not deviate significantly from the actual OS (Fig. 8b).