HCC is a highly malignant tumor, often diagnosed at an advanced stage, making the selection of treatment options very challenging and the prognosis poor[19]. At the same time, HCC is a highly heterogeneous disease. Due to significant differences in the tumor microenvironment (TME), even the same treatment methods can lead to varying therapeutic outcomes[20]. Therefore, it is necessary to identify a new biomarker for the early detection of HCC, which could improve patients' prognosis. In this study, we analyzed cells at the single-cell level from HCC patients. We downloaded single-cell RNA sequencing data from 20 HCC patients from the GEO database (GSE235057). Data quality control was performed on the 20 samples. The cells were annotated into eight types, including T Cells, NK Cells, Plasma Cells, B Cells, Macrophages, Hepatocytes, and Fibroblasts. We found that NK cells were highly present in all 20 HCC samples. We then carried out further analysis on the NK cells, aiming to identify new biomarkers based on NK cells that can effectively predict the prognosis and immunotherapy response of HCC.
In this study, we obtained 1396 High_NK_Score genes from scRNA analysis. Next, we identified 250 NK-related genes from weighted gene co-expression network analysis (WGCNA). Eventually, we identified 227 intersecting genes for further analysis. KEGG enrichment analysis was performed to identify potential pathways. The results showed that the intersected genes were primarily enriched in the 'Natural Killer Cell-mediated cytotoxicity' pathway, making them suitable for further analysis. To identify genes related to prognosis and construct a model, we performed univariate Cox analysis. We used 101 types of machine learning combination algorithms to calculate the C-index for each model and select the most suitable one. Lasso + StepCox was regarded as the best model. We then conducted multivariate Cox analysis to further explore the prognostic genes. Finally, we identified the following six prognostic genes: KLRB1, STAT4, EHD1, PIK3IP1, PSMA3, and TNFAIP8.
Based on the 6 prognostic genes, we constructed a new prognostic risk score model. The risk score model demonstrated stable and good performance in predicting the prognosis of HCC patients; specifically, high-risk patients in both the TCGA training cohort and the ICGC-LIRI cohort exhibited worse clinical outcomes. Additionally, based on the risk score and survival time, we constructed a Kaplan-Meier curve to show the different survival outcomes between the high-risk and low-risk groups. Finally, a ROC curve was utilized to assess the stability and accuracy of our model. Surprisingly, we found that our model performed well in both the TCGA and ICGC-LIRI cohorts.
To determine the correlation between model genes and NK cell infiltration, we conducted extensive exploration. We explored the correlation between the levels of model genes and NK cell infiltration scores. Our findings showed that the model genes were positively correlated with NK cell infiltration. We then analyzed the correlation between clinical signatures and the risk score. We found that a high risk score was associated with higher T and N stages. Our nomogram demonstrated good accuracy in predicting the outcomes for HCC patients. We explored the distribution of model genes in NK cells at the single-cell level. KLRB1 has shown the highest content in NK cells of HCC and has shown a protective role in prognosis[21, 22, 23, 24]. The results of GSEA indicated that the model genes were enriched in five potential pathways. The results from Monocle 3 analysis showed that KLRB1 significantly affects the growth trajectory and may be a potential therapeutic target[24, 25].
We explored the relationship between risk and immune response. We found that the Natural Killer (NK) cell score was higher in the low-risk group in both the training and validation cohorts. The results showed that low-risk HCC patients had a higher NK cell infiltration score than high-risk patients, indicating better clinical outcomes[26, 27, 28, 29]. The higher NK cell Infiltration score means a lower risk and better clinical outcomes. At the same time, we found that a lower risk score was associated with higher IMMUNEScore, STROMALScore, and ESTIMATEScore. Previous studies have shown that higher IMMUNEScore, STROMALScore, and ESTIMATEScore indicate higher tumor purity and a greater likelihood of benefiting from immunotherapy[30, 31, 32, 33]. We can infer that the low-risk patients of our study had a better clinical outcome and immunotherapy response.
The high sensitivity to all types of immune checkpoint drugs and significant differences between different risk groups were examined, suggesting that they may have therapeutic effects on HCC. Using a composite score of the six genes, including KLRB1, STAT4, EHD1, PIK3IP1, PSMA3, and TNFAIP8, we can identify patients who are not sensitive to immunotherapy drugs and develop more precise treatment plans. Then we identified 12 kinds of common immunotherapy drugs that were sensitive to HCC patients.
Finally, we downloaded online immunohistochemistry results for the model genes to validate our model. From the images, we observed that KLRB1 and STAT4 were more overexpressed in normal tissue, while EHD1, PSMA3, and TNFAIP8 were more overexpressed in HCC tissue. These observations demonstrate that our model has good predictive performance for the prognosis of HCC and is suitable for clinical application.