Comprehensive Analysis of METTLs (METTL1/13/18/21A/23/25/2A/2B/5/6/9) and Associated mRNA Risk Signature in Hepatocellular Carcinoma

Currently, 80%–90% of liver cancers are hepatocellular carcinomas (HCC). HCC patients develop insidiously and have an inferior prognosis. The methyltransferase-like (METTL) family principal members are strongly associated with epigenetic and tumor progression. The present study mainly analyzed the value of METTLs (METTL1/13/18/21A/23/25/2A/2B/5/6/9) and associated mRNA risk signature for HCC. METTLs expression is upregulated in HCC and is a poor prognostic factor in HCC. METTLs were upregulated in patients older than 60 and associated with grade. Except for METTL25, the remaining 10 genes were associated with the HCC stage, invasion depth (T). In addition, METTLs showed an overall alteration rate of 50%. Except for METTL13/2A/25/9, the expression of the other seven genes was significantly associated with overall survival, disease-specific survival, and progression-free survival. Multivariate studies have shown that METTL21A/6 can be an independent prognostic marker in HCC. A total of 664 mRNAs were selected based on Pearson correlation coefficient (R > 0.5), unsupervised consensus clustering, weighted coexpression network analysis, and univariate Cox analysis. These mRNAs were significantly associated with METTLs and were poor prognostic factors in HCC patients. The least absolute shrinkage and selection operator (lasso) was used to construct the best METTLs associated with mRNA risk signature. The mRNA risk signature was significantly associated with age, stage, and t grade. The mRNA high-risk group had higher TP53 and RB1 mutations. This study constructed a nomogram with the mRNA risk profile and clinicopathological features, which could better predict the OS of individuals with HCC. We also analyzed associations between METTLs and mRNA risk signatures in epithelial-mesenchymal transition, immune checkpoints, immune cell infiltration, tumor mutational burden, microsatellite instability, cancer stem cells, tumor pathways, and drug sensitivity. In addition, this study constructed a protein interaction network network including METTLs and mRNA risk signature genes related to tumor microenvironment remodeling based on single-cell sequencing. In conclusion, this study provides a theoretical basis for the mechanism, biomarker screening, and treatment of HCC.


Introduction
Hepatocellular carcinoma (HCC) is one of the digestive tract tumors and an aggressive type of tumor [1].Although HCC incidence ranks sixth worldwide, HCC mortality ranks third worldwide [2].Currently, the diagnosis of HCC is limited, resulting in the fact that most patients with HCC are advanced stage when diagnosed [3].In addition, HCC treatment approaches are also limited, and patients with advanced HCC who are treated also develop tumor cell metastasis, ultimately leading to patient death [4].Therefore, finding effective biomarkers and therapeutic targets for HCC is crucial.
Among the 33 METTL family members, 11 METTLs (METTL1/13/18/21A/23/25/2A/2B/5/6/9) with upregulated expression in HCC and poor prognosis for HCC patients were enrolled in this study.In addition, this study employed multiple algorithms to construct the mRNA risk signature associated with these 11 genes separately.Ultimately, this study comprehensively analyzed the correlation between METTLs and mRNA risk signatures in terms of prognosis, clinicopathological features, mutations, tumor microenvironment (TME), epithelial-mesenchymal transition (EMT), immune cell infiltration, tumor mutational burden (TMB), microsatellite instability (MSI), cancer stem cells (RNAss), cancer pathways, and drug sensitivity.In conclusion, this study provides a theoretical basis for the prognosis, biomarker screening, mechanism, and drug screening of HCC.

