Risk modeling of single-cell transcriptomes reveals the heterogeneity of immune infiltration in hepatocellular carcinoma

Abstract

Hepatocellular carcinoma (HCC) is one of the most common primary hepatic malignancies. E2F transcription factors play an important role in the tumorigenesis and progression of HCC, mainly through the RB/E2F pathway. Prognostic models for HCC based on gene signatures have been developed rapidly in recent years; however, their discriminating ability at the single-cell level remains elusive, which could reflect the underlying mechanisms driving the sample bifurcation. In this study, we constructed and validated a predictive model based on E2F expression, successfully stratifying patients with HCC into two groups with different survival risks. Then we used a single-cell dataset to test the discriminating ability of the predictive model on infiltrating T cells, demonstrating remarkable cellular heterogeneity as well as altered cell fates. We identified distinct cell subpopulations with diverse molecular characteristics. We also found that the distribution of cell subpopulations varied considerably across onset stages among patients, providing a fundamental basis for patient-oriented precision evaluation. Moreover, single-sample gene set enrichment analysis revealed that subsets of CD8⁺ T cells with significantly different cell adhesion levels could be associated with different patterns of tumor cell dissemination. Therefore, our findings linked the conventional prognostic gene signature to the immune microenvironment and cellular heterogeneity at the single-cell level, thus providing deeper insights into the understanding of HCC tumorigenesis.

Keywords: cellular heterogeneity; hepatocellular carcinoma; immune microenvironment; overall survival; prognosis.

Figure 1

E2F transcription factors are essential regulators in HCC.A and B, gene set enrichment analysis showed the five most enriched biological processes for the differentially expressed genes in the TCGA cohort (A) and the combined GEO cohort (B). C and D, differential expression of E2F-target genes in the TCGA cohort (C) and the integrated GEO cohort (D).

Figure 2

Performance of the predictive gene signature.A–C, Kaplan–Meier survival curves for comparison of the overall survival rates between patients in the low-risk group and the high-risk group for the TCGA training cohort (A), the TCGA validation cohort (B), and the external validation cohort from an independent GEO dataset (C). D, the nomogram predicting the 0.5-, 1-, and 1.5-year overall survival probability of the patients with HCC was created by integrating the gene signature with common clinical features. E, the calibration plot showing the predicted and the actual observed survival rates for patients with HCC at the 0.5-, 1-, and 1.5-year time points. F, survival curves of the nomogram in predicting overall survival probabilities for patients with HCC.

Figure 3

Deconvolution of bulk sequencing data showed the immune cell infiltration landscape in patients with HCC from the high-risk and low-risk groups.A, distribution in fractions of diverse tumor infiltrating immune cells in high-risk and low-risk patients with HCC. B, comparisons of the proportions of different tumor infiltrating immune cells between the high-risk and low-risk patients. Only the comparisons with statistical significance were visualized.

Figure 4

Discriminating T-cell populations based on single-cell RNA sequencing profiles.A, Uniform Manifold Approximation and Projection (UMAP) plot displaying 22 T-cell subpopulations identified by the predictive gene model based on the existing 11 clusters. Each dot represents one single cell. The existing 11 clusters: C01_CD8, naive CD8⁺ T cells; C02_CD8, effector CD8⁺ T cells; C03_CD8, mucosal-associated invariant T cells (MAIT); C04_CD8, exhausted CD8⁺ T cells; C05_CD8, intermediate state CD8⁺ T cells between effector and exhausted CD8⁺ T cells; C06_CD4, naive CD4⁺ T cells; C07_CD4, peripheral T regulatory cells (Tregs); C08_CD4, tumor Tregs; C09_CD4, mixed state cells; C10_CD4, exhausted CD4⁺ T cells; C11_CD4, cytotoxic CD4⁺ T cells. B, the proportions of class I and class II subpopulations in each stage of the patients with HCC. C, the proportions of each T-cell type between the class I and class II populations across peripheral blood, tumor, and adjacent normal tissues. D, UMAP plot of all T cells from peripheral blood and normal and tumor tissues, with all cells colored according to the gene expression of GZMK, CTLA4, and IL18R1, respectively. The two cell subpopulations between which the respective gene was differentially expressed were marked with red circles. E, UMAP plot of all T cells, with only tumor cells colored according the gene expression of GZMK, CTLA4, and IL18R1, respectively. Nontumor cells were colored gray. The two cell subpopulations between which the respective gene was differentially expressed were marked with red circles. F, differential expression of GZMK, CTLA4, and IL18R1 in all cell samples of CD4⁺ and CD8⁺ T cells. G, differential expression of GZMK, CTLA4, and IL18R1 in tumor cell samples of CD4⁺ and CD8⁺ T cells.

Figure 5

Different developmental cell fates between the class I and class II cells.A, faceted pseudotime plot according to T-cell types for the class I and class II subpopulations. B, heatmap plot for the expression of the 20 most dynamic genes along pseudotime for class I and class II subpopulations. C, the dynamic expression patterns of the marker genes CXCL13 and LAYN in class I and class II subpopulations. D, Gene Ontology results showed the biological processes that associated with the significantly dynamic genes from the branch-dependent expression analysis. Each node represents an enriched Gene Ontology term. E, single sample gene set enrichment analysis results for comparing the expression of the dynamic genes during branch evolution for the “Positive Regulation of Cell Adhesion” term and the “Positive Regulation of Immune System Process” term.