Comprehensive molecular classification predicted microenvironment profiles and therapy response for HCC

Abstract

Background and aims: Tumor microenvironment (TME) heterogeneity leads to a discrepancy in survival prognosis and clinical treatment response for patients with HCC. The clinical applications of documented molecular subtypes are constrained by several issues.

Approach and results: We integrated 3 single-cell data sets to describe the TME landscape and identified 6 prognosis-related cell subclusters. Unsupervised clustering of subcluster-specific markers was performed to generate transcriptomic subtypes. The predictive value of these molecular subtypes for prognosis and treatment response was explored in multiple external HCC cohorts and the Xiangya HCC cohort. TME features were estimated using single-cell immune repertoire sequencing, mass cytometry, and multiplex immunofluorescence. The prognosis-related score was constructed based on a machine-learning algorithm. Comprehensive single-cell analysis described TME heterogeneity in HCC. The 5 transcriptomic subtypes possessed different clinical prognoses, stemness characteristics, immune landscapes, and therapeutic responses. Class 1 exhibited an inflamed phenotype with better clinical outcomes, while classes 2 and 4 were characterized by a lack of T-cell infiltration. Classes 5 and 3 indicated an inhibitory tumor immune microenvironment. Analysis of multiple therapeutic cohorts suggested that classes 5 and 3 were sensitive to immune checkpoint blockade and targeted therapy, whereas classes 1 and 2 were more responsive to transcatheter arterial chemoembolization treatment. Class 4 displayed resistance to all conventional HCC therapies. Four potential therapeutic agents and 4 targets were further identified for high prognosis-related score patients with HCC.

Conclusions: Our study generated a clinically valid molecular classification to guide precision medicine in patients with HCC.

Graphical abstract

FIGURE 1

Single-cell transcriptomic analysis identified distinct cell subclusters. (A) t-SNE visualization of 248,478 cells annotated in 7 major cell types and 40 subclusters. Cells from different subclusters are marked by color code. (B) The tumor cells were divided into 5 CNV subclusters by applying the k-means clustering based on inferred single-cell CNV profiles. (C) Heatmap of marker gene expression in each cell subcluster. (D) Barplot displaying the number of cells, tissue type (adjacent normal or tumor), and the fraction of cells in each cell subcluster from left to right. Abbreviations: CNV, copy number variation; t-SNE, distributed stochastic neighbor embedding.

FIGURE 2

Identification of prognosis-related subclusters and transcriptomic subtypes in HCC. (A) Flowchart of screening prognosis-related subclusters and establishing transcriptomic classes in HCC. (B) Univariate Cox regression of subcluster percentages and subcluster-specific scores in the training cohort. (C) Consensus clustering identified 5 transcriptomic subtypes for the training cohort. (D) Kaplan-Meier curves of 5 transcriptomic subtypes in the training cohort. (E) Distribution of prognosis-related subclusters in 5 subtypes. (F) Multilass AUCs of the internal training set and validation set based on 8 individual machine-learning algorithms. (G) Kaplan-Meier curves comparing the OS of 5 subtypes in the microarray validation cohort. (H) Kaplan-Meier curves showed the OS of 5 subtypes in the RNA-seq validation cohort. (I) Multiplex immunofluorescence staining of 2 CNV_2 SSMs (TFF2 and AGR2). (J) Multiplex immunofluorescence staining of 2 FB_3 SSMs (THBS4 and ACTG2). Abbreviations: CNV, copy number variation; OS, overall survival; SSM, subcluster-specific marker; SSS, subcluster-specific score.

FIGURE 3

