Upregulation of multiple key molecules is correlated with poor prognosis and immune infiltrates in hepatocellular carcinoma by bulk and single-cell RNA-seq

Abstract

Background: Recent discoveries in hepatocellular carcinoma (HCC) unveil key molecules. However, due to liver cancer's high heterogeneity, predicting patient prognosis is challenging. This study aims to construct a model for predicting HCC prognosis using multiple key genes.

Methods: TCGA provided RNA expression and clinical data, differentially analyzed by DESeq2, edgeR, and Limma. The hub gene was pinpointed via CytoHubba's degree algorithm in Cytoscape. GO and KEGG analyses illuminated potential pathways. Single-cell sequencing detailed key gene expression in diverse cell types. The LASSO regression model predicted patient prognosis.

Result: In the RNA-seq analysis using three R packages, we identified 762 differentially expressed genes, with Cytoscape revealing ten key genes showing significant prognostic value (P < 0.05). GO and KEGG analyses highlighted key biological processes and pathways. IHC confirmed higher expression in cancer tissues. Reduced immune cell infiltration was observed in HCC tissues, and immune checkpoint analysis showed a strong correlation between PD1, CTLA4, and hub genes. Single-cell sequencing indicated higher expression of key genes in immune cells than hepatocytes. Cox analysis validated the riskScore as a reliable, independent prognostic marker (HR = 4.498, 95% CI: 2.526-8.007).

Conclusions: The results from differential analysis using three R packages are robust, revealing genes closely linked to immune cell infiltration in the tumor microenvironment. Additionally, a validated prognostic model for liver cancer was established based on key genes.

Keywords: hepatocellular carcinoma; immune infiltration; immunotherapy; prognostic signature; tumor microenvironment.

Figure 1

A flow chart of the manuscript.

Figure 2

Difference of genomic landscape between normal and LIHC tissues. (A) Hierarchical clustering visualizing the intersections of DEGs with KM analysis. (B–D) Gene Ontology functional enrichment analyses for differentially expressed genes. (B) Biological process. (C) Molecular function. (D) Cellular component. (E) KEGG pathway enrichment analyses for differentially expressed genes. All enriched pathways were significant. The color depth represented enriched adjusted p-value.

Figure 3

Identification of key molecules in LIHC. (A) Venn diagram visualizing the intersections of DEGs with KM analysis. (B) Construction of a protein–protein interaction (PPI) network among differentially expressed genes. (C) The relationship among the ten key molecules at the protein level. (D) Volcano plot constructed with the cut-off criterion p < 0.05 and |logFC|≥ 1. Red, up-regulated genes; Green, down-regulated genes. The circle represented each gene and the identified key molecules were marked. (E) Principal component analysis for the key molecules revealed two completely disjoint populations, suggesting these key molecules could well distinguish LIHC samples from normal samples. Blue, normal samples; Red, tumor samples.

Figure 4

Immune cell infiltration and correlation analysis. (A) Differences in 28 TME infiltration cells between normal liver and LIHC tissues (^*P < 0.05; ^**P < 0.01; ^***P < 0.001). (B) The correlation between each key molecule and each immune checkpoint. Red, positive; Purple, negative. (C) Immunotherapy efficacy of 10 key genes.

Figure 5

Single-cell RNA-seq analysis. (A) The different cells were annotated. (B, C) Single-cell analysis was used to monitor the expression of 10 key genes in different immune cells.

Figure 6

Construction of riskScore signature. (A) Least absolute shrinkage and selection operator (LASSO) coefficient profiles of the ten key molecules. (B) Penalty plot for the LASSO model for the 10 prognostic genes with error bars denoting the standard errors. (C) The optimal cut-off point to dichotomize riskScore into low and high groups was determined by the function surv_cutpoint. The optimal cut-off point was 3.58. (D) Proportion of deaths in high and low risk groups as riskScore values increased. Hierarchical clustering of seven key genes between low and high risk groups. Red, up-regulated; Blue, down-regulated. (E) Survival analyses for low and high riskScore groups using Kaplan-Meier curves. (P < 0.0001, Log-rank test) (F) Predictive efficacy of riskScore on prognosis.

Figure 7

Prognostic value of the riskScore gene signature. (A) Forest plot showing the riskScore was an independent prognostic biomarker using multivariate analyses. (B) The nomogram was constructed to predict the probability of patient mortality. (C) The Predictive efficacy of nomogram score on prognosis. (D) The calibration plot of nomograms between predicted and observed 3-year and 5-year outcomes. The 45-degree line represented the ideal prediction. (E) The GSEA enrichment reveals two significantly activated signaling pathways, including the cell cycle pathway.

Figure 8

Association of riskScore with immunity. (A) Correlation of riskScore with immune cells. (B) Tumor stem cell relevance. (C) Differences in TMB scores in high and low risk groups. (D, E) Correlation of Mismatch Repair gene and immune checkpoints with riskScore. (F, G) Assessing differences in risk scores between immune efficacy groups.

Figure 9

Drug sensitivity and model validation. (A) The sensitivity of various drugs was assessed between the high and low risk groups. (B) GSE14520 and IMvigor210 were used to verify the accuracy of the model.