Prognostic Analysis of Risk Model and m6A-related lncRNAs
Figure 1 shows a coexpression network map between m6A and lncRNAs in HNSCC. Blue represents lncRNA, and red represents coexpressed m6A. We developed a signature for predicting the prognosis of HNSCC by performing Lasso Cox regression analysis on 21 m6A-related lncRNAs based on the TCGA database. Next, we obtained twelve genes (SAP30L-AS1, AC022098.1, LINC01475, AC090587.2, AC008115.3, AC015911.3, AL122035.2, AC010226.1, AL513190.1, ZNF32-AS1, AL035587.1 and AL031716.1) to build the risk model, and the risk score was calculated by the coefficients of these genes(Fig. 2A). Based on the median cutoff value of the risk score, patients with HNSCC were separated into low-risk or high-risk groups. As shown in Figs. 2B, C, Kaplan–Meier survival curves indicated that high-risk patients had significantly worse OS than low-risk patients as the risk score increased(P < 0.01). The heat map (Fig. 2E) shows that the expression of all twelve genes decreased with increasing risk score. Risk score and survival status distributions are plotted in Fig. 2D. The ROC curve showed that the risk score had strong predictive ability, with an AUC of 0.715 in the training set (Fig. 2F). Our results demonstrated that the risk model may serve as an important indicator for evaluating the prognosis of HNSCC patients.
Stratification Analysis of the m6A-related lncRNAs
We attempted to explore the association between clinicopathological features and the risk score. The results revealed that HNSCC patients with cluster 1, low immune scores, and T3-4 stage disease (Figs. 3B–D) had higher risk scores, while the risk score was not associated with age or clinical stage (Fig. 3A, 3E). We performed a stratification analysis to confirm whether the risk score retains its ability to predict OS in various subgroups. Our results showed that higher-risk HNSCC patients had worse OS in patients aged ≤ 65 or > 65 years, patients with G1-2 or G3-4 and patients with M0 (Fig. 3F-I, N). Moreover, the m6A-related lncRNA risk signature could distinguish the prognostic difference between HNSCC patients with N0 disease and patients with stage I-II disease (Fig. 3J-M).
Effects of the risk score and clinicopathological variables on the prognosis of HNSCC patients
The expression of the twelve selected m6A-related lncRNAs and clinicopathological variables in the high- and low-risk groups were shown in the heat map. We found Significant differences between the two groups in terms of immune score (P < 0.05) and cluster (P < 0.05) (Fig. 4A). Age, risk score and stage were found linked with OS significantly by univariate analyses (P < 0.05, Fig. 4B). The same results were observed by Multivariate analyses(P < 0.05, Fig. 4C). These results suggest that the risk signature is a risk factor for HNSCC patients and can be independent prognosis biomarkers for HNSCC patients.
Identification of the relationship between m6A-related lncRNAs and 22 immune cells
Comprehensive analysis of the results of correlation analyses (Fig. 5A-O) showed that eosinophils, activated dendritic cells, activated mast cells, M0 macrophages, M2 macrophages, resting NK cells and resting CD4 memory T cells specifically had positive correlations with the gene signature, and naive B cells, plasma cells, CD8 T cells, activated CD4 memory T cells, regulatory T cells (Tregs), follicular helper T cells, gamma delta T cells and resting mast cells showed negative correlations with the gene signature.
Gene Set Enrichment Analysis
In view of the correlation between cluster and HNSCC prognosis, we conducted GSEA between cluster 1 and cluster 2. GSEA revealed that the majority of the novel m6A-related lncRNA prognostic signatures were positively correlated in the cluster 2 group, including Huntington’s disease, Parkinson’s disease, pyrimidine metabolism, Alzheimer‘s disease, oxidative phosphorylation, spliceosome, proteasome, RNA degradation, RNA polymerase and nucleotide excision repair (Fig. 6A-J).
Construction of m6A-related lncRNA clusters
21 m6A-related lncRNAs were mapped to the expression profile of HNSCC samples to perform consistent clustering by CCP tool, so that we can separate HNSCC samples into different tumor clusters by their different immune phenotypes. We set 9 as the maximum number of clusters (Fig. 7A). The most stable results were dividing the different tumor clusters into two clusters showed by CCP analysis, naming cluster 1 and cluster 2 (Fig. 7B, 7C). We used 487 HNSCC samples from TCGA dataset as a training group to investigate the prognosis of m6A-related lncRNA clusters. All patients in the training group were divided into cluster 1 and cluster 2 according to m6A-related lncRNAs. The immune score, stromal score and tumor purity (estimate score) was calculated by the ESTIMATE algorithm according to the gene expression profile data of 487 HNSCC samples. Then, we analyzed the disparities of stromal score, immune score, and tumor purity (estimate score) in different m6A-related lncRNA clusters. Cluster 1 had a lower immune score and higher stromal score than cluster 2 (Fig. 7E, P = 0.025; Fig. 7D, P = 0.0019). Similarly, cluster 2 also had a lower tumor purity with no statistical significance (Fig. 7F, P = 0.66). The cluster 1 had a shorter overall survival (OS) than cluster 2 showed by Kaplan-Meier survival curves (Fig. 7G, P = 0.005). The differentiation ratio of each of 22 tumor immune cells between cluster 1 and cluster 2 was showed by violin plot. The significance test was accessed by Wilcoxon rank-sum (Fig. 7H).
Comparison of composition and immune cell infiltration of the TME (tumor microenvironment) in different m6A-related lncRNA clusters
The immune cell composition of 487 HNSCC samples was estimated by the CIBERSORT method. Furthermore, we normalized the results to relative proportions by cell type. Differences were observed in 7 types of immune cells. In cluster 1, the relative levels of immune cell types were higher than those in cluster 2, such as M0 macrophages (P = 0.0067), activated mast cells (P = 0.029), and neutrophils (P = 0.017), unlike naive B cells (P = 0.0026), CD8 T cells (P = 0.012), follicular helper T cells (P = 0.003) and regulatory T cells (Tregs) (P = 0.00012), which were lower in cluster 1 (Fig. 8).