Integrated Analysis of N1-Methyladenosine Methylation Regulators-Related lncRNAs in Hepatocellular Carcinoma

Simple Summary The relationship between m1A-related lncRNAs and HCC is still unclear. In this study, five m1A-related lncRNAs (AL031985.3, NRAV, WAC-AS1, AC026412.3, and AC099850.4) were identified and used to develop a prognostic signature. The prognostic signature was an independent risk factor related to OS in HCC patients. Synergistic effects on patient survival were observed after combining with TP53 or TMB. In addition, we also screened small molecules which could be potential drugs for HCC patients. Our results suggested that five m1A-related lncRNAs generated a prognostic signature that could be a promising prognostic prediction approach and therapeutic response assessment tool for HCC patients. To the best of our knowledge, this is the first study of m1A-related lncRNAs in HCC. Abstract N1-methyladenosine (m1A) and long non-coding RNAs (lncRNAs) play significant roles in tumor progression in hepatocellular carcinoma (HCC). However, their association with HCC is still unclear. In this study, lncRNAs related to m1A were extracted from the mRNA expression matrix in The Cancer Genome Atlas (TCGA) database. Five m1A-related lncRNAs (AL031985.3, NRAV, WAC-AS1, AC026412.3, and AC099850.4) were identified based on lasso Cox regression and they generated a prognostic signature of HCC. The prognostic signature was identified as an independent prognosis factor in HCC patients. Moreover, the prognostic signature achieved better performance than TP53 mutation status or tumor mutational burden (TMB) scores in the stratification of patient survival. The immune landscape indicated that most immune checkpoint genes and immune cells were distributed differently between both risk groups. A higher IC50 of chemotherapeutics (sorafenib, nilotinib, sunitinib, and gefitinib) was observed in the high-risk group, and a lower IC50 of gemcitabine in the low-risk group, suggesting the potential of the prognostic signature in chemosensitivity. In addition, fifty-five potential small molecular drugs were found based on drug sensitivity and NRAV expression. Together, five m1A-related lncRNAs generated a prognostic signature that could be a promising prognostic prediction approach and therapeutic response assessment tool for HCC patients.


Introduction
Hepatocellular carcinoma (HCC) is the most common hepatic cancer, accounting for about 90% of primary liver cancers. It ranks sixth in incidence and fourth in cancer-related deaths worldwide [1]. However, HCC ranks fifth in incidence and second in cancer-related deaths in China [2]. Currently, the mechanism of HCC pathogenesis is unclear. Chronic hepatic B/C virus infection, alcoholic consumption, and non-alcoholic steatohepatitis are

