We used UCSCXenaTools (18) to download bulkRNA sequence counts (transcriptomics analysis of pooled cell populations, tissue sections or biopsies) from the TCGA-LIHC cohort, extracting tumor data for 368 patients (Supplementary Table S1). In addition, RNA sequence counts for 232 patients in the LIRI-JP cohort were downloaded from the International Cancer Genome Consortium (ICGC) database (Supplementary Table S2). Finally, data were extracted from the GSE111389 and GSE156625 datasets in the GEO database (https://www.ncbi.nlm.nih.gov/gds/?term=GSE111389[Accession] & https://www.ncbi.nlm.nih.gov/gds/?term=GSE156625[Accession]). GSE111389 contains high-throughput sequencing of tumor-infiltrating lymphocytes (TILs) in HCC, including CD8⁺ T cells expressing no, intermediate or high levels of PD-1 (19). All high-throughput sequencing counts were normalized using the TPM method. Samples with low-frequency counts were excluded. Genes were excluded if their expression levels were below the lower end of the interquartile range (IQR) of expression across all genes in the sample.

HJURP expression was compared between HCC and normal tissues using GEPIA2 (http://gepia2.cancer-pku.cn/). The survival analysis module in this software was used to perform univariable survival analysis. Then a multivariable Cox regression model was built using the “survival” package in R 4.1.0 (20) in order to examine whether HJURP expression was an independent prognostic factor.

The Immuno-Oncology-Biological-Research (IOBR) package in R integrates 6 commonly used algorithms (MCPcounter, TIMER, xCell, CIBERSORT, EPIC and quanTiseq) to separately analyze tumor-infiltrating immune cells (TILs) in the TME, and it draws on 255 gene signatures related to tumors and the TME (21). We calculated scores for TILs and for the TME signature in the TCGA-LIHC cohort, then we compared the scores between 92 samples whose HJURP expression was below the lower end of the IQR and 92 samples whose HJURP expression exceeded the upper end of the IQR.

Using data from the Human Protein Atlas (https://www.proteinatlas.org/), we compared HJURP expression between normal liver tissue and peripheral blood mononuclear cells (PBMCs) at the single-cell level. We also downloaded single-cell transcriptome sequencing data for 57 HCC specimens from the GSE156625 dataset (22), which we subjected to a standard preprocess workflow in Seurat 4.0.4 (23). We filtered out genes that were expressed in fewer than three cells, and we removed cells whose mitochondrial gene expression accounted for >5% of total gene expression as well as cells with “nFeature_RNA” below 200 and more than 5000. With the remaining 78,430 cells, we reduced the dimensionality of the data using principal component analysis (PCA). Based on Elbow analysis, we selected the first 15 principal components in order to cluster cells (Supplementary Figure S1). Marker genes for each cluster were identified using the Wilcoxon test (adjusted p < 0.05), then each cluster was annotated using the marker list (22). We visualized the clustering results using t-distributed stochastic neighbor embedding (TSNE). We used the Wilcoxon test to identify cell clusters strongly expressing HJURP. Gene Set Enrichment Analysis (GSEA) was performed on the cluster of HJURP-expressing malignant hepatocytes.

HJURP expression was compared across CD8⁺ T cells expressing negative, intermediate or high levels of PD-1. The TILs came from six HCC specimens in the GSE111389 dataset and had been sorted using fluorescence-activating cell sorting (FACS).

Among the TCGA samples, we extracted the subset showing HJURP expression at levels higher than 75% of all samples, we excluded low-quality samples and outliers (Supplementary Figure S2A), then we performed weighted gene co-expression network analysis (WGCNA) using the WGCNA package in R (24). A scale-free network with power = 7 was constructed based on SFT.R.sq and mean connectivity (Supplementary Figure S2B). Module eigengenes (MEs) were identified, and those significantly related to HJURP expression (p < 0.05) were analyzed for potential correlation with the TME signature score. For the resulting TME-related module eigengenes, we calculated their gene significance (GS), defined as the absolute value of the correlation between the gene and the trait; eigengene connectivity (KME), defined as the degree of connection between a gene and other genes; and module membership (MM), defined as the correlation between module eigengenes and the gene expression profile. Finally, we screened for hub genes based on criteria of GS > 0.4 and KME >0.8, and we used Cytoscape (25) to visualize co-expression networks with vital edges in each module.

We analyzed hub genes for their enrichment in Gene Ontology (GO) biological processes using the clusterprofiler package (26). We defined enrichment as an adjusted p < 0.05, where adjustment was performed using the Benjamini-Hochberg method. Results were visualized using enrichplot (27).

We assessed the prognostic potential of hub genes using univariable Cox regression. Then we established the following risk score model using multivariable Cox regression. Variables were selected for Cox regression using forward selection based on the likelihood ratio (forward LR). The risk score was developed using the TCGA-LIHC cohort, and externally validated using the ICGC-LIRI-JP cohort. The final risk score was r i s k   s c o r e = ( H J U R P × 0.09523 ) + ( P P T 1 × 0.04176 ) + ( P M L × − 0.04464 ) + ( C L E C 7 A × − 0.1061 )

The prognostic potential of the risk score and other prognostic parameters was assessed based on the area under time-dependent receiver operating characteristic curves (AUC). The four genes in the risk score were used to build a random survival forest model, whose tree number was set to 1000.

In the end, we developed a prognostic model based on the risk score, age, sex, grade, prothrombin time, and tumor stage. The ability of the model to predict survival was assessed using calibration curves and AUCs.

Data were analyzed statistically using R 4.1.0 (http://www.r-project.org) and SPSS 21 (IBM, Armonk, NY, United States). Differences between groups were assessed for significance using the Student-Newman-Keuls test, or the Wilcoxon test as appropriate. Potential correlations were assessed using Pearson coefficients. Overall survival was compared between groups using the Kaplan-Meier method and log-rank test.

Prognostic factors were assessed in terms of hazard ratios (HRs) and their 95% confidence intervals (CIs) based on univariable and multivariable Cox regression within the “survival” package in R. Time-dependent receiver operating characteristic curves and their AUCs were calculated using the “timeROC” package in R. Statistical significance was set at p < 0.05, and all p values were two-tailed.

The GEPIA2 database showed that HJURP expression was significantly higher in HCC than normal tissues (p < 0.05; Figure 1A), and higher expression correlated with lower overall and progression-free survival (p < 0.001; Figures 1B,C). These results suggest that increased expression of HJURP may be a risk factor for poor prognosis.

A multivariable Cox regression model was defined using data on HCC patients from the TCGA-LIHC cohort, HJURP expression and other independent prognostic factors, including age, gender, Child-Pugh liver function classification, and tumor stage (Figure 1D). The HR for HJURP was 1.23 (p < 0.01), which confirmed that HJURP is an independent prognostic factor for HCC.

To further explore whether HJURP is associated with the TME in HCC, we used IOBR to define TME signatures. We found that tumor-related signature scores were significantly higher in samples expressing high HJURP than in samples expressing low HJURP. These signatures included genes involved in cell cycle regulation, cell cycle, DNA damage repair (DDR), mismatch repair, and homogeneous recombination, nature metabolism hypoxia, molecular cancer m6A, exosome, positive regulation of exosomal secretion (Figure 2A). Contrary to the above signatures, ferroptosis scores was significantly lower in samples expressing high HJURP (Figure 2A). Scores for TME-related signatures for myeloid-derived suppressor cells (MDSCs), immune checkpoints, and CD8⁺ T cell depletion were also significantly higher in samples overexpressing HJURP than in those not overexpressing it (Figure 2B). These results link high HJURP expression with high scores for TME-related signatures, which indicates a more immunosuppressive TME.

Next we used six algorithms (MCPcounter, quanTiseq, xCell, CIBERSORT, EPIC and TIMER) to estimate immune cell infiltration in HCC. Infiltration by B cells was significantly higher in samples expressing high HJURP than in samples expressing low HJURP, except for naive B cells as calculated by CIBERSORT and plasma cells as calculated by CIBERSORT and xCell (Figure 3).

Infiltration by follicular helper T cells and regulatory T cells, as calculated by CIBERSORT, was significantly higher in samples expressing high HJURP than in samples expressing low HJURP. Infiltration by CD8⁺ T cells, as calculated by MCPcounter, quanTIseq, and TIMER, were significantly higher in samples expressing high HJURP. Some T cell subtypes showed a tendency toward higher infiltration in samples expressing high HJURP, as calculated by xCell, but this tendency was not as strong as the difference calculated by other algorithms (Figure 3).

Infiltration by dendritic cells (DCs), as calculated by TIMER, MCPcounter, CIBERSORT, and quanTIseq, was significantly higher in samples expressing high HJURP. Different algorithms gave different results for infiltration by macrophages. Infiltration was significantly higher in samples expressing high HJURP in the case of macrophages, as calculated by TIMER; M0 macrophages, as calculated by CIBERSORT; M2 macrophages, as calculated by quanTIseq; and M2 macrophages, as calculated by CIBERSORT. Conversely, infiltration was significnatly lower in samples expressing high HJURP in the case of macrophages, as calculated by xCell and EPIC; and M2 macrophages, as calculated by xCell (Figure 3).

Infiltration by neutrophils did not differ significantly between samples showing low or high HJURP expression, except when analyzed using TIMER. The various algorithms gave conflicting results for monocytes and natural killer cells (Figure 3).

In general, the various algorithms gave consistent results for all types of B cells, some types of T cells, and dendritic cells: their infiltration was higher in samples expressing high HJURP (Figure 3).

Single-cell data from the Human Protein Atlas showed that HJURP was highly expressed in T cells (c-15), erythrocytes(c-16), and B cells (c-6) in normal liver (Figure 4A). Figure 4B shows the expression of marker genes of various cell clusters in the normal liver. In PBMCs, HJURP was highly expressed in plasmablasts, Treg cells, exhausted memory B cells, memory CD4 Th1 T cells, CD8⁺ memory effector T cells and other cells (Figure 4C).

To further explore which cell clusters highly express HJURP in the HCC TME, we processed single-cell RNA sequencing data from HCC tissues using the standard Seurat procedure and obtained 15 cell clusters (Figure 5A). The expression of marker genes in each cell cluster are shown in Supplementary Table S3. HJURP was highly expressed in CD8⁺ T cells, dendritic cells, and hepatocytes (Figures 5B,C; Supplementary Table S4), which verified the correlation between HJURP and tumor-infiltrating immune cells. HJURP was highly expressed in the hepatocytes of HCC tissues, most of which are malignant, but not in hepatocytes of normal liver. These results suggest that the observed overexpression of HJURP in tumor tissues may reflect overexpression in malignant hepatocytes.

We classified the hepatocytes into four clusters and annotated them using a marker list (Figure 6A). HJURP expression was higher in the cluster Hepatocytes4 than in other clusters (Figure 6B). We used GSEA to analyze biological processes occurring in each hepatocyte cluster (Figure 7). The cluster Hepatocytes1 was found to be related mainly in metabolism and immune-related biological processes; Hepatocytes2 and Hepatocytes3, metabolism and processes related to the genetic central dogma; and Hepatocytes4, immune processes and processes related to oxidative respiratory chain regulation. Most biological processes in Hepatocytes4 were related to immune responses, in contrast to other clusters, so we focused on Hepatocytes4 as an important immune-related cluster. GSEA enrichment scores were negative for non-specific and specific immune response biological processes in this cluster, suggesting downregulation of immune-related biological processes (Figure 8).

In order to further explore the relationship between immune molecular characteristics of CD8⁺ T cells and expression of HJURP, we compared the expression of HJURP in CD8⁺ T cells expressing negative, intermediate or high PD-1. HJURP expression was significantly higher in CD8⁺ T cells expressing high levels of PD-1 (p < 0.05; Figure 9). These results implicate HJURP in immune checkpoints.

A gene co-expression network based on data from samples strongly expressing HJURP was built up and classified into 95 module eigengenes (MEs) (Supplementary Figure S2C). HJURP-related module eigengenes were defined as those correlating with HJURP expression, based on a correlation p < 0.05. Among HJURP-related module eigengenes, MEcyan and MEgreenyellow were found to be associated with TME signatures (Supplementary Figure S3), and their expression correlated positively with hepatic tissue inflammation and liver fibrosis (p < 0.01; Supplementary Figure S4). There was also a tendency for MEcyan to be associated with survival time, though the correlation did not achieve significance. Scatter plots showed good correlations of GS with MM and KME (Supplementary Figure S5).

Relationships among MEcyan, MEgreenyellow, and TME signatures are shown in Supplementary Figure S6. A total of 27 hub genes were screened out from MEcyan, but none from MEgreenyellow. Scatter plots and Pearson correlation analyses confirmed correlations between HJURP and all hub genes except TCHH (Supplementary Figure S7). Therefore we excluded TCHH as well as two genes encoding long non-coding RNAs (LINC01094, L3MBTL4−AS1), leaving 24 hub genes for subsequent analysis. Cytoscape was used to depict parts of the network based on MEcyan and MEgreenyellow (Supplementary Figure S8).

GO enrichment analysis of the 24 hub genes confirmed that they had an immune-related function. A cluster tree plot of the top 30 GO bioprocesses showed the main pathway clusters: cellular detection biotic interferon−gamma, neuron death oxidative stress, immune assembly phagocytosis bioprocesses, gland migration morphogenesis duct, and interleukin−1 cysteine−type production endopeptidase (Figures 10A, 11A). Associations were observed between the genes and bioprocesses (Figures 10B, 11B), and the bioprocesses formed a network whose edges connected overlapping gene sets (Figure 10C). Most hub genes identified through single-cell RNA sequencing of HCC tissues were associated with immune cells (Supplementary Figure S9).

Given the prognostic potential of MEcyan, we performed univariable Cox regression on each hub gene. All hub genes met the assumption of proportional hazards (Supplementary Figure S10A). We tested HJURP and all hub genes in multivariable Cox regression models. Variables were selected using likelihood ratio forward selection. We built a significant Cox regression model using four genes: HJURP, PPT1, PML, and CLEC7A. With this model, we defined a risk score (Supplementary Figure S10B). Figures 12A,B show the distribution of risk scores and patient survival in the TCGA and ICGC cohorts. Low-risk groups showed significantly longer overall survival than high-risk groups (p < 0.01; Figures 12C,D). The prognostic risk score also proved to be an independent prognostic factor (p < 0.01; Figure 13).

Using the risk score and TCGA data, we obtained AUCs of 0.702 for 1-year survival, 0.648 for 3-year survival, and 0.661 for 5-year survival (Figure 14A). Next we validated the risk score externally using ICGC, obtaining AUCs of 0.738 for 1-year survival, 0.764 for 3-year survival and 0.783 for 5-year survival (Figure 14B). These results suggest that the risk score predicts survival of HCC patients moderately well. The risk score gave comparably good AUCs as other prognostic parameters in many cases, and it outperformed the other parameters for predicting long-term survival. In addition, the AUCs for the risk score appeared to be more stable than those for other prognostic parameters (Figures 14C,D).

We use the TCGA dataset as the training set to construct a random survival forest model. The random survival forest model constructed by HJURP, PPT1, PML, and CLEC7A has an out of bag (OOB) error of 0.41. Supplementary Figure S11A shows the relationship between OOB errors and the number of trees. We calculated the predicted values of the model for the training set (TCGA), plotted the TimeDependentROC curve and calculated the AUC, in which the 1-year survival AUC was 0.836; the 3-year survival AUC was 0.842; and the 5-year survival AUC was 0.892 (Supplementary Figure S11B). We calculated the predicted values of the model for the test set (LIHC). TimeDependentROC curves were used to assess the predictive performance of the model in the test set. The 1-year survival AUC in the test set was 0.726; the 3-year survival AUC was 0.748; and the 5-year survival AUC was 0.731 (Supplementary Figure S11C). The random survival forest model performed well in the training set, with AUC values above 0.8, outperforming the performance of the risk score based on Cox regression. However, its performance on the test set is significantly poorer than that on the training set, and its performance is slightly poorer than that of the risk score on the test set. Because only four features were used to construct the model, the robustness of the model was decreased, leading to overfitting of the model. Therefore, in addition to the risk score constructed based on Cox regression, random survival forest can also be used as another method to construct prognostic predictors. The other parameters of the random survival forest modfel are shown in Supplementary Table S5.

Finally, we established a prognostic model based on risk score, age, sex, grade, prothrombin time and stage, and we constructed the corresponding nomogram (Figure 15A). Testing of the model against the TCGA cohort gave a C-index of 0.6631 and moderate fit for the calibration curve (Figure 15B). The prognostic model performed better than other prognostic parameters at predicting survival, particularly long-term survival (Supplementary Figure S11).