Acquisition of mRNAs Closely Associated with METTLs.
First, Pearson correlation analysis in the R software "limma" package was used to screen the related mRNAs of METTLs in this study (Pearson r > 0.5 and P <0:05).Unsupervised consensus clustering is an algorithm for k-means machine learning [19].The present study subsequently employed unsupervised consensus clustering to analyze the TCGA cohort based on the expression of METTLs.The R software "ConsensusClusterPlus" package was used for the clustering.And OS curves were used to evaluate the prognosis of different clusters for HCC patients.The R software "limma" and "ggplot2" packages were adopted to screen the differential mRNAs between different clusters (|log FC| > 0.585, P <0:05).This study further used R software "limma" and the "WGCNA" package to construct a coexpression network of 11 genes in the METTL family based on weighted gene coexpression network analysis (WGCNA) [20].The present study intersected the relevant mRNAs of 11 genes in the METTL family, the upregulated mRNAs in the subcohort, and the coexpressed mRNAs in the WGCNA.The intersection genes were further screened in this study by univariate Cox analysis.This study defines mRNAs with poor prognosis for HCC as genes closely related to METTLs.Finally, "clusterProfiler" and "org.HS.Eg. db" packages were adopted to perform Kyoto Encyclopedia of Genes and Genomes (KEGG) and GO enrichment analysis of closely related mRNAs of METTLs.
2.4.Construction of an mRNA Risk Signature.The present study constructed an associated risk signature based on the closely related mRNAs of METTLs.This study employed R software "glmnet" and "survival" packages for lasso regression and multivariate Cox analysis to construct the most suitable risk signature.The risk score for patients was calculated from the normalized expression level of each gene and the corresponding regression coefficient with the formula: Risk score = ∑coefi * exp.The R software "survival," "survminer," and "survminer" packages were used to draw the OS curves and time-dependent ROC curves for the mRNA risk signature.The R software "Rtsne" and "ggplot2" packages were used to draw principal component analysis (PCA) and t-SNE scatter plots to distinguish risk patients among the risk features.The data above are from the TCGA cohort (tumor: n = 370).In addition, the validation set was obtained from the ICGC database (tumor: n = 240).