Correlation between transcriptomic subtypes and cancer stemness. (A) t-SNE visualization of CNV subclusters in tumor cells. (B) Hierarchical clustering of NMF programs identified 5 prominent metaprograms in epithelial subclusters. (C) The expression of HCC stem cell markers CD133, EpCAM, and KRT19 in 5 CNV subclusters. (D) Violin plot displaying the stemness scores of cancer stem cell (Palmer, 2012) and Stem.Sig (Zhang, 2022) across CNV subclusters. (E) t-SNE plot showing the regulon activity of TP53 and MYC. (F) Differentiation trajectory of CNV subclusters (left) displaying pseudotime (middle) and CytoTRACE score (right). Colored dots indicate cells annotated in (A). (G) The proportion of CNV subclusters in distinct transcriptomic classes. (H) Stemness scores of cancer stem cells (Palmer, 2012) and Stem.Sig (Zhang, 2022) for 5 transcriptomic classes. (I) Specific activated or suppressed regulons for 5 transcriptomic classes and the enrichment pathways of activated regulons of class 5 (right). (J) The hub PPI network of stemness-related regulons in class 5. (K) The activities of PITX2 and MYBL2 regulons in CNV subclusters (left) and HCC transcriptomic classes (right). (L) Multiplex immunofluorescence staining of stemness regulons including TP53, EPCAM, PITX2, and MYBL2 in tumor cells of HCC (marked Glypican3). Scale bars=50 µm. (M) Sphere-forming assays of Hep3B and HCCLM3 cells after PITX2 and MYBL2 knockdown. Representative images are shown, and the number of spheroids is quantified in bar graphs. Scale bars=200 µm. **p<0.01, ***p<0.001. Abbreviations: CNV, copy number variation; NMF, non-negative matrix factorization; PPI, protein-protein interaction; t-SNE, t-distributed stochastic neighbor embedding.

FIGURE 4

Immune characteristics of HCC transcriptomic subtypes in the training cohort and the Xiangya real-world cohort. (A) Comparing the infiltration level of 22 immune cell types in 5 distinct transcriptomic classes. (B) Violin plot of immune score and T-cell exclusion score in 5 transcriptomic classes. (C) Violin plot of the expression of inhibitory immune checkpoint genes PDCD1, CD274, CTLA4, and TIGIT in 5 transcriptomic classes. (D) Activities of different cancer immunity cycle steps in transcriptomic subtypes. (E) Heatmap showing the overall status, BCLC stage, cirrhosis, vascular invasion, ICB response, targeted therapy response, TACE response, and the expression of 15 immune checkpoint genes in the transcriptomic classes of the Xiangya HCC cohort. (F) Overall response rate of 5 transcriptomic classes after ICB therapy in the Xiangya real-world cohort (chi-square test, p=7.6e−06). (G) Kaplan-Meier survival curves displaying the OS (log-rank test, p=0.0025) and PFS (log-rank, p=0.11) of 5 transcriptomic subtypes in the Xiangya real-world cohort. (H) Arterial-phase MRI image displayed the differences in individual immunotherapy responses. The red arrow points to the tumor area. (I) Multiplex immunofluorescence staining in different transcriptomic classes with magnifications of T-helper cells (CD3+; CD8− and FOXP3−), cytotoxic T lymphocytes (CTL; CD3+, CD8+, and FOXP3−), and regulatory T cells (Treg; CD3+, CD8−, and FOXP3+). Scale bar=100 µm. (J) The heatmap showed the association between the immune cell infiltration percentage and the transcriptomic classes among the Xiangya HCC cohort. Annotation of immune cells based on canonical markers. (K) Gating strategy for identifying CD45 immune cell populations. (L) Dimplot of transcriptomic subtypes and cell clusters based on CyTOF. (M) A heatmap showing the differential expression of 29 immune markers in the 9 cell clusters. (N) The distribution of immune checkpoint CTLA4 in different cell clusters. (O) The proportion of CD45 immune cells in the 5 HCC transcriptomic subtypes (ANOVA, p=0.0062). (P) Frequency diagram of immune cell subclusters in each patient with different transcriptomic subtypes. Abbreviations: BCLC, Barcelona Clinic Liver Cancer; ICB, immune checkpoint blockade; OS, overall survival; TACE, transarterial chemoembolization.

FIGURE 5