Data Collection and Processing
The flowchart of our study is shown in Figure 1. TCGA database (https://gdc.cancer. gov, accessed on 1 March 2021) was used to download the RNA-seq expression profile (FPKM values), somatic mutation data, and clinicopathological and survival information. In addition, ten m1A methylation regulators were obtained from the published literature [13,21,22], including four methyltransferases (TRMT6, TRMT61A, TRMT61B, and TRMT10C), two demethylases (ALKBH1 and ALKBH3), and four RNA-binding proteins (YTHDF1, YTHDF2, YTHDF3, and YTHDC1). The lncRNA expression matrixes and these 10 m1A methylation regulators were extracted from the mRNA expression profile using the "limma" R package, and further, 126 m1A-related lncRNAs were screened with a criterion of the absolute value of the correlation coefficient > 0.50 and p < 0.001. Univariate Cox regression analysis was performed with p < 0.05 as a cut-off value to investigate the prognosis values of 126 m1A methylation regulators-related lncRNAs. In addition, the Pan-Cancer Atlas Hub (UCSC Xena, http://xena.ucsc.edu, accessed on 1 March 2021) was used to obtain DNA and RNA stemness scores for subsequent analysis.

Unsupervised Clustering Analysis of m1A-Related lncRNAs
Based on the extracted m1A-related lncRNAs, unsupervised clustering analysis was performed to classify patients through the "ConsensusClusterPlus" R package. Gene set variation analysis (GSVA) was performed using the "c2.cp.kegg.v7.4.symbols" gene set obtained from the MSigDB database to compare enriched functional differences between the two clusters. In addition, the distribution of 23 immune cells between both clusters was also identified using the R packages "GSVA" and "GSEABase". The optimal number of clusters at which the magnitude of the cophenetic correlation coefficient starts to decrease is the k value.

m1AScore Construction
All patients (n = 370) were randomly divided into two cohorts, including the training (n = 186) and testing (n = 184) cohorts, to develop and validate prognostic risk models. The lasso Cox regression analysis is a method to improve prediction accuracy and interpretability of statistical models and realize variable selection and regularization. Our study used lasso regression to screen the most valuable prognostic predictors in the m1A-related lncRNAs and used it to develop an m1AScore. m1AScore = (β lncRNA1 × expression level of lncRNA1) + (β lncRNA2 × expression level of lncRNA2) + · · · + (β lncRNAn × expression level of lncRNAn). The prognostic signature was developed based on the median cut-off of m1AScore, and patients were allocated into high-or low-risk groups.

Predictive Performance of the Prognostic Signature
Clinicopathological risk factors related to overall survival (OS) were investigated among the training, testing, and entire cohorts using univariate and multivariate analyses. In addition, subgroup analysis was used to detect the predictive ability of the prognostic signature in patients with different characteristics. A time-dependent receiver operating characteristic (ROC) curve was applied to investigate the predictive ability on 2-year survival using the "timeROC" R package.
Gene Set Enrichment Analysis (GSEA) was conducted using the gene set "c2.cp.kegg.v7.4.symbols" to determine significantly enriched functional pathways between the high-and low-risk groups, using the R packages "org.Hs.eg.db" and "clusterProfiler". The top five enriched pathways in both groups were displayed.

Generating a Nomogram
After integrating with clinicopathological parameters and m1AScore, a predictive nomogram was built to assess 1-year, 3-year, and 5-year OS using the "rms" R package. Calibration curves were used to identify the predictive accuracy of the established nomogram. In addition, ROC curve and decision curve analyses were used to compare predictive performance between the different clinical parameters and the prognostic signature.

Correlation between Single-Nucleotide Variants and the Prognostic Signature
The distribution of somatic variation between the high-and low-risk groups was investigated, and the top 20 driver genes with the highest mutational frequencies were displayed using the R package "maftool". Subsequently, tumor mutational burden (TMB) was calculated by counting the total non-synonymous mutations, and the median TMB score was used as a cut-off value to differentiate high-and low-TMB patients. Prognostic differences among the high-TMB + high-risk (H-TMB + H-Risk) group, the high-TMB + low-risk (H-TMB + L-Risk) group, the low-TMB + high-risk (L-TMB + H-Risk) group, and the low-TMB + low-risk (L-TMB + L-Risk) group were analyzed. Furthermore, after combining TP53 mutation status with risk levels, all patients were divided into four classes. Similarly, survival analysis among the TP53-mutation + high-risk (TP53-M + H-Risk), TP53mutation + low-risk (TP53-M + L-Risk), TP53-wild type + high-risk (TP53-W + H-Risk), and TP53-wild type + low-risk (TP53-W + L-Risk) groups was performed.

Assessment of the Prognostic Signature in Immune Landscapes
Correlation analysis was used to investigate the interrelationship between 23 immune cells and m1AScore using Spearman's method. In addition, the expression patterns of 29 immune checkpoint genes between both risk groups were also compared, according to the previous study [23].

Therapy Response Assessment of the Prognostic Signature
The "pRRophetic" R package was applied to investigate drug sensitivity between the various risk groups based on the half-maximal inhibitory concentration (IC 50 ) [24]. In addition, HCC patients' immunophenoscores (IPS) were obtained from The Cancer Im- munome Database (TCIA, https://tcia.at/home, accessed on 2 March 2021) and compared between the two risk groups. IPS is a representative gene score related to immunogenicity, comprised of four determinants (MHC molecules, immunomodulators, effector cells, and suppressor cells) [25]; higher IPS indicates improved immunogenicity. The correlation between drug sensitivity and identified m1A-related lncRNAs was analyzed using the CellMiner database (version 2021.2, database 2.7) (https://discover.nci.nih.gov/cellminer/, accessed on 2 March 2021) [26] and Pearson's correlation coefficient analysis. All drugs used for correlation analysis were approved by the Food and Drug Administration or identified by clinical trials.

Specimen Collection
Ten pairs of fresh tumor tissues and their adjacent paratumor tissues from ten surgically resected HCC patients were obtained from Zhongshan Hospital, Fudan University, and stored at -80 • C until RNA extraction. Informed consent forms were collected from all patients. The Ethics Committee of the Zhongshan Hospital, Fudan University, approved the protocol of this study.

Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)
Total RNA was extracted using Trizol reagent (Invitrogen) and reverse transcribed to cDNA using PrimeScript RT Reagent Kit (Takara). SYBR Premix Ex Taq (Takara) was used to perform qRT-PCR according to the manufacturer's instructions. Our primers were designed with reference to the previous study [27]. We import the sequence into Primer3 for primer design. The primer sequences used in this study are listed in Supplementary  Table S1. GAPDH was the endogenous control.

