Identification of Immune Cell Correlations and Patient Consensus Clustering
We performed correlation analysis on immune cells in LIHC samples. In Fig. 1A, T cells CD4 naive and Neutrophil cells with the highest positive correlation because the red sector has the largest area. NK cells activated and Macrophages M0 has the highest negative correlation. Then we identified three subtypes, Cluster-A, B, and C based on the results of CIBERSORT (k=3, Fig. 1B, C). The survival curve (Fig. 1D) shows that Cluster-B had the best survival, while Cluster-C had the most dismal prognosis (p=0.032).
Identification of Immune-Related Subtype
To further understand the difference between different subtypes, we made a heatmap to show the infiltration degree of different immune cells among each subtype (Fig. 1E). We found that Cluster-B had the most abundant enrichment of immune cells, involving CD8 T cells and NK cells which are associated with tumor killing. Others included T cells follicular helper and plasma cells. CD4 memory resting T cells were significantly higher in Cluster-A, and Macrophages M0 cells were dramatically increased in Cluster-C. Then the result of ssGSEA shows that the infiltration degree of most immune cells was higher in Cluster-B too (Fig. 1F). Finally, by the result of ESTIMATE, we found no differences in the Stromal Scores of the three subtypes (Fig. 1G). Cluster-B had the highest Immune Score and ESTIMATE Score, and Cluster-C had the lowest score, which was consistent with the prognostic analysis. Cluster-C tumors have the highest purity and the worst prognosis, while Cluster-B tumors have the lowest purity and the best prognosis. Therefore, we believe that Cluster-A and Cluster-C fall into the cold tumor, while Cluster-B belongs to the hot tumor.
Response of Three Clusters to Immunotherapy and Chemotherapy
ICIs such as PD-1 and CTLA-4 are of great promise in LIHC therapy. We validated the expression of 13 immune checkpoint-related genes in 3 clusters. The results show that PD1 related (Fig. 2A), CTLA4 related (Fig. 2B), and other checkpoint genes (Fig. 2C) were highly expressed in Cluster-B.
In Cluster-A and Cluster-C, the expression of most immune checkpoint-related genes was not significantly different, but all were lower than in Cluster-B. Previous studies reported that cold tumors were insensitive to treatment with ICIs. In our opinion, the results of immune checkpoint genes validate our judgment on the three clusters of subtypes. We evaluated the response of chemotherapy drugs. Then, we found 9 drugs with potential predominance in Cluster-A (Fig. 2D), 53 in Cluster-B (Fig. 2E), and 9 in Cluster-C (Fig. 2F).
Functional Enrichment Analysis and PPI Network
Cluster B belongs to the hot tumor, while Cluster-A and Cluster-C are cold tumors. To further explore the genetic differences between the two tumors, we identified 42 DEGs in Cluster-B vs Cluster-A or Cluster-C. The results in GO (Fig. 3A) and KEGG (Fig. 3B) enrichment analysis suggest that the differences between hot and cold tumors may involve PD-L1 expression and PD-1 checkpoint pathway in cancer, and other immune pathways. In addition, there are proliferation, differentiation, and activation of various immune cells. Finally, we mapped the protein interaction network of 42 DEGs (Fig. 3C).
Establish the Prognostic Model
Considering the immune heterogeneity of the 3 Cluster, we have established a prognostic model based on this heterogeneity. We incorporated 42 DEGs into the Lasso regression analysis to reduce the number of predicted variables. Then, five hub genes were selected by the LASSO-COX regression analysis (Fig. 3D). The formula of Riskscore is as follows: . Exp is the expression value of the hub gene (Table S1). Fig. 3E and 3F show the distribution of Riskscore and the survival status of LIHC patients. The K-M survival curve shows that the low-risk group has a better prognosis (p<0.001) (Fig. 3G). Fig. 3H exhibits the distribution of various clinical features and hub gene expression levels in two risk groups. The expression of EOMES, GZMH, NKG7, and CD8A are higher in the low-risk group. The expression of ANKRD22 is higher in the high-risk group. Furthermore, immunohistochemical staining revealed differences in partial hub genes between LIHC and normal tissues (Fig. 3I). Fig. 3J shows higher mortality in the high-risk group, and Fig. 3K indicates that there is a significant difference in Riskscore between dead and alive patients (P<0.001).
We included age, gender, stage, and Riskscore in univariate (Fig. 3L) and multivariate (Fig. 3M) regression analysis. Results reveal that Riskscore is an independent prognostic factor for LIHC patients.
To assess the predictive accuracy of our models, we employed the ROC curve (Fig. 3N). The area under the curve (AUC) value (0.68) of Riskscore is higher than other clinicopathological characteristics, indicating that Riskscore has the best predictive ability.
Immune Infiltration and Sensitivity to Immunotherapy
As shown in Fig. 4A, the low-risk group had higher tumor infiltration of diverse anti-cancer immune cells, including CD8+ T cells, macrophages, Th1 cells, NK cells, and dendritic cells. However, M2 macrophages, which might have a suppressive TME, are more accumulated in the high-risk group. Fig. 4B shows that 21 immune cells were more strongly infiltrated in the low-risk group. The results in Fig. 4C indicate that the low-risk group had more abundant stromal components and immune infiltration. The higher tumor purity in high-risk groups may lead to a poor prognosis. The PD1 related (Fig. 4D), CTLA4 related (Fig. 4E), other immune checkpoint genes (Fig. 4F), and T cell activation agonists (Fig. 4G) are apparently higher in the low-risk group. In consideration of the positive correlation between the expression of the immune checkpoint genes and the response to immunotherapy, the low-risk group may benefit much more from immunotherapy. In conclusion, the low-risk group has stronger antitumor immunity, thereby displaying better immunotherapy response and has a more favorable prognosis.
Somatic mutation and stem cell index in two groups, and Nomogram construction
To further investigate the heterogeneity of the two risk groups, we analyzed the somatic mutation distribution and stem cell index. We found that in the TCGA-LIHC cohort, the high-risk group had a higher somatic mutation rate. (Fig. 5A, B). The mutation landscape showed the top 30 mutant genes. Then we compared the cancer stem cell index of the two groups and observed a higher stem cell index in the high-risk group (Fig. 5C). The Riskscore had a weaker positive correlation with the stem cell index (Fig. 5D). Fig. 5E is a nomogram predicting the overall survival of LIHC, combined with Riskscore and clinical features.