The total transcriptome sequencing (RNA-seq) data and clinical information of patients with HCC were acquired and downloaded from the International Cancer Genome Consortium (ICGC) portal (https://dcc.icgc.org/projects/LIRI-JP) and The Cancer Genome Atlas Liver Hepatocellular Carcinoma (TCGA-LIHC) dataset (https://tcga-data.nci.nih.gov/tcga/), respectively. Furthermore, the 115-NRGs list was obtained from the National Center for Biotechnology Information, United States National Library of Medicine (https://www.ncbi.nlm.nih.gov/gene/).

The univariate cox analysis was used to screen the NRGs for prognostic value. The cut-off p-value was set at 0.05, and the selected survival-related genes were used for further study. Next, we utilized a LASSO regression model (Alhamzawi and Ali, 2018) for predicting the prognosis of patients with HCC, using the “glmnet” R package and integrating survival time, survival status, and gene expression data. Then, we used the 10-fold cross-validation method to obtain the most optimized model. Finally, eight critical NRGs were included in the LASSO regression model, and the minimum criteria determined the penalty parameter (λ). The risk score was calculated using the formula: risk score = (gene A expression × a) + (gene B expression × b) … + (gene N expression × n). Then, a, b, and n represented the regression coefficients. Based on the risk score median, patients with HCC were divided into the high-risk group and the low-risk group.

To validate the model’s effectiveness, the survival analysis between two risk groups was carried out using the “survminer” R package and the log-rank t test. Next, we used the “pROC” R package to analyze the time-dependent ROC based on survival time, survival status, and risk score of patients with HCC. In the meantime, we analyzed the relationship between different risk scores and follow-up time, the survival status, and the changes in eight critical gene expressions. The results were presented by using online bioinformatic analysis tools such as Sangerbox 3.0 (http://vip.sangerbox.com/home.html).

The ICGC-LIRI-JP dataset was regarded as the internal validation cohort, and the TCGA-LIHC dataset was regarded as the external validation cohort to validate the accuracy and availability of our model. Finally, a clinical trial was carried out in our center, and the obtained HCC samples (clinical cohort) were used for further validation.

The PCA and t-SNE analysis were performed by using “Stats” and “Rtsne” R packages, respectively, to visualize the sample distribution between the high-risk group and the low-risk group. The results were then presented by using the “ggplot2” R package and the website iDEP.95 (http://bioinformatics.sdstate.edu/idep/) (Ge et al., 2018).

The DEGs screening was performed by using the “limma” R package. Briefly, we first carried out a multiple linear regression by using the function “lmFit” and then computed the moderated t-statistics, moderated F-statistics, and log-odds of differential expression by empirical Bayes moderation of the standard errors towards a common value. The DEGs were screened by the threshold p value < 0.05, FDR < 0.05 and fold change > 1.5 for the ICGC-LIRI-JP and TCGA-LIHC datasets, while p value < 0.05 and fold change > 1.5 for the clinical cohort.

To visualize the distribution of DEGs, the map and clusters of DEGs were created using the dimension reduction algorithm t-SNE, and this was created by the website iDEP.95 (http://bioinformatics.sdstate.edu/idep/). In the meantime, the top 20 upregulated and downregulated DEGs heat maps were presented using the “ggplot2” R package.

To understand the potential mechanisms and signal pathways which DEGs participated in, the gene ontology (GO) enrichment analyses and circle plots were performed using the “clusterProfiler” R package (Mi et al., 2019). The minimum number and maximum number of genes in the cluster were 5 and 5,000, respectively. The significantly different GO terms and signal pathways were defined as the threshold p value < 0.05 and FDR < 0.1. The results were visualized by using the “ggplot2” R package.

To understand the potential relationships of DEGs in protein level, the authors constructed a PPI network using the Search Tool for the Retrieval of Interacting Genes (STRING) online database (http://version10.string-db.org/) (Szklarczyk et al., 2015). The interactions with an interaction score <0.9 and the genes that had no direct/indirect interactions with ADRB2 would be hidden.

By integrating risk score, age, sex, race, and TNM stage into the Cox regression model, we evaluated the significance of these factors in predicting the OS of patients with HCC in the TCGA-LIHC dataset. Also, a novel prognostic nomogram was developed to offer a reliable and quantifiable method for predicting patients with HCC survival. The results were presented by using online bioinformatic analysis tools like Sangerbox 3.0 (http://vip.sangerbox.com/home.html).

The connectivity map database (CMAP; https://clue.io/, data version: 1.1.1.2) was used to explore the potential targeted drugs for HCC treatments, which had a high-risk score based on the prognostic model. As a collection of genome-wide transcriptional expression data from cultured human cells treated with bioactive small molecules, the CMAP database can assist researchers in discovering the functional connections between genes, drugs, and diseases through the transitory feature of common gene-expression changes. We could acquire the candidate drugs which resulted in opposite gene changes in HCC by inputting upregulated genes of DEGs and downregulated genes of DEGs into the CMAP. The top six potential drugs in HCC cell lines were listed (ranked by the correlation score).

To further explore the interactions between the candidate drugs and genes, the STITCH (version 5.0) database was used, and the correlated genes were shown in the network (http://stitch.embl.de/) (Szklarczyk et al., 2016).

This was a prospective observational study. The study complied with the Helsinki Declaration and the Consolidated Standards of Reporting Trials (CONSORT) statement, and this was approved by the Renji Hospital Ethics Committee (KY2020-185). The written informed consents were obtained from all patients or authorized family members. The inclusion criteria were as follows: 1) age ≥ 18, 2) primary HCC, and 3) received HCC excision surgery. The patients were excluded if they 1) suffered from multiple metastases, 2) combined with other type of cancers, or 3) the clinical data were missing. Patients with HCC who met the criteria were recruited, and HCC samples were collected in the operation room immediately and stored in −80°C refrigerator. All samples of HCC were confirmed by the pathological diagnosis after surgery.

The HCC samples were sent for RNA-seq to explore expression of NRGs and functional enrichment. The HCC tissue (2 cm × 2 cm) was immediately put into liquid nitrogen for preservation after excision from the patients. Then, the HCC tissue was ground and lysed in TRIzol reagent (Invitrogen, United States ), and the total RNA was extracted for mRNA sequencing. Concentrations and RNA integrity were verified before the library preparation. Library preparation and sequencing were performed by the Biomarker Technologies Corporation, Beijing, China. Sequencing was performed on a HiSeq2500 instrument (Illumina, United States ) with 150 bp paired-end reads. By calculating risk scores, patients were separated into the high-risk group and the low-risk group. Sequentially, DEGs and functions enrichment were proceeded.

Statistical analyses were completed using IBM SPSS Statistics 23.0 (SPSS Inc., Armonk, NY, United States). Quantitative data are presented as the mean ± SD, and categorical variables are presented as frequency (n) or proportion (%). Differences between the two groups were analyzed with an independent samples/paired Student’s t-test. Categorical variables were compared using the χ ² test with the Yates correction or Fisher’s exact test. Survival curves were created using Kaplan–Meier survival analysis with a log-rank t-test. All the statistical tests were two sided with p values < 0.05 being considered statistically significant.

As shown in Figure 1, 203 patients with HCC were identified from the ICGC-LIRI-JP dataset, and they were used as the internal validation cohort. By performing the univariable Cox survival analysis of 115 NRGs in the ICGC-LIRI-JP dataset, nine critical genes (CHRNA3, GABRR2, GRM2, CHRNG, GRIA2, GRM6, GRIN2B, ADRA2C, and GRID2), which were significantly correlated with the OS of the patients with HCC, were identified (Figure 2A). Next, these genes were included in the LASSO regression analysis to establish the prognostic model by integrating survival time, survival status, and gene expression data. Finally, 8 genes were successfully included into the model, and the risk score formula is as follows: risk score = 1.389 × CHRNA3 + 1.065 × GABRR2 + 0.560 × GRM2 + 1.683 × CHRNG + 0.400 × GRIA2 + 1.608 × GRM6–0.606 × GRIN2B + 0.006 × ADRA2C (Figures 2B,C). Every patient got a risk score by integrating the expression level of each gene into the formula. The results in Figure 2D showed that patients with higher risk scores displayed more deaths or shorter survival years. Furthermore, the heat map analysis also indicated that patients with high-risk score had generally high expression levels of the eight critically predictive NRGs, except GRIN2B (Figure 2D). Moreover, the survival curve was consistent with the heat map analysis by showing that patients with higher risk score had a worse OS (p < 0.001, HR = 3.15, Figure 2E). The subsequent ROC analysis suggested that this risk score model could effectively predict patients with HCC survival, with the area under the curve (AUC) of a 4-year survival reaching to 0.81 (Figure 2F).

The 203 patients from the ICGC-LIRI-JP dataset were then divided into two groups based on the median of all their risk scores. Furthermore, the DEGs were screened. As shown in Figure 3A; Supplementary Table S2, 412 downregulated DEGs and 290 upregulated DEGs were identified. In addition, the heat map of the top 20 upregulated and downregulated DEGs between the two groups is shown in Figure 3B. Next, by visualizing the DEGs into the different clusters (Figure 3C), it suggested that the DEGs were mainly involved in the immune system process, lipid metabolic process, organic acid metabolic process, oxoacid metabolic process, etc. Subsequently, GO pathway enrichment was performed for upregulated and downregulated DEGs, respectively. The results revealed that metabolic processes, especially organic acid and lipid metabolic processes, cell motility, and leukocyte activation, were enriched (Figures 3D,E; Supplementary Figure S1).

The TCGA-LIHC dataset was used for the external validation to validate the risk score model’s effectiveness. Three hundred forty seven patients with HCC were identified from the TCGA-LIHC dataset, and each of them got a score calculated with the risk score formula. They were then separated into two groups based on the median of all their risk scores. As shown in Figures 4A,B, patients from the high-risk score group had a worse prognosis compared to those from the low-risk score group (p < 0.001, HR = 1.61). Additionally, the prognostic model exhibited a promising predictive capability, and the AUC of a 6 months survival reached to 0.74 (Figure 4C).

A volcano plot of DEGs showed that 2,320 genes were downregulated and 445 genes were upregulated between the two groups (Figure 4D; Supplementary Table S3). Moreover, the heat map of the top 20 upregulated and downregulated DEGs, respectively, was shown in Figure 4E. DEGs visualization and GO enrichment analysis revealed that similar pathways were enriched compared to the ICGC-LIRI-JP dataset (Figures 4F–H and Supplementary Figure S2), indicating that NRGs might play key roles in acid substances and lipid metabolic processes, cell adhesion, and migration in HCC.

A total of 447 mutual DEGs which were both in the ICGC-LIRI-JP dataset and the TCGA-LIHC dataset were screened (Figure 5A). The PPI network revealed that complex connections among mutual DEGs, with some critical proteins accounting for tumor progression [i.e., Src, matrix metalloproteinase (MMP) family, and chemokine family], are being identified (Figure 5B) (Kim et al., 2009; Kessenbrock et al., 2010). GO enrichment analysis also suggested that mutual DEGs mainly participate in acid substance and lipid metabolic processes, which is consistent with the above results in Figures 3, 4 (Figures 5C–E).

We performed a clinical trial to further validate the model’s effectiveness. From April 2021 to May 2021, nine patients with HCC who met the criteria were recruited, and the collected HCC samples were sent for RNA-seq. The dataset for RNA-seq was supplied as the Supplementary Material. By calculating the risk score, patients were separated into two groups. The clinical characteristics of these patients are shown in Table 1 and Figure 6A. Patients in the high-risk score group exhibited more advanced tumor stage and Child-Pugh stage, larger tumor size, more vascular invasion, and portal vein tumor thrombus (PVTT). PCA and t-SNE analyses also suggested that the gene distribution pattern was significantly different in patients of the two groups (Figure 6B). Nine hundred fifteen DEGs were identified between the two groups and were shown in the volcano plot (Figure 6C; Supplementary Table S4). Furthermore, the heat map of the top 20 upregulated and downregulated DEGs was shown in Figure 6D, respectively. In signaling pathway prediction, top-enriched GO pathways were mainly involved in the cellular metabolic process, response to stress, and immune-related processes, which were consistent with the results both in the ICGC-LIRI-JP dataset and the TCGA-LIHC dataset. Organelle organization, DNA-related pathways, and glycolipid metabolic processes were enriched as well (Figures 6E,F). These results support the above hypothesis that NRGs may play critical roles in cellular metabolic processes (especially organic acid, inorganic acid, and lipid metabolism) and in immune response in HCC.

The PPI network was constructed to show the relationship of DEGs. As shown in Supplementary Figure S3, several critical proteins in tumor development and progression, such as PTEN, MAPK8, CUL1, and RAC3, were enriched.

By integrating risk score, age, sex, race, and TNM stage, we developed a novel prognostic nomogram to establish a reliable and quantifiable method for predicting patients with HCC survival. According to Figure 7A, risk score and TNM stage significantly affected the OS of patients with HCC (both p < 0.05). In addition, the calibration plots of the nomogram also suggested a reliable prediction effect based on the model (Figure 7B).

Potential targeted drugs for HCC treatment were predicted by using the CMAP dataset. The DEGs of the clinical HCC samples were put into the dataset, and the top six candidate drugs which could reverse changes in DEGs were shown in Table 2. Interestingly, all six candidate drugs were tubulin inhibitors, and studies have revealed the importance of tubulin in tumor progression (Gudimchuk and McIntosh, 2021). Two of them were benzimidazole-related, which were mainly used for worm infections. Others were vinca-alkaloid-related, which has been widely used for tumor treatment (Jordan and Wilson, 2004). In addition, interactions between the top 6 drug candidates and proteins were presented in Figure 8. These results suggest that benzimidazole-related drugs and vinca-alkaloid-related drugs may be potential targeted drugs for patients with HCC with high-risk scores.