The RNA-seq data of LIHC patients from The Cancer Genome Atlas (TCGA, http://gdc.cancer.gov/) database and three datasets, including GSE14323 (HCC, n = 55; normal, n = 60), GSE14520 (HCC, n = 225; normal, n = 220) and GSE25097 (HCC, n = 268; normal, n = 289), from Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database was downloaded to analysis the AURKA expression in LIHC patients. Gene expression levels were normalized by Robust Multi-Array Average (RMA).

Differential gene expression analysis of mRNAs was performed based on TCGA database. Limma package of R software (R version 3.6.2) was used to conduct differential gene expression analysis. The adjusted p-value was analyzed to correct for false positive results in TCGA or GTEx. “Adjusted p < 0.05 and Log (Fold Change) >1 or Log (Fold Change) <−1” were regarded as the thresholds of differential expression of mRNAs. The study flowchart was presented in Figure 1.

STRING (version 11.0; http://string.embl.de/) is an open-access biological database that predicts comprehensive interactions of genes at the protein level from multiple organisms including Homo sapiens (Mering et al., 2003). To screen AURKA-related genes, we used STRING to explore AURKA-related genes as well as conduct PPI network analysis on AURKA-related genes. The protein–protein interactions (PPI) with medium confidence >0.4 were regarded as significant.

To further confirm the underlying function of potential targets, the data were analyzed by functional enrichment. Gene Ontology (GO) is a widely-used tool for annotating genes with functions, especially molecular function (MF), biological pathways (BP), and cellular components (CC). Kyoto Encyclopedia of Genes and Genomes (KEGG) Enrichment Analysis is a practical resource for analytical study of gene functions and associated high-level genome functional information. ClusterProfiler package (version: 3.18.0) in R was employed to analyze the GO function of potential targets and enrich the KEGG pathway. p < 0.05 was set as the cut-off criterion.

Gene Expression Profiling Interactive Analysis (GEPIA) is an interactive web application based on The Cancer Genome Atlas (TCGA) and Genotype-Tissue Expression databases (Tang et al., 2017). We used the Gene Expression Profiling Interactive Analysis to screen the prognostic value of AURKA and AURKA-related genes. The identification of prognostic genes was based on the following criteria: 1) significant differences in gene expression levels between normal liver samples and liver tumor samples and 2) significant association of the gene with both overall survival (OS) and disease-free survival (DFS) of patients with liver hepatocellular carcinoma (LIHC). p-values <0.05 denoted statistically significant differences.

Primary screening of the LIHC data from TCGA database was performed for missing data. After deleting the samples with missing data, The available LIHC data from the TCGA database (340 samples) was divided into two subsets: a training set (n = 240) and a validation set (n = 100) randomly using the random function of Microsoft Excel program. The training set was used to train a predictive model and the test set was applied for validation. Another database International Cancer Genome Consortium (ICGC) (dcc.icgc.org) was also used as validation set.

