Ubiquitin-related molecular classification and risk stratification of hepatocellular carcinoma.

The roles of ubiquitin-related genes in hepatocellular carcinoma (HCC) have not been thoroughly investigated. This study aimed to systematically examine ubiquitin-related genes and identify subtypes and stratify prognosis of HCC by using ubiquitin-related signatures. Survival, biological processes, tumor microenvironment (TME), and genomic alterations of the HCC subtypes were investigated. Patients with HCC were classified into two subtypes (clusters 1 and 2) with distinct survival outcomes, pathways, and genomic alterations. Cluster 2 had better prognosis than did cluster 1. Hepatitis B, hepatitis C, Janus tyrosine kinase (JAK)-signal transducer and activator of transcription (STAT) pathway, and natural killer cell-mediated cytotoxicity were enriched in cluster 1. Moreover, cluster 2 had a higher immune score and immune cell infiltrations, whereas cluster 1 had a lower immune score and immune infiltrations. Additionally, mutations, amplifications, and deletions among the phosphatidylinositol 3-kinase (PI3K)-AKT, p53, and receptor tyrosine kinase (RTK)-RAS pathways more frequently occurred in cluster 1, while those among the Hippo, MYC, and Notch signaling pathways were found in cluster 2. Finally, a prognostic signature, consisting of eight ubiquitin-related genes, was established and validated. In brief, our study established a new classification and developed a prognostic signature for HCC.

Introduction

Hepatocellular carcinoma (HCC) is one of the most common malignancies and the second leading cause of cancer-related deaths worldwide. In spite of the advances in treatment and diagnostic methods, the clinical outcome of patients with HCC remains poor. Therefore, screening populations at high risk of developing HCC, discovering new therapeutics targets, and improving prognosis are urgently needed.

Ubiquitin is a small protein and serves as a post-translational protein modifier by marking proteins for degradation. Ubiquitin transfer cascades (known as ubiquitination), including ubiquitin-activating enzymes (E1s), ubiquitin-conjugating enzymes (E2s), and ubiquitin-protein ligases (E3s), constitute a complex network to modify protein substrates. This ubiquitination process can be reversed by deubiquitinating enzymes (DUBs), that is, removing ubiquitin from modified proteins. In addition, many proteins with ubiquitin-binding domains (UBDs) or ubiquitin-like domains (ULDs) also play important roles in the ubiquitin systems. These ubiquitin system-related genes regulate a number of biological processes such as protein degradation, DNA damage repair, signal transduction, and cell cycle.6, 7, 8 Ubiquitin chain dysregulation and ubiquitin-related protein malfunctions are involved in the development of various diseases such as cancers, metabolic diseases, and neurodegenerative diseases.,

Over the past few decades, multiple studies have revealed that ubiquitin-related genes are aberrantly expressed in cancers and regulate many cancer-related genes such as tumor suppressor genes (TSGs: VHL, PTEN, p27, p53, and RB) and oncogenes (EGFR and MYC).9, 10, 11, 12 These deregulated ubiquitin-related genes can cause aberrant activation or inactivation of cancer-associated pathways and play important roles in carcinogenesis., In HCC, several ubiquitin-related proteins interplay with cancer-related proteins and act as oncogenic proteins. For instance, the E3 ligase TRIM25 promotes HCC cell survival by targeting Keap1 for ubiquitination and degradation, activating Nrf2 signaling and reducing reactive oxygen species levels during endoplasmic reticulum stress. TRAF6, another E3 ligase, was reported to interact with histone deacetylase 3 to increase gene expression levels and the protein stability of MYC. Targeting ubiquitin-related genes has been a promising strategy for anticancer drug development. However, the roles of ubiquitin-related genes in HCC have not been thoroughly investigated. A better understanding of the mechanisms and roles of ubiquitin-related genes could provide a basis for understanding HCC and the development of novel therapeutics. Therefore, this study aimed to systematically study ubiquitin-related genes in HCC. Given the high heterogeneity in HCC, our study also intended to identify HCC subtypes from tumor samples and stratify the risk and prognosis of patients with HCC by using ubiquitin-related signatures.

Results

Identification of HCC subtypes based on ubiquitin-related genes

