4.1 Expression profile of m6A regulatory genes in LUAD
In this study, we obtained a TCGA dataset involving 473 LUAD tissues and 41 adjacent normal tissues. The expression of 23 m6A regulatory genes in LUAD was abstracted from the TCGA database, including 9 m6A writers, 13 readers and 2 erasers. We found that there was a difference in the expression of writers METTL3, METTL14, METTL16, WTAP, ZC3H13, RBM15, RBM15B, and KIAA1429; readers YTHDC2, YTHDF1, YTHDF2, YTHDF3, HNRNPC, LRPPRC, HNRNPA2B1, IGFBP2, IGFBP3, and RBMX; and erasers FTO and ALKBH5 in LUAD and adjacent normal tissues (Fig. 1A-B). Furthermore, these differentially expressed m6A regulatory genes were associated with prognostic risk in patients with LUAD (Fig. 1C). Currently, it has been reported that the regulation of lncRNAs and m6A plays a significant role in LUAD. Thus, we extracted 2002 lncRNAs related to the expression of m6A regulatory genes from the TCGA database through code. The correlation between them is shown in Fig. 1D.
4.2 The m6A-lncRNAs related to prognosis in LUAD
Based on the aforementioned relationship between m6A regulators and lncRNAs in LUAD, lncRNAs play an important role in the tumorigenesis of lung cancer. Meanwhile, numerous lncRNAs have been demonstrated to be independent prognostic factors in lung cancer. Therefore, we continued to explore which m6A-related lncRNA was of significance in the prognosis of LUAD. The results of univariate Cox regression analysis showed that 36 of 2002 lncRNAs were significantly correlated with prognosis in LUAD (Fig. 2A). The clinical baseline information of
370 HCC cases in TCGA database is presented at (Table1) Meanwhile, a heatmap was used to show the expression of these prognosis-related m6A-lncRNAs in LUAD (Fig. 2B). Finally, through the TCGA database, we identified expression differences in 36 lncRNAs in almost all LUAD and adjacent normal tissues (Fig. 2C).
Table 1
Clinical Characteristics of the LUAD Cases in TCGA Database.
Characteristics | | Total TCGA (N = 514) |
Age (years) | ≦ 65 > 65 | 227 240 |
Gender | Male Female | 222 264 |
Stage | Stage I–II Stage III–IV | 374 104 |
T | T1-2 T3-4 | 423 60 |
M | M0 M1 MX | 333 24 125 |
N | N0 N1-3 | 312 162 |
Pharmaceutical Therapy | Yes No Unknown | 172 314 46 |
Radiation Therapy | Yes No Unknown | 100 386 46 |
Pack_years_smoked | > 60 <=60 Unknown | 56 142 154 |
4.3 Clustering study based on m6A-related lncRNAs with prognostic value
To obtain more m6A-related lncRNAs in LUAD, consensus clustering was used for cluster analysis of m6A-related lncRNAs. The results show that k = 2, which reveals that the interference between subgroups is minimal and that the distinction is significant (Fig. 3A). All of these results were based on the CDF curve of the consistency matrix (Fig. 3B-C). Together, a total of 473 patients with LUAD were classified into two clusters. The results of the heatmap show that the expression of 36 m6A-related lncRNAs in Cluster 2 was higher than that in Cluster 1. Furthermore, we compared the clinicopathological characteristics between the two subgroups and found that there was a large proportion in Cluster 2 with more stages and associated with the late M stage. However, no significant differences in N and T stage were found between the two groups. Moreover, we observed that patients who experienced treatment were more concentrated in Cluster 2, but the significance of radiotherapy and chemotherapy alone or combined with radiotherapy and chemotherapy was not reflected between the two groups. Meanwhile, age, sex and smoking history were meaningless (Fig. 3E). Most importantly, the OS of LUAD in Cluster 2 was significantly shorter than that in Cluster 1 (Fig. 3D).
4.4 Potential regulatory mechanisms of the two clusters
Many studies have shown that immune checkpoint blockade is significant in tumor immunotherapy and that m6A methylation is closely related to it, but the potential mechanism is still unclear. In this study, the immune checkpoint genes PD-L1, CTLA4, HAVCR2 and B7-H3, which are effective targets of immunotherapy, were selected. First, we determined the expression of these immune checkpoint genes in the two clusters and found that PD-L1 and B7-H3 were significant in Cluster 1 and Cluster 2 (p < 0.001) and normal and tumor tissues (p < 0.001) (Fig. 4A-H). Then, we explored the correlation between PD-L1, B7-H3 and m6A-related lncRNAs and found that most of them were positively correlated with PD-L1 and B7-H3 (P < 0.05) (Fig. 4J-K). These results implied the effect of m6A-related lncRNAs in regulating immunomodulators and provided a new strategy for LUAD immunotherapy. Furthermore, the immune function of m6A-related lncRNAs in LUAD was explored. T cells, NK cells and macrophages play an important role in the immune microenvironment. We recognize that M1 macrophages (p < 0.001), resting NK cells (p < 0.001), activated memory CD4 T cells (p < 0.001), and CD8 T cells (p < 0.05) were more highly expressed in Cluster 2 than in Cluster 1. However, activated mast cells (p < 0.05), NK cells (p < 0.001), resting mast cells (p < 0.001), regulatory T cells (p < 0.01), and plasma cells (p < 0.001) were more highly expressed in Cluster 1 than in Cluster 2 (Fig. 5A-I). These results suggest that m6A-related lncRNAs may regulate the tumor immune microenvironment in LUAD.
The biological mechanism of the heterogeneity of the two clusters was investigated by GSEA. The results showed that more tumor-related signaling pathways were enriched in Cluster 2, including the CELL CYCLE (NES: 2.5257452, FDR p value < 0.05), ERBB (NES: 1.7973282, FDR p value < 0.05), MTOR (NES: 1.7134794, FDR p value < 0.05), p53 (NES: 2.302949, FDR p value < 0.05) and WNT (NES: 1.6371493, FDR p value < 0.05) signaling pathways and the IN CANCER pathway (NES: 1.8742931, FDR p value < 0.05) (sFigure 1 A-G). All these differences between the two clusters indicated the role of m6A-related lncRNAs in the tumorigenicity of LUAD.
4.5 Construction and verification of a clinical prognosis model for LUAD based on m6A-related lncRNAs
We identified 36 differentially expressed m6A-related lncRNAs that were significantly associated with OS by univariate Cox regression analysis. Then, the 18 most significant m6A-related lncRNAs associated with prognosis were identified by LASSO regression analysis (Fig. 6A-C). Finally, through regression coefficients and expressions of the 18 m6A-related lncRNAs, a risk model for LUAD was constructed, as follows: risk score = (-0.24842) * expression (AC008763.1) +(-1.28861) * expression (AC010999.2) + 0.339972 * expression (AC008494.3) + 0.141112 * expression (AL606489.1) + (-0.27363) * expression (TSPOAP1-AS1) + 1.424991*expression (FRMD6-AS1) + (-1.45406) * expression(AL031600.2) + 0.033235* expression (LINC01876) +(-0.19802) * expression (AP002026.1)+(0.022912)*expression(LINC02587)+(-0.10636)*expression(AC087752.3)+ -0.25166)*expression(AC090617.5) + 0.094322*expression(LINC02178)+(-0.10619)*expression(AC026355.2)+(-1.20335)*expression(AC0185). Patients with LUAD were divided into high- and low-risk groups based on the median value of the risk score. The predictability of the risk model was evaluated by receiver operating characteristic (ROC) curve analysis, and our analysis showed that the areas under the ROC curve (AUCs) for OS in the training and testing cohorts were 0.779 and 0.758, respectively (Fig. 6D-E). We found that the higher the risk score was, the worse the survival status of patients in the training and testing cohorts (Fig. 6F-I). The heatmap indicates that 9 lncRNAs were highly expressed in the high-risk group in both the training and testing cohorts, and 9 lncRNAs were highly expressed in the low-risk group in both the training and testing cohorts (Fig. 6J-K). More significantly, we found that m6A-related lncRNAs expressed in high-risk groups were indeed highly expressed in LUAD compared with adjacent normal tissues, and each had a poor prognosis (sFigure 2, sFigure 3). However, m6A-related lncRNAs expressed in low-risk groups were highly expressed in adjacent normal tissues compared with LUAD tissues and had a better prognosis (sFigure 2, sFigure 3). Thus, our results suggested that the model based on the m6A-related lncRNA risk score can better predict the prognosis of patients.
4.6 The risk model of m6A-related lncRNAs as an independent prognostic factor for LUAD
We found that the risk model based on m6A-related lncRNAs could be used as an independent prognostic factor by univariate and multivariate regression analyses (Fig. 7A-D). The results of these analyses showed that OS was obviously related to stage (P < 0.001) and riskScore (P < 0.001) in the testing cohort and related to stage (P < 0.001) and riskScore (P < 0.001) in the training cohort (Fig. 7E-F). Then, through survival analysis, we revealed that high-risk patients with LUAD exhibited worse OS than low-risk patients in the testing (p < 0.001) and training cohorts (p < 0.001). Furthermore, the high- and low-risk groups showed heterogeneity in different clinical subgroups. We compared the OS of the two groups in different clinical subgroups, indicating that the survival rates of the two clusters had significant differences in the incidences of stage I + stage II (P < 0.001), stage Ⅲ + stage Ⅳ (P = 0.007), M0 (P < 0.001), M1b (P < 0.039), T1-2 (P < 0.001), T3-4 (P < 0.001), N0 (P < 0.001), age ≤ 65 (P < 0.001), age ≤ 30 (P < 0.001), and age > 60 (P = 0.015). All these results showed that the prognosis of high-risk patients was worse than that of low-risk patients.
4.7 Prognostic risk model was associated with tumorigenesis mechanism
We explored the relationship between the risk prognosis model and tumor immunity. The results showed that the expression of B7-H3 was higher in the high-risk group than in the low-risk group (Fig. 8B). However, the expression of PD-L1 was not significantly different between the two risk groups (Fig. 8A). Then, we found that the risk model was associated with immune cells, and the results showed that the higher the risk score was, the greater the number of M0 macrophages (R = 0.15, p = 0.0019), CD4 memory activated T cells (R = 0.16, p = 0.00054), M1 macrophages (R = 0.21, p = 0.0000056), and gamma T cells (R = 0.13, p = 0.0081) present. However, M2 macrophages (R=-0.094, p = 0.48), plasma cells (R=-0.11, p = 0.017), and monocytes (R=-0.15, p = 0.0023) were more common when the risk score was low (Fig. 8C-J). Furthermore, through KEGG enrichment analysis, it was found that the high-risk group was mainly concentrated in the tumor-related cell cycle and p53 signaling pathways, which are closely related to tumor immunity. (Fig. 8K-M). Collectively, the risk score established based on m6A-related lncRNAs is of great significance to immunotherapy and tumorigenesis in patients with LUAD.