Multi-omics analysis and validation of the tumor microenvironment of hepatocellular carcinoma under RNA modification patterns

1.   Introduction

Hepatocellular carcinoma (HCC), the most common type of liver cancer, is a malignancy with significant heterogeneity and strong aggressiveness ^[1]. It results from an interplay between multiple factors, such as viral infections and dietary factors ^[2,3], and commonly presents insidiously ^[4]. In this context, the 5-year survival rate of HCC patients worldwide shows a falling trend ^[5], despite the advances in treatment ^[4]. Patients suffering from HCC generally have a poor prognosis, which is a critical difficulty in current clinical treatment. Research has shown that genetic and transcriptional changes are responsible for the molecular diversity in HCC patients, while the heterogeneity resulting from this diversity leads to treatment failure to a large extent ^[6]. To improve the prognosis of HCC patients, it is necessary to identify effective and reliable signature genes that may contribute to clinical treatment.

RNA modifications are diverse in human cells, and they are important in regulating the biological functions of human bodies ^[7]. As research deepens, mRNA modifications and their fine-regulation mechanisms have received increasing attention. Previous research has established a strong correlation between RNA modifications and the occurrence and treatment of different human diseases. RNA modifications have been found to impact the development and progression of type 2 diabetes ^[8]. Utilizing RNA modifications, researchers can specifically target vascular endothelial growth factor (VEGF) to effectively treat choroidal neovascularization ^[9,10]. Additionally, dysregulation of RNA modifications has been associated with an elevated risk of developing cancer ^[11,12]. Presently, multiple common mRNA modification patterns have been identified, including m1A (N1-methyladenosine), m5C (5-methylcytosine), m6A (N6-methyladenosine) and m7G (N7-methylguanosine). It has been established that the genes modified by m1A, m5C, m6A and m7G can promote the occurrence, metastasis and progression of various tumors ^{[13,14,15,16]}. For instance, Wang et al. ^[17] proved that alteration of the m1A-modified genes tRNA methyltransferase 6 (TRMT6) and TRMT61A promotes the occurrence of HCC. Other studies reported that the m5C-modified gene Aly/REF export factor (ALYREF) exhibited increased expression in HCC tissue, which was associated with a poorer prognosis and a higher degree of malignancy in HCC patients ^[18,19]. m6A methyltransferase like 3 (METTL3) was reported with up-regulated expression in multiple cancers including HCC, and it was also found that METTL3 could promote the malignant progression of HCC by mediating the expression of miR-589-5p ^[20]. The study of Dai et al. ^[21] revealed that the overexpression of METTL1 and WD repeat domain 4 (WDR4) potentiated the progression of intrahepatic cholangiocarcinoma (ICC) through regulating m7G, and the up-regulation of METTL1 and WDR4 predicted a poor prognosis in patients.

A growing number of studies have proven the critical role of RNA modifications in occurrence, development and prognosis of HCC. Based on previous publications, the present study focused on m1A/m5C/m6A-modified genes that are closely linked with HCC and m7G-modified genes that are critical for cancer to explore the functional and immune status of HCC patients with different RNA modifications and gene patterns. Additionally, effective and reliable signature genes were identified to establish a novel prognostic model, termed RNAM-S, for HCC, and the expression of the signature genes was validated by in vitro cellular experiment. Moreover, PLOD2, which had significantly increased expression in HCC, was analyzed for its interplay with immune cells using single-cell communication analysis. Finally, the signature genes were subjected to pan-cancer analysis to confirm their application value and association with immunity. The flowchart of this study is shown in Figure 1.

Figure 1.  Process flowchart.

DownLoad: Full-Size Img PowerPoint

2.   Materials and methods

2.1.   Data source and preliminary processing