The expression profiles of 777 ubiquitin-related genes from The Cancer Genome Atlas (TCGA) were used for consensus clustering analysis of HCC. The optimal cluster was obtained when the k value was 2 according to the cumulative distribution function (CDF) curves (Figures 1A and 1B). The 371 HCC patients were classified into two subtypes as follows: cluster 1 (n = 250) and cluster 2 (n = 121; Figure 1C). A survival analysis revealed that the overall survival (OS) time and progression-free survival (PFS) time in cluster 1 were shorter than those in cluster 2 (p < 0.05; Figures 1D and 1E). The principal component analysis (PCA) revealed that the samples from cluster 1/2 were well separated from each other (Figure 1F). The results suggested significant differences between cluster 1 and cluster 2.

Figure 1

Identification of two subtypes of HCC in TCGA cohort

(A) Consensus clustering cumulative distribution function (CDF) for k = 2–9. (B) Relative change in area under the CDF curve for k = 2–9. (C) The consensus score matrix of HCC samples when k = 2 (1, cluster 1; 2, cluster 2). (D) Kaplan-Meier curves of overall survival (OS) of clusters 1 and 2. (E) Kaplan-Meier curves of progression-free survival of clusters 1 and 2. (F) Principal component analysis of the gene-expression profiles in the TGGA HCC cohort.

Clinicopathological features, biological processes, and pathways in clusters 1 and 2

The expression levels of the ubiquitin-related genes in clusters 1 and 2 are shown in the heatmap (Figure S1). The clinicopathological features of the two subtypes were compared. We found that the proportion of patients with lower tumor stage (stage I/II) was significantly higher in cluster 2 (p = 0.0015; Table 1). Moreover, the proportions of T1 and T2 were significantly higher in cluster 2 (p = 0.0471; Table 1), while HCC samples with advanced tumor stages (stage III/IV) and higher tumor sizes (T3/T4) were included in cluster 1. There were significantly more patients with hepatitis B and C infections (p = 0.0104) in cluster 1. No significant differences in the distribution of age, sex, histological grade, and N stage were found between the two subtypes (p > 0.05; Table 1).

Table 1

The demographic and clinicopathological characteristics of two clusters in the TCGA HCC cohort

Variables	Group	Cluster 1 (n = 250)	Cluster 2 (n = 121)	p value	Method
Age (mean ± SE)		59.6 ± 0.9	59.2 ± 1.1	0.8296	t-test
Sex	female	87 (34.8%)	34 (28.1%)	0.1968	χ2 test
	male	163 (65.2%)	87 (71.9%)
Clinical stage	I	100 (40.0%)	71 (58.7%)	0.0015	Fisher’s exact test
	II	60 (24.0%)	26 (21.5%)
	III	67 (26.8%)	18 (14.9%)
	IV	2 (0.8%)	3 (2.5%)
	NA	21 (8.4%)	3 (2.5%)
T stage	T1	110 (44.0%)	71 (58.7%)	0.0471	Fisher’s exact test
	T2	68 (27.2%)	26 (21.5%)
	T3	61 (24.4%)	19 (15.7%)
	T4	10 (4.0%)	3 (2.5%)
	TX	1 (0.4%)	2 (1.7%)
N stage	N0	163 (65.2%)	89 (73.6%)	0.1371	Fisher’s exact test
	N1	4 (1.6%)	0 (0.0%)
	NX/NA	83 (33.2%)	32 (26.4%)
M stage	M0	172 (68.8%)	94 (77.7%)	0.0206	Fisher’s exact test
	M1	1 (0.4%)	3 (2.5%)
	MX	77 (30.8%)	24 (19.8%)
HBV/HCV/HBV+HCV	yes	122 (48.8%)	42 (34.7%)	0.0104	χ2 test
	none	128 (51.2%)	79 (65.3%)
Histological grade	G1	36 (14.4%)	19 (15.7%)	0.8140	Fisher’s exact test
	G2	117 (46.8%)	60 (49.6%)
	G3	86 (34.4%)	36 (29.8%)
	G4	7 (2.8%)	5 (4.1%)
	NA	4 (1.6%)	1 (0.8%)