Identified prognostic genes were submitted to the multivariate Cox regression model to calculate each prognostic gene’s coefficient and risk scores. We used X-tile plot (version 3.6.1, http://www.tissuearray.org/rimmlab) to determine the optimum cutoff of AURKA-related gene signature risk score. X-tile which could calculate the best cut-point of sub-populations, is a software developing by team Rimm Laboratory from Yale University (Camp et al., 2004). A train set and two validation sets were divided into three sub-groups using the same best cut-off value. Risk score analysis, including risk score distribution, survival status, and gene expression heatmap among sub-groups were performed. Kaplan–Meier curves were analyzed using Kruskal-Wallis test and visualized through GraphPad Prism (version 8.0). Receiver operating characteristic (ROC) curve analysis was performed to assess the predictive value of AURKA-related gene signature. The clinical characteristics of three cohorts were shown in Table 1. The correlation among clinicopathological characteristics and risk groups were analyzed by the chi-square test. For sample sizes of less than 40 or theoretical frequencies (T) of less than 1, Fisher exact probability method was used.

To identify the somatic mutations of patients with LIHC, we downloaded single-nucleotide polymorphism (SNV) data and clinical follow-up information from TCGA database. The downloaded single-nucleotide polymorphism data were organized in the multiple alignment (MAF) format and visualized using the “maftools” package in R software. The horizontal histogram showed the genes with the highest frequency of mutation.

To explore the associations between different subgroups and immune cells infiltration, we employed Tumor Immune Estimation Resource (TIMER), which is a useful resource for comprehensive analysis of tumor-infiltrating immune cells (Li et al., 2017). The infiltration of six type of immune cell, including B cell, Macrophage, Myeloid dendritic cell, Neutrophil, T cell CD4⁺ and T cell CD8⁺, were calculated. SIGLEC15, IDO1, CD274, HAVCR2, PDCD1, CTLA4, LAG3 and PDCD1LG2 were selected to be immune checkpoints and the expression values of these eight immune checkpoints among sub-groups were explored. Differences between the three groups were assessed using the Kruskal-Wallis test. p < 0.05 was considered statistically significant. All the above analysis methods and R package were implemented by R foundation for statistical computing (2020) version 4.0.3 and software packages ggplot2 and pheatmap.

A nomogram model was constructed to predict the probability of survival at 3 and 5 years for liver cancer patients. Briefly, the prognostic value of the clinicopathological characteristics for OS was estimated through univariate and multivariate Cox regression analyses in both training set and validation set. The performance of the risk model was validated by internal validation and external validation. Internal validation was performed by bootstrap Cox proportional regression analysis based on 1,000 bootstrap samples. Validation set was conducted based on another HCC patients from the TCGA database. Those parameters with p-values <0.05 in both training set and validation set were identified as potential prognostic factors, which were were included in multivariate Cox regression model and visualized using R package “rms.” The Calibration curves were plotted to analyze the diagnostic performance of the nomogram. The ROC curve were conducted to determine the clinical value of the nomogram.

The results of the DEG analysis are shown in Figure 2. In TCGA database, a total of 421 HCC samples were selected for the DEG analysis. The analysis showed that the AURKA gene was upregulated (log₂ fold change >1.5 and adjusted p-value <0.05), which indicated that the mRNA expression levels of AURKA differ significantly between normal liver tissue and liver cancer tissue. The mRNA expression levels of AURKA were also determined using three Gene Expression Omnibus data sets (GSE14323, n = 115; GSE14520, n = 445; and GSE25097, n = 557). The results also showed that AURKA was significantly highly expressed in HCC.

According to the predictive results of the Search Tool for the Retrieval of Interacting Genes (STRING) database, another ten genes were identified as AURKA-related genes with significant interaction, namely targeting protein for Xklp2 (TPX2), cyclin dependent kinase 1 (CDK1), polo like kinase 1 (PLK1), DLG associated protein 5 (DLGAP5), cell division cycle 20 (CDC20), baculoviral IAP repeat containing 5 (BIRC5), transforming acidic coiled-coil containing protein 3 (TACC3), centromere protein A (CENPA), cyclin B2 (CCNB2), and ubiquitin conjugating enzyme E2 C (UBE2C). The protein-protein interaction (PPI) network of AURKA and AURKA-related genes was constructed and visualized using the online STRING database (Figure 3A).

The Gene Ontology (GO) enrichment analysis was composed of three parts: GO biological process (GO-BP), GO cellular component (GO-CC), and GO molecular function (GO-MF). In GO-CC (Figure 3B), these genes were enriched in condensed nuclear chromosome, centromeric region, mitotic spindle, and spindle pole. In GO-BP (Figure 3C), AURKA and its related genes were significantly enriched in mitotic nuclear envelope disassembly, positive regulation of ubiquitin protein ligase activity, and anaphase-promoting complex-dependent catabolic process. In GO-MF (Figure 3D), genes were mainly enriched in protein serine/threonine kinase activity, protein heterodimerization activity, and protein kinase binding. For the Kyoto Encyclopedia of Genes and Genomes pathway analysis, nine pathways (Figure 3E) were observed, namely progesterone-mediated oocyte maturation, oocyte meiosis, cell cycle, TP53 signaling pathway, forkhead box O (FOXO) signaling pathway, ubiquitin-mediated proteolysis, cellular senescence, viral carcinogenesis, and human T-lymphotropic virus type I (HTLV−I) infection.

The results of the identification of prognostic genes are shown in Figure 4. According to the screening strategy and criteria described above, two AURKA-related genes, namely CCNB2 (OS Log-rank p > 0.05) and UBE2C (OS Log-rank p > 0.05), were excluded. Eight other AURKA-related genes were significantly highly expressed in HCC tissue compared with normal liver tissue and correlated with prognosis.

The AURKA and the other eight AURKA-related genes were used to construct a risk model. The risk score was calculated through multivariate regression analyses. The cutoff value of the risk score was identified using the X-tile software. According to the results, the HCC samples were divided into three subgroups. The cutoff values of the risk score were 0.17 and 0.95. HCC samples with a risk score <0.17 and >0.95 classified into the low- and high-risk groups, respectively. The remaining samples were assigned to the moderate-risk group.

The details of the risk model showed in Figure 5, which revealed that higher risk scores were associated with higher expression levels of the nine-gene signature. Furthermore, higher risk scores also indicated worse OS. These results were similar in both the training and other two validation sets. Collectively, these results suggest that the risk model had potential value in predicting the prognosis of HCC. Furthermore, as the ROC curves showed, the area under the ROC curve (AUC) of the training set and two validation sets were higher than 0.5, which indicated that AURKA-related gene signature risk model had important value in predicting prognosis (Figures 5D–F).

The clinicopathological characteristics among the risk groups in training set were shown in Table 2. Age, T stage, pathological TNM (pTNM) stage, grade, recurrence, and survival endpoint were significantly different between the three risk groups (p < 0.05).

Figure 6 illustrates the somatic landscape of the three risk subgroups. Information on the mutation status of each gene in each sample was shown in the waterfall plot, where different colors with specific annotations at the bottom indicated the various types of mutation. The barplot above the legend exhibited the number of mutations. The results showed that CTNNB1 was the most commonly mutated gene in the low-risk group. Tumor protein p53 (TP53) was the most frequently mutated gene in both the moderate- and high-risk groups. Hence, we further grouped the HCC samples into two groups based on the TP53 mutation status. A total of 101 and 241 HCC samples were assigned to the TP53 mutant- and wild-type groups, respectively. The survival analyses of both the TP53 mutant- and wild type cohorts yielded similar results. Higher risk scores were associated with worse prognostic outcome (p < 0.05).

We further investigated differences in the degree of immune infiltration in various risk groups. Using the TIMER database, we evaluated the immune cell infiltration in samples from the three aforementioned groups. As shown in Figure 7, we found that the level of immune cell infiltration (macrophages, myeloid dendritic cells, neutrophils, and CD4T cells) was significantly different between the three groups. To further investigate the levels of immune cell infiltration on the gene level, the following eight immune-related genes were selected: sialic acid binding Ig like lectin 15 (SIGLEC15), T cell immunoreceptor with Ig and ITIM domains (TIGIT), CD274, hepatitis A virus cellular receptor 2 (HAVCR2), programmed cell death 1 (PDCD1), cytotoxic T-lymphocyte associated protein 4 (CTLA4), lymphocyte activating 3 (LAG3), and programmed cell death 1 ligand 2 (PDCD1LG2). We compared the expression levels of immune-related genes among the three risk groups. The results are shown in Figure 8. Except for CTLA4, LAG3, and PDCD1, the expression levels of the other five immune-related genes differed significantly between the three groups. Overall, the high-risk group showed significantly higher levels of immune gene expression and immune cell infiltration compared with the other groups.

As shown in Table 3, we performed the univariate and multivariate analyses using the SPSS software (version 23.0; IBM Corporation, Armonk, NY, United States) to identify independent prognostic factors predicting OS in patients with HCC. In the training data set, the risk score, T stage, M stage, pTNM stage, and recurrence were identified as independent prognostic factors. The results of the testing and training data sets were similar. Overall, the risk score, T stage, and pTNM stage were identified as independent prognostic factors in both data sets.

Hence, we used age, T stage, N stage, M stage, pTNM stage, recurrence, and risk group as estimated factors in the construction of our model. A nomogram was constructed to estimate the probabilities for three- and 5-year survival. Calibration curves analysis showed that the new nomogram model had good predictive accuracy (Figure 9B). Furthermore, the new nomogram model achieved an area under the curve (AUC) of 0.78 at 5 years, which was better than that of a model without the gene signature (AUC = 0.73) or pTNM stage (AUC = 0.70) (Figure 9C).