Statistical Analysis
The t-test was used to compare differences among continuous variables, and the Chi-square test or Fisher's exact test was applied to compare differences in categorical variables. The log-rank test and Cox regression analysis were used to detect the prognostic value of these variables, while PCA was used to identify the outstanding performance of clustering information and prognostic signature. The correlation between the m1AScore and the stemness index of the tissue samples containing DNA methylation-based stemness scores (DNAss) and mRNA expression-based stemness scores (RNAss) was analyzed by the Spearman's correlation test. All data were analyzed using R (version 4.1.1), and a two-tailed p-value < 0.05 indicated a significant difference.

Identification of Two Clusters Based on the Co-Expressed lncRNAs
The expression matrixes of ten m1A methylation regulators and all lncRNAs were obtained from the mRNA profile. On correlation analysis, 126 co-expressed lncRNAs were identified as m1A-related lncRNAs. Forty-three co-expressed lncRNAs were classified as prognosis-related genes using univariate Cox regression analysis (Supplementary Table S2). The network-linked lncRNAs and m1A methylation regulators are displayed in Supplementary Figure S1A. Most co-expressed lncRNAs were related to YTHDC1 expression. The expression patterns of these prognosis-related lncRNAs were significantly different between the normal and tumor tissues (Supplementary Figure S1B).
In all tumor samples, two clusters were identified through unsupervised clustering analysis based on the expression matrix of screened lncRNAs (Supplementary Figure S1C). Patients were distinguished visually through PCA analysis (Supplementary Figure S1D). Survival analysis demonstrated that patients in cluster 2 suffered a worse OS (Figure 2A). A heatmap containing clinical parameters, cluster types, and lncRNA expression data was displayed to compare patient characteristics ( Figure 2B). There were significant differences in T stage and pathological grade between cluster 1 and cluster 2. GSVA analysis indicated that pathways such as RNA_DEGRADATION, SPLICEOSOME, NUCLEOTIDE_EXCISION_ REPAIR, and CELL_CYCLE, were enriched in cluster 2; while OLFACTORY_TRANSDUC TION, FOLATE_BIOSYNTHESIS, ARGININE_AND_PROLINE_METABOLISM, and ARACHIDONIC_ACID_METABOLISM were enriched in cluster 1 ( Figure 2C). The immune cell infiltrations were compared between both clusters, and the infiltration levels of activated CD4 T cells and type2 T helper cells were higher in cluster 2 than cluster 1. In contrast, activated CD8 T cells, CD56+ natural killer cells, eosinophils, MDSCs, mast cells, natural killer cells, neutrophils, plasmacytoid dendritic cells, and type1 T helper cells were higher in cluster 1 than cluster 2 ( Figure 2D). analysis based on the expression matrix of screened lncRNAs (Supplementary Figure  S1C). Patients were distinguished visually through PCA analysis (Supplementary Figure  S1D). Survival analysis demonstrated that patients in cluster 2 suffered a worse OS ( Figure  2A). A heatmap containing clinical parameters, cluster types, and lncRNA expression data was displayed to compare patient characteristics ( Figure 2B). There were significant differences in T stage and pathological grade between cluster 1 and cluster 2. GSVA analysis indicated that pathways such as RNA_DEGRADATION, SPLICEOSOME, NUCLEO-TIDE_EXCISION_REPAIR, and CELL_CYCLE, were enriched in cluster 2; while OLFAC-TORY_TRANSDUCTION, FOLATE_BIOSYNTHESIS, ARGININE_AND_PROLINE_ME-TABOLISM, and ARACHIDONIC_ACID_METABOLISM were enriched in cluster 1 (Figure 2C). The immune cell infiltrations were compared between both clusters, and the infiltration levels of activated CD4 T cells and type2 T helper cells were higher in cluster 2 than cluster 1. In contrast, activated CD8 T cells, CD56+ natural killer cells, eosinophils, MDSCs, mast cells, natural killer cells, neutrophils, plasmacytoid dendritic cells, and type1 T helper cells were higher in cluster 1 than cluster 2 ( Figure 2D).