NA, not available; TX, unknown T stage; MX, unknown M stage; NX, unknown N stage; HBV, hepatitis B virus; HCV, hepatitis C virus.

A gene set enrichment analysis (GSEA) was performed to identify pathways enriched in the two subtypes. Under the threshold of an adjusted p value of <0.05, there are 18 significantly enriched pathways in cluster 1 and there were no significant pathways involved in cluster 2. In detail, the results showed that the cytosolic DNA-sensing pathway, hepatitis B, hepatitis C, the Janus tyrosine kinase (JAK)-signal transducer and activator of transcription (STAT) pathway, natural killer cell-mediated cytotoxicity, the neuroactive ligand-receptor interaction, the non-obese diabetic (NOD)-like receptor pathway, pentose and glucuronate interconversions, the retinoic acid-inducible gene (RIG)-I-like receptor pathway, and the Toll-like receptor pathway were highly enriched in cluster 1 (normalized enrichment score > 1, adjusted p < 0.05; Figure S2).

Tumor microenvironment (TME) characterization in clusters 1 and 2

GSEA results indicated that some immune-related pathways, such as natural killer cell-mediated cytotoxicity, the NOD-like/RIG-I-like/Toll-like receptor pathways, and JAK-STAT pathway were significantly activated in cluster 1. Subsequently, the TME (tumor purity and infiltrating stromal/immune cells in tumor samples) in clusters 1 and 2 were examined. We found that the immune score was significantly higher in cluster 2 (p < 0.05; Figure 2A). However, we found no significant differences in stromal score, ESTIMATE score, and tumor purity between the two subtypes (Figures 2B–2D). Given the significant differences in immune score between the two subtypes, immune cell infiltration was further investigated to characterize the immune landscapes of clusters 1 and 2 (Figure 2E). Cluster 2 had a significantly higher abundance of B cells, cytotoxic cells, T cells, dendritic cells, and neutrophils (p < 0.05; Figure 3A). The expression levels of eight immune checkpoint genes (ICGs: PD1, PD-L1, PD-L2, CTLA4, VTCN1, TIM3, LAG3, and TIGIT) in the two subtypes were further investigated. The results indicated that cluster 1 exhibited higher expression levels for the eight ICGs (except for the LAG3 gene) than cluster 2 (p < 0.05; Figure 3B).

Figure 2

Tumor microenvironment characterization in clusters 1 and 2

(A–D) Comparison of immune score, stromal score, ESTIMATE score, and tumor purity in the two subtypes. (E) Heatmap describing the abundance of 24 immune cell types in the two subtypes.

Figure 3

Immune cells abundance and expression levels of immune checkpoint genes in clusters 1 and 2

(A) Comparison of the abundance of 24 immune cell types in clusters 1 and 2. (B) Expression levels of eight immune checkpoint genes in clusters 1 and 2. Level of gene expression is reported as log2-transformed count.

Comparison of mutations and CNVs between clusters 1 and 2

To identify potential drug targets to reverse the poorer survival in cluster 1, we explored the significant mutations and copy number variations (CNVs) between clusters 1 and 2. The top 30 most frequently mutated genes in clusters 1 and 2 are presented in Figures 4A and 4B. Cluster 1 had a significantly higher mutation frequency of TP53 (TSG in the p53 pathway) than did cluster 2 (35% versus 20%). The mutation frequency of CTNNB1 (oncogene in the Wnt pathway) was equal in the two subtypes (both 25%). The mutation frequencies of 10 critical oncogenic pathways in clusters 1 and 2 were summarized (Figure 4C). Mutations in the phosphatidylinositol 3-kinase (PI3K)-AKT, p53, and receptor tyrosine kinase (RTK)-RAS pathways more frequently occurred in cluster 1 (Figure 4C). The Hippo, MYC, and Notch pathways displayed higher mutation frequencies in cluster 2 (Figure 4C). The mutation frequency of the cell-cycle pathway was predominant in both subtypes (Figure 4C). In addition, differences in somatic CNV in clusters 1 and 2 were evaluated using GISTIC 2.0. The CNV analysis revealed that amplifications of 11q13.3 (CCND1 [oncogene in the cell-cycle pathway]), 19q12 (CCNE1 [oncogene in the cell-cycle pathway]), and 7q31.2 (MET [oncogene in the RTK-RAS pathway]), and deletions of 13q14.2 (RB1 [TSG in the cell-cycle pathway]), 9p21.3 (CDKN2A [TSG in the cell-cycle pathway]), 10q23.31 (PTEN [TSG in the PI3K-AKT pathway]), and 9q34.13 (TSC1 [TSG in the PI3K-AKT pathway]) were identified in cluster 1 (Figures 5A and 5B). Moreover, amplifications of 11q13.3 (CCND1 [oncogene in the cell-cycle pathway]), 6p21.1 (VEGFA [oncogene]), and 5p15.33 (TERT [oncogene]), and deletions of 9p21.3 (CDKN2A and CDKN2B [TSGs in the cell-cycle pathway]), 13q14.2 (RB1 [TSG in the cell-cycle pathway]), and 4q21.3 (FAT4 [TSG in the Hippo pathway]) were identified in cluster 2 (Figures 5C and 5D).