Mutational Analysis and Construction of a Nomogram
Prediction Model for Risk Signature.The SangerBox 3.0 tool was used (https://doi.org/10.1002/imt2.36).Mutation differences in significant HCC genes between high-and low-risk  The R software "timeROC," "ggdca," "survminer," and "survival" packages were used for analysis and visualization.
2.6.Integrated Analysis of METTLs and mRNA Risk Signatures.The R software "pRRophetic" package was used and based on the drug IC50 values to screen sensitive drugs to METTLs and mRNA risk profiles [21].Spearman correlation analysis of METTLs and mRNA risk signatures with EMT core genes, immune checkpoints, immune cells, tumor-associated macrophage (TAM) surface genes, TMB, MSI, RNAss, and IC50 were assessed using R software "reshape2" and "RColorBrewer" packages.Gene set variant analysis (GSVA) was employed to evaluate the biological functions of METTLs/ mRNA risk signatures based on the gene sets of KEGG terms and hallmark terms.

2.7.
Processing and Analysis of Single-Cell Data.The HCC single-cell dataset of this study (GSE146115, n = 4) was obtained from the GEO database.R software "Seurat," "dplyr," and "singleR" packages were used for quality control processing and subpopulation annotation of single-cell sequencing data in this study [22,23].Quality control conditions were set to retain genes expressed in at least three intracellular compartments, cells expressing at least 200 genes, cells with the number of detected genes between 200 and 2,500, and cells with less than 5% expression of mitochondrial genes.Lognormalize was used for normalization in this study.Cell communication was analyzed using the "sqjin/CellChat" package.Based on Genemania database construction (http://genemania.org/)[24] to construct a TME remodeling associated protein interaction network (PPI).
2.8.Statistical Analysis.R software 4.0.3 was used to perform all statistical analyses.Differences between the two groups were analyzed using the Wilcoxon rank sum test or paired t-test.

Results
Ultimately, the related mRNAs of METTLs, the upregulated mRNAs in the subpopulation, and the coexpressed mRNAs in WGCNA were intersected.The intersection RNAs were subjected to univariate Cox analysis.We obtained a total of 664 mRNAs (Figure 4(j)).These mRNAs showed a close positive correlation with METTLs.GO and KEGG analyses were performed on these mRNAs.The biological functions of these RNAs were shown to be related to DNA replication, chromosome splicing, RNA transport, and cell cycle (Figures 4(k) and 4(l)).

Construction of an mRNA Risk Signature Associated with
METTLs.In this study, based on 664 mRNAs associated with METTLs, the best mRNA risk signature was constructed by multivariate Cox regression analysis using lasso (Figures S1(a (g-and h) soft thresholding and dendrograms, (i) heatmap of correlation between module eigengenes and METTLs, (j) intersection genes and results of univariate Cox analysis, (k) GO analysis, and (l) KEGG analysis.12 Analytical Cellular Pathology best risk signature.To verify the reliability of these 11 mRNAs, the Heatmap of coexpression between the construct METTLs and these 11 mRNAs was performed using Spearman's correlation test.As shown in Figure S2(a), METTLs were significantly and positively correlated with these 11 mRNAs (r > 0.3, P <0:05).Moreover, these 11 mRNAs were all upregulated in HCC (Figure S2(b)).Moreover, in the OS curve, the high expression of these 11 mRNAs showed a significant decrease in survival (P <0:05) (Figure 5(a)).In the DSS curves, the survival of the high-expression group of mRNAs except TMEM69 decreased significantly (P <0:05) (Figure 5(b)).In the PFI curves, the high expression of mRNAs except PSRC1 and TMEM69 showed a significant decrease in survival (P <0:05) (Figure 5 6(a)).The area under the curve (AUC) at 1, 3, and 5 years in the time-dependent ROC curve was 0.783, 0.727, and 0.719, respectively (Figure 6(b)).There were more deaths in the high-risk group (Figure 6(c)).PCA and t-SNE scatter plot results indicated that HCC patients at different risks were well able to be separated into two clusters (Figure 6(d)).The same was true for the mRNA validation set (ICGC) (Figure 6(e)-6(h)).Moreover, combined with clinicopathological features, univariate/multivariate Cox analysis indicated that the mRNA risk signature was an independent poor prognostic factor for HCC patients (Figure 6(i)).
3.6.Construction of an mRNA Risk Signature Nomogram Prediction Model.In this study, a nomogram prediction model was constructed to evaluate the predictive value of the mRNA risk signature for the OS of HCC.In this study, a nomogram was constructed combining clinicopathological features and mRNA risk features to predict the survival rates of HCC patients with 1-, 3-, and 5-year OS (Figure 7(a)).The calibration curves indicated that the nomogram prediction model was able to predict the 1-, 3-, and 5-year OS of HCC patients with better accuracy (Figure 7(b)).The time ROC curves indicated that the 3-and 5-year nomogram prediction models outperformed the mRNA risk signature (Figure 7(c)).DCA indicated that the mRNA risk signature was more accurate (Figure 7(d)).

Correlation of mRNA Risk Signature with Clinicopathological
Features and Genetic Alterations.This study evaluated the correlation of mRNA risk signature with clinicopathological features and genetic alterations.The results showed significant differences between T, grade, and stage with high-and low-risk groups of mRNAs (Figure 8(a)).In addition, mRNA scores were significantly different from age, grade, stage, and t (Figure 8(b)).In this study, we further compared the mutational differences between high-and low-risk groups in 20 genes (TP53, CTNNB1, ALB, AXIN2, KEAP1, BAP1, NFE2L2, LZTR1, RB1, PIK3CA, KRAS, IL6ST, CDKN2A, ARID2, ARID1A, ACVR2A, NRAS, HISR1H1C, PTEN, and ERRFI1) that are predominantly mutated in HCC [25].The results showed that TP53 and RB1 mutations were significantly higher in the high-risk group (P <0:05) (Figure 8(c)).

Multiangle Analysis of METTLs and mRNA Risk Signatures.
This study further used multiple angles to assess the role of METTLs and mRNA risk signatures in HCC, including EMT, immune checkpoints, immune cell infiltration, TAM markers, TMB, MSI, and RNAss.The risk network plots of METTLs and mRNA risk signature were first constructed in this study (Figure 9(a)).The EMT results indicated that METTLs and mRNA risk scores were positively correlated with most of the EMT core genes, among which they were all significantly correlated with MMP9 (P <0:05) (Figure 9(b)).The results of immune checkpoint analysis indicated that METTLs and mRNA risk score were significantly and positively correlated with most immune checkpoints (Figure 9(c)).Immune cell infiltration results indicated that METTLs and mRNA risk scores were positively correlated with multiple immune cell infiltration, among which all genes and mRNA risk scores except METTL25 were significantly correlated with M0 macrophages (P <0:05) (Figure 9(d)).We further analyzed the correlations between the METTLs and mRNA risk scores and TAM markers.The results showed that the METTL9 and mRNA risk scores were significantly (P <0:05) correlated with all TAM markers (Figure 9(e)).In addition, the METTL1/18/23/2A/2B/5/6/9 and mRNA risk scores were significantly positively correlated with TMB values (P <0:05) (Figure 9(f)).There was a significant positive correlation between METTL23/5 and mRNA risk score and MSI (P <0:05) (Figure 9(g)).The METTL1/13/18/21A/23/2A/ 2B/5/6 and mRNA risk scores showed significant positive correlations with RNAss (P <0:05) (Figure 9(h)).In conclusion, METTLs and mRNA risk signature play a crucial role in HCC and may serve as novel therapeutic targets in tumor therapy.

Discussion
HCC is a malignant tumor, and the pathogenesis remains unclear to date [26].HCC has a high mortality rate, mainly because most patients are advanced when diagnosed or have developed local metastases [27].There are many therapeutic approaches for HCC, including early diagnosis, surgery, drugs, immunotherapy, and so forth.However, current treatments do not achieve the desired efficacy for patients with advanced HCC stages [28].Therefore, it is crucial to find out effective biomarkers and therapeutic targets for HCC.In the present study, METTLs (METTL1/13/18/21A/23/25/2A/2B/ 5/6/9), among 33 METTL family members are poor prognostic factors in HCC patients and upregulated in.In addition, METTLs have a better ability to discriminate between tumor and nontumor.Therefore, we focused on these 11 genes as research objects and comprehensively analyzed their values for HCC.
In the present study, OS results indicated that the METTLs high-expression group had shorter; the DSS results indicated a shorter survival time in the high-expression group of the remaining eight genes, except for METTL13/ 23/2A; the PFS results indicated a shorter survival time in the high-expression group of the remaining eight genes except for METTL13/25/9.Univariate Cox analysis indicated that both METTLs were poor prognostic factors in HCC.Combined with multiple HCC clinicopathological features, multivariate Cox analysis indicated that METTL21A/6 was an independent prognostic marker for HCC.Therefore, overall METTLs are poor prognostic for HCC.Previous studies have shown that the probability of having HCC increases progressively with age as well [29].Changes in the stage or grade of HCC can also be accompanied by alterations in specific regulated genes [30].Invasion depth (T) is also one of the significant factors contributing to the poor prognosis of HCC patients [31].In the clinical relevance analysis of this study, METTLs were significantly upregulated in patients older than 60 years.For the HCC grade, the expression of METTLs was significantly different from the grade, and the higher grade increased the expression of METTLs.For stage, except for METTL25, the expression of the remaining 10 genes increased with stage.For T, except for METTL13/25/ 2B/9, the expression of the remaining seven genes was significantly different from that of T.Moreover, 50% of the total genetic alterations in METTLs have been reported.Currently, METTL21A/23/25/2A/2B/9 are not studied in HCC.However, in other tumors, METTL2a is a potential N3 methylcytidine (M3C)-the related gene that promotes breast cancer development [32].The rise of METTL9 is associated with peritoneal metastasis of gastric cancer [33].In multiple previous studies, METTL1/13/18/5/6 were able to promote HCC progression as oncogenes.Chen et al. [10] showed that the knockdown of METTL1 can reduce m7G tRNA modification, reducing the translation of mRNA and thereby inhibiting HCC development.Decreasing METTL1 also reduces HCC recurrence [34].Li et al. [35] showed that upregulation of METTL13 could promote HCC growth and metastasis.Downregulation of METTL18 can inhibit the proliferation and metastasis of HCC [17].Peng et al. [36] showed that METTL5 is associated with 18S rRNA m6A modification, and downregulating METTL5 can inhibit HCC progression.Bolatkan et al. [12] showed that downregulating METTL6 could reduce metastasis of HCC.Taken together, METTLs have an essential role in HCC prognosis and progression.
In this study, 11 mRNA profiles of METTLs-related risk features were constructed by multiple algorithms.All these 11 mRNAs were closely associated with METTLs.These 11 mRNAs were upregulated in HCC and were poor prognostic factors for HCC.The mRNA risk signature could well stratify HCC patients into high-and low-risk groups.In the present study, we found that the mRNA risk signature was an independent poor prognostic factor for HCC patients.In clinical correlation analysis, the mRNA risk signature was associated with age, grade, stage, and T in HCC patients.In addition, the nomogram prediction model we constructed had high accuracy for OS 1, 3, and 5 years in patients with HCC.This nomogram prediction model had higher 3-and 5-year accuracies than the mRNA risk signature.However, DCA indicated more robust accuracy of mRNA risk signature.In the mutation analysis, we found that TP53 and RB1 were mutated with increased frequency in the high-risk group.It is well known that TP53 and RB1 mutations are one of the  most common inheritances in HCC, and the more TP53 mutations, the worse the prognosis of HCC patients [37,38].Therefore, the mRNA risk signature we constructed has an essential effect on the prognosis and development of HCC.Previous studies have shown that seven mRNAs in the mRNA risk signature can act as oncogenes to promote HCC progression, including PPM1G [39], PHOSPHO2 [40], YBX1 [41], PSRC1 [42], UCK2 [43], EZH2 [44], and CDCA8 [45].However, MEX3A, TAF3, TMEM69, and DYNC1LI1 have not been experimentally studied in HCC.However, in other tumors, they are also able to act as oncogenes to promote tumor progression [46][47][48][49].These studies further illustrate that the mRNA risk signature we constructed is an important factor that promotes HCC occurrence.Further integrative analysis in this study found that METTLs and mRNA risk signatures were associated with EMT, multiple immune checkpoints, multiple immune cell infiltration, TMB, MSI, RNAs, and cancer pathways.In terms of drug sensitivity, 27 drugs sensitive to METTLs and mRNA risk profiles were screened in this study.In conclusion, this study provides a theoretical basis for HCC prognosis, progression, and drug screening.
Tumor production causes alterations or accumulation of surrounding immune or nonimmune cells [50].These immune and nonimmune cells, some noncellular components, and the environment that tumor cells share is called the TME [51].The TME has a dual role in that it can both inhibit tumor cell growth and allow tumor cells to evade surveillance and then promote tumor cell metastasis [52].The TME also influences HCC progression and prognosis [53].Previous studies have demonstrated that METTLs can influence the TME and immune therapy.For example, Liu et al. [54] showed that METTL1 is a key immune regulator in the TME, and targeting METTL1 can inhibit MDSC recruitment and enhance anti-PD-1 efficacy.Xu et al. [55] found that knocking down METTL5 can downregulate the expression of PD-L1, thereby inhibiting immune evasion in HCC.Zheng et al. [56] demonstrated that a model constructed with METTL9 can predict the immune therapy and prognosis of HCC.However, the impact of the remaining eight METTLs on the TME and immune therapy requires further investigation.In the present study, four cell populations were identified by single-cell data, including liver parenchymal cells, monocytes, NK cells, and T cells.In the present study, based on the communication between these four types of cells, we predicted the receptors and ligands between these four types of cells.Previous studies show that T cells and NK cells show that ligand imbalance leads to immune escape against tumor cells [57,58].Furthermore, targeting monocytes can enhance the therapeutic effect of HCC [59].Therefore, a PPI network including METTLs, 11 mRNA risk signature proteins, four cellular receptors and ligands, and

28
Analytical Cellular Pathology potential connexins was established in this study.This PPI network was associated with the communication between these four types of cells.In addition, the METTLs and mRNA risk signature genes were significantly associated with the receptors and ligands of these four cells.Therefore, this PPI network is associated with TME remodeling.In summary, this study performed a multiangle comprehensive analysis of METTLs (METTL1/13/18/21A/23/25/ 2A/2B/5/6/9) and their mRNA risk signature, and these results provide a theoretical basis for HCC prognosis, biomarker screening, mechanism, and drug screening.

3. 10 .
Screening for Sensitive Drugs Shared with METTLs and mRNA Risk Signatures.This study evaluated the correlation of METTLs and mRNA risk profiles with drug IC50 and screened for common-acting sensitive drugs.We selected drugs with significant associations, including 31 drugs in common with METTL1/13/18/21A/23, 48 drugs in common with METTL25/2A/2B/5/6, and 95 drugs in common with METTL9 and the mRNA risk signature (Figure11(a)).All three had 27 drugs in common with significant associations with METTLs and mRNA risk scores (Figure11(b)).Moreover, the IC50 values of these 27 drugs were all significantly inversely

FIGURE 7 :
FIGURE 7: Construction of the mRNA risk signature nomogram prediction model.(a) Nomogram prediction model, (b) calibration curve, (c) the time ROC curves for the particular periods show the AUC values for the individual parameters, and (d) decision curve analysis (DCA).

FIGURE 8 :ðhÞFIGURE 9 :
FIGURE 8: Correlation of mRNA risk signature with clinicopathological features and mutation difference analysis.(a) Differences between mRNA high-and low-risk groups and clinicopathological characteristics, (b) differences between mRNA risk scores and clinicopathological characteristics, and (c) mutational differences in mRNA high-and low-risk groups of the 20 mostly mutated genes in HCC.

4 CorrelationFIGURE 10 :
FIGURE 10: Correlation analysis of METTLs and mRNA risk signature with potential pathways.(a) KEGG terms calculated by GSVA and (b) GSVA calculated hallmark terms.

1 0FIGURE 11 :
FIGURE 11: Screening for sensitive drugs shared with METTLs and mRNA risk signatures.(a and b) Screening for drugs in common with METTLs and mRNA risk profiles; (c) Heat map of correlation between IC50 values of 27 drugs and METTLs and mRNA risk scores.

4
Analytical Cellular Pathology groups were compared, and waterfall plots were plotted to present the results.The nomogram prediction model was based on multivariate Cox regression analysis and was used to predict 1-, 2-, and 3-year OS in patients with HCC.The calibration curve was used to predict survival in terms of the OS of HCC patients.Nomogram prediction model parameters included mRNA risk score, age, grade score, T, and stage.Time ROC curve and decision curve analysis (DCA) were employed to evaluate the accuracy of the nomogram prediction model and other parameters.
(c)).In conclusion, our resulting 11 mRNAs were reliable.Subsequently, we constructed mRNA risk signatures.The risk score was calculated as follows: Risk score (mRNA) = (0005 * PPM1G exp) + (0007 * MEX3A exp) + (0119 * PHOSPHO2 exp) + (0001 * YBX1 exp) + (0036 * UCK2 exp) + (0068 * TAF3 exp) + (0002 * EZH2 exp) + (0031 * PSRC1 exp) + (0009 * TMEM69 exp) + (0003 * CDCA8 exp) + (004 * DYNC1LI1 exp).For the mRNA training set (TCGA) risk signature, the high-risk group showed significantly worse OS than the lowrisk group (P <0:05) (Figure Analysis of mRNA risk signature.(a) Overall survival difference between high-and low-risk group in the mRNA training set, (b) time-dependent ROC curve in the mRNA training group, (c) high-and low-risk score median values and survival status distribution in the mRNA training set, (d) principal component analysis (PCA) and t-SNE scatter plots in the mRNA training set, (e) overall survival difference between high-and low-risk groups in the mRNA validation set, (f ) time-dependent ROC curve in the mRNA validation set, (g) high-and low-risk score median values and survival status distribution in the mRNA validation set, (h) principal component analysis and t-SNE scatter plots in the mRNA validation set, and (i) univariate/multivariate Cox analysis of the training and validation sets.