3.1. Acquisition and analysis of m6A-related lncRNA data of HCC patients
Figure 1 was the work chart of this study. Firstly, we identified 203 m6A-related lncRNAs based on 23 m6A regulators (Figure 2A). Univariate Cox regression analysis was conducted to identify m6A-related lncRNA with prognostic roles. And 43 m6A-related lncRNAs which were significantly related to the OS of HCC patients were identified (Figure 2B). Figure 2C and 2D were the boxplot and heatmap to show the expression level of 43 OS-related lncRNAs. The result indicted that all the expression of 43 lncRNAs were significantly upregulated in tumor tissue than normal tissue (p< 0.05).
3.2 Construction of a m6A-related lncRNAs signature
The LOSSO regression was used in the expression profiles of 43 m6A-related lncRNAs associated with OS and 8 m6A-related lncRNAs were further selected (Figure 3A). A stepwise multivariate Cox regression analysis was performed and finally, a 5 m6A-related lncRNAs risk signature was constructed (Figure 3B). The risk score was computed with the formula: risk score = (0.302121 × expression of AC145207) + (0.34025 × expression of AL031985) + (0.198432 × expression of SNHG4) + (0.051602× expression of AC099850) + (1.312568×expression of AC026412) (Table 1). Subsequently, we calculated the risk score for each patient and classified them into a high-risk group (n = 44) and a low-risk group (n = 321) according to the cut-off point of 1.91, which was obtained by the R package “survminer” (Figure 3C).
The result of PCA indicated that the clusters of the low-risk group were independent of the high-risk group without obvious intersection in m6A-related lncRNAs expression level (Figure 3D). ROC analysis indicated that the value of m6A-related signature in predicting OS (AUC= 0.754) was higher than those of tumor stage (AUC= 0.664), age (AUC= 0.535), tumor grade (AUC= 0.505) and gender (AUC= 0.513) (Figure 3E). Further analysis showed poorer prognosis and more deaths for the patients with higher risk scores (Figure 3F). A significantly shorter OS was shown in the high-risk group by Kaplan-Meier survival analysis (p < 0.001) (Figure 3G).
3.3 Relationship between m6A-related lncRNAs signature and clinicopathological variables
Figure 4A summarized the clinical characteristics of patients in different risk groups. Patients in low-risk group were older than patients in the low-risk group, with higher tumor grades and more advanced tumor stage (p< 0.05). Besides, all of 5 m6A-related lncRNAs in the risk signature were augmented in the high-risk group. The univariate Cox regression analysis revealed that the risk score was associated with a poor OS indicator (HR=1.402, 95% CI= 1.285- 1.535, p< 0.001, Figure 4B). Besides, multivariate Cox regression analysis indicated that the risk score was an independent risk predictor (HR= 1.400, 95% CI=1.270–1.542, p < 0.001, Figure 4C). In the subgroups classified by age, gender, grade, or tumor stage, patients in the low-risk group continued to have longer OS compared with those in the low-risk group (all p< 0.001, Figure 5A-5H).
3.4 M6A-related lncRNAs signature and cancer somatic variants
The distributions of HCC somatic variants between the high- and low-risk groups were represented in Figure 6A-B. The top 20 driver genes with the highest alteration frequency were shown by “maftools”. Compared with the low-risk group, the alteration frequency of TP53 (p< 0.001) was significantly higher in the high-risk group, while the alteration frequency of CTNNB1 (p= 0.022) was significantly lower. However, the TMB between high and low-risk groups did not meet the significant difference (p= 0.5, Figure 6C). Patients with low TMB had better OS outcomes than those with high TMB (p< 0.001, Figure 6D). Stratified survival analysis in subgroups showed significant differences in OS according to the risk score among the TMB subgroups (p< 0.001, Figure 6E). The result showed that the m6A-related lncRNAs signature could be a survival predictor independent of the TMB.
3.5 Estimation of the immune cell infiltration and immune checkpoints expression
We applied the CIBERSORT method to analyze the enrichment situation of 22 immune cells. The samples with p < 0.05 were used for further analysis (Figure 7A). High-risk group had lower infiltration level of naive B cells (p= 0.010) and regulatory T cells (p= 0.004, Figure 7B). Furthermore, the result of Spearman correlation analysis revealed that activated NK cell was negatively correlated with risk score (r= -0.28, p= 0.038, Figure 7C).
Next, we compared the expression of immune checkpoints expression between these two subgroups. The results indicated that immune checkpoints such as CTLA4, PDCD1, LAG3, TIGIT, and HAVCR2 were upregulated in the high-risk group (all p< 0.05, Figure 7D-7H). Although the expression of PDL1 (Figure 7I) was not significantly different between the two subgroups, Spearman correlation analysis indicated that all these immune checkpoints expression (CTLA4, PDCD1, PDL1, LAG3, TIGIT, and HAVCR2) were positively correlated with risk score (Figure 8A-8F). These results revealed a more suppressed immune status in the high-risk group.
3.6 Analysis of the therapeutic responses and underlying regulatory mechanism
Firstly, we used the TIDE score to predict the efficiency of immunotherapy. As Figure 9A showed, the TIDE score in the high-risk group was significantly lower than the low-risk group (p< 0.001). This result suggested that patients in the high-risk group had a superior response to immunotherapy.
Besides, we analyzed the association between risk subgroups and the efficacy of common chemotherapeutics and targeted therapy in treating HCC (Figure 9B-9E). These results verified that patients in the high-risk group had a lower IC50 of chemotherapeutics such as cisplatin (p = 0.028), doxorubicin (p < 0.001), mitomycin (p < 0.001), while with a higher IC50 for sorafenib (p < 0.001). Thus, the results indicated that the m6A-related lncRNA risk signature could be a potential predictor for drug sensitivity.
The result of GSEA demonstrated that the patients with the high-risk score were enriched in some cancer and proliferation pathways: cell cycle, mismatch repair, oocyte meiosis, pathway in cancer, and RNA degradation. Interestingly, patients with low-risk scores were mainly enriched in drug metabolism, such as peroxisome, drug metabolism cytochrome, and metabolism of xenobiotics by cytochrome P450 (Figure 9F). These results could explain why the low-risk group was resistant to conventional chemotherapy.