Figure 4

Comparison of mutations in clusters 1 and 2

(A and B) The top 30 most frequently mutated genes in two HCC subtypes. (C) The mutation frequencies of ten critical oncogenic pathways in two HCC subtypes.

Figure 5

Comparison of CNVs in clusters 1 and 2

(A–D) GISTIC 2.0 amplifications and deletions in clusters 1 (A and B) and 2 (C and D). Chromosomal locations of peaks of significantly recurring focal amplifications (red) and deletions (blue) are displayed. The q values, representing the statistical significance, are displayed along the bottom. Regions with q values < 0.25 (green lines) were considered significantly altered. The locations of the peak regions of maximal copy-number change and the known cancer-related genes within those peaks are indicated to the right of each panel.

Establishment and evaluation of a ubiquitin-related signature

From 777 ubiquitin-related genes, 336 genes with prognostic value (p < 0.01) were identified by univariate analyses. From the 336 genes, 22 genes were selected by the least absolute shrinkage and selection operator (LASSO) regression analysis (Figure S3A). Finally, we identified eight genes to establish the prognostic signature by the stepwise multivariate regression analysis (Table S1; Figure S3B). The risk score was calculated using the following formula: risk score = (0.1553 × ExpUBE2S) + (−0.3136 × ExpSOCS2) + (0.2454 × ExpRNF2) + (0.3879 × ExpHECTD3) + (0.2298 × ExpATG10) + (0.0634 × ExpBRSK2) + (0.1460 × ExpRNF133) + (0.3940 × ExpTRIM6-TRIM34). According to the median risk score, the patients in the TCGA HCC cohort (training cohort) were assigned to a high- or low-risk subgroup. The survival analysis result indicated that patients in the high-risk subgroup displayed worse OS (p < 0.0001; Figure 6A). The area under the receiver-operating characteristic (ROC) curve (AUC) for OS was 0.826 at 1 year and 0.748 at 3 years (Figure 6B). The robustness and effectiveness of the ubiquitin-related signature were evaluated on the validation cohort. Similarly, the patients in the high-risk subgroup had poorer OS than did those in the low-risk group (p = 0.041; Figure 6C). The 1- and 3-year AUCs were 0.701 and 0.616, respectively, in the International Cancer Genomics Consortium (ICGC) cohort (Figure 6D). We further investigated the correlation between the two risk subgroups and clinicopathological features in the TCGA cohort. Significant differences were found between the high- and low-risk subgroups, which were marked with ∗ in the heatmap (Figure 6E). The high-risk subgroup correlated with advanced T stage (p < 0.01), higher tumor stage (p < 0.001), and higher histological grade (p < 0.01).

Figure 6

Construction and evaluation of a ubiquitin-related signature

(A and C) Kaplan-Meier OS curves for patients assigned to high- and low-risk groups based on the risk score in the TCGA (A) and ICGC (C) cohorts. (B and D) Time-dependent ROC curves were performed in the TCGA (B) and ICGC (D) cohorts. (E) The heatmap of the eight ubiquitin-related genes in low- and high-risk groups. The distribution of clinicopathological features was compared between the low- and high-risk groups. ∗p < 0.05 and ∗∗∗p < 0.001.

