Bulk and single-cell RNA sequencing identify prognostic signatures related to FGFBP2⁺ NK cell in hepatocellular carcinoma

Data collection and preprocessing

RNA sequencing data and clinical follow-up data of HCC samples were downloaded from The Cancer Genome Atlas (TCGA) database through Genomic Data Commons (GDC) Application Programming Interface (API). The ICGC-LIRI-JP dataset was collected from the HCCDB database. Then, after deleting samples that did not have clinical follow-up data or status, the Ensembl was transformed to a gene symbol, and the expression average value was taken for multiple gene symbols. Finally, 370 tumor samples and 50 control samples were acquired in TCGA-LIHC cohort, utilizing as the training set. The ICGC-LIRI-JP cohort contained 212 HCC samples, which was served as the validation set.

The scRNA-seq data of GSE162616 dataset, including three HCC samples and three healthy liver samples, was derived from the Gene Expression Omnibus (GEO) database. For filtering the scRNA-seq data, each gene was set to be expressed in a minimum of three cells, and each cell expressed at least 200 genes. Then, the cells with nFeature_RNA >300, nCount_RNA >3,000, and mitochondrial gene expression percent.mito <25% were reserved. The NormalizeData function was applied for log conversion, and the FindVariableFeatures function was used to identify highly variable genes. Next, the ScaleData function was utilized to normalize the expression values of all genes, the RunPCA function was employed for principal component analysis (PCA), and the “harmony” R package was used to remove batch effects between different samples (Zulibiya et al., 2023).

Cell clustering analysis

Cell clustering analysis of GSE162616 dataset was conducted using the “Seurat” R package (Du et al., 2024). Firstly, UMAP was conducted on the top 15 principal components (PCs) for dimensionality reduction. Then, the cells were clustered using the FindNeighbors and FindClusters functions, with cluster resolution of 0.1 for all cells and 0.3 for NK cells. Finally, cell types were annotated based on the marker genes provided by CellMarker2.0 database (Lu et al., 2024).

Single-sample gene set enrichment analysis

The correlation between the prognosis of HCC and each type of NK cells was examined based on the single-sample gene set enrichment analysis (ssGSEA) score calculated utilizing “GSVA” R package for the tumor samples in TCGA-LIHC cohort (Wang et al., 2023). The surv_cutpoint function was applied to find the optimal cut-off point, and HCC patients were divided into low- and high-score groups.

Pseudo-time trajectory analysis

In order to invesitgate the role of FGFBP2⁺ NK cell in the progression of HCC, Monocle2 was used to construct pseudo-time trajectory (Tao et al., 2024). The cds object was established by newCellDataSet function, and the genes expressed in less than 10 cells were filtered out. The differentially expressed genes (DEGs) between HCC samples and healthy samples were identified by differentialGeneTest function, and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis was conducted on these DEGs using the “ClusterProlifer” R package (Zhang et al., 2024). Subsequently, the reduceDimension function was employed for dimensionality reduction (max_components = 2, method = “DDRTree”), and the orderCells function was applied to order the cells and construct the pseudo-time trajectory of FGFBP2⁺ NK cell.

Cellular communication analysis

The “CellChat” R package was employed to conduct cell–cell interaction analysis in HCC (Lian et al., 2024). The number of interactions and interaction weights/strength was analyzed, as well as the key ligand–receptor pairs between FGFBP2⁺ NK cell and other cell subpopulations were visualized by a bubble plot.

Establishment and verification of RiskScore model

Univariate Cox regression analysis (p < 0.05) was performed on the marker genes of FGFBP2⁺ NK cell. The number of genes in the model was reduced by LASSO Cox regression analysis was conducted using the “glmnet” R package (Zeng & Chen, 2024). By 10-fold cross validation, we selected the optimal lambda value as the result of LASSO regression for subsequent analysis. Further, stepwise regression analysis was performed, and the prognostic gene signature related to FGFBP2⁺ NK cell in HCC was identified to establish RiskScore model. The RiskScore of each patient in TCGA-LIHC cohort was obtained according to the following formula (Li et al., 2024a): ${\displaystyle \text{RiskScore}=\sum \mathrm{\beta }\text{i*ExPi}}$

β is the coefficient of gene in Cox regression model, and i is the gene expression.

Z-score was utilized for standardization, and based on the threshold of RiskScore = 0, high- and low-risk groups were classified. To evaluate the prognostic prediction performance of RiskScore model, receiver operating characteristic (ROC) analysis was conducted using the “timeROC” R package (Lin et al., 2023). Additionally, the robustness of RiskScore model was validated in the ICGC-LIRI-JP cohort.

Immune cell infiltration and immunotherapy response analysis