Single-cell RNA-seq and scTCR/BCR-seq revealed the immune repertoire of transcriptomic subtypes. (A) UMAP visualization annotated with 10 major cell types in classes 1, 3, 4, and 5. (B) The proportion of 10 major cell types in different transcriptomic subtypes. (C) The expression of immunosuppressive markers (FOLR2, SPP1, and SELENOP) and MHC class II genes (HLA-DRB1, HLA-DQA1, and HLA-DQB1) in myeloid cells across identified transcriptomic subtypes. (D) UMAP visualization of T-cell phenotypes and steady-state RNA velocity of T-cell phenotypes. (E) RNA velocity and expression of FOXP3 and CTLA4. (F) Violin plot displaying the expression of canonical markers in distinct T-cell phenotypes. (G) UMAP plot annotated with TCR clonotypes. (H) The proportion of 10 TCR clonotypes in different transcriptomic subtypes. (I) TCR expansion frequency in distinct T-cell phenotypes. (J) Resident, exhausted, cytotoxicity, and costimulatory scores in T cells of different transcriptomic subtypes. (K) TCR expanding frequency in transcriptomic subtypes. (L) The percentage of T-cell phenotypes in transcriptomic subtypes. (M) The expression pattern of precursor-exhausted CD8 T cells (Ly108+CD69−) and terminal exhausted CD8 T cells (Ly108−CD69+) in different transcriptomic subtypes. (N) Multiplex immunofluorescence staining of CD3/CD8/FOXP3/PD1/PD-L1 in different transcriptomic subtypes. Scale bars=50 µm. Abbreviations: BCR, B cell receptor; MHC, major histocompatibility complex; TCR, T cell receptor; UMAP, uniform manifold approximation and projection.

FIGURE 6

Construction of PRS based on machine learning. (A) Flowchart of PRS construction. (B) Survival curves of different PRS groups (top) and time-dependent ROC curves (bottom) of PRS in the combined training cohort, microarray validation cohort, RNA-seq validation cohort, and Xiangya HCC cohort. (C) Sankey diagram displaying the correlation of transcriptomic subtypes, PRS groups, and OS status. (D) ROC curves of PRS and clinical characteristics for predicting OS in the TCGA-LIHC, CHCC-HBV, and LIRI-JP cohorts. (E) Forest plots showing the HR of PRS and clinical characteristics using both univariate and multivariate Cox regression in the TCGA-LIHC, CHCC-HBV, and LIRI-JP cohorts. (F) Correlation between PRS and activities of specific hallmark pathways, PROGENy pathways, stemness-related regulons, and cancer stemness scores. (G) The association between PRS and clinical characteristics in the TCGA-LIHC cohort. (H) The association between PRS and clinical characteristics in the CHCC-HBV cohort. The statistical significance of the difference was determined using the Kruskal-Wallis test. (I) The correlation analysis of PRS and immune checkpoints. (J) The expression of PDCD1 and CTLA4 in high PRS and low PRS groups. (K) The proportion of CD45+ cells in high PRS and low PRS groups. (L) PRS in immunotherapy responsive and nonresponsive groups of the Xiangya HCC cohort. (M) Multiplex immunofluorescence staining of PD1/PD-L1/CTLA4 in high and low PRS groups. Scale bars=20 µm. Abbreviations: OS, overall survival; PRS, prognosis-related score.

FIGURE 7

Identification of potential drugs, therapeutic targets, and pathways for HCC transcriptomic subtypes and PRS. (A) The response rate of 5 transcriptomic classes after sorafenib treatment (chi-square test, p=7.8e−7). (B) The response rate of 5 transcriptomic classes after TACE treatment (chi-square test, p=7e−6). (C) Violin plot showing predicted IC50 value of sorafenib in 5 transcriptomic classes based on CGP2016 (Kruskal-Wallis, p=1.9e−14) and CTRP2 database (Kruskal-Wallis, p<2.2e−16). (D) Spearman correlation analysis of PRS and 6 CGP2016-predicted compounds (top) and 7 CTRP-predicted compounds (bottom). (E) Violin plot displaying the differences in IC50 values of 6 CGP2016-predicted compounds and 7 CTRP-predicted compounds between high PRS and low PRS groups. (F) The CMap score of screened compounds in multiple HCC cell lines. (G) Venn diagram for intersecting the druggable targets based on mRNA expression and CERES score. (H) Volcano plot of Spearman correlation analysis between PRS and mRNA expression (left) or CERES score (right) of drug targets. Red and blue dots indicate significant correlations (p<0.05 and Spearman r >0.5 or <−0.5). (I) Scatter plot of PRS with mRNA expression (top) and CERES score (bottom) of 4 intersected targets. (J) PROGENy pathway activities in 5 transcriptomic subtypes. Abbreviations: PRS, prognosis-related score; TACE, transarterial chemoembolization.