Establishment of a nomogram

Univariate analyses were performed to examine the prognostic values of several clinicopathological features. Consequently, tumor stage, T stage, and M stage correlated with OS (p < 0.05; Figure S4A). Tumor and T stages also correlated with PFS (p < 0.05; Figure S4B). The multivariate regression analyses suggested that the risk score was an independent prognostic indicator of OS and PFS in the TCGA cohort (p < 0.001; Figures S4C and S4D). In addition, the ROC analysis revealed that the sensitivity and specificity of the ubiquitin-related signature were greater than those of the other clinicopathological features (Figure 7A).

Figure 7

Integration of the ubiquitin-related signature and clinicopathologic features

(A) ROC curves show the sensitivity and specificity of the ubiquitin-related signature and clinicopathologic features in predicting the OS of HCC patients. (B) Nomogram constructed to predict the 1-, 2-, and 3-year OS in the TCGA cohort. (C and D) Calibration curves of the nomogram for predicting the probability of OS at 1 and 3 years.

A nomogram was constructed, which included risk score, tumor stage, and T/M stage, to predict the probability of OS in HCC patients from the TCGA cohort (Figure 7B). The calibration curves at 1 and 3 years showed good consistency between the actual OS and predicted OS (Figures 7C and 7D).

Discussion

HCC is an aggressive and heterogeneous tumor with a high incidence and a short survival time. Subtype characterization and identification is critical for risk and treatment stratification of HCC patients. In the present study, ubiquitin-related genes were investigated to identify different molecular classifications in HCC samples and to stratify the risk levels of patients with HCC.

Subtype characterization and identification

On the basis of the expression profiles of the ubiquitin-related genes, patients with HCC were classified into two clusters with distinct survival outcomes, pathological features, pathways, TME, and genomic alterations. In detail, cluster 2 had a significantly longer survival time, including OS and PFS, and a higher proportion of patients with early-stage tumors. Hepatitis B and C virus infections are the most common risk factors in the progression of HCC. Of all HCC cases, 75% were associated with hepatitis infection. In our study, we found that cluster 1 with worse prognosis had a higher proportion of patients with hepatitis B and C virus infections, which is consistent with the GSEA results. The GSEA revealed that hepatitis B and hepatitis C virus infection-related pathways were highly enriched in cluster 1. Moreover, the other pathways identified in the GSEA, such as the JAK-STAT pathway and NOD-like/RIG-I-like/Toll-like receptor pathways were frequently reported to play key roles in hepatitis infection or immune processes in HCC.17, 18, 19, 20, 21, 22 For instance, hepatitis B virus infection upregulates Toll-like receptor 2 to promote the invasion of hepatitis B virus-related HCC cells. Therefore, we supposed that clusters 1 and 2 had distinct immune microenvironments. Furthermore, TME characterization in clusters 1 and 2 were investigated to confirm our hypothesis. As a result, cluster 2 had a significantly higher immune score and greater abundance of B cells, cytotoxic cells, T cells, dendritic cells, and neutrophils, while cluster 1 had relatively lower immune cell infiltration and lower immune score. Briefly, cluster 1 had lower immune infiltration and a poor prognosis, while cluster 2 had higher immune infiltration and a favorable prognosis. These findings were consistent with those of a previous study that showed that abundant immune infiltration was associated with better survival in HCC. In addition, cluster 1 exhibited higher expression of ICGs, such as PD-L1, CTLA4, and PD-1 than did cluster 2. In a latest meta-analysis, high PD-L1 expression level in tumor tissue was reported to correlate with shorter OS, poor tumor differentiation, hepatitis, and tumor-infiltrating lymphocytes in HCC. Higher immune checkpoint expression could decrease immune cell infiltration and inhibit the immune response in cancer. Our results (a negative association between immune checkpoint expression and immune cell infiltration) were consistent with these previous studies. The combined effects of lower immune infiltration, high expression levels of ICGs, advanced tumor stage, large tumor size, and hepatitis B and C virus infections may be responsible for the worse survival in cluster 1.

