Integrated omics landscape of hepatocellular carcinoma suggests proteomic subtypes for precision therapy

Abstract

Patients with hepatocellular carcinoma (HCC) at the same clinical stage can have extremely different prognoses, and molecular subtyping provides an opportunity for individualized precision treatment. In this study, genomic, transcriptomic, proteomic, and phosphoproteomic profiling of primary tumor tissues and paired para-tumor tissues from HCC patients (N = 160) are integrated. Proteomic profiling identifies three HCC subtypes with different clinical prognosis, which are validated in three publicly available external validation sets. A simplified panel of nine proteins associated with metabolic reprogramming is further identified as a potential subtype-specific biomarker for clinical application. Multi-omics analysis further reveals that three proteomic subtypes have significant differences in genetic alterations, microenvironment dysregulation, kinase-substrate regulatory networks, and therapeutic responses. Patient-derived cell-based drug tests (N = 26) show personalized responses for sorafenib in three proteomic subtypes, which can be predicted by a machine-learning response prediction model. Overall, this study provides a valuable resource for better understanding of HCC subtypes for precision clinical therapy.

Keywords: Sorafenib; hepatocellular carcinoma; machine learning; multi-omics; precision therapy; proteomic subtypes; response prediction model.

Graphical abstract

Figure 1

Proteomic characterization identified three HCC subtypes (A) Consensus clustering of 152 HCC tumors. The associations of HCC proteomic subtypes with clinical characteristics are annotated in the upper panel (chi-squared test, ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001). The heatmap depicts the relative abundance of signature proteins (log₂-transformed). Each column represents a patient sample, and rows indicate proteins. (B) Kaplan-Meier (KM) curves of OS and RFS for each proteomic subtype (log-rank test). (C–E) KM curves of OS for each proteomic subtype in Jing et al.’s cohort (C), Gao et al.’s cohort (D), and Ng et al.’s cohort (E). The p values are calculated by log-rank test. (F) ssGSEA reveals the pathways that are significantly enriched in the three respective proteomic subtypes. The specific enriched signaling pathways of each subtype are analyzed and summarized into four clusters. C1, upregulated in SI only; C2, upregulated in both SI and SII; C3, upregulated in SII and SIII; C4, upregulated in SIII only. (G) Signal pathway changing trend in the four clusters (C1–C4). See also Figures S1 and S2.

Figure 2

The robustness and universality of HCC proteomic subtypes and their simplified discriminating panel (A) Workflow of cross-validation of proteomic signatures in three cohorts (Gao et al.’s cohort: N = 159; Jiang et al.’s cohort: N = 101; present cohort: N = 152). (B) Validation of Jiang et al.’s and Gao et al.’s signatures in our cohort. The left panel shows the KM curves of OS and RFS according to the proteomic subtypes (log-rank test). The right alluvial plot shows the comparison between these subtypes and the original subtypes. (C) Validation of our signatures in Jiang et al.’s and Gao et al.’s cohorts. The confusion matrices are provided. (D) Validation of Jiang et al.’s and Gao et al.’s signatures in each other’s cohorts. The confusion matrices are provided. (E) Concordance rate of three proteomic signatures in three cohorts. (F) Workflow for developing the simplified panel for discriminating proteomic subtypes. (G) Receiver-operating characteristic accuracy, sensitivity, and specificity of simplified panel for discriminating proteomic subtypes in the validation set. (H) Alluvial plot shows the comparison between proteomic subtypes identified by simplified panel and the original subtypes in the validation set. The predictor represents subtypes from the SP9, while the response represents the original subtypes from the full panel. See also Figure S3.

Figure 3

Proteogenomic and immune landscape of three HCC proteomic subtypes (A) Illustration of 152 paired HCC cases used in the individual omics experiments. The omics experiments are colored blue, and the tumor tissues and non-tumor tissues were detected in pairs. (B) WES-based genomic landscape of the three HCC proteomic subtypes. The top panel represents the mutation profile, the middle panel represents TMB and TNB levels, and the bottom panel represents the CNV. (C–E) Feature importance ranking by the random forest algorithm for SI (C), SII (D), and SIII (E). The feature with the highest rank of importance score indicates the highest association (either positive or negative) with the proteomic subtype. (F) KM curves for OS of HCC patients with CTNNB1 mutation or wild type (log-rank test). (G) Heatmap shows the immune cell populations of HCC patients belonging to different proteomic subtypes. (H) Principal component analysis plot of immune scores of immune cell populations based on proteomic data in three proteomic subtypes. (I) Proteome-based immune scores in three proteomic subtypes (two-tailed Wilcoxon test). Box plots show median (central line), upper and lower quartiles (box limits), and 1.5× interquartile range (whiskers). (J) Proteomics-based immune scores of anti-tumor immunity and pro-tumor immune suppression in three proteomic subtypes (two-tailed Wilcoxon test). Box plots show median (central line), upper and lower quartiles (box limits), and 1.5× interquartile range (whiskers). (K) Correlation between anti-tumor immunity and pro-tumor immune suppression based on proteome in three proteomic subtypes. Each point represents a patient sample, the lines represent the fitted curves of correlation in each subtype, and the shaded area represents 95% confidence interval. Pearson’s correlation coefficient (r) and p values are present in the table. The p values were calculated using Pearson’s correlation method. See also Figure S4.

Figure 4