Generating an m1AScore for Prognostic Prediction
A total of 370 HCC patients were allocated randomly to the training (n = 186) and testing (n = 184) cohorts, and no statistical differences were found between them (Supplementary  Table S3). Furthermore, based on the forty-three co-expressed m1A-related lncRNAs, five m1A-related lncRNAs were screened as the most valuable predictors for prognostic assessment to generate an m1AScore using lasso Cox regression (Supplementary Figure S1E,F); the correlation coefficients are shown in Table 1.  Figure 3A-C shows the distribution of the relative m1A-related lncRNAs expression and m1AScore among the training, testing, and entire cohorts, respectively. Patients were distinguished visually through PCA analysis based on the prognostic signature (Supplementary Figure S1G). The Kaplan-Meier curves showed that patients with low risk displayed superior OS compared with high risk in the training, testing, and entire cohorts (log-rank test: p < 0.001, p = 0.003, and p < 0.001, respectively) ( Figure 3D-F). The univariate and multivariate regression analysis also indicated that m1AScore was an independent prognostic indicator in both the training and testing cohorts (p < 0.001) (Supplementary Table S4). In addition, the time-dependent ROC curve analysis showed that the area under curves (AUCs) for 2-year OS were 0.752, 0.716, and 0.718 in the training, testing, and entire cohorts, respectively ( Figure

Function Analysis and Nomogram Construction
A Sankey diagram was performed to analyze the interrelationship between the cluster types and the prognostic signature. It showed that most patients in cluster 1 belonged

Function Analysis and Nomogram Construction
A Sankey diagram was performed to analyze the interrelationship between the cluster types and the prognostic signature. It showed that most patients in cluster 1 belonged to the low-risk group with favorable outcomes (Supplementary Figure S3A), and the m1AScore in cluster 2 was higher than in cluster 1 (Supplementary Figure S3B). GSEA function analysis revealed that patients in the high-risk group were enriched in pathways such as CELL_CYCLE, while metabolism-related pathways were involved in the low-risk group patients. Thus, these findings supported the GSVA results in the above clustering analysis ( Figure 4A,B).   A nomogram was established to assess the 1-, 3-, and 5-year OS ( Figure 4C), and the calibration curve revealed the nomogram's ideal consistency for predicting 1-, 3-, and 5-year survival probability ( Figure 4D). The ROC and decision curves validated the nomogram's predictive performance ( Figure 4E,F). Overall, these results suggested the predictive accuracy of the generated nomogram.

Survival Stratification Based on the Prognostic Signature and Single-Nucleotide Variant
Somatic alterations were found to affect the expression of oncogenes and tumor suppressor genes associated with tumor progression [28]. In this study, mutational frequencies of the top 20 driver genes were compared between the groups, and the results showed that the driver genes' alteration frequencies in the high-risk group were significantly higher than the low-risk group, especially in TP53 (43% vs. 14%) ( Figure 5A,B). Because TP53 mutation is associated with advanced tumor biological features and poor prognosis, high mutation frequencies of TP53 and others could be the potential reasons of dismal prognosis in the high-risk group.
More and more studies have suggested that TMB could be an independent indicator for prognosis and immunotherapy [29,30]. Therefore, the combination of TMB and m1AScore could be useful for prognostic stratification and immunotherapeutic response assessment. The correlation between m1AScore and TMB was insignificant ( Figure 5C). However, their combination could further stratify patient survival ( Figure 5D). The OS of patients with H-TMB + H-Risk or L-TMB + H-Risk group was worse than H-TMB + L-Risk or L-TMB + L-Risk group, suggesting that prognostic signature could differentiate the prognosis in the H-TMB or L-TMB patients. These findings demonstrated synergistic effects of TMB and m1AScore in prognosis prediction.
Similarly, we further detect the synergistic effects of m1AScore and TP53 in prognostic assessment. The prognostic signature could stratify outcomes in patients with TP53-M or TP53-W, while TMB scores failed to differentiate patient survival in high-or low-risk patients ( Figure 5E), suggesting promising prognostic values of the prognostic signature.

The Immune Landscape of m1AScore
To explore the differences in immune landscape between the two groups, we also performed correlation analysis to investigate the interrelationship between the m1AScore and twenty-three immune cells. Activated CD4 T cells, activated dendritic cells, immature dendritic cells, plasmacytoid dendritic cells, regulatory T cells, type17 T helper cells, and type2 T helper cells were found to be positively associated with m1AScore. In contrast, activated CD8 T cells and eosinophils were negatively correlated to m1AScore (Supplementary Figure S3E). Immune checkpoint genes are closely related to the immunotherapy of malignant tumors. In this study, the expression levels of 20 immune checkpoint molecules were significantly higher in the high-risk group than the low-risk group, while FGL1 expression was lower in the high-risk group ( Figure 5F). The correlation between m1AScore and DNAss or RNAss was also investigated, and it was found that there were no strong correlations between them (Supplementary Figure S3C More and more studies have suggested that TMB could be an independent indicator for prognosis and immunotherapy [29,30]. Therefore, the combination of TMB and m1AScore could be useful for prognostic stratification and immunotherapeutic response assessment. The correlation between m1AScore and TMB was insignificant ( Figure 5C).

Therapeutic Response Assessment and Drug Sensitivity
The IC50 is a marker of response to chemotherapeutic drugs of tumor cells. The sensitivities of sorafenib, lapatinib, nilotinib, sunitinib, and gefitinib were enhanced in the high-risk group compared to the low-risk group ( Figure 6A-E). In contrast, the sensitivity of gemcitabine was improved in the low-risk group, suggesting the potential therapeutic value of these drugs in various groups. Moreover, IPS-CTLA(−)-PD1(−), IPS-CTLA(−)-PD1(+), and IPS-CTLA(+)-PD1(−) showed significant differences between both groups ( Figure 6F, p < 0.05), suggesting that the prognostic signature could be a promising predictor for assessing therapeutic drug responses. The RNA expression data and activity of one m1A-related lncRNA, NRAV, were extracted from the CellMiner database. Detailed information of significantly related drugs is listed in Supplementary Table S5. The top eight drugs with the highest correlation between NRAV expression and drug sensitivity are displayed in Figure 6G.

Validation of the Expression Patterns of Five Screened lncRNAs
The expression patterns showed that AL031985.3, NRAV, WAC-AS1, AC026412.3, and AC099850.4 were expressed more in the tumor tissues than normal tissues ( Figure 7A). PCR also validated similar expression patterns using six tumor samples and paired normal samples ( Figure 7B). The Kaplan-Meier curves demonstrated that low expression of AL031985.3, NRAV, WAC-AS1, AC026412.3, and AC099850.4 is associated with better survival in the TCGA database ( Figure 7C). These findings indicated the potential of these genes for diagnosis and prognostic assessment. The RNA expression data and activity of one m1A-related lncRNA, NRAV, were extracted from the CellMiner database. Detailed information of significantly related drugs is listed in Supplementary Table S5. The top eight drugs with the highest correlation between NRAV expression and drug sensitivity are displayed in Figure 6G.

Validation of the Expression Patterns of Five Screened lncRNAs
The expression patterns showed that AL031985.3, NRAV, WAC-AS1, AC026412.3, and AC099850.4 were expressed more in the tumor tissues than normal tissues ( Figure 7A). PCR also validated similar expression patterns using six tumor samples and paired normal samples ( Figure 7B). The Kaplan-Meier curves demonstrated that low expression of AL031985.3, NRAV, WAC-AS1, AC026412.3, and AC099850.4 is associated with better survival in the TCGA database ( Figure 7C). These findings indicated the potential of these genes for diagnosis and prognostic assessment.

Discussion
Dynamic RNA methylation modifications are associated with tumor progression and regulators of m1A methylation are involved in cell apoptosis, death, and carcinogenesis [11,16,17]. A recent study demonstrated that m1A regulatory genes may play a critical role in HCC progression and could be used as biomarkers for diagnosis and prognostic assessment [13]. Previous studies indicated that RNA methylation modifications could mediate lncRNA expression [31,32], and lncRNAs affect RNA methylation modification regulators [33,34]. Wang et al. reported that the prognostic signature based on six m6A/m5C/m1Arelated lncRNAs were associated with survival, immune microenvironment, TMB, and immunotherapy in head and neck squamous cell carcinoma patients [19]. However, the role of m1A methylation regulators-related lncRNAs in HCC remains unclear. To the best of our knowledge, this is the first study of m1A-related lncRNAs in HCC.
In this study, two clusters were identified based on m1A lncRNA expression, and patients in cluster 2 showed worse survival than cluster 1. An m1A-related lncRNA risk model was generated through lasso Cox regression to improve the performance of prognosis prediction. Survival analysis identified that the prognostic signature could discriminate patient outcomes among the training, testing, and entire cohorts. Univariate and multivariate Cox regression analysis demonstrated that m1AScore might be a valuable predictor for HCC patients independent of age, gender, grade, and tumor stage. A nomogram was conducted by integrating clinical parameters and the prognostic signature, and the calibration curves revealed good consistency for 1-, 3-, and 5-year survival probability.

Discussion
Dynamic RNA methylation modifications are associated with tumor progression and regulators of m1A methylation are involved in cell apoptosis, death, and carcinogenesis [11,16,17]. A recent study demonstrated that m1A regulatory genes may play a critical role in HCC progression and could be used as biomarkers for diagnosis and prognostic assessment [13]. Previous studies indicated that RNA methylation modifications could mediate lncRNA expression [31,32], and lncRNAs affect RNA methylation modification regulators [33,34]. Wang et al. reported that the prognostic signature based on six m6A/m5C/m1A-related lncRNAs were associated with survival, immune microenvironment, TMB, and immunotherapy in head and neck squamous cell carcinoma patients [19]. However, the role of m1A methylation regulators-related lncRNAs in HCC remains unclear. To the best of our knowledge, this is the first study of m1A-related lncRNAs in HCC.
In this study, two clusters were identified based on m1A lncRNA expression, and patients in cluster 2 showed worse survival than cluster 1. An m1A-related lncRNA risk model was generated through lasso Cox regression to improve the performance of prognosis prediction. Survival analysis identified that the prognostic signature could discriminate patient outcomes among the training, testing, and entire cohorts. Univariate and multivariate Cox regression analysis demonstrated that m1AScore might be a valuable predictor for HCC patients independent of age, gender, grade, and tumor stage. A nomogram was conducted by integrating clinical parameters and the prognostic signature, and the calibration curves revealed good consistency for 1-, 3-, and 5-year survival probability. These findings thus suggested the reliable performance of established prognostic signatures in HCC patients.
HCC development and progression are associated with gene mutations [39]. Survival analysis demonstrated that the H-TMB + H-Risk group showed the worst survival and the L-TMB + L-Risk group the best prognosis, suggesting a synergistic effect after combining TMB and the prognostic signature. Moreover, the prognostic signature could discriminate OS in the H-TMB or L-TMB patients, while TMB status failed to stratify survival in the patients with low-risk; similar results were observed in the combination of TP53 status and the prognostic signature. These findings indicated that the prognostic signature is better than TMB or TP53 status in predicting patient OS. Moreover, combining the prognostic signature and TMB or TP53 status can improve the prognostic values.
Exploration of response rates of chemotherapeutic drugs to tumor cells is valuable for drug screening. Drug sensitivity analysis demonstrated that sorafenib, nilotinib, sunitinib, and gefitinib had higher sensitivity in the high-risk group, while gemcitabine displayed enhanced sensitivity in the low-risk group. Sorafenib is one of the first-line therapeutic strategies for advanced HCC, and gemcitabine is an effective chemotherapeutic drug for advanced HCC in clinical practice [40,41]. Nilotinib, an orally available receptor tyrosine kinase inhibitor, can induce autophagy in HCC through AMPK activation [42]. Sunitinib is a multi-targeted receptor tyrosine kinase inhibitor, similar to sorafenib. The PRODIGE 16 study showed that TACE plus sunitinib as first-line therapy was feasible for HCC patients when surgical resection was not suitable [43]. Gefitinib, a selective EGFR tyrosine kinase inhibitor, blocks EGFR activity and has an antitumor effect on HCC development in DENexposed rats [44]. The sensitivity of screened lncRNAs to different small molecule drugs was also investigated. Sapitinib was highly correlated with NRAV expression and could induce apoptosis and suppress phospho-EGFR and its downstream pathways [45]. These findings suggested the prospect of targeting NRAV for HCC treatment.
Nevertheless, there are several limitations to this study. First, the m1AScore was generated based on the TCGA database, and although the testing cohort was used to verify its values, it is necessary to examine it in an independent database. Second, the expression patterns of five m1A-lncRNAs were identified by qRT-PCR; however, their function should be studied in vivo and in vitro.

Conclusions
In this study, five m1A-related lncRNAs were identified that generated prognostic signatures and can be used as an independent prognostic risk factor for HCC patients. The prognostic signature revealed its potential performance for diagnosis, prognostic prediction, and assessment of therapeutic responses.  Figure S3. (A) The relationship among the cluster, risk score and survival state displayed by a Sankey diagram. (B) The distribution levels of m1AScore between cluster 1 and cluster 2. Correlation analysis between the m1AScore and DNAss (C) or RNAss (D). (E) Correlation analysis between the m1AScore and 23 immune cells. Table S1. The primers sequences used in the qRT-PCR experiments. Table S2. Prognostic values of co-expressed m1A-related lncRNAs. Table S3. Baseline characteristics of patients in the training, testing and entire cohorts. Table S4. Univariate and Multivariate Cox regression analysis based on different clinical characteristics and overall survival in HCC patients. Table S5. Correlation analysis between NRAV expression and drug sensitivity.
Author Contributions: Conceptualization, K.Z. and J.Z.; methodology, D.S. and X.W.; software, D.S. and X.W.; validation, D.S. and X.W.; formal analysis, D.S. and X.W.; data curation, all authors; writing-original draft preparation, all authors; writing-review and editing, all authors; All authors have read and agreed to the published version of the manuscript.