We further investigated the genomic alterations and intended to detect potential drug targets. At the gene level, cluster 1 had a higher mutation frequency of TP53, which leads to p53 pathway activation, than did cluster 2. The mutation frequency of CTNNB1, which leads to Wnt pathway activation, was equal in two subtypes. At the pathway level, the mutation frequency of the cell-cycle pathway is predominant in both subtypes. Mutations in the PI3K-AKT, p53, and RTK-RAS pathways more frequently occurred in cluster 1. The Hippo, MYC, and Notch pathways displayed higher mutation frequencies in cluster 2. The amplifications of oncogenes such as CCND1, CCNE1, and MET and the deletions of TSGs such as RB1, CDKN2A, and PTEN were identified in cluster 1. Therefore, we speculated that the hyperactivated p53, cell-cycle, RTK-RAS, or PI3K-AKT pathway might be responsible for the worse survival in cluster 1. In addition, the amplifications of oncogenes such as CCND1 and VEGFA were identified in cluster 2. All in all, our study demonstrated that certain genomic alterations in these critical pathways may be related to prognosis in clusters 1 and 2. These genes with aberrant mutations, amplifications, or deletions could be considered as effective therapeutics targets for different HCC subtypes. In fact, some corresponding inhibitors that target these key pathways or genes showed promising outcomes in trials, such as MET. Targeting molecular agents that control multiple signaling pathways are also under development.

Risk stratification

A prognostic signature was developed using eight ubiquitin-related genes, which could significantly distinguish the high- and low-risk patients from TCGA and the ICGC cohorts. Among the eight genes, SOCS2 was a protective factor and the rest (UBE2S, RNF2, HECTD3, ATG10, BRSK2, RNF133, and TRIM6-TRIM34) were risk factors. Some of these genes have been reported in HCC. A study demonstrated that UBE2S, as a member of the E2s, is overexpressed in HCC and promotes the progression of HCC cells by enhancing the ubiquitination of p53. RNF2, as an E3 ligase, promotes HCC cells growth and metastasis by targeting SIK1 for degradation. Another study showed that ATG10 rs10514231 might affect the expression of ATG10 and was significantly associated with HCC susceptibility. SOCS2 are reported to inhibit the migration, invasion, and metastasis of HCC cells. These genes may represent promising therapeutic strategies for HCC treatment. The other four genes have not been reported in HCC but in other tumors. Future studies are needed to examine the functional impact of these four genes in the carcinogenesis and development of HCC.

Finally, we integrated risk score, tumor stage, and T/M stage to construct a nomogram, and we found good agreement between the actual OS and predicted OS at 1 and 3 years. These results indicated that the ubiquitin-related signature combined with tumor and T/M stages might be a promising prognostic tool for HCC patients.

Strength and limitations

To the best of my knowledge, our study first identified HCC subtypes in tumor samples and stratified the risk and survival of HCC patients according to ubiquitin-related signatures. The differences in signaling pathways, molecular mechanisms, TME, and genomic alterations between the two subtypes identified suggest that the therapy sensitivity of the two subtypes will be distinct and should be targeted under specific therapeutic strategies. Moreover, a new prognostic signature consisting of eight ubiquitin-related genes was identified and validated. This signature can be used as a screening tool for patients at high risk of developing HCC. The identification of high-risk patients is crucial for early intervention and survival.

This study has a few limitations. First, this was a retrospective study, and all HCC patients were collected from public databases. Second, our HCC samples were smaller than those used in other studies that have integrated multiple datasets and may generate more comprehensive outcomes. Third, a large sample size of patients with HCC from our own hospital is needed for further prospective external validation, and future functional studies are essential to elucidate the precise roles of ubiquitin-related genes in the development of HCC.

Conclusions

Our study established a new classification for HCC based on the expression profiles of ubiquitin-related genes. In addition, we developed a prognostic signature using eight ubiquitin-related genes to screen patients at high risk of developing HCC.

Materials and methods

Data acquisition and processing