Phosphoproteomic profile and kinase-substrate regulatory network of three HCC proteomic subtypes (A) Summary of the identification of phosphoproteins and kinases. (B) Pathway alterations of phosphoproteins in three proteomic subtypes. (C) Enriched functions of three proteomic subtypes by ssGSEA. (D) Kinase activation enrichment of differentially abundant phosphosites among three proteomic subtypes. Each column represents a patient sample, and each row indicates a kinase. (E) Kinase regulation-pathway network. Rhombus indicates signaling pathways, and roundness indicates kinases differentially expressed in phosphoproteomic data. (F) Kinase-substrate regulation networks in three proteomic subtypes. The edges represent Pearson’s correlation coefficient between kinases and the corresponding phosphosubstrates. (G) Distribution of Pearson’s correlation coefficients of kinase-substrate networks in (F) (two-tailed Wilcoxon test). Box plots show median (central line), upper and lower quartiles (box limits), and 1.5× interquartile range (whiskers). (H) Correlation between kinase abundance and kinase activity. r represents Pearson’s correlation coefficient. The p values were calculated using the Pearson’s correlation method. (I and J) KM curves of OS (I) and RFS (J) for BCKDK activity (log-rank test). See also Figure S5.

Figure 5

Integrated multi-omics analysis of three HCC proteomic subtypes (A) Correlation between WES, transcriptome, proteome, and phosphoproteome. The overlap patients of every two individual omics were analyzed. (B) Spearman’s correlation of CNV to mRNA and protein. Each dot represents a transcript/protein. Attenuated proteins are represented in red using a Gaussian mixture model with two mixture components. (C) Enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of the attenuated proteins. (D) Scatterplots depicting log₂ (T/N) of protein (x axis) and phosphoprotein (y axis) abundance. The red dots indicate positive correlations, and the blue dots indicate negative correlations. The p values were calculated using the Pearson’s correlation method. (E) Enriched KEGG pathways for negative-correlated proteins in (D). (F) Hierarchical clustering analysis map of significantly changed phosphosite-to-protein correlations among three proteomic subtypes. Pearson’s correlation coefficients between matched pairs of phosphosite abundances versus protein abundances were calculated. (G) Functional enrichment for significant phosphosite-to-protein correlations in each cluster. (H) Phosphosite-to-protein co-varying MCODE complexes/subnetworks of significantly changed phosphosite-to-protein correlations among three proteomic subtypes. Top five phosphosite-to-protein co-varying MCODE complexes/subnetworks are shown in different colors. (I–K) Top specific MCODEs in SI (I), SII (J), and SIII (K) positivity are combined with upstream kinases. See also Figure S6.

Figure 6

Drug prediction and key drug target screening for three HCC subtypes (A) Phosphosubstrates of kinases with clinical available drugs and fold change at proteomic and phosphoproteomic levels for kinases and substrates, respectively. The top section displays the abundance of drug-targeting proteins across the three subtypes, with each row representing a drug-targeting protein and “k” labeling representing drug-targeting kinase. The middle section shows the phosphorylation site abundance of substrate for differentially abundance kinases among the three subtypes, with each row representing the phosphorylation sites of substrate. The bottom section displays the abundance of the protein where the phosphorylation site is located, with each row representing a protein. The color gradient represents the abundance of drug-targeting proteins (kinases), phosphorylation sites, and substrate proteins, with green indicating low expression and red indicating high expression. Red labels indicate statistically significant differences in abundance among the three proteomic subtypes. (B) Kinase activity of FDA-approved drug targets in three proteomic subtypes (two-tailed Wilcoxon test). (C) Pathways based on the selected phosphosubstrates and kinases, with relevant drugs shown by targets. (D) Prognostic risk scores of each target from FDA-approved HCC clinical drugs. The x axis indicates log₂-based hazard ratio for each target (log-rank test); y axis indicates log₂-based T/N fold change for each target (two-tailed Wilcoxon test). (E) KM curves of OS and RFS for RAF1. The p values were calculated by log-rank test. (F) Abundance of RAF1 among three proteomic subtypes (two-tailed Wilcoxon test). (G) Correlation between substrate RAF1 and upstream kinase activity among three proteomic subtypes (two-tailed Wilcoxon test). See also Figure S6.

Figure 7

Subtype-specific drug sensitivities and machine-learning-based efficacy prediction model for sorafenib (A) Workflow of pharmacological tests using PDC models. (B) Immunofluorescence of PDCs. Blue represents DAPI staining of the cell nucleus, green represents the distribution of GPC3 on the cell membrane, and red represents the localization of α-fetoprotein (AFP) and albumin (ALB) in the cytoplasm. Scale bar, 20 μm. (C) Sorafenib sensitivity results of the PDC models in three proteomic subtypes (two-tailed Wilcoxon test). Box plots show median (central line), upper and lower quartiles (box limits), and 1.5× interquartile range (whiskers). (D) Sorafenib sensitivity results of the PDC models under different concentrations in three proteomic subtypes (two-tailed Wilcoxon test). Box plots show median (central line), upper and lower quartiles (box limits), and 1.5× interquartile range (whiskers). (E) Proportion of patients who achieved IC₅₀ response in three proteomic subtypes. IC₅₀, half-maximal inhibitory concentration. (F) Correlations between elastic-net-predicted and observed AUC in the validation set. Correlations were calculated by Pearson’s correlation method. (G) Regulation networks of selected protein features in sensitivity prediction model for sorafenib (Pearson’s correlation). (H) Correlations between enriched pathways and observed AUC in PDC models. (I) Enrichment of pathways associated with sorafenib sensitivity in three proteomic subtypes (two-tailed Wilcoxon test). See also Figure S7.