Immune cell infiltration was assessed between high- and low-risk groups in TCGA-LIHC cohort. StromalScore, ImmuneScore, and ESTIMATEScore were calculated by the “estimate” R package (Dong et al., 2023). Tumor Immune Estimation Resource (TIMER) was applied to assess the infiltration of six immune cells (Xiao et al., 2020). Microenvironment cell populations-counter (MCP-counter) algorithm was utilized to assess the infiltration of 10 immune cells (Chen et al., 2022).

The response of different risk groups to immunotherapy was predicted by TIDE algorithm (Li et al., 2024b). Exclusion, Dysfunction, and TIDE scores of different risk groups were calculated in TCGA-LIHC cohort. Moreover, we analyzed the relationship of RiskScore and 9 immune checkpoint genes.

Correlation analysis between RiskScore and drug sensitivity

The IC₅₀ values of drugs for HCC patients in TCGA-LIHC cohort were predicted using the “oncoPredict” R package (Maeser, Gruener & Huang, 2021). Then, the association between RiskScore and drug sensitivity was analyzed (p < 0.05 and |cor| > 0.3).

Cell culture and transfection

Human hepatic astrocyte cell line LX-2 (BNCC337957) and HCC cell line Huh7 (BNCC337690) were acquired from the BeNa Culture Collection (BNCC) Biotechnology Co. (Beijing, China). Then, LX-2 and Huh7 cell lines were separately cultivated in RPMI-1640 (BNCC338360) and DMEM (BNCC363314) contained 10% fetal bovine serum (FBS), and were all incubated in the condition of 5% CO₂ and 37 °C. Hereafter, according to the protocol of Lipofectamine 2000 (Invitrogen, Waltham, MA, USA), Huh7 cells were transfected with the small interfering (si) RNA of FTL (si-FTL: 5′-TCCCAGATTCGTCAGAATTATTC-3′, Sangon, China), (si) RNA of PTP4A2 (si-PTP4A2: 5′-GAGGTTCTATGTGCCATAATTAA-3′, Sangon, China) and negative control (si-NC). We have performed short tandem repeat (STR) identification on the cells, and the mycoplasma detection results turned out to be negative.

Quantitative real-time PCR

Total RNA of LX-2 and Huh7 cells were acquired employing the Trizol reagent (B610409, Sangon, China). Then, the complementary DNA (cDNA) was synthesized through reverse transcription applying the RevertAid First Strand cDNA Synthesis Kit (B300538, Sangon, China). Thereafter, quantitative real-time PCR (qRT-PCR) amplification was performed using the SYBR Green (B110031, Sangon, China) on the basis of manufacturer’s instructions. The qPCR conditions were: 94 °C for 30 s first, then 40 cycles of 94 °C for 5 s and 60 °C for 30 s. The primer pairs of this study were shown in Table S1. GAPDH was utilized as an internal control to normalize the relative mRNA expressions of each gene by 2^−ΔΔCT method (Xu et al., 2020b).

Proliferation assay

Huh-7 cells that had been transfected were maintained for 48 h until they reached the logarithmic growth phase, after which they were moved into 96-well plates. The EdU Cell Proliferation Assay Kit (RiboBio, Guangzhou, China) was utilized to evaluate cell proliferation. In accordance with the established protocol, the cells underwent staining, followed by examination and imaging using a fluorescence microscope (Nikon, Toyko, Japan). EdU-positive cell counts were conducted using ImageJ software.

Wound healing assay

The effect of FTL and PTP4A2 silencing on the migration of HCC cells Huh7 was measured by wound healing assay (Song et al., 2021). The transfected Huh7 cells (1 × 10⁵) were inoculated into a 6-well plate and grown overnight. Afterwards, the wound was created using a sterile 20 µL pipette tip and Huh7 cells were sustainably cultured in serum-free medium for 48 h. Finally, the pictures of wound areas at 0 hour (h) and 48 h were obtained under a microscope (ECLIPSE Ei, Nikon, Tokyo, Japan), and the wound closure (%) of Huh7 cells was estimated with the ImageJ2 software.

Transwell assay

Transwell assay was conducted to examine the influence of FTL and PTP4A2 silencing on the invasion of HCC cells Huh7 (Huang et al., 2021b). The diluted Matrigel (BD Biosciences, Franklin Lakes, NJ, USA) was pre-coated into the Transwell chamber (8.0 µm, Corning, Corning, NY, USA). Next, the transfected Huh7 cells (1 × 10⁵) were starved in the upper chamber encompassing 200 µL non-serum medium. 500 µL DMEM containing 10% FBS was filled into the lower chamber. After 48 h of incubation, the invaded Huh7 cells were fixed by 4% paraformaldehyde (YTB1299, bjbalb, Beijing, China) and dyed by 0.1% crystal violet (YT913, bjbalb, Beijing, China). Further, the number of invaded Huh7 cells was counted under the same microscope as above.