The data used in this study were mainly extracted from three sources. First, data on RNA modification-regulating genes (m1A, m5C and m6A) significantly associated with the occurrence and development of hepatocellular carcinoma (HCC) were obtained from a review published by Xu et al. ^[22], and data on m7G-related regulatory genes associated with human cancer were obtained from a study by Tomikawa ^[23]. Second, the mRNA and single-cell gene sequencing (scRNA-seq) data of HCC samples were extracted from TCGA and GEO. The mRNA datasets downloaded from TCGA (TCGA, https://portal.gdc.cancer.gov/) and GEO (GEO, https://www.ncbi.nlm.nih.gov/geo/) were designated as cohorts 1 (TCGA-LIHC, including 374 HCC tissue samples and 50 normal liver tissue samples) and 2 (GSE14520, including 247 HCC tissue samples and 241 normal liver tissue samples), respectively. Expression data (fragments per kilobase of exon model per million mapped fragments [FPKM] format) and the corresponding clinicopathological and survival data were extracted from the TCGA-LIHC cohort using the "TCGAbiolinks" package (R software [version 4.1.1]). The two datasets were preliminarily processed. Data in the FPKM format were converted to those in the TPM format. All data were standardized and filtered (samples with missing information and those with a survival time of 0 day were excluded). The scRNA-seq data in the GSE210679 dataset were filtered using the "Seurat" package. Third, a pan-cancer dataset was downloaded from the UCSC database (https://xenabrowser.net/). All data and resources used are publicly available, and ethical review was not required for this study.

2.2.   Preliminary study of changes in RNA methylation in HCC

Metascape (https://metascape.org/gp/index.html) allows the analysis of multi-integrated data. Enrichment analysis is the core function of Metascape (analysis of biological processes, pathways and protein localisation) ^[24]. The biological functions of RNA modification-regulating genes can be identified and visualized in the form of gene clusters using Metascape, which may help to understand the regulatory process of RNA modification. Bradner et al ^[25] reported that transcriptional disorder caused by gene alteration is the root cause of cancer. Genetic changes were analyzed to examine changes in RNA modification (m1A/m5C/m6A/m7G) within HCC. Somatic mutations are often found in tumor cells ^[26]. Therefore, based on the mutation data of the HCC cohort, our team studied the mutations and copy number variations (CNV) of RNA modification regulatory genes in HCC, including the frequency of CNV occurrence and the location of CNV changes on chromosomes. In addition, we used the Wilcoxon rank-sum test to identify differentially expressed RNA modification regulatory genes between normal and HCC tissues, and analyzed the relationship between genetic variations and the expression of these genes in HCC.

2.3.   Identification of different modification patterns and their relationship with the clinical characteristics and prognosis of patients with HCC

Recursive consensus clustering (RCC) was performed using an optimal number of clusters obtained via n iterations ('ConsensusClusterPlus' R package), and patients with HCC were divided into different groups according to the expression of RNA modification-regulating genes. This method is more robust than one in which a single iteration is used. Furthermore, in order to examine the clinical value of different modification patterns, this study compared the survival outcomes of HCC patients under different modification patterns. The R packages "survival" and "survminer" were used to plot Kaplan-Meier (K-M) survival curves and evaluate the differences in overall survival (OS) among different modification patterns. The R package "pheatmap" was used to generate heatmaps and explore the relationship between different modification patterns and clinical features of HCC patients, including age, gender and tumor stage.

2.4.   GSVA and TME analysis under different modification patterns

Gene set variation analysis (GSVA) was used to evaluate changes in the biological functions of genes in different modification patterns and analyze the significant biological processes of genes in each modification pattern. Tumor growth is affected by two-way communication between cells and their microenvironment ^[27]. The influence of the tumor microenvironment (TME) on tumors is dynamic. TME comprises various cells (stromal and immune cells), systems (vascular system) and tumor components. Therefore, changes in the composition of TME have different effects on tumor growth ^[28]. To understand the composition of the TME under different modification patterns, we utilized the R package "ESTIMATE" to evaluate the stromal score, immune score and estimate score of HCC patients. By comparing the differences in TME among different modification patterns, we further determined the status of 23 immune cell subtypes in the TME. Using single-sample gene set enrichment analysis (ssGSEA), we calculated the immune infiltration levels and differences in infiltration in the TME of HCC patients. The Wilcoxon rank-sum test was used for the comparison.

2.5.   Identification and functional annotation of DEGs between samples with different modification patterns

We used the R package "limma" to identify differentially expressed genes (DEGs) with significant expression changes between different modification patterns. The significance criteria were set as |log2FC| > 1 and FDR < 0.05. These DEGs play a crucial role in the process of cellular state changes. To investigate the biological processes and signaling pathways in which these DEGs are involved in influencing cellular state processes, we performed gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis using the R package "clusterProfiler". This helped us determine the functional and regulatory roles of the DEGs in these processes.

2.6.   Construction of a prognostic RNAM-S model related to different RNA modification patterns

A prognostic model, named the risk nomogram of clinical characteristics combined with a risk score (RNAM-S), was constructed based on RNA modification patterns. The model was used to calculate risk scores to predict the prognosis of patients with HCC. Univariate Cox analysis was performed to identify prognostic DEGs in patients with different modification patterns. Based on their expression, these DEGs were clustered via RCC according to the optimal cluster number. To verify the stability and effectiveness of the prognostic model, both TCGA-LIHC and GSE14520 cohorts were randomly divided into two groups (the split ratio was 1:1, one group was used to construct the RNAM-S model, whereas the other group was used to verify it). The model was constructed based on prognostic DEGs using the least absolute shrinkage and selection operator (LASSO) and multivariate Cox analyses. To avoid the risk of over-fitting, genes that were most significantly associated with the prognosis of HCC were used to construct the RNAM-S model. Based on the model, risk scores were calculated using the following formula: Risk score = ${\displaystyle \sum {}_{\text{x-}1}^{\text{n}}}(\text{Coefx}\ast RNA\text{expx})$. In the abovementioned equation, "Coefx" is the risk coefficient of each gene, and "RNAexpx" is the expression level of each gene.

In order to assess the predictive performance of RNAM-S, patients with HCC were divided into two groups (high- and low-risk groups) based on the median gene expression. Principal component analysis (PCA) was performed to verify the categorization of patients. The survival of patients in the two risk groups was compared, and receiver operating characteristic (ROC) curves were plotted to verify the predictive accuracy of the prognostic model. The same methods were used to validate RNAM-S. Eventually, the predictive accuracy of RNAM-S among patients with different clinicopathological characteristics was determined via hierarchical analysis.

2.7.   Cell culture

Human HCC cells (HUH-7 and MHCC-97H) and normal liver cells (LO-2) were cultured in high-glucose DMEM (Gibco) supplemented with 10% fetal bovine serum (Procell) in a 5% CO2 incubator (Heal force) at a constant temperature of 37℃.

2.8.   Western blotting

Total protein was extracted from cells (HUH-7, MHCC-97H and LO-2) lysed in RIPA buffer (Solarbio) and quantified using a BCA kit (Solarbio). The expression of KIF20A (1:1000, A15377, ABclonal), PLOD2 (1:1000, A6946, ABclonal), NR1I2 (1:1000, A17038, ABclonal) and NR2F1 (1:1000, A16437, ABclonal) in the three cell lines was detected via western blotting (WB).

2.9.   Real-time reverse transcription polymerase chain reaction

Real-time reverse transcription polymerase chain reaction (qPCR) was performed to examine the mRNA expression of KIF20A, PLOD2, NR1I2 and NR2F1 in HUH-7, MHCC-97H and LO-2 cells. Total RNA was extracted from cells using the Promega kit and quantified. cDNA was synthesized using a reverse transcription kit. Primers for KIF20A, PLOD2, NR1I2 and NR2F1 were synthesized by Bioengineering Co., Ltd. (Shanghai), and the sequences are shown in Table S1. The PCR conditions were as follows: 35 cycles of pre-denaturation (95℃, 2 min), denaturation (95℃, 30 s), annealing (58℃, 30 s) and extension (72℃, 30 s). Beta-actin was used as an internal control, and the relative expression levels of KIF20A, PLOD2, NR1I2 and NR2F1 were evaluated using the 2^-△△Ct method.

2.10.   Analysis of immune status, cancer stem cell index, mutations and drug sensitivity in different risk groups

To evaluate the abundance of tumor-infiltrating immune cells (TIICs) in the two risk groups, the infiltration scores of 22 types of immune cells were calculated using the CIBERSORT algorithm, and the immune status of the groups was visualized. The relationship between the four prognostic genes involved in the RNAM-S model and immune cell infiltration was analyzed. In addition, immune checkpoints with differential expression between the two risk groups were identified and visualized on a box plot. The cancer stem cell (CSC) index helps to identify the degree of tumor differentiation, and mutation analysis helps to identify gene mutations and evaluate their frequency ^[29]. Combining the two analyses may help to identify therapeutic targets. The sensitivity of patients to clinical drugs was assessed by evaluating the half-maximal inhibitory concentration (IC50) of drugs. The aforementioned statistical methods were conducted using the Wilcoxon rank-sum test and Spearman correlation analysis.

2.11.   Construction and validation of a risk nomogram and gene nomogram

Nomograms are widely used to predict the risk and prognosis of tumors ^[30,31]. In this study, a RNAM-S of patients and a gene nomogram based on four prognostic genes was constructed to solve the practical problem that patients with HCC cannot obtain their single or multiple genes during detection in clinical practice.

2.12.   Communication among immune cells in the TME of HCC

First, we utilized dimension reduction algorithms to identify various cell types within the HCC TME. Next, we analyzed the cell type within the HCC TME that is associated with PLOD2, a feature gene significantly upregulated in poor prognosis RNAM-S. With the "CellChat" package, which identifies intercellular communication, we obtained the features and potential pathways of cell communication within the HCC scRNA-seq dataset (GSE210679). To discover more critical intercellular communication within the HCC TME, our team analyzed receptor-ligand pairs associated with PLOD2 and obtained the corresponding signaling pathways. Finally, we presented the expression distribution of the signaling genes using violin plots.

2.13.   Correlation between the expression of four genes and prognosis

To examine the role of the four prognostic genes (KIF20A, PLOD2, NR1I2 and NR2F1) in various cancers, their expression was examined in pan-cancer. The relationship between the expression of these genes and prognosis was examined using data from a pan-cancer dataset from the UCSC database. We also investigated the relationship between the expression of these genes and the OS of pan-cancer patients using univariate Cox regression analysis.

2.14.   Immune regulation, immune checkpoint expression and immune infiltration in pan-cancer

Immune regulatory processes associated with the four genes were analyzed in pan-cancer. The expression data of each characteristic gene and 150 immune modulators were extracted from the pan-cancer dataset, and Spearman coefficients were evaluated to examine the correlation between the expression of the four genes and that of immune modulators. Additionally, the correlation between the four genes and two types of immune checkpoints was analyzed. Finally, four characteristic genes were observed and the correlation among pan-cancer related immune cells was calculated by using six immune cell evaluation algorithms.

3.   Results

3.1.   Overview of changes in RNA methylation-regulating genes in HCC

The role of genes regulating RNA modification (m1A, m5C, m6A and m7G) in HCC was examined. Figure 2A shows the modification sites, molecular structures and related "writer", "eraser" and "reader" proteins of m1A, m5C, m6A and m7G. Enrichment analysis of RNA modification-regulating genes was performed using Metascape. Figure 2B shows the biological processes of the genes, which are mainly involved in the translation, degradation and modification of RNAs. To the best of our knowledge, this study is the first to report mutations in RNA modification-regulating genes in HCC. Gene mutations were found in approximately 10.51% of samples (n = 371 HCC samples) (Figure S1). Additionally, CNVs were widespread in regulatory genes. AGO2 and ALYREF had significant copy number amplifications, followed by METTL3, NCBP2 and EIF4E1B, whereas EIF4G3 and EIF4A1 had significant copy number deletions (Figure 2C). Figure 2D shows the location of CNVs on the chromosome of regulatory genes. Furthermore, the differential expression of RNA modification-regulating genes was examined between normal and HCC tissues. The expression of all genes, except EIF4E1B and NCBP2L, was different between normal and HCC tissues. NUDT10 and EIF4E3 were significantly downregulated and other regulatory genes were significantly upregulated in HCC samples (Figure 2E). Based on the above-mentioned results, a correlation was observed between CNVs and gene expression. For example, AGO2 with copy number amplification was significantly upregulated in HCC tissues, and EIF4E3 with copy number deletion was significantly downregulated in HCC tissues. However, the expression of EIF4E1B with copy number amplification was not different between HCC and normal tissues. Therefore, whether CNVs can regulate gene expression remains unclear. Overall, the results suggest that CNVs influence the expression of RNA modification-regulating genes, which play a potential role in HCC.

Figure 2.  Changes in the expression of RNA modification-regulating genes in HCC. (A) Modification sites, molecular structures and related "writer", "eraser" and "reader" proteins of the four types of RNA modifications. (B) Enrichment analysis of RNA modification-regulating genes (those with similar functions are in the same cluster, and cluster annotations are denoted by various colors). (C) CNV analysis (the pink dots indicate amplification, whereas the green dots indicate deletion). (D) RNA modification regulates the position of CNVs on chromosomes. (E) Differential expression of RNA modification-regulating genes between normal and HCC tissues (* * *, P < 0.001; * *, P < 0.01; *, P < 0.05).

DownLoad: Full-Size Img PowerPoint

3.2.   Identification of different modification patterns and their correlation with clinical prognosis

The expression data of TCGA-LIHC and GSE14520 cohorts were integrated to identify RNA modification patterns in HCC. Cox regression and survival analysis were performed to examine the relationship between the expression of RNA modification-regulating genes and the prognosis of patients with HCC. Figure S2 demonstrates the genes that regulate RNA modification with significant prognostic value in HCC. Based on the expression of these prognostic genes, patients with HCC were divided into two clusters with significant differences in modification patterns, according to the optimal number of clusters obtained via n iterations (Figures 3A, 3B; 329 patients in cluster A and 291 patients in cluster B). Patients in cluster A had a better prognosis (Figure 3C, P = 0.024). Subsequently, the relationship between the two modification patterns and the clinical characteristics of patients was examined. Figure 3D demonstrates significant differences in gene expression and the clinical characteristics of patients between the two clusters. Gene expression was upregulated in cluster B. GSVA was performed to examine the biological functions of genes in the two clusters. Genes in cluster A were significantly enriched in metabolism-related (sugar and amino acid metabolism) and immune pathways (Figure 3E), whereas those in cluster B were mainly enriched in pathways associated with RNA modification and cell cycle. In addition, TME scores (including stromal, immune and total scores) were higher in cluster A than in cluster B (Figure 3F). ssGSEA revealed that the infiltration levels of 12 types of immune cells were significantly different between clusters A and B. The infiltration levels of activated CD4 T cells and type 2 helper T cells were higher in cluster B, whereas the levels of the other 10 types of immune cells were higher in cluster A (P < 0.01, Figure 3G).

Figure 3.  Identification and correlation analysis of two modification patterns. (A) The optimal number of clusters obtained via n iterations (k = 2). (B) Principal component analysis showed significant differences between the two modification patterns. (C) Survival was better in cluster A than in cluster B. (D) Relationship between the two modification patterns and the clinical characteristics of patients with HCC. (E) GSVA of the two modification patterns. (F) TME scores were higher in cluster A. (G) The infiltration levels of 12 types of immune cells were different between the two clusters.

DownLoad: Full-Size Img PowerPoint

3.3.   Identification of gene patterns based on DEGs

A total of 1140 DEGs with RNA modification patterns (Table S2) was identified between clusters A and B. These DEGs were found to be significantly involved in pathways associated with chromosomal regulation, DNA, transcription factors, metabolism, cell proliferation and cancer (Figures 4A, B). A predictive model, named RNAM-S, was constructed to estimate the prognosis of patients with HCC. A total of 862 DEGs associated with the prognosis (OS) of HCC were selected via univariate Cox analysis. Patients were divided into two clusters based on the expression of these prognostic DEGs (Figure 4C). OS was worse in gene cluster A compared to gene cluster B (Figure 4D, P < 0.001). The expression of RNA modification-regulating genes was significantly different between the two clusters, with the expression being higher in gene cluster A (Figure 4E). Compared with genes in cluster B, those in cluster A were strongly associated with advanced TNM (tumor, node and metastasis) stages (stages Ⅲ and Ⅳ) (Figure 4F). Additionally, GSVA revealed that gene cluster B was significantly enriched in metabolism-related pathways, whereas gene cluster A was mainly enriched in cell proliferation- and cancer-related pathways (Figure 4G).

Figure 4.  Identification and correlation analysis of two gene patterns. (A) GO enrichment analysis of DEGs. (B) KEGG enrichment analysis of DEGs. (C) The optimal number of clusters obtained via n iterations (k = 2). (D) Survival was better in cluster B than in cluster A. (E) Differential expression of RNA modification-regulating genes between two gene patterns. (F) The relationship between the two gene patterns and the clinical characteristics of patients with HCC. (G) GSVA of two gene patterns.

DownLoad: Full-Size Img PowerPoint

3.4.   Construction and validation of the prognostic RNAM-S model

The prognostic RNAM-S model was constructed based on DEGs identified in clusters A and B. Figures 5A–C demonstrate the relationship among risk scores, prognosis, modification patterns and gene patterns. In the modeling group, LASSO and multivariate Cox analysis were performed to select prognosis-related DEGs, and four genes (kinesin family member 20A (KIF20A), Procollagen-Lysine, 2-Oxoglutarate 5-Dioxygenase 2 (PLOD2), nuclear receptor subfamily 1 group Ⅰ member 2 (NR1I2) and Nuclear Receptor Subfamily 2 Group F Member 1 (NR2F1)) most significantly related to the prognosis of HCC were eventually identified and used for the construction of RNAM-S. Risk scores were calculated as follows:

Figure 5.  Construction and validation of the prognostic RNAM-S model. (A) Sankey diagram of modification patterns, gene patterns, RNAM-S and clinical survival statuses. (B) RNAM-S distribution in two modification patterns. (C) RNAM-S distribution in two gene patterns. (D) PCA. (E) KM curves for high and low RNAM-S patients. (F) ROC curves for 1-, 3- and 5-year OS. (G-H) Distributions of RNAM-S and survival status in high and low RNAM-S cohorts. (I) The heatmap of four characteristics genes. (J-K) Univariate and multivariate Cox regression analyses.

DownLoad: Full-Size Img PowerPoint

Risk score = (0.2539 * KIF20A) + (0.3209 * PLOD2) + (–0.1964 * NR1I2) + (0.1805 * NR2F1). Patients were divided into two groups (high- and low-risk groups) according to the median risk score. Figure 5D demonstrates that the median risk score appropriately divided the patients into two groups. Survival analysis showed that the prognosis was worse in the high-risk group than in the low-risk group (Figure 5E). ROC curves revealed that the area under curve (AUC) values of the risk scores for predicting survival at 1, 3 and 5 years were 0.804, 0.719 and 0.719, respectively, indicating adequate predictive performance of RNAM-S (all AUC values were > 0.7, Figure 5F). Increases in risk scores decreased the survival of patients (Figures 5G–H). Additionally, the expression of NR1I2 was high in the low-risk group, whereas the expression of KIF20A, PLOD2 and NR2F1 was high in the high-risk group (Figure 5I).

To verify the predictive performance of RNAM-S, the risk scores of the entire datasets (TCGA and GSE14520 cohorts) and validation group were calculated. Patients in each dataset were divided into two groups (high- and low-risk groups) based on their median risk scores. The results revealed that RNAM-S exhibited a stable and effective predictive performance in different datasets (Figures S3A–F and S4A–F). RNAM-S was identified as an independent prognostic factor for HCC via Cox analysis (P < 0.001) (Figures 5J, K). In addition, it efficiently predicted the prognosis of patients with different clinical characteristics (Figure S5A–F).

3.5.   Analysis of immune status, CSC, mutations and drug sensitivity in the two risk groups

Scatter plots (Figure 6A) were plotted to demonstrate the relationship between risk scores and the abundance of immune cells. Risk scores showed positive correlation with the abundance of activated mast cells and activated CD4 memory T cells and negative correlation with the abundance of resting CD4 memory T cells. Significant differences (P < 0.05) were observed in the expression of 40 immune checkpoints between the two risk groups (Figure 6B). Additionally, the CSC index was positively correlated with risk scores. With an increase in risk scores, the higher the similarity between tumor and stem cells in HCC, the greater the possibility of tumor metastasis (Figure 6C). The mutation frequency of TTN was highest in the low-risk group (24%), whereas that of TP53 was highest in the high-risk group (35%) (Figures 6D, E). Furthermore, the sensitivity of the two risk groups to commonly used antineoplastic drugs was evaluated. The IC50 values of commonly used chemotherapeutic drugs (cisplatin, cytarabine, doxorubicin, gemcitabine, methotrexate and paclitaxel) were lower in the high-risk group, whereas those of drugs targeting epidermal growth factor receptor (EGFR) (Mortisanib, Afatinib, Gefitinib, Lapatinib and Remfitinib) were lower in the low-risk group. In addition, the IC50 values of camptothecin and shikonin were lower in the high-risk group (Figures 7A–P).

Figure 6.  Analysis of immune status, CSC index and mutations in the two risk groups. (A) Relationship between risk scores and the abundance of immune cells. (B) Expression of 40 immune checkpoints in two risk groups. (C) Relationship between the CSC index and risk scores. (D, E) Somatic mutations in the two risk groups.

DownLoad: Full-Size Img PowerPoint

Figure 7.  Drug sensitivity analysis. (A–F) Relationship between risk scores and chemosensitivity. (G) Relationship between risk scores and sensitivity to drugs targeting EGFR. (O–P) Relationship between risk scores and sensitivity to biological drugs.

DownLoad: Full-Size Img PowerPoint

3.6.   Risk and gene nomograms

To assess the clinical applicability of RNAM-S in predicting the OS rate of patients with HCC, a risk nomogram integrating clinical characteristics and risk scores (RNAM-S) and a gene nomogram based on the four prognostic genes (KIF20A, PLOD2, NR1I2 and NR2F1) was constructed for predicting the 1-, 3- and 5-year OS rate of patients with HCC (Figures 8A, D). Calibration curves were plotted to verify the predictive accuracy of the risk and gene nomograms. The calibration curve of the risk nomogram demonstrated that the predicted diagonal line in the risk nomogram was close to and parallel to that of an ideal model (grey diagonal line), indicating that the accuracy of the risk nomogram is satisfactory (Figure 8B). Furthermore, the predictive accuracy of the risk nomogram was compared with that of RNAM-S and the currently used clinical features. The predictive accuracy of the nomogram was superior to that of clinical staging, indicating the excellent predictive performance of the nomogram. In addition, calibration curves (Figures 8E, F) validated the adequate predictive accuracy of the gene nomogram. The predicted line (red) was close to that of an ideal model (grey), and the predicted survival probability (red) was similar to the actual survival probability (blue).

Figure 8.  Construction and validation of prognostic RNAM-S. (A) The risk nomogram for 1-, 3- and 5-year OS of the patients of HCC. (B, C) The calibration curves of the risk nomogram. (D) The gene nomogram for 1-, 3- and 5-year OS of the patients of HCC. (E, F) The calibration curves of the gene nomogram.

DownLoad: Full-Size Img PowerPoint

3.7.   Verification of the expression of the four prognostic genes

WB and qRT-PCR were used to verify the expression of the four prognostic genes (KIF20A, PLOD2, NR1I2 and NR2F1) in vitro. WB revealed that the expression of the low-risk gene NR1I2 was higher in LO-2 cells (human normal liver cells) than in MHCC-97H and HUH-7 cells (human HCC cells), whereas the expression of the high-risk genes KIF20A, PLOD2 and NR2F1 was higher in MHCC-97H and HUH-7 cells than in LO-2 cells. KIF20A and PLOD2 were significantly overexpressed in HUH-7 and MHCC-97H cells. NR2F1 was significantly overexpressed in HUH-7 cells and considerably overexpressed in MHCC-97H cells. The results of qRT-PCR and WB were consistent, indicating that the four prognostic genes are significantly dysregulated in HCC.

3.8.   Analysis of single-cell sequencing data

An algorithm for dimensionality reduction was used to identify various types of cell populations in the TME of HCC. The TME of HCC contained six types of cell populations: T cells, NK cells, B cells, hepatocytes, monocytes and endothelial cells (Figure S6). Cell populations associated with PLOD2 were further analyzed. PLOD2 was significantly overexpressed in prognostic RNAM-S and associated with poor prognosis of HCC. PLOD2 expression correlated with the abundance of hepatocytes and endothelial cells within the TME of HCC. The "CellChat" package was used to assess intercellular communication and related pathways using data from the scRNA-seq dataset GSE210679. Hepatocytes, endothelial cells, monocytes and NK cells were found to be closely related to other cell populations (Figures 11D, E). To examine the critical role of intercellular communication in TME, PLOD2-related receptor and ligand pairs and the related signaling pathways were identified. The CD74–CXCR4 and CD74–CD44 complexes and MIF signaling pathway were found to play a key role in the TME of HCC (Figure 11F). Furthermore, the effects of these ligand–receptor pairs on immune cells were examined. The CD74–CXCR4 complex mainly transmits the signals of hepatocytes and endothelial cells, which interact with T cells, B cells and monocytes. The CD74–CD44 complex mainly transmits the signals of hepatocytes, endothelial cells and T cells, which interact with each other (Figures 11G, H). The expression and distribution of the signaling genes were visualized on a violin plot. As shown in Figure 11I, CD74 is expressed in all six cell populations; CXCR4 is mainly expressed by B cells, endothelial cells, monocytes, NK cells and T cells and CD44 is mainly expressed in B cells, monocyte, NK cells and T cells.

Figure 9.  Verification of the expression of the four prognostic genes via WB.

DownLoad: Full-Size Img PowerPoint

Figure 10.  Verification of the expression of the four prognostic genes via qRT-PCR: (A) KIF20A; (B) PLOD2; (C) NR1L2; (D) NR2F1.

DownLoad: Full-Size Img PowerPoint

Figure 11.  Construction and validation of prognostic RNAM-S. (A-C) The expression of PLOD2 in different immune cells. (D, E) The number of ligand receptor interactions detected between different cell types. (F) Ligand receptor contribution to the overall signaling pathway. The CD74-CXCR4 ligand receptor pair contributed the most, followed by the CD74-CD44 ligand receptor pair. (G, H) Receptor ligand pair interactions between immune cells. (I) The distribution and expression level of signal genes involved in the signal pathway networks.

DownLoad: Full-Size Img PowerPoint

3.9.   Expression of the four prognostic genes in pan-cancer and their relationship with prognosis

The expression of the four prognostic genes (KIF20A, NR1I2, NR2F1 and PLOD2) was evaluated in various cancers. KIF20A expression was higher in all cancer tissues compared to the corresponding normal tissues (Figure S7A). NR1I2 expression was high in 19 types of cancers but low in 8 types of cancers, including LIHC (Figure S7B). NR2F1 expression was high in 7 types of cancers, including LIHC, but low in 22 types of cancers (Figure S7C). PLOD2 expression was high in 20 types of cancers, including LIHC, but low in PRAD, ACC and KICH (Figure S7D). Furthermore, the relationship between the expression of these genes and OS was examined in pan-cancer. High KIF20A expression was significantly correlated with the prognosis of 17 types of cancers. However, KIF20A expression was low and the prognosis was poor in THYM (Figure S7E). Low NR1I2 expression was significantly associated with the poor prognosis of LIHC, READ and ACC, whereas high NR1I2 expression was associated with the prognosis of 5 types of cancers (Figure S7F). NR2F1 expression was associated with the prognosis of 8 types of cancers. High NR2F1 expression was correlated with the prognosis of 5 types of cancers. Patients with GBMLGG, HNSC and KIRC with low expression of NR2F1 had a poor prognosis (Figure S7G). High PLOD2 expression was significantly associated with the poor prognosis of 17 types of cancers (Figure S7H).

3.10.   Immune analysis of the four prognostic genes in pan-cancer

Understanding the effects of the four prognostic genes on the immune system is important for determining the cancer types that may benefit from immunotherapy based on these genes. The expression of the four prognostic genes was positively correlated with that of most immunomodulators (Figures S8A–D) and immune checkpoints (Figures S8E–H) in various cancer types. The expression of KIF20A was positively correlated with that of CD276 and HMGB1 in most cancer types. The expression of CD276 and VEGFA was also positively correlated with that of PLOD2 in most cancers. Additionally, the four genes showed specificity in terms of immune infiltration (Figures S9A–D). These results suggest that the four prognostic genes are significantly associated with immune checkpoints and immune cell infiltration in pan-cancer.

4.   Discussion

HCC is a major health issue, and its high invasiveness and difficulty of treatment make it a challenge that urgently needs to be addressed within the current medical field. Therefore, research on the prognosis of HCC has received increasing attention. A multicenter study has shown that traditional pathological factors such as TNM staging, alpha-fetoprotein (AFP) and tumor size are related to the prognosis of HCC ^[32]. However, these clinical factors lack early predictive and diagnostic accuracy. Based on this, many molecular biomarkers have been used in the research of HCC prognosis and treatment. For example, another study has shown that aldoketo reductase 1B10 (AKR1B10) is an adverse prognostic indicator for HCC, and its predictive effect is superior to AFP ^[33]. RNA modification plays an important role in the development of liver cancer ^[34]. To date, most studies have focused on a single type of RNA modification, whereas only a few studies have investigated multi-type RNA modifications (m1A, m5C and m6A) ^[35]. Moreover, the relationship between the overall regulation of multi-type RNA modifications and the TME of HCC remains unknown. Therefore, in this study, we analyzed RNA modification-regulating genes (m1A, m5C, m6A) that are significantly associated with the occurrence and development of HCC and m7G modification-regulating genes that are associated with human cancers. We analyzed changes in various types of RNA modifications in HCC at the gene and transcription levels and examined the relationship between RNA modifications and the TME of HCC.

Based on the expression of RNA modification-regulating genes in HCC samples, two modification patterns were identified, which were characterized by different pathological characteristics, prognoses, TMEs and immune infiltration levels. Compared with genes in cluster B, those in cluster A were significantly enriched in metabolic and immune-related pathways. Additionally, TME scores were significantly higher in cluster A. The infiltration levels of various immune cells were higher in cluster A, which resulted in a better prognosis. The DEGs identified between clusters A and B were mainly enriched in metabolic and cell proliferation-related pathways. Based on these DEGs, two gene patterns were further identified. The characteristics of the two gene patterns and modification patterns were similar. High expression of RNA modification-regulating genes in both gene and modification patterns resulted in a worse OS rate. These findings suggest that multi-type RNA modifications can be used to predict the prognosis and treatment response of patients with HCC. Therefore, a prognostic model named RNAM-S, based on the DEGs related to modification patterns, was constructed and its predictive ability was verified in multiple aspects. Additionally, the expression of the four prognostic genes (KIF20A, NR1I2, NR2F1 and PLOD2) involved in RNAM-S was examined in normal liver cells and two types of HCC cells through in vitro experiments. Based on the median risk score, patients were divided into high- and low-risk groups. Significant differences were observed in immune status, CSC index, mutations and drug sensitivity between the two groups. Furthermore, risk and gene nomograms were constructed to assess the clinical applicability of RNAM-S. The prognostic genes with significantly high expression in RNAM-S and a significant correlation with poor prognosis of HCC were associated with abundant hepatocytes and endothelial cells in the TME of HCC. Receptor and ligand pairs related to PLOD2, namely, CD74–CXCR4 and CD74–CD44, and the MIF signaling pathway, play a key role in the TME of HCC. Finally, the role of the four genes involved in RNAM-S was examined in multiple cancer types. KIF20A, NR1I2, NR2F1 and PLOD2 were differentially expressed in pan-cancer and were associated with the prognosis and immune response of patients. The prognostic RNAM-S model has excellent predictive ability and can be used to predict the prognosis of several types of cancer. It can help understand the molecular mechanisms underlying the development of not only HCC but also several other types of cancer to a certain extent and facilitate the development of new strategies for immunotherapy.

KIF20A's function in various tumor cells and its involvement in regulating tumorigenesis have been observed. Via the JAK/STAT3 signaling pathway, KIF20A promotes the growth of colorectal cancer cells ^[36]. Mediation of the KDELR2-KIF20A axis promotes Golgi-mediated MMP secretion to drive bladder cancer growth and metastasis ^[37]. KIF20A activates androgen receptor autocrine to promote prostate cancer progression and is a promising target for cancer immunotherapy ^[38]. A regulatory role of KIF20A in HCC has also been observed. KIF20A knockdown inhibited the proliferation of HCC cells and promoted the chemosensitivity of HCC cells to sorafenib and cisplatin ^[39]. In addition, the downregulation of KIF20A expression exacerbated chromosomal instability in HCC and inhibited HCC cell viability in vitro and in vivo ^[40]. Pregnane X receptor (PXR, NR1I2) increases apoptosis of tumor cells, inhibits the proliferation of tumor cells and protects normal cells from cancer. In colon cancer, NR1I2 promotes apoptosis and inhibits the proliferation of human colon cancer cells by arresting the colon cancer cell cycle at G0/G1 phase, which is mediated by upregulated p21 and downregulated E2F/RB signaling ^[41]. The ligand RIF activates NR1I2, which arrests the cell cycle at the G2/M phase and inhibits the growth and proliferation of cells in breast cancer and cervical cancer ^[42]. However, because NR1I2 promotes drug metabolism mediated by MDR1/CYP3A4 and facilitates the repair of DNA damage, it is tremendously tolerant to genotoxic drugs. Thus, NR1I2 plays a role in drug resistance during cancer treatment ^[43]. NR2F1, a biomarker of tumor dormancy, inhibits metastasis and proliferation of cancer cells in several malignancies. Administration of NR2F1 agonists to primary HNSCC induced tumor dormancy and inhibited metastasis ^[44]. However, NR2F1 can also promote tumor invasion and metastasis in breast cancer ^[45]. NR2F1 showed different expression trends in the two HCC cells in the present study. Overall, different cell lines used and different cancer types may lead to the conflicting findings described above. Nevertheless, NR2F1 with RNA modifying function is indeed involved in cancer onset and progression. PLOD2 is a key gene that mediates the formation of stable collagen crosslinks. Accelerated collagen secretion and extracellular matrix (ECM) remodeling in cancer cells can lead to cancer cell invasion and metastasis ^[46]. Numerous studies have shown that expression of PLOD2 is upregulated in cancers and promotes cancer cell growth. In addition, its expression is also significantly upregulated in osteosarcoma, which can promote migration, invasion and angiogenesis in osteosarcoma in both in vivo and in vitro experiments ^[47]. A study of gastric cancer revealed that PLOD2 was substantially linked to peritoneal dissemination; PLOD2 had a considerably higher expression at both mRNA and protein levels, and silencing PLOD2 substantially attenuated cell invasion and migration in vitro ^[48]. Hypoxia induced PLOD2 expression, and its expression level was significantly elevated in HCC and associated with shorter disease-free survival ^[49]. In summary, all four feature genes in the RNAM-S model play a crucial role in cancer onset and progression and can act as RNA-modification regulatory genes. Therefore, this study further explores their effects on the HCC TME under genetic and transcriptional alterations and constructs a reliable prognostic model for survival prediction among HCC patients.

TME contains various types of cells. As an important part of TME, immune cells participate in various immune activities of tumors and may help to understand tumor immune responses ^[50]. In this study, the modification and gene patterns with more significant immune scores and immune cell infiltration resulted in better OS. Additionally, in both modification and gene patterns, higher risk scores were associated with lower immune scores and immunosuppression in cluster B, whereas lower risk scores were associated with enhanced immune activation and higher immune scores in cluster A. These findings highlight the key role of TME in HCC. T cells play an important role in immune defense in HCC ^[51,52]. The occurrence of liver diseases, including liver tumors, has a strong impact on the immune function of the liver and has a significant impact on the progression of liver cancer ^[53]. After comprehensive analysis, Lurje found that activating CD8+ T cells through the tolerance or immunogenicity direction of mediated T-cell reaction can effectively induce anti-tumor defense capacity in the liver ^[54]. Macrophages are also a member of immune cells in TME that participate in various immune reactions. This study revealed that the main ligand and receptor pairs related to PLOD2, which is highly expressed in RNAM-S but is significantly associated with the poor prognosis of HCC, affect the interaction between T cells and monocytes, which may be the main reason for the poor prognosis. Therefore, the treatment and prognosis of HCC are closely related to the interaction and mediation of immune cells in TME. This phenomenon should be considered while developing therapeutic strategies for HCC in the future.

Cell-to-cell communication is crucial for biological signal transduction and plays an important role in complex diseases. Through single-cell data analysis, it is advantageous to reveal the functional characteristics of individual cells in the TME and the potential molecular mechanisms of intercellular crosstalk ^[55]. In our study, we found close relationships between hepatocytes, endothelial cells, monocytes, NK cells and other cell populations. Among them, the CD74-CXCR4 and CD74-CD44 receptors, as well as the macrophage migration inhibitory factor (MIF) signaling pathway, play important roles in the HCC TME. MIF can regulate the polarization of tumor-associated macrophages (TAM) in tumors; macrophages usually differentiate into pro-tumor M2 TAM, promoting tumor cell proliferation and metastasis while inhibiting anti-tumor immune responses ^[56]. MIF is an immunomodulator in the TME, inducing the formation of a tumor inhibitory immune microenvironment ^[56].

The expression of the four prognostic genes differed in several cancers and was significantly correlated with the prognosis and immunity of patients. Therefore, KIF20A, NR1I2, NR2F1 and PLOD2 could potentially serve as biomarkers and should be thoroughly investigated in future studies. However, this study does have some limitations. The mRNA and scRNA-seq datasets used in the study are published datasets. The expression of the four prognostic genes used to construct RNAM-S was assessed through in vitro experiments and pan-cancer analysis. However, it should be validated via clinical and in vivo experiments. In summary, our study demonstrates that RNA modifications and their mediated effects have significant impacts on HCC. It is worthwhile to utilize more methods to further investigate the role of RNA modifications in tumors, such as imaging techniques ^[57,58] and metabolomics ^[59]. Through biological imaging analysis and tumor cell metabolic fingerprint analysis, the dynamic changes of RNA modifications in tumors can be unraveled. With the rapid development of life sciences and computer sciences, these techniques have been widely applied in tumor research. Therefore, this will be our future direction for further research.

5.   Conclusions

The comprehensive analysis of various RNA modification-regulating genes in HCC revealed that RNA modifications play an extensive regulatory role in HCC. The overall mediation of RNA modifications affects the TME, clinical characteristics and prognosis of HCC. The prognostic RNAM-S model constructed based on RNA modification-regulating genes can help predict the OS of HCC. The risk and gene nomograms constructed in this study validate the predictive accuracy of RNAM-S. The findings of this study highlight the clinical significance of RNA modifications and may guide the development of individualized treatment strategies for HCC.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

Yuanqian Yao performed the research and wrote the paper. Xiaohua Hong designed the research study. Jianlin Lv and Guangyao Wang analyzed the data. This work was supported by The National Natural Science Foundation of China (Grant No. 81860839) and The Natural Science Foundation of Guangxi (Grant No. 2020GXNSFAA238020, 2022GXNSFAA035446).

Availability of data and materials

The datasets used and analyzed during the current study are available from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/, TCGA-LIHC) and Gene expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/).

Conflict of interest

The authors declare there is no conflict of interest.