Combined analysis of bulk and single-cell RNA sequencing reveals novel natural killer cell-related prognostic biomarkers for predicting immunotherapeutic response in hepatocellular carcinoma

Introduction: Natural killer (NK) cells play an irreplaceable and important role as a subtype of innate immune cells in the contemporary setting of antitumor immunity.

Methods: We chose a total of 1,196 samples for this analysis from the public dataset’s six separate cohorts. To identify 42 NK cell marker genes, we first carried out a thorough study of single-cell RNA sequencing data from the GSE149614 cohort of hepatocellular carcinoma (HCC).

Results: Using the NK cell marker genes in the TCGA cohort, we next created a seven-gene prognostic signature, separating the patients into two categories with distinct survival patterns. This signature’s prognostic prediction ability was well verified across several validation cohorts. Patients with high scores had higher TIDE scores but lower immune cell infiltration percentages. Importantly, low-scoring patients had superior immunotherapy response and prognosis than high-scoring patients in an independent immunotherapy cohort (IMvigor210). Finally, we used CD56 and TUBA1B antibodies for immunohistochemical labeling of HCC tissue sections, and we discovered a lower number of CD56+ cells in the HCC tissue sections with high TUBA1B expression.

Discussion: In summary, our research created a unique prognostic profile based on NK cell marker genes that may accurately predict how well immunotherapy would work for HCC patients.

1 Introduction

It is generally recognized that a diverse tumor microenvironment (TME) surrounds tumor cells in hepatocellular carcinoma (HCC) (1). The TME, which has a very complicated makeup, is crucial to the development and growth of tumors. Additionally, the interaction between immune cells and tumor cells in the TME not only influences a patient’s prognosis but may also alter a patient’s response to immunotherapy (2). The importance of innate immune cells has not gotten enough attention in the contemporary setting of antitumor immunity, which has mostly focused on adaptive T-cell responses, like CD4+ CD25+ Foxp3 regulatory T cells (Tregs), cytotoxic T lymphocytes (CTLs) et all (3, 4). By specifically identifying and eliminating tumor cells and encouraging adaptive T-cell immunity responses, natural killer (NK) cells are a subtype of innate immune cells that can reduce the proliferative and invasive potential of tumor cells at an early tumor stage (5). The balance of inhibitory and activating receptors that can interact with ligands on target cells determines how well NK cell function. NK cells can collaborate with T cells to control the spread of cancer and play a crucial part in the development of antitumor immunity. Cancer risk is increased by decreased NK cell activity in peripheral blood (6). Additionally, more tumor-infiltrating NK cells are strongly linked to improved prognosis across a variety of tumor types (7). Given the function of NK cells in immunity, earlier research has focused on their molecular features in cancer and infectious disorders (8, 9), but little is known about their complete molecular analysis in HCC.

Unprecedented chances to unveil the molecular properties of various immune cell populations in TME have been made possible by the advent of single-cell RNA sequencing (scRNA-seq) technology and related data processing methodologies (10). The prognosis and immunotherapeutic response in cancer patients may be accurately predicted by examining gene expression patterns based on molecular characterization of immune cells acquired from scRNA-seq data, according to previous research (11, 12). In this work, we first carried out a thorough examination of scRNA-seq data in HCC to characterize the molecular properties of tumor-infiltrating NK cells and to identify NK cell flag genes. Then, using bulk RNA-seq analysis, NK cell marker gene-related signatures for predicting the prognosis of HCC were created. Additionally, the link between NKCMGS and HCC immunotherapy response was examined, and the predictive ability of NKCMGS was verified in three separate cohorts from the ICGC and the Gene Expression Omnibus (GEO) database. Multiple datasets from TCGA, GEO and ICGC cohort were analyzed in our study for constructed NK cell-related genetic signature. We obtained a more parsimonious gene signature over existing studies, which contains seven genes, and provided better prediction for immunotherapeutic effect and drug sensitivity.

2 Methods

2.1 Data collection and pre-processing