Statistical analysis

The bioinformatic analysis was conducted using R programing language (version 4.1.0). The differences between different groups were compared by the Wilcoxon rank-sum test. Kaplan–Meier (K-M) analysis of overall survival (OS) was conducted by log-rank test. The Spearman method was utilized for correlation analysis. All experimental data of independent triplicates were expressed as mean ± standard deviation,, and statistical analysis was carried out by GraphPad Prism8.0. For data that did not conform to a normal distribution, non-parametric tests such as the Mann–Whitney U test or Kruskal-Wallis test were applied. p < 0.05 signified statistical significance.

Single-cell atlas of HCC revealed the reduction proportion of NK cell

The scRNA-seq data of GSE162616 dataset was analyzed to delineate the single-cell atlas of HCC. After cell filtration, standardization, dimensionality reduction, and clustering, a total of 47,550 cells were obtained and divided into 11 cell clusters (Fig. S1). Using the marker genes, these cells were further annotated as seven cell types (Fig. 1A), comprising NK cell (GZMB, GZMH, FGFBP2, GNLY, CD160), T cell (CD2, CD3D, CD3E, IL7R), macrophage (AIF1, MS4A7, LYZ), plasma cell (DERL3, IGHG1, MZB1), CD8⁺ T cell (MKI67, STMN1), B cell (CD79A, MS4A1), and hepatocyte (ALB, APOA1, KRT18) (Fig. 1B). Additionally, the proportion of each cell type in different samples was analyzed, showing that the number of NK cell in HCC samples was markedly decreased compared with healthy samples (Figs. 1C and 1D). These results indicated that the proliferation and activity of NK cell may be inhibited in HCC microenvironment.

FGFBP2⁺ NK cell was associated with the prognosis of HCC patients

NK cell as an important immune cell can recognize and kill tumor cells. Hence, we extracted the NK cell from HCC tissues for re-clustering, obtaining three cell subpopulations (Fig. 2A). The top20 highly expressed genes in each NK cell subpopulation were displayed by a heatmap (Fig. 2B). Using the marker genes, these cell subpopulations were annotated as FGFBP2⁺ NK cell (FGFBP2, FCGR3A, GZMB, GZMH), CD160⁺ NK cell (GZMK, CD160), and IL7R⁺ NK cell (IL7R, SELL, LMNA, CD44) (Fig. 2C). Compared with healthy samples, the proportion of FGFBP2⁺ NK cell and IL7R⁺ NK cell was elevated while CD160⁺ NK cell was decreased in HCC samples (Fig. 2D). Furthermore, K-M survival curve suggested that FGFBP2⁺ NK cell exhibited an association with the prognosis of HCC patients (p =0.047), with higher survival probability in high ssGSEA score group in comparison to the low ssGSEA score group (Fig. 2E). Whereas, CD160⁺ NK cell (p =0.082) and IL7R⁺ NK cell (p =0.1) was not closely related to the prognosis of HCC patients (Fig. S2).

Pseudo-time trajectory of FGFBP2⁺ NK cell from normal to HCC was constructed

For further exploring the role of FGFBP2⁺ NK cell in HCC progression, the pseudo-time trajectory of FGFBP2⁺ NK cell from normal to HCC was constructed by Monocle2, and the branch with more healthy samples was set as the starting point, and the branch with more HCC samples was set as the ending point (Figs. 3A and 3B). Moreover, DEGs analysis between HCC samples and healthy samples revealed two differential expression gene clusters, and the gene expression in Cluster1 was gradually up-regulated with the increase of pseudo-time (Fig. 3C). KEGG enrichment analysis demonstrated that the DEGs in Cluster1 were primarily involved in the endocytosis, protein processing in nuclear factor (NF)-kappa B signaling pathway, endoplasmic reticulum, apoptosis, antigen processing and presentation, NK cell mediated cytotoxicity (Fig. 3D). These pathways might play a crucial role in the process of FGFBP2⁺ NK cells transitioning from the normal state to be involved in HCC.

FGFBP2⁺ NK cell exhibited extensive intercellular communication in HCC

Cellular communication analysis was performed to explore the potential interactions between different cell types in HCC. It was observed that the number and strength of ligand–receptor interactions was complicated (Fig. 4A). FGFBP2⁺ NK cell as signal sender showed extensive communication with Macrophage, B cell, IL7R⁺ NK cell, T cell, CD160⁺ NK cell, and CD8⁺ T cell (Fig. 4B). By extracting the critical ligand–receptor pairs, we found that FGFBP2⁺ NK cell communicated with other cell subpopulations via tumor necrosis factor (TNF)-TNFRSF1B and macrophage migration inhibitory factor (MIF)-(CD74+CD44) (Fig. 4C), while other cell subpopulations communicated with FGFBP2⁺ NK cell via TNFSF14-TNFRSF14 and MIF-(CD74+CXCR4) (Fig. 4D).