In total, 600 patients with HCC were selected from two datasets, including 371 American patients (a training cohort) from TCGA and 229 Japanese patients (an external validation cohort) from the ICGC. The ICGC Liver Cancer-RIKEN, Japan (LIRI-JP) dataset, including transcriptome data and clinical information, was downloaded from the ICGC website (https://dcc.icgc.org). Multi-omics data from the TCGA HCC cohort, including transcriptome data, somatic mutation data, CNV data, and clinical information, was downloaded from the TCGA website (https://portal.gdc.cancer.gov/). By using the maftools R package, somatic mutation data were analyzed, summarized, visualized. For CNV data, GISTIC 2.0 was performed to determine genes with significant amplification or deletion, as based on q values <0.25.

List of ubiquitin-related genes

A list of 807 ubiquitin-related human genes was collected from the iUUCD 2.0 database (http://iuucd.biocuckoo.org/), including 8 E1s, 38 E2s, 501 E3s, 97 DUBs, 121 UBDs, and 42 ULDs. We extracted 777 ubiquitin-related genes with available mRNA expression profiles from the TCGA HCC dataset.

Screening for HCC subtypes

On the basis of the expression profiles of ubiquitin-related genes, molecular subtypes were clustered and identified using the ConsensusClusterPlus R package. The clustering was performed using the following settings: 50 iterations, 80% resampling rate, and Euclidean distance. The optimal cluster number was determined by constructing CDF curves. A PCA was performed to compare the gene-expression patterns among different HCC subtypes. A GSEA using the clusterprofiler R package was performed to investigate the pathways correlated with the different HCC subtypes. The Kyoto Encyclopedia of Genes and Genomes (KEGG) gene sets retrieved from the Molecular Signatures Database were used for the GSEA, and the significance threshold was set at an adjusted p value of <0.05.

TME analysis

To evaluate the heterogeneity in the TME among different HCC subtypes, we inferred the immune and stromal cell components. Immune and stromal scores that can represent the infiltration of TME cells were calculated using ESTIMATE algorithm. Specially, stromal score represents the infiltration of stroma cells in neoplastic samples, immune score represents the presence of immune cells in neoplastic samples, and ESTIMATE score infers tumor purity, which is the proportion of cancer cells in the admixture. In addition, the immune cell infiltration was estimated by performing a single-sample GSEA (ssGSEA) using the GSVA R package. Marker genes for 24 types of immune cells, including 11 innate and 13 adaptive immune cells, were acquired from a published study. Based on the gene-expression profiles and marker genes, infiltrating immune cells were quantified by ssGSEA in individual tumor samples.

Establishment and evaluation of a ubiquitin-related signature and nomogram

Univariate Cox, LASSO-penalized, and stepwise multivariate Cox regression analyses were performed sequentially to establish a ubiquitin-related prognostic signature. Univariate regression analyses were used to identify ubiquitin-related genes with prognostic value. If p < 0.01, the corresponding genes were considered as prognostic genes. LASSO-penalized regression analysis was then performed to screen prognostic genes. Finally, stepwise multivariate regression analysis was conducted to further screen prognostic genes by using the lowest Akaike information criterions (AIC) value and establish the ubiquitin-related signature. The risk score was calculated as follows: risk score = Σ(Ci × Expi), where “C” is the coefficient of gene derived from the multivariate Cox analysis and “Exp” is the gene expression level. We drew ROC curves to compare the prediction efficiency of the signature and several clinicopathological features. Univariate and multivariate Cox regression analyses were performed to investigate prognostic values for the signature and several clinicopathological features. A nomogram was constructed to predict the probability of OS in HCC patients using the rms R package. Calibration curves were drawn to evaluate the effectiveness of the nomogram using the rms R package.

Statistical analyses

A chi-square test was performed to compare the distribution of clinicopathological features, including age, sex, T/N/M pathological stage, clinical stage, and histological grade, between the different groups. Differences between two groups were tested using the Wilcoxon rank test for non-normally distributed variables, and the unpaired t test for normally distributed variables. The patients were categorized into low- and high-risk groups according to median risk score. Survival analyses were performed to compare the survival of the patients in the different subtypes or in the high- and low-risk groups. All statistical analyses were performed using the R software (version 3.5.2) and GraphPad Prism software programs (version 7). All statistical results with a p value of <0.05 were considered significant.