Introduction
Liver cancer is one of the most prevalent malignancies with increasing incidence worldwide.1,2 Hepatocellular carcinoma (HCC) accounts for the majority of liver cancer, with an overall 5-year survival rate of less than 20%. Due to the insidious feature of HCC, most patients are diagnosed at advanced stages, and thus lose the opportunity for radical treatments.3 Systemic treatments with targeted therapy dramatically alter the situation of treatment for extrahepatic metastases and vascular invasion.4 Nevertheless, the improvement of survival is limited. Recent reports have suggested that immunotherapy based on immune checkpoint inhibitors is a promising therapeutic strategy for patients with advanced HCC, while its efficiency remains less than 20%.5 Therefore, advances in the treatment of HCC are urgently needed, which requires comprehensive investigation on the pathogenesis of HCC.
N6-methyladenosine (m6A), the most abundant mRNA modification in eukaryotic cells, has been found to play important roles in various pathophysiological processes, including the occurrence, development, and metastasis of cancers, including HCC.6–8 Moreover, it is a dynamic reversible modification implicated in post-transcriptional gene regulation, referring to the RNA splicing, transporting, translation, stability, and the high-level structure.9 Notably, m6A modification is highly conserved and can be catalyzed or removed by specific enzymes, including the m6A methyltransferases (“writers”) and the demethylases (“erasers”). It is recognized by m6A-binding proteins (“readers”).10 Increasing evidence reveals that m6A is a novel regulatory mechanism regulating gene expression and cell fate decisions. The potential mechanisms are briefly summarized in Figure S1.
With the advent of next-generation sequencing, the landscape of m6A in the transcriptome has been determined, and the characteristics and patterns of m6A modification on mRNAs have also been delineated.11 Interestingly, researchers have found that the m6A modification is not only ubiquitous in mRNA but also prevalent in noncoding RNAs (ncRNAs), including microRNAs (miRNAs), long noncoding RNAs (lncRNAs), and circular RNAs (circRNAs).11–13 The m6A modification of ncRNAs potentially affects their functions. As such, the interaction between m6A modification and ncRNAs provides a novel approach to explore the pathophysiological mechanisms implicated in diseases, including cancer.14
LncRNAs are a large class of ncRNAs with a length of above 200 nucleotides. Generally, they lack protein-coding capacity and primarily regulate gene expression.15 With the identification of their multiple functions, lncRNAs have been supposed as key players in various physiological processes and diseases by regulating gene expression at both transcriptional and post-transcriptional levels.16
Recent studies have revealed that lncRNAs can also be extensively m6A-modified, and the potential regulatory mechanisms for lncRNA functions have also been identified.13,17–22 First, the occurrence of m6A modifications might provide binding sites for the m6A “readers” or trigger RNA-binding protein entry by altering the structure of local RNA. Second, m6A modifications potentially modulate the RNA-DNA triple helix structure and influence the binding between lncRNAs and corresponding DNA sites. Nonetheless, the crosstalk and mechanisms between m6A modification and lncRNAs remain to be further elaborated.
In the present study, we comprehensively investigated the role of m6A-related lncRNAs in HCC. First, the m6A-related lncRNAs were identified in HCC based on the differentially expressed m6A-related genes between tumor and normal tissues in The Cancer Genome Atlas (TCGA) database. Subsequently, the m6A-related lncRNAs with prognostic significance were screened out and used for consensus clustering, survival analysis, and pathway enrichment evaluation. A novel risk score model based on the m6A-related lncRNAs was constructed by the least absolute shrinkage and selection operator (LASSO) Cox regression. Additionally, multivariate Cox analysis was used to assess the independent prognostic significance in the training and validation groups. A prognostic nomogram was established for the overall survival (OS) of HCC patients. Finally, we evaluated the correlation between the risk model and the immune checkpoint expression, immune subtype, tumor stemness, and drug susceptibility. Our findings provide novel perspectives on the interaction between m6A and lncRNAs in HCC, and a novel prognostic indicator for HCC patients.
Methods
Datasets and patients
The level 3 RNA-seq transcriptome profiles of 374 HCC samples and 50 normal samples were obtained on March 1, 2021, from the TCGA-LIHC database (https://portal.gdc.cancer.gov/ ). The corresponding clinical and pathological characteristics, including age, sex, tumor grade, TNM stage, follow-up time, and survival status, were also collected. Duplicated sample information from the same patient was combined as averages. Finally, the data from 370 HCC patients were included. The transcriptome data of TCGA-LIHC were further distinguished with mRNA and lncRNAs, and 19,604 mRNA and 14,086 lncRNAs in HCC were obtained.
TCGA is a public database, and all cases involved in the database have been consented to use for TCGA and obtained ethical approval, which is available online to be downloaded and analyzed by individual researchers. This work was based on the open-source data and associated identification information of all enrolled cases had been blinded; hence, the study protocol was exempted from ethical approval. The study was also performed and reported in accordance with the Helsinki Declaration. The data analysis procedure is shown in Figure 1.
Identification of the m6A-related lncRNAs with prognostic significance in HCC
First, the m6A-related genes were summarized by reviewing the previous literature, from which 36 genes were obtained (summarized in Table S1). Then, the comparison analyses between tumors and normal controls were conducted using the Wilcoxon test with R package “limma” to identify the differentially expressed m6A-related genes (m6A-DEGs) in HCC. Subsequently, the lncRNAs associated with m6A-DEGs were screened through Pearson’s correlation analyses (coefficient >0.5, and p<0.001), through which 259 lncRNAs were obtained. Finally, univariate Cox regression analyses using the R package “survival” were conducted to identify m6A-related lncRNAs with significant prognostic values for the OS of HCC patients. A p-value <0.001 was considered a prognostic candidate for m6A-related lncRNA.
Consensus clustering analysis
Based on the consensus expression of prognostic m6A-related lncRNAs in tumors, the 370 HCC cases were clustered into two subgroups using the R package “ConsensusClusterPlus”. A survival analysis based on the Kaplan-Meier method was conducted using the R package “survival” to compare the differential OS between clusters. The clinical characteristics were also compared by the Chi-square test. Gene set enrichment analysis (GSEA) was conducted to evaluate the enriched pathways within different subgroups. Finally, the expression of several immune checkpoint-related genes was evaluated between the two subgroups.
Construction of the risk score model based on the prognostic m6A-related lncRNAs in HCC
Based on the results of univariate Cox analyses, the LASSO Cox regression was utilized to eliminate the highly correlated lncRNAs using the R package “glment”. Ultimately, only six m6A-related lncRNAs were identified and entered to construct the novel risk score model. With the coefficients and expression levels of lncRNAs, the risk score for each HCC patient was calculated by the following formula:
Risk score=∑i=1n Coefi∗xi(Coefi=coefficient, xi=lncRNA expression).
Then, all patients were randomly divided into two cohorts (186 cases in the training group, and 184 cases in the validation group). In each cohort, the HCC patients were further divided into low- and high-risk groups by the median of the risk score. Thereafter, the survival analyses using the R package “survival” for OS were conducted. Receiver operating characteristic (ROC) curves were also adopted to evaluate the prognostic accuracy. Besides, the survival analyses were further evaluated in different subgroups divided by the clinical factors.
To evaluate the independent prognostic significance of the risk score model in HCC, we further conducted both univariate and multivariate Cox regression analyses. Then, a nomogram for OS prediction was established. The prediction performance was assessed by the concordance index (C index) and calibration curves.
Correlations between risk score model and immune-related genes and immune subtypes in HCC
To explore the potential mechanisms of m6A-related lncRNAs in HCC, we further analyzed the relationship between the scoring model and tumor immunity. First, we obtained the immune-related genes with negative regulation functions from the Tracking Tumor Immunophenotype (TIP) database (http://biocc.hrbmu.edu.cn/TIP/index.jsp ).23 The difference of these negative immunoregulatory genes between the high- and low-risk group were compared using the Wilcoxon test.
Besides, the profile of immune subtypes for HCC in the TGGA cohort was obtained from the UCSC-Xena (https://xenabrowser.net/datapages/ , version: 2018-04-03). The Wilcoxon test was used to compare the risk score between different immune subtypes.
Statistical analysis
All statistical analyses were conducted using R (https://www.r-project.org/ , version 4.0.3). A p-value <0.05 (two-sided) was considered statistically significant.
Results
Identification of the differentially expressed m6A-related genes in HCC
By comparing the expression of 36 m6A-related genes in 374 HCC tumors and 50 normal tissues, 24 genes with significant differential expression were obtained (Fig. 2A, details in Table S2). Out of these, 22 genes were upregulated (METTL3, METTL16, WTAP, VIRMA, RBM15B, PCIF1, YTHDF1, YTHDF3, HNRNPC, FMR1, LRPPRC, HNRNPA2B1, RBMX, NXF1, TRMT112, ZCCHC4, CPSF6, CBLL1, SETD2, SRSF3, SRSF10, and PRRC2A; all p<0.05), while two genes (IGF2BP1 and IGF2BP3, all p<0.05) were down-regulated in the tumor. The expression profiles, correlations, and potential regulatory function and pathways of these m6A-related genes in HCC are summarized in Figure S2.
Identification of the m6A-related lncRNAs with prognostic significance in HCC
With annotation of the TCGA-LIHC transcriptome profile, 14,086 lncRNAs expressed in HCC were identified. After the correlation analyses of m6A-DEGs and lncRNAs, 259 m6A-related lncRNAs were screened out (Fig. 2B, details in Table S3). Further univariate Cox analyses revealed that 29 m6A-related lncRNAs were related to the OS of HCC (all p<0.001), and all significant m6A-related lncRNAs were the risk factors for the prognosis (all hazard ratio >1) (Fig. 2C). The expression profiles of each m6A-related lncRNAs in the HCC cohort are shown in the heatmap (Fig. 2D), and all lncRNAs appeared up-regulated in HCC tumors, which indicated the important roles of these m6A-related lncRNAs in HCC. Furthermore, correlation analyses revealed that nearly all m6A-related lncRNAs in HCC were significantly correlated (Fig. 2E).
Consensus clustering of HCC patients based on the prognostic m6A-related lncRNAs
To evaluate the significance of m6A-related lncRNAs in the development of HCC, consensus clustering analysis was performed to separate the HCC tumors into various clusters based on the expression profile characteristics of these prognostic m6A-related lncRNAs. The cumulative distribution function (CDF) with different clustering strategies from k=2–9 and the relative changes of the area under CDF curves are illustrated in Figure 3A–B. The sample distribution determined by the clustering strategies from k=2–9 is shown in Figure 3C. Considering the CDF increase and consistent expression of m6A-related lncRNAs in HCC, k=2 was adopted as the clustering criterion (Fig. 3D). There were a total of 310 and 60 cases in clusters 1 and 2, respectively.
Prognostic and functional differences between HCC clusters
Survival analysis revealed that HCC patients in cluster 2 had worse OS than patients in cluster 1 (p<0.001) (Fig. 3E). Furthermore, the Chi-square tests confirmed the significant difference of several clinical factors between the two clusters (Fig. 3F). The proportions of female patients (p<0.05) and advanced pathological grade (grade 3–4) (p<0.001), T stage (T3–4) (p<0.01), and clinical stage (stage III–IV) (p<0.001) were higher in cluster 2 than in cluster 1. Besides, the expression of all enrolled m6A-related lncRNAs was upregulated in cluster 2. These results suggested that the clustering strategy based on the m6A-related lncRNAs had important clinical significance.
Additionally, we evaluated the potential functional differences between the two clusters divided by the m6A-related lncRNAs. GSEA indicated that the pathways of the spliceosome, cell cycle, RNA degradation, nucleotide excision repair, and ubiquitin-mediated proteolysis were enriched in cluster 2, while the pathways of PPAR signaling, retinol metabolism, complement and coagulation cascades, drug metabolism cytochrome P450, and fatty acid metabolism were enriched in cluster 1 (Fig. 3G). The differentially enriched pathways in different clusters potentially account for the distinct expression patterns of m6A-related lncRNAs.
Based on the role of m6A modification in the regulation of the tumor microenvironment (TME), we attempted to investigate whether there were differences in the immune checkpoint genes’ expression between the two clusters. Results showed that the expression of programmed cell death 1 (PD-1), cytotoxic T lymphocyte-associated antigen-4 (CTLA-4), lymphocyte-activation gene 3 (LAG3), T-cell immunoglobulin and mucin domain 3 (TIM3), and T cell immunoreceptor with immunoglobulin and ITIM domain (TIGIT) were all higher in cluster 2 (all p<0.001), while no significance was only found with the PD-1 ligand 1 (PD-L1) (p>0.05) (Fig. 3H). This indicated that differentially expressed m6A-related lncRNAs might play a role in the regulation of the expression of immune checkpoint genes, which could help to guide the treatment of immune checkpoint inhibitors for HCC patients and provide potential for efficacy prediction.
Development of the novel prognostic risk score model based on the m6A-related lncRNAs
Based on the 29 prognostic m6A-related lncRNAs from univariate Cox analyses, we applied LASSO Cox regression analysis for further screening of independent prognostic factors for OS of HCC (Fig. 4A–B). Six critical lncRNAs (KDM4A-AS1, NRAV, AC074117.1, SNHG3, AC025176.1, and AL031985.3) were finally identified. Then, the risk score based on the six m6A-related lncRNAs for each HCC patient both in the training (n=186) and validation (n=184) group was calculated according to their expression levels and coefficients determined by LASSO Cox regression. In each group, HCC patients were further divided into high- and low-risk subgroups by the median of risk scores. Survival analyses revealed that patients with high risk had worse OS than those with low risk in the training group (p<0.001; Fig. 4C), and the prognostic difference was also noted in the validation group (p=0.012; Fig. 4D). The time-dependent ROCs were further plotted. In the training group, the area under the curve (AUC) for 1-, 3-, and 5-year OS was 0.794, 0.771, and 0.694, respectively (Fig. 4E). In the validation group, the AUC for 1-, 3-, and 5-year OS was 0.661, 0.630, and 0.692, respectively (Fig. 4F). Besides, 42% of high-risk patients died during the follow-up period, whereas only 24% of low-risk patients died in the training group (Fig. 5A). Similar results were found in the validation group. Further, 43% of patients died in the high-risk subgroup, while only 31% died in the low-risk subgroup (Fig. 5B). The risk plots of both the training and validation groups showed the risk score distribution, survival status, and expression of the six lncRNAs of each HCC patient (Fig. 5C–D). These findings indicated that the risk score model based on the m6A-related lncRNAs could effectively differentiate and accurately predict the OS of HCC patients.
Independent prognostic significance of the novel risk score model based on the m6A-related lncRNAs
Both univariate and multivariate Cox analyses were performed to identify the independent prognostic factors for the OS of HCC patients. In the training group, only the clinical stage and the novel risk score were significant in the univariate [stage, p<0.001, HR=1.772 (1.275–2.462); risk score, p<0.001, HR=1.181 (1.113–1.254)] (Fig. 6A) and multivariate Cox analyses [stage, p=0.009, HR=1.573 (1.118–2.212); risk score, p<0.001, HR=1.166 (1.092–1.246)] (Fig. 6B). Similar results were also observed in the validation group, and only the clinical stage and risk score were found with significance both in the univariate [stage, p<0.001, HR=1.624 (1.244–2.119); risk score, p<0.001, HR=1.218 (1.091–1.361)] (Fig. 6C) and multivariate Cox analyses [stage, p<0.001, HR=1.585 (1.205–2.086); risk score, p=0.004, HR=1.190 (1.056–1.341)] (Fig. 6D). Furthermore, we conducted various subgroup analyses within different clinical factors to evaluate the prognostic significance of our risk score model in all HCC patients. We found that patients in the high-risk group showed worse OS both with age ≤60 years (p<0.001; Fig. 6E) and >60 years (p<0.001; Fig. 6F), female sex (p=0.003; Fig. 6G) and male sex (p<0.001; Fig. 6H), grade 1–2 (p<0.001; Fig. 6I) and 3–4 (p=0.010; Fig. 6J), and stage I–II (p<0.001; Fig. 6K) and III–IV (p=0.009; Fig. 6L). The findings showed that when patients were divided based on the risk score, a significant difference was noted in OS between the high- and low-risk patients regardless of the clinical factor subgroups.
Besides, the distribution of different clinical factors was compared between high- and low-risk groups of all HCC patients. Results showed significant differences with pathological grade, clinical stage, and cluster between these two risk groups. There were more patients with advanced pathological grade (3–4) (p<0.001) and stage (III–IV) (p<0.05) in the high-risk group, and all patients in cluster 2 were classified into the high-risk group (p<0.001) (Fig. 7A). Higher average risk scores in grade 3–4, stage III–IV, and cluster 2 (all p<0.001) were also noted (Fig. 7B).
A novel prognostic nomogram was constructed based on the clinical stage and risk score (Fig. 7C), which could provide an easy-to-use method for predicting the 1-, 3-, and 5-year OS of HCC. The nomogram exhibited good efficacy in estimating the OS with a high C-index of 0.674. The calibration curves for the 1-, 3-, and 5-year OS rates were also largely overlapped with the standard lines (Fig. 7D).
Correlations between the risk score and expression of immune-related genes and immune subtypes of HCC
Considering the important role of the immunosuppressive microenvironment in tumors, we sought to find out whether the risk score based on m6A-related lncRNAs was correlated with the expression of immune-negative regulatory genes. The immune-negative regulatory genes were first extracted from the list of immune-related genes obtained from the TIP database, and a total of 54 genes were identified. Then, the expression of these genes was compared between the high and low-risk patients in the training and validation groups, respectively. Consequently, 20 and 26 immune-negative regulatory genes displayed significant differences in the training and validation groups, respectively (Fig. S3A–B). Only 19 differentially expressed genes were found in both groups, including PDCD1, CTLA4, HAVCR2 (TIM3), ARG1, CCL28, CXCL8, DNMT1, EZH2, HAVCR1, ICAM1, LAIR1, LGALS9, MICA, MICB, SMC3, TGFB1, TNFRSF14, VEGFA, and VTCN1. PD-1, CTLA-4, and TIM3 expression in patients with different risk scores, and the correlations between these genes and risk score are illustrated both in the training (Fig. 8A) and validation (Fig. 8B) groups. Higher expression levels of PD-1, CTLA-4, and TIM3 were found in the high-risk group and positive correlations were also identified between these immune-related genes and risk scores.
Furthermore, the distribution of risk scores was compared in different immune subtypes of HCC. The subtype of C1, C2, C3, C4, and C6 is the wound healing type, interferon-gamma dominant type, inflammatory type, lymphocyte depleted type, and transforming growth factor-beta dominant type, respectively. The results showed that the risk scores of patients belonging to C1 and C2 immune subtypes were significantly higher than those belonging to C3 and C4 immune subtypes (Fig. 8C). In the high-risk group, the proportions of C1 and C2 subtypes were higher, whereas the C3 subtype was significantly less (Fig. 8D).
Discussion
HCC is one of the leading causes of cancer-related death.1 A high recurrence rate restricts its prognosis critically, while no effective prevention methods have been reported. With the advancement of systematic therapies, increasing drugs have emerged after sorafenib, but they are still unsatisfying in improving survival.3 Immunotherapy based on immune checkpoints has also been recently used in HCC.5 However, its efficiency is limited to less than 20%. Therefore, there is a desperate and urgent need for comprehensive research on the mechanisms of development, progression, recurrence, metastasis, and resistance of HCC. This will help reveal key regulatory pathways or networks in cancer, as well as promote the development and improvement of corresponding therapies. Besides, this provides the potential to break the current treatment dilemma.
The m6A modification of RNA provides a new mechanism of gene expression regulation, which is mainly involved in the regulation of RNA splicing, stability, transporting, and translation.24,25 Dysregulation of m6A-related regulators, including different “writers”, “readers”, and “erasers”, are involved in all stages of tumor development and progression.9,10 In the present study, the expressions of 24 m6A-related genes were found dysregulated in HCC tumors, which were also evaluated by other researchers. Zhang et al.26 discovered that 11 out of 13 m6A-related genes were significantly increased in HCC. Most of the differentially expressed genes in their study were also included in our results. However, because of much more m6A-related genes enrolled by reviewing the latest literature, we obtained additional differentially expressed genes. This indicates that research on the modification of m6A is rapidly progressing and increasing discoveries are being revealed.
In addition to the modifications on mRNA, increasing evidence has confirmed a close correlation between m6A modification and ncRNAs, specifically for lncRNAs.11,20–22 Previous studies have demonstrated that m6A methylation potentially promote the lncRNA X-inactive specific transcript (XIST)-mediated transcriptional repression.27 The lncRNA of metastasis-associated lung adenocarcinoma transcript 1 (MALAT1) could be modified by METTL3 and facilitate the aggravation of renal fibrogenesis.28 Besides, the m6A methylation could also modify lncRNAs to participate in carcinogenesis and metastasis of various tumors, including head and neck cancer, esophageal cancer, lung cancer, gastric cancer, pancreatic cancer, colorectal cancer, and HCC.29 The m6A level of XIST was found to be increased by METTL14, which could inhibit the invasion of colorectal cancer.30 However, LINC00958 was found to be upregulated by METTL3 and promoted the tumor migration and invasion in HCC, and METTL16 could interact directly with MALAT1 to facilitate the tumor proliferation.19,21 Hence, these findings demonstrated significant crosstalk between lncRNAs and m6A in tumorigenesis and development, which provided novel mechanisms and potential targets for cancer regulation and treatment.
Based on the significant roles of m6A in cancer and the close interactions between m6A and lncRNAs, we investigated the potential regulatory mechanisms and prognostic values of m6A-related lncRNAs in HCC on the whole with the transcriptome data from TCGA. After rigorous selection, six critical m6A-related lncRNAs (NRAV, SNHG3, KDM4A-AS1, AC074117.1, AC025176.1, and AL031985.3) were identified and used to construct a novel risk score model, which could distinguish the prognosis of patients with high accuracy. Notably, the risk model was also identified as an independent prognostic factor for HCC. Then, a prognostic nomogram providing high accuracy for 1-, 3-, and 5-year OS prediction of HCC was established. This extended the accessibility of m6A-related lncRNAs for clinical use.
Unexpectedly, significant correlations were observed between m6A-related lncRNAs and immune-related genes and immune subtypes of HCC. Considering the essential support of the TME in tumorigenesis and progression, the malignant behavior of the tumor is also regulated by other cells and factors in the microenvironment. The interaction between cancer cells, immune cells, and the immunosuppressive microenvironment regulates the tumor progression at all stages. Therefore, the changes of tumor biological behaviors mediated by m6A modification may play a part via immune-related mechanisms.31 In the present study, the expression of immune checkpoints, including PD-1, PD-L1, CTLA-4, LAG3, TIM3, and TIGIT, were found with significant differences between the two clusters divided based on the 29 m6A-related lncRNAs. The expression levels of all the immune checkpoint genes were higher in cluster 2, which had higher expression of m6A-related lncRNAs. When the risk score model based on the six critical m6A-related lncRNAs was further constructed, the expression profiles of immune checkpoint genes with PD-1, CTLA-4, and TIM3 were highly expressed in the high-risk group. These findings indicated the potential role of m6A in the regulation of immune checkpoints expression in the TME. Furthermore, the distribution of immune subtypes in the high- and low-risk groups was significantly different. The proportions of C1 (wound healing type) and C2 (interferon-gamma dominant type) were higher in the high-risk group, whereas the C3 proportion (inflammatory type) was higher in the low-risk group of HCC. These findings provided potential guidance for the selection of immune checkpoint inhibitor-based therapy for HCC patients.
As for the six critical m6A-related lncRNAs, the roles of NRAV and SNHG3 in cancer have been investigated in some studies. Xu et al.32 discovered that NRAV was upregulated in HCC tumors and could serve as a poor prognostic factor for HCC patients, which was consistent with our findings. However, the NRAV was identified as an m6A-related lncRNA of HCC in our study, whereas Xu et al.32 reported an immune-related lncRNA. This indicated the complex roles of lncRNA in tumors. SNHG3 has been identified as a tumor-promoting lncRNA in various cancers.33–36 Wu et al.33 found that SNHG3 was upregulated in HCC and could promote cell proliferation, migration, and invasion. Zhang et al.37 found that SNHG3 could induce epithelial-mesenchymal transition and sorafenib resistance in HCC. Notably, a recent study revealed that SNHG3 could epigenetically regulate its neighboring MED18 gene methylation by binding to EZH2 and inhibiting its expression in gastric cancer.38 Besides, Zhang et al.34 found that the increased expression of SNHG3 induced by the platinum treatment could regulate the m6A level by targeting METTL3 in esophageal cancer. These findings indicated that the important role of SNHG3 in epigenetic regulation of cancer progression and resistance. Yet, the role of SNHG3 in the regulation of m6A in HCC remains unreported. As such, additional studies on these newly discovered lncRNAs are essential.
This work has several limitations worth mentioning. First, the data used and the results obtained in this study are based on the transcriptome data from TCGA, and therefore should be further verified by local experiments. Secondly, the training and validation groups were divided from the TCGA cohort, which required an external cohort for further validation. Thirdly, the transcriptome data of HCC in TCGA are mainly derived from Americans, which will cause inevitable selection bias.
In conclusion, a novel risk score model was developed based on the m6A-related lncRNAs with significant potential for HCC prognosis prediction. Our findings provide new perspectives for the potential mechanisms of m6A-related lncRNAs in regulating the immune microenvironment of HCC.
Supporting information
Fig. S1
Schematic diagram of the potential mechanisms of RNA N6-methyladenosine (m6A) methylation.
(TIF)
Fig. S2
Expression profile and function enrichment of m6A-DEGs in HCC.
(A) Heatmap of 24 differentially expressed m6A-related genes in HCC tumors and normal controls. (B) Correlation among these m6A-related genes. (C) Function enrichment of differentially expressed m6A-related genes in HCC. (D) Pathway analysis of differentially expressed m6A-related genes in HCC. HCC, hepatocellular carcinoma; m6A, N6-methyladenosine; m6A-DEGs, differentially expressed m6A-related genes.
(TIF)
Fig. S3
Expression profile of immune-negative regulatory genes in HCC patients with different risk scores.
(A) Differentially expressed immune-negative regulatory genes (n=20) in the training group. (B) Differentially expressed immune-negative regulatory genes (n=26) in the validation group. *p<0.05, **p<0.01, ***p<0.001. HCC, hepatocellular carcinoma; TCGA, The Cancer Genome Atlas.
(TIF)
Table S1
m6A-related genes.
(DOCX)
Table S2
Coefficients of m6A-related lncRNAs in the risk score model.
(DOCX)
Table S3
Clinical characteristics of HCC patients in the training and validation groups.
(DOCX)
Abbreviations
- AUC:
area under the curve
- C-index:
concordance index
- CDF:
cumulative distribution function
- circRNAs:
circular RNAs
- CTLA-4:
cytotoxic T lymphocyte-associated antigen-4
- GSEA:
gene set enrichment analysis
- HCC:
hepatocellular carcinoma
- HR:
hazard ratio
- LASSO:
least absolute shrinkage and selection operator
- LAG3:
lymphocyte-activation gene 3
- lncRNAs:
long noncoding RNAs
- m6A:
N6-methyladenosine
- m6A-DEGs:
differentially expressed m6A-related genes
- MALAT1:
metastasis-associated lung adenocarcinoma transcript 1
- miRNAs:
microRNAs
- ncRNAs:
noncoding RNAs
- OS:
overall survival
- PD-1:
programmed cell death 1
- PD-L1:
PD-1 ligand 1
- ROC:
receiver operating characteristic
- TCGA:
The Cancer Genome Atlas
- TIGIT:
T cell immunoreceptor with immunoglobulin and ITIM domain
- TIM3:
T-cell immunoglobulin and mucin domain 3
- TIP:
Tracking Tumor Immunophenotype
- TME:
tumor microenvironment
- XIST:
X-inactive specific transcript
Declarations
Acknowledgement
The results shown here are in whole or part based upon data generated by the TCGA Research Network: https://www.cancer.gov/tcga. We also thank Home for Researchers editorial team (www.home-for-researchers.com) for providing assistance in preparing the manuscript.
Data sharing statement
The datasets presented in this study can be found in online repositories. The names of the repositories can be found in the article.
Funding
The work was supported in part by grants from the Guangdong Natural Science Foundation (Nos. 2015A030313038, 2015A030312013), Science and Technology Program of Guangzhou City (No. 201607010024), and Guangdong Key Laboratory of Liver Disease Research (No. 2017B030314027).
Conflict of interest
The authors have no conflict of interests related to this publication.
Authors’ contributions
Study concept and design (TD, YY, GW), acquisition of data (JL, LY, HY), analysis and interpretation of data (TD, MD, WL, HL), drafting of the manuscript (TD, JL, LY), critical revision of the manuscript for important intellectual content (WL, HL, YY, GW), and administrative, technical, or material support, study supervision (YY, GW).