RiskScore model was constructed and verified based on 8-genes prognostic signature

Firstly, the marker genes in FGFBP2⁺ NK cell were subjected to univariate Cox regression analysis (p < 0.05). LASSO and stepwise regression analysis was further performed to reduce the gene number, and the model was optimal when the lambda was 0.0286 (Figs. 5A and 5B). Then, eight prognostic signatures related to FGFBP2⁺ NK cell in HCC were identified, including six “risk” genes (UBE2F, AHSA1, PTP4A2, CDKN2D, FTL, RGS2) and two “protective” genes (KLF2, GZMH) (Fig. 5C). Next, we established a RiskScore model of “RiskScore = −0.352*GZMH−0.303*KLF2+0.153*FTL+0.312*PTP4A2+0.442*UBE2F+ 0.234*CDKN2D+0.141*RGS2+0.427*AHSA1”.

Furthermore, according to the threshold of RiskScore = 0, all samples in TCGA-LIHC cohort were separated into high- and low-risk groups (Fig. 5D). The performance of RiskScore model was evaluated by the area under ROC curve (AUC). It was showed that 1-, 3-, 5-years AUC values of the RiskScore model were 0.78, 0.77, 0.76 (Fig. 5E), which manifested that the RiskScore model exhibited good prognostic prediction performance. K-M curve analysis suggested that the OS rate of high-risk group was lower (Fig. 5E), showing that HCC patients with high RiskScore may have a poor prognosis. The robustness of RiskScore model was verified in ICGC-LIRI-JP cohort, and the results of which were similar to the TCGA-LIHC cohort (Figs. 5F and 5G). These outcomes demonstrated the reliability of the RiskScore model in the prognostic predicting of HCC.

RiskScore model exhibited potential in predicting immunotherapy response for HCC

We first elucidated the association between RiskScore and immune cell infiltration. ESTIMATE algorithm results showed that compared to low-risk group, high-risk group had lower StromalScore, ImmuneScore, and ESTIMATEScore (Fig. 6A). Based on the TIMER database, the infiltration levels of neutrophil, B cell, macrophage, and dendritic cell (DC) in high-risk group were significantly higher, while CD8 T cell was notably lower in comparison to the low-risk group (Fig. 6B). MCP-counter algorithm suggested that the infiltration levels of CD8 T cell, cytotoxic lymphocytes, NK cell, and endothelial cell were markedly lower, yet monocytic lineage was higher in high-risk group (Fig. 6C). In addition, the immunotherapy response between different risk groups was predicted in TCGA-LIHC cohort. In comparison to low-risk group, high-risk group exhibited higher exclusion score and TIDE score (Fig. 6D), demonstrating that high-risk HCC patients may be less likely to benefit from taking immunotherapy. RiskScore was positively correlated with several immune checkpoint genes, including CD44, CD276, CD80, LGALS9, and CTLA4 (Fig. 6E), showing that HCC patients with higher RiskScore may be more possibly to experience immune escape. Furthermore, we examined the association between drug sensitivity and RiskScore, and screened 16 drugs (such as Doramapimod, Nutlin.3a, Selumetinib, Sepantronium bromide, Tozasertib) that showed significant correlation with RiskScore (p < 0.05 and —cor— > 0.3) (Fig. 6F).

In vitro HCC cell-based model to validate key genes

The qRT-PCR analysis demonstrated that in comparison to the human hepatic astrocyte cells LX-2, the relative mRNA expression levels of GZMH was remarkably lower, yet FTL, PTP4A2, UBE2F, CDKN2D, RGS2, and AHSA1 were notably higher in HCC cells Huh7 (Fig. 7A). Since the pro-carcinogenic role of FTL in HCC has been reported in related studies, and the functional mechanism of PTP4A2 in HCC lacks systematic studies, we chose these two genes for in vitro functional experiments to further validate their roles in HCC cell migration, invasion and proliferation. Thereafter, we verified the knockdown efficiency of these two genes in Huh-7 cells (Fig. S3). We observed a significant decrease in proliferation of HCC cell lines after silencing FTL and PTP4A2 (Fig. 7B). In addition, wound healing assay displayed that FTL and PTP4A2 silencing could decrease the wound closure rate of Huh7 cells (Fig. 7C). Transwell assay suggested that the number of invading Huh7 cells was reduced by the silencing of FTL and PTP4A2 (Fig. 7D). Consequently, these results indicated the crucial involvement of FTL and PTP4A2 in HCC progression.