A total of 1196 samples, 31396 cells, including 10 HCC samples with single-cell RNA-sequencing (scRNA-seq) data from the GSE149614 cohort, 342 HCC samples from the Cancer Genome Atlas (TCGA) cohort (https://xenabrowser.net/), 230 HCC samples from the International Cancer Genome Consortium (ICGC) cohort (https://dcc.icgc.org/), 221 HCC samples from the GSE14520 cohort, 95 HCC samples from the GSE76427 cohort, and 298 samples treated with immunotherapy from the IMvigor210 cohort (http://research-pub.Gene.com/imvigor210corebiologies/), were enrolled in this study. In the GSE149614 dataset, with each gene expressed in at least three cells and each cell expressing more than 250 genes, single cells were initially screened for scRNA-seq data. The percentage of mitochondria and rRNA was then calculated by the Seurat package (13). Further testing of the single cells involved caused each one to express at least 5000 genes with a UMI > 100. The mitochondrial content was no more than 30%. In the end, 31396 cells were still present for identifying the NK cell marker genes of HCC. To find survival-related genes and create prognostic signatures, the bulk transcriptome data (FPKM normalized) and clinical details of HCC patients in the TCGA were employed. For external validation, three separate datasets were acquired: ICGC, GSE14520, and GSE76427.

2.2 Identification of NK cells and their hub genes

The GSE149614 dataset contains scRNA-seq data from 10 HCC samples, which we again examined. Following log normalization of the expressed genes, uniform flow-form approximation and projection techniques were used to reduce nonlinear dimensionality. We used the FindNeighbors and FindClusters () algorithms to arrange individual cells into 17 separate subgroups at dim=50 and resolution=0.1. Three marker genes, including CD3D, CD3E, and NKG7, were identified in NK cells. Using the FindAllMarkers program with logFC=0.5, minpct=0.25, and adjusted p-values less than 0.05, marker genes were found for each NK subpopulation. A univariate Cox regression analysis with P less than 0.05 was then used to further identify the genes among these NK marker genes that are associated with prognosis. We used the LASO-Cox regression to compress the number of genes and created a risk profile based on the outcomes of the multivariate Cox model using the equation: Risk score =∑iCoefficient (Genei)*Expression (Genei). Depending on their risk assessments, patients were separated into high- and low-risk groups. The receiver operating characteristic curve (ROC) analysis and Kaplan-Meier survival analysis were used to examine the risk profile’s ability to predict survival outcomes. The validation cohort underwent a similar examination. GSEA was used to examine KEGG, GO, and HALLMARK elements that had drastically changed across the different categories.

2.3 Analysis of the immune landscape

Based on the gene expression patterns of HCC patients, the stromal and immune scores were computed using the ESTIMATE software (14). The abundance ratio of immune cells was evaluated using the CIBERSORT (15), MCPcounter (16), and TIMER (17) databases to learn more about the TME.

2.4 Response to immune checkpoint blockade (ICB) and analysis of the sensitivity to potential therapeutic drugs

The TIDE database (http://tide.dfci.harvard.edu/), which estimates how often immunotherapy for HCC patients will be effective, was first performed. We also retrieved the matched clinical and transcriptome data from the IMvigor210 cohort of patients who were receiving anti-PD-L1 medication. Additionally, we evaluated multiple immune checkpoint gene expression changes such as PD1, PD-L1, CTLA4, and PD-L2 in different subgroups. Finally, we searched the Cellminer database for delicate medications that successfully address this risk profile (18). If a drug’s adjusted P-value was less than 0.001 and its Pearson correlation coefficient was larger than 0.3, it was categorized as tumor-sensitive. The discrepancies in the half-maximal inhibitory doses (IC50) of several classes of tumor-sensitive medications were then studied.

2.5 Evaluation of TUBA1B expression in clinical samples

Using 30 samples of HCC and related paracancerous tissues that had undergone standard pathological evaluation in our pathology department, a validation cohort was developed. The clinicopathological information for all patients were shown in Table 1. Following the tissue wax blocks’ serial sectioning, sample sections were gathered and stored for subsequent use in a 4°C freezer. We next used TUBA1B (Abcam, ab108629) and CD56 (Abcam, ab75813) antibodies in IHC experiments on formalin-fixed de-paraffinized slices and captured pictures using microscopy as previously reported (19).

TABLE 1

Table 1 Characteristics of patients and tumor samples studied (n=30).

2.6 Statistical analysis

R software was used to conduct all statistical analyses (v4.1.2). Pearson correlation was used to calculate the correlation analysis. For comparisons between the two groups, the chi-square test and grouped t-test were used, respectively. Kaplan-Meier survival analysis and a Log-rank test were used to assess survival differences between groups. Using the RMS software, a nomogram was produced following the signature. Statistics were deemed significant when the P value was less than 0.05.

3 Results

3.1 Identification of NK cells in the scRNA-seq samples

Figure S1 displays the full outcomes of data preparation. After log-normalization and dimensionality reduction, 17 clusters were found. The TSNE plots showing the distribution of the 17 clusters are displayed in Figure 1A. As shown in Figure 1B and S2A, based on the expression of three marker genes (CD3D, CD3E, NKG7, CXCR3 and IL2RB), two NK cell subsets were discovered (Figure 1C). The fact that neither of the two NK cell clusters expressed the epithelial cell-specific gene (PECAM1) proves that NKs were correctly identified (Figure 1D). Further analysis of CD19 and CD14 expression in 17 clusters was performed for ruled out the interference of other cell types (Figure S2B). The expression of the top 10 DEGs in the two clusters is shown in Figure 1E. The 2 NK cell clusters contained 42 DEGs (marker genes recognized as NK clusters). The percentage of the two clusters in each sample was shown in Figure 1F. Furthermore, we computed the ssGSEA scores for the marker genes of each NK cluster (the top 10 DEGs of the NK clusters) in the TCGA cohort to examine the connection between NK clusters and prognosis. The samples in the high ssGSEA score group in the NK 0 cluster had a better prognosis than those in the low ssGSEA score group, while the opposite finding was seen in the NK 1 cluster (Figure 1G). Finally, to further analyze the function and mechanism of NKs marker genes in HCC, we performed molecular subtype identification analysis of the TCGA dataset by non-negative matrix decomposition (NMF) algorithm. The HCC samples were split into two subclasses based on the NKs marker genes after it was found that two clusters were the ideal number (Figures 2A, S3). Significant disparities in patient survival existed between the two subgroups (Figure 2B). Additionally, the two subgroups’ TME features were contrasted. Figures 2C–E demonstrates that as compared to samples from cluster 2, samples from cluster 1 had higher immunological, stromal, and ESTIMATE scores. HCC patients in Cluster 1 had a larger percentage of immune cell infiltration in their TME, as shown by the results of the TIMER (Figure 2F), MCPcounter (Figure 2G), and CIBIS. ORT (Figure 2H).

FIGURE 1

Figure 1 Identification of NK cells in the scRNA-seq samples. (A) The TSNE plots showed the distribution of the 17 clusters. (B) The expression of three marker genes (CD3D, CD3E, and NKG7) in the 17 clusters. (C) Two NK cell subsets were discovered. (D) The expression of PECAM1 gene. (E) The expression of the top 10 DEGs in the two clusters. (F) The percentage of the two clusters in each sample. (G) The samples in the high ssGSEA score group in the NK 0 cluster had a better prognosis than those in the low ssGSEA score group, while the opposite finding was seen in the NK 1 cluster.

FIGURE 2

Figure 2 Subtypes identification by NMF algorithm based on NKs marker genes. (A) Samples were split into two subclasses. (B) Significant disparities in patient survival existed between the two subgroups. (C-E) Difference of immunological, stromal, and ESTIMATE scores in different molecular subtypes. The difference in the percentage of immune cell infiltration in different molecular subtypes was analyzed by the TIMER (F), MCPcounter (G), and CIBISORT (H). ns, not significant; *p < 0.05; **p <0.01; ***p < 0.001.

3.2 Screening for NK-associated hub genes

Using univariate Cox regression analysis, the prognostic value of these DEGs was evaluated to create a signature, with 10 genes displaying prognostic values (Figure 3A). To reduce the number of genes, Lasso-Cox regression analysis was used (Figure 3B). Seven genes were left with a lambda value of 0.0139 (Figure 3C). After multivariate Cox regression analysis, we finally included these seven genes (CREM, PFN1, KLRB1, TUBA1B, APOC1, ACTG1, and HSPA1A) in the signature. Following is the final seven-gene signature formula: Score = (0.1574×CREM) + (0.4582×PFN1) - (0.3804×KLRB1)+(0.1094×TUBA1B)-(0.0821×APOC1)+(0.3267×ACTG1) + (0.1718×HSPA1A). After each sample’s score was determined, the groups of high- and low-risk individuals were created (Figure 3D). The association between score and clinical characteristics was first evaluated and found that higher scores were associated with HBV infection, advanced TNM stage, later grade, later T stage, and recurrence (Figure S4). In both the TCGA, ICGC, and the GEO cohort, Kaplan-Meier survival analyses showed that high-risk patients had considerably worse survival outcomes than low-risk patients (Figures 3E, 4A). The TCGA cohort’s AUC values for the model for 1- to 3-year survival range from 0.71 to 0.77 (Figures 3F), whereas those for the ICGC and GEO cohorts ranged from 0.57 to 0.72 (Figures 4B).

FIGURE 3

Figure 3 Screening for NK-associated hub genes. (A) Using univariate Cox regression analysis. (B) Lasso-Cox regression analysis. (C) Seven genes were left with a lambda value of 0.0139. (D) After each sample’s score was determined, the groups of high- and low-risk individuals were created. (E) Kaplan-Meier survival analysis. (F) ROC analysis.

FIGURE 4

Figure 4 Validation of the signature in the ICGC and GEO cohorts. (A) Kaplan-Meier survival analysis. (B) ROC analysis.

3.3 Mutation and functional enrichment analysis of the signature

The results of Gene Set Enrichment Analysis (GSEA) revealed that the majority of the impacted HALLMARK, GO, and KEGG components were engaged in DNA synthesis and replication, mitosis, chromosome segregation, and other biological processes related to the cell cycle (Figure S5). Studies on genetic modification that focused on significantly altered genes showed that the mutation rates in the two groups were very different from one another (Figure S6A). After tumor mutational burden (TMB) values for each HCC patient were analyzed, we found that patients in the high-score group with greater TMB values had the lowest overall survival rates, while the opposite results were found in patients in the low-score group with lower TMB values (Figure S6B).

3.4 Correlation analysis between the signature and immunity

Figure 5A showed that samples in the low-scoring group had higher immune, stromal, and ESTIMATE scores compared to samples in the high-scoring group. As shown in the results of TIMER (Figure 5B), MCPcounter (Figure 5C), and CIBISORT (Figure 5D), a greater percentage of immune cell infiltration was found in the TME of HCC patients in the low-scoring group.

FIGURE 5

Figure 5 Correlation Analysis between the signature and immunity. (A) samples in the low-scoring group had higher immune, stromal, and ESTIMATE scores compared to samples in the high-scoring group. The difference in the percentage of immune cell infiltration in different molecular subtypes was analyzed by the TIMER (B), MCPcounter (C), and CIBISORT (D). ns, not significant; *p < 0.05; **p <0.01; ***p < 0.001.

3.5 The signature’s response to PD-L1 blockade immunotherapy

The Tumor Immune Dysfunction and Exclusion (TIDE) analysis showed that, although the exclusion scores had the opposite impact, the TIDE and dysfunction scores were significantly greater in the group with higher risk scoring than in the group with lower risk scoring (Figure 6A). When the projected immunotherapy response rate was included, the proportion of “respond” was lower in the high-risk group (Figure 6B). Then, using data from the IMvigor210 cohort, we assessed the predictive efficacy of immune checkpoint treatment risk factors. In comparison to the high-scoring group, patients in the low-scoring group saw notable clinical benefits and considerably longer overall life (Figure 6C). Figure 6D showed that patients with progressing disease (PD) or stable disease (SD) had greater risk ratings than those who had a complete response (CR)/partial response (PR). Finally, we discovered that the levels of the genes PD-1, PD-L1, PD-L2, CTLA4, CD4, CXCR4, LAG3, and LL6 were higher in patients with lower scores than in patients with higher scores (Figure 6E), suggesting that these ICIs may be more beneficial for patients with lower scores.

FIGURE 6

Figure 6 The signature’s response to PD-L1 blockade immunotherapy. (A) The TIDE analysis. (B) When the projected immunotherapy response rate was included, the proportion of “respond” was lower in the high-risk group. (C) In comparison to the high-scoring group, patients in the low-scoring group saw notable clinical benefits and considerably longer overall life in the IMvigor210 cohort. (D) Patients with PD/SD had greater risk ratings than those who had a complete response CR/PR. (E) The expression levels of the ICIs genes in the two groups. ns, not significant; **p <0.01; ***p < 0.001.

3.6 Construction of a nomogram model and exploration of potential drug sensitivity

As shown in Figure S7A, we created a nomogram incorporating clinical features and the signature to maximize the predictive performance of risk characteristics. The calibration plots demonstrated that the nomogram was capable of accurately forecasting the final survival rate (Figure S7B). In addition, we identified 7 drugs with tumor sensitivity (Figure S8A). As shown in Figure S8B, we also found that the IC50 for Cladribine, Fludarabine, and Clofarabine was lower in patients with higher scores (20–22).

3.7 TUBA1B expression in HCC

We first initially investigated the expression and prognostic value of these seven genes in HCC in the GEPIA database (23). We found that only TUBA1B showed differential expression in HCC and normal tissues (Figure S9A), although several genes including KLRB1, TUBA1B, APOC1, ACTG1, and HSPA1A had high prognostic values (Figure S9B). We then focused our main attention on TUBA1B. Using the Human Protein Atlas database (HPA) (24), we found significant variability in the protein expression of TUBA1B in normal and HCC tissues (Figure 7A). This phenomenon was confirmed in clinical HCC and normal tissue sections (Figures 7B, C). Last but not least, we found differential expression of CD56 in these HCC tissues (Figure 7D) and a lower number of CD56-positive cells in the HCC tissue sections with high TUBA1B expression (Figure 7E).

FIGURE 7

Figure 7 TUBA1B expression in HCC. (A) TUBA1B expression explored in HPA database. (B-C) TUBA1B expression explored by IHC in clinical samples. (D) CD56 expression explored by IHC in HCC samples. (E) A lower number of CD56-positive cells in the HCC tissue sections with high TUBA1B expression.

4 Discussion

Researchers are learning more about the variety and heterogeneity of TME as well as the molecular properties of tumor-infiltrating immune cells thanks to the quick development of scRNA-seq technology (25). However, the majority of recent research has concentrated on adaptive immune cells, and the function of innate immune cells has not received enough attention, which may have a significant impact on the prognosis and effectiveness of immunotherapy in patients with tumors (26). If neighboring cells exhibit surface markers linked to oncogenic transformation, NK cells can quickly destroy many of them by improving antibody and T-cell responses (27). The prognosis of patients with various tumors, including lung adenocarcinoma (28), gastric cancer (29, 30), liver cancer (31), melanoma (32), and colorectal cancer (33), is highly correlated with the number of tumor-infiltrating NK cells. The overall survival of HCC following hepatectomy is significantly impacted by the low frequency of NK cells relative to myeloid and other lymphocytes seen in HCC tumor tissue (34). The number of NK cells inside the TME, which has a favorable correlation with patient survival, also influences how well patients respond to sorafenib therapy (35). In the current work, we aimed to investigate NK cell marker genes in HCC by bulk and scRNA-seq analysis as well as to create a transcriptional signature based on NK cell marker genes to evaluate NK cell infiltration in TME. By boosting NK scores, we discovered a substantial classification of HCC patients’ prognosis that was well verified across three separate cohort datasets. Additionally, we discovered that immunotherapy response rates were much greater for patients with low NK scores than for patients with high NK scores, indicating that immune checkpoint blockade treatment is better suitable for individuals with low NK scores.

Immune cells that invade tumors and contribute considerably to tumor growth might have a negative impact on a patient’s prognosis if they have HCC (36). By using the TIMER, MCPcounter, and CIBERSORT algorithms to estimate and compare the abundance of immune cell infiltration between high and low NK score populations, we discovered higher levels of immune cell infiltration, particularly T and B cells, in low NK score tumors, indicating that low NK score tumors are referred to as “hot tumors” with increased antitumor activity (37). The greater survival rate of patients with low NK scores may be partially explained by the strong immune cell infiltration, which may promote tumor cell attenuation to avoid immune monitoring and impede tumor development.

Taking into account that variations in immune cell infiltration between different NK score subgroups affect the effectiveness of immunotherapy that patients receive, we first examined the clinical response to immunotherapy in HCC patients using the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm (38). We found significantly higher TIDE scores in the higher NK score group than in the lower NK score group, and a lower proportion of “respond” in the high NK score group. Subsequently, we validated the predictive power of our NK score with an immunotherapy cohort (Imvigor210). We found that patients with progressing disease (PD) or stable disease (SD) had greater risk ratings than those who had a complete response (CR)/partial response (PR). In light of the possibility that complex TME can cause HCC cells to develop resistance to immune checkpoint inhibitors (ICIs), which could affect the efficacy of immunotherapy, we also looked at the differences in the expression levels of various immune checkpoint genes between high- and low-NK scores subgroups. It has been demonstrated that patients with lower scores had larger prevalences of the genes PD-1, PD-L1, PD-L2, CTLA4, CD4, CXCR4, LAG3, and LL6 than people with higher scores. In conclusion, NK scores may be a valid biomarker for predicting response to immunotherapy, and patients with low NK scores are more likely to benefit from it.

Finally, utilizing the CellMiner database, we discovered seven medicines that are tumor sensitive. Since individuals with higher NK scores had lower IC50s for Cladribine, Fludarabine, and Clofarabine, it is obvious that these patients are more susceptible to these medications. Cladribine and Clofarabine are nucleoside analogs that are frequently used to treat hematologic cancers and target B and T cells (39, 40). Cladribine has been successfully utilized as a first-line therapy for hairy cell leukemia for some time now (41). Unfortunately, when used to treat multiple sclerosis, cladribine can lead to acute, specific liver harm in individuals (42). Additionally, clofarabine is employed as an anticancer therapy for several solid tumors, including bladder (43), stomach (44), and breast malignancies (45). Since fludarabine dramatically reduces the release of HBV progenitor DNA, it has been used to treat HBV infection and enhance the prognosis of HCC that is related to HBV (46). Combining fludarabine with fusion proteins comprising the poliovirus receptor (PVR) and the programmed death-1 (PD-1) extracellular structural domain improves long-term tumor-specific immunosurveillance and CD8+ T cell-mediated anticancer effects (47). However, additional research is required to confirm if these drugs can eventually increase tumor cell death by targeting NK cells.

Despite the encouraging findings, there are several limitations to this study. First, the candidate genes involved in our study were limited to NK cell marker genes, and the tumor immune microenvironment is highly spatially heterogeneous; second, a sizable multicenter randomized controlled trial will be needed in the future to evaluate this signature. Finally, all mechanistic analyses in our study were descriptive. Future studies must explore the potential mechanisms between NK marker gene expression and HCC prognosis.

5 Conclusion

In summary, a prognostic seven-gene signature built on NK cell marker genes was discovered and proven to have a strong performance in predicting immunotherapy response in HCC patients. It may be used as a prognostic biomarker to aid in the selection of suitable individuals who would benefit from immunotherapy and support therapeutic decision-making about customized prediction.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found within the article/Supplementary Materials.

Ethics statement

The studies involving human participants were reviewed and approved by Ethics Committees of Zhengzhou University. The patients/participants provided their written informed consent to participate in this study.

Author contributions

KZ designed this work and analyzed the data; EY performed experiments and analyses; KZ helped for providing tumor samples; KZ wrote the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work is supported by PhD research startup foundation of the Third Affiliated Hospital of Zhengzhou University (2021079).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2023.1142126/full#supplementary-material

Abbreviations

NK, Natural killer; HCC, hepatocellular carcinoma; TME, tumor microenvironment; scRNA-seq, single-cell RNA sequencing; GEO, Gene Expression Omnibus; TCGA, the Cancer Genome Atlas; ICGC, International Cancer Genome Consortium; ICB, immune checkpoint blockade; IC50, half-maximal inhibitory doses; NMF, non-negative matrix decomposition; GSEA, Gene Set Enrichment Analysis; TMB, tumor mutational burden; TIDE, Tumor Immune Dysfunction and Exclusion; PD, progressing disease; SD, stable disease; CR, complete response; PR, partial response; HPA, Human Protein Atlas; ICIs, immune checkpoint inhibitors.

References

1. Lu Y, Yang A, Quan C, Pan Y, Zhang H, Li Y, et al. A single-cell atlas of the multicellular ecosystem of primary and metastatic hepatocellular carcinoma. Nat Commun (2022) 13(1):4594. doi: 10.1038/s41467-022-32283-3

PubMed Abstract | CrossRef Full Text | Google Scholar

2. He Q, Yang J, Jin Y. Immune infiltration and clinical significance analyses of the coagulation-related genes in hepatocellular carcinoma. Briefings Bioinf (2022) 23(4). doi: 10.1093/bib/bbac291

CrossRef Full Text | Google Scholar

3. Granito A, Muratori L, Lalanne C, Quarneti C, Ferri S, Guidi M, et al. Hepatocellular carcinoma in viral and autoimmune liver diseases: Role of CD4+ CD25+ Foxp3+ regulatory T cells in the immune microenvironment. World J Gastroenterol (2021) 27(22):2994–3009. doi: 10.3748/wjg.v27.i22.2994

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Saeidi V, Doudican N, Carucci JA. Understanding the squamous cell carcinoma immune microenvironment. Front Immunol (2023) 14:1084873. doi: 10.3389/fimmu.2023.1084873

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Sajid M, Liu L, Sun C. The dynamic role of NK cells in liver cancers: Role in HCC and HBV associated HCC and its therapeutic implications. Front Immunol (2022) 13:887186. doi: 10.3389/fimmu.2022.887186

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Wu SY, Fu T, Jiang YZ, Shao ZM. Natural killer cells in cancer biology and therapy. Mol Cancer (2020) 19(1):120. doi: 10.1158/1557-3125.HIPPO19-IA12

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Liu S, Galat V, Galat Y, Lee YKA, Wainwright D, Wu J. NK cell-based cancer immunotherapy: from basic biology to clinical development. J Hematol Oncol (2021) 14(1):7. doi: 10.1186/s13045-020-01014-w

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Björkström NK, Strunz B, Ljunggren HG. Natural killer cells in antiviral immunity. Nat Rev Immunol (2022) 22(2):112–23. doi: 10.1038/s41577-021-00558-3

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Terrén I, Orrantia A, Vitallé J, Zenarruzabeitia O, Borrego F. NK cell metabolism and tumor microenvironment. Front Immunol (2019) 10:2278. doi: 10.3389/fimmu.2019.02278

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Li X, Sun Z, Peng G, Xiao Y, Guo J, Wu B, et al. Single-cell RNA sequencing reveals a pro-invasive cancer-associated fibroblast subgroup associated with poor clinical outcomes in patients with gastric cancer. Theranostics (2022) 12(2):620–38. doi: 10.7150/thno.60540

PubMed Abstract | CrossRef Full Text | Google Scholar

11. He J, Xiong X, Yang H, Li D, Liu X, Li S, et al. Defined tumor antigen-specific T cells potentiate personalized TCR-T cell therapy and prediction of immunotherapy response. Cell Res (2022) 32(6):530–42. doi: 10.1038/s41422-022-00627-9

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Zhang Z, Wang ZX, Chen YX, Wu HX, Yin L, Zhao Q, et al. Integrated analysis of single-cell and bulk RNA sequencing data reveals a pan-cancer stemness signature predicting immunotherapy response. Genome Med (2022) 14(1):45. doi: 10.1186/s13073-022-01050-w

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Butler A, Hoffman P, Smibert P, Papalexi E, Satija R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol (2018) 36(5):411–20. doi: 10.1038/nbt.4096

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun (2013) 4:2612. doi: 10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Gentles AJ, Newman AM, Liu CL, Bratman SV, Feng W, Kim D, et al. The prognostic landscape of genes and infiltrating immune cells across human cancers. Nat Med (2015) 21(8):938–45. doi: 10.1038/nm.3909

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Zhu K, Xiaoqiang L, Deng W, Wang G, Fu B. Development and validation of a novel lipid metabolism-related gene prognostic signature and candidate drugs for patients with bladder cancer. Lipids Health Dis (2021) 20(1):146. doi: 10.1186/s12944-021-01554-1

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS, et al. TIMER: A web server for comprehensive analysis of tumor-infiltrating immune cells. Cancer Res (2017) 77(21):e108–10. doi: 10.1158/1538-7445.AM2017-108

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Reinhold WC, Sunshine M, Liu H, Varma S, Kohn KW, Morris J, et al. CellMiner: a web-based suite of genomic and pharmacologic tools to explore transcript and drug patterns in the NCI-60 cell line set. Cancer Res (2012) 72(14):3499–511. doi: 10.1158/0008-5472.CAN-12-1370

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Zhang G. Regulatory T-cells-related signature for identifying a prognostic subtype of hepatocellular carcinoma with an exhausted tumor microenvironment. Front Immunol (2022) 13:975762. doi: 10.3389/fimmu.2022.975762

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Bayas A, Christ M, Faissner S, Klehmet J, Pul R, Skripuletz T, et al. Disease-modifying therapies for relapsing/active secondary progressive multiple sclerosis - a review of population-specific evidence from randomized clinical trials. Ther Adv neurological Disord (2023) 16:17562864221146836. doi: 10.1177/17562864221146836

CrossRef Full Text | Google Scholar

21. Murthy GSG, Kim S, Estrada-Merly N, Abid MB, Aljurf M, Assal A, et al. Association between the choice of the conditioning regimen and outcomes of allogeneic hematopoietic cell transplantation for myelofibrosis. Haematologica (2023). doi: 10.3324/haematol.2022.281958

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Li Y, Zhao J, Xue Z, Tsang C, Qiao X, Dong L, et al. Aptamer nucleotide analog drug conjugates in the targeting therapy of cancers. Front Cell Dev Biol (2022) 10:1053984. doi: 10.3389/fcell.2022.1053984

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Tang Z, Li C, Kang B, Gao G, Li C, Zhang Z. GEPIA: a web server for cancer and normal gene expression profiling and interactive analyses. Nucleic Acids Res (2017) 45(W1):W98–w102. doi: 10.1093/nar/gkx247

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Colwill K, Gräslund S. A roadmap to generate renewable protein binders to the human proteome. Nat Methods (2011) 8(7):551–8. doi: 10.1038/nmeth.1607

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Wang X, Miao J, Wang S, Shen R, Zhang S, Tian Y, et al. Single-cell RNA-seq reveals the genesis and heterogeneity of tumor microenvironment in pancreatic undifferentiated carcinoma with osteoclast-like giant-cells. Mol Cancer (2022) 21(1):133. doi: 10.1186/s12943-021-01484-7

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Liu Y, Zhang Q, Xing B, Luo N, Gao R, Yu K, et al. Immune phenotypic linkage between colorectal cancer and liver metastasis. Cancer Cell (2022) 40(4):424–437.e425. doi: 10.1016/j.ccell.2022.02.013

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Shimasaki N, Jain A, Campana D. NK cells for cancer immunotherapy. Nat Rev Drug Discovery (2020) 19(3):200–18. doi: 10.1038/s41573-019-0052-1

CrossRef Full Text | Google Scholar

28. Song P, Li W, Guo L, Ying J, Gao S, He J. Identification and validation of a novel signature based on NK cell marker genes to predict prognosis and immunotherapy response in lung adenocarcinoma by integrated analysis of single-cell and bulk RNA-sequencing. Front Immunol (2022) 13:850745. doi: 10.3389/fimmu.2022.850745

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Li H, Wang C, Lan L, Behrens A, Tomaschko M, Ruiz J, et al. High expression of vinculin predicts poor prognosis and distant metastasis and associates with influencing tumor-associated NK cell infiltration and epithelial-mesenchymal transition in gastric cancer. Aging (2021) 13(4):5197–225. doi: 10.18632/aging.202440

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Xie MZ, Tang YP, Hu BL, Li KZ, Li JL, Liang XQ. Percentage of natural killer (NK) cells in peripheral blood is associated with prognosis in patients with gastric cancer: A retrospective study from a single center. Med Sci monitor Int Med J Exp Clin Res (2021) 27:e927464. doi: 10.12659/MSM.927464

CrossRef Full Text | Google Scholar

31. Sun C, Xu J, Huang Q, Huang M, Wen H, Zhang C, et al. High NKG2A expression contributes to NK cell exhaustion and predicts a poor prognosis of patients with liver cancer. Oncoimmunology (2017) 6(1):e1264562. doi: 10.1080/2162402X.2016.1264562

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Marin ND, Krasnick BA, Becker-Hapak M, Conant L, Goedegebuure SP, Berrien-Elliott MM, et al. Memory-like differentiation enhances NK cell responses to melanoma. Clin Cancer Res an Off J Am Assoc Cancer Res (2021) 27(17):4859–69. doi: 10.1158/1078-0432.CCR-21-0851

CrossRef Full Text | Google Scholar

33. Väyrynen JP, Haruki K, Lau MC, Väyrynen SA, Ugai T, Akimoto N, et al. Spatial organization and prognostic significance of NK and NKT-like cells via multimarker analysis of the colorectal cancer microenvironment. Cancer Immunol Res (2022) 10(2):215–27. doi: 10.1158/2326-6066.CIR-21-0772

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Wu M, Mei F, Liu W, Jiang J. Comprehensive characterization of tumor infiltrating natural killer cells and clinical significance in hepatocellular carcinoma based on gene expression profiles. Biomedicine pharmacotherapy = Biomedecine pharmacotherapie (2020) 121:109637. doi: 10.1016/j.biopha.2019.109637

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Cariani E, Pilli M, Barili V, Porro E, Biasini E, Olivani A, et al. Natural killer cells phenotypic characterization as an outcome predictor of HCV-linked HCC after curative treatments. Oncoimmunology (2016) 5(8):e1154249. doi: 10.1080/2162402X.2016.1154249

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Sheng J, Zhang J, Wang L, Tano V, Tang J, Wang X, et al. Topological analysis of hepatocellular carcinoma tumour microenvironment based on imaging mass cytometry reveals cellular neighbourhood regulated reversely by macrophages with different ontogeny. Gut (2022) 71(6):1176–91. doi: 10.1136/gutjnl-2021-324339

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Donne R, Lujambio A. The liver cancer immune microenvironment: Therapeutic implications for hepatocellular carcinoma. Hepatol (Baltimore Md) (2022). doi: 10.1002/hep.32740

CrossRef Full Text | Google Scholar

38. Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med (2018) 24(10):1550–8. doi: 10.1038/s41591-018-0136-1

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Hermann R, Krajcsi P, Fluck M, Seithel-Keuth A, Bytyqi A, Galazka A, et al. Review of transporter substrate, inhibitor, and inducer characteristics of cladribine. Clin Pharmacokinet (2021) 60(12):1509–35. doi: 10.1007/s40262-021-01065-3

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Bonate PL, Arthaud L, Cantrell WR Jr., Stephenson K, Secrist JA 3rd, Weitman S. Discovery and development of clofarabine: a nucleoside analogue for treating cancer. Nat Rev Drug Discovery (2006) 5(10):855–63. doi: 10.1038/nrd2055

CrossRef Full Text | Google Scholar

41. Chihara D, Arons E, Stetler-Stevenson M, Yuan CM, Wang HW, Zhou H, et al. Randomized phase II study of first-line cladribine with concurrent or delayed rituximab in patients with hairy cell leukemia. J Clin Oncol Off J Am Soc Clin Oncol (2020) 38(14):1527–38. doi: 10.1200/JCO.19.02250

CrossRef Full Text | Google Scholar

42. Saraceno L, Pirro F, Stigliano R, Agostoni EC, Protti A. Acute idiosyncratic liver injury after cladribine treatment for multiple sclerosis: first case report and review on associated hepatic disorders. Multiple sclerosis (Houndmills Basingstoke England) (2022) 28(13):2142–5. doi: 10.1177/13524585221125360

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Ertl IE, Lemberger U, Ilijazi D, Hassler MR, Bruchbacher A, Brettner R, et al. Molecular and pharmacological bladder cancer therapy screening: Discovery of clofarabine as a highly active compound. Eur Urol (2022) 82(3):261–70. doi: 10.1016/j.eururo.2022.03.009

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Khalafi S, Zhu S, Khurana R, Lohse I, Giordano S, Corso S, et al. A novel strategy for combination of clofarabine and pictilisib is synergistic in gastric cancer. Trans Oncol (2022) 15(1):101260. doi: 10.1016/j.tranon.2021.101260

CrossRef Full Text | Google Scholar

45. Lubecka K, Kaufman-Szymczyk A, Cebula-Obrzut B, Smolewski P, Szemraj J, Fabianowska-Majewska K. Novel clofarabine-based combinations with polyphenols epigenetically reactivate retinoic acid receptor beta, inhibit cell growth, and induce apoptosis of breast cancer cells. Int J Mol Sci (2018) 19(12):3970. doi: 10.3390/ijms19123970

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Yang J, König A, Park S, Jo E, Sung PS, Yoon SK, et al. A new high-content screening assay of the entire hepatitis b virus life cycle identifies novel antivirals. JHEP Rep Innovation Hepatol (2021) 3(4):100296. doi: 10.1016/j.jhepr.2021.100296

CrossRef Full Text | Google Scholar

47. Zhang H, Zhang Y, Dong J, Zuo S, Meng G, Wu J, et al. Recombinant adenovirus expressing the fusion protein PD1PVR improves CD8(+) T cell-mediated antitumor efficacy with long-term tumor-specific immune surveillance in hepatocellular carcinoma. Cell Oncol (Dordrecht) (2021) 44(6):1243–55. doi: 10.1007/s13402-021-00633-w

CrossRef Full Text | Google Scholar