NK cell–associated LncRNA signature predicts prognosis and immune landscape in hepatocellular carcinoma
Xiangyang Li, Jinze Li, Yuxuan Li, Yuhao Fan, Wei Tan, Fan Shi, Ya Zhu, Cheng Gong
TL;DR
This study creates a model using lncRNAs linked to NK cells to predict liver cancer outcomes and understand immune interactions.
Contribution
A novel NK cell-associated lncRNA model is developed for HCC prognosis and immune landscape analysis.
Findings
A four-lncRNA model robustly predicts HCC patient outcomes.
NK cells synergize with CD8+ T cells in the HCC tumor microenvironment.
New HCC molecular subgroups linked to immune profiles and drug sensitivity are identified.
Abstract
The purpose of this investigation was to develop and validate a risk model based on natural killer (NK) cell-associated lncRNAs to predict outcomes in individuals with hepatocellular carcinoma (HCC). To further explore the role of NK cells in the HCC tumor microenvironment, we leveraged single-cell RNA sequencing data and the TCGA-LIHC dataset to identify NK cell-associated lncRNAs. Using Cox regression and LASSO techniques, we pinpointed four key lncRNAs as prognostic markers for the model. The model demonstrated robust predictive power across the training set, validation set, and entire dataset. Additionally, we identified a synergistic interaction between NK cells and other immune cells, particularly CD8 + T cells, in HCC. Moreover, we uncovered novel molecular subgroups of HCC and their associations with the immune microenvironment and drug sensitivity. To further validate these…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
Figure 10
Figure 11
Figure 12
Figure 13
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8
Figure 9- —Hubei Province Natural Science Foundation of China
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsImmune Cell Function and Interaction · Cancer-related molecular mechanisms research · Cancer Immunotherapy and Biomarkers
Introduction
Primary liver cancer, a leading cause of cancer-related deaths, ranks among the top three in 46 countries and the top five in 90 countries worldwide. In 2020, approximately 905,700 people were diagnosed with liver cancer, and 830,200 individuals died from the disease [1]. Hepatocellular carcinoma (HCC), which constitutes the majority of primary liver cancers, remains particularly challenging to treat due to late detection and resistance to current therapies [2]. The recent emergence of immune checkpoint inhibitors (ICIs), particularly the combination therapy of Atezolizumab and the monoclonal antibody Bevacizumab, has set a new standard in HCC treatment, ushering in a new era for advanced-stage HCC management. Currently, there are at least 30 phase III immunotherapy trials for HCC at various stages, reflecting a dynamic and ongoing effort in the field [3].
Despite the significant promise that immunotherapy holds for individuals with HCC, the tumor immune microenvironment (TIME) remains a double-edged sword in cancer treatment. Accurate molecular characterization of TIME is essential for optimizing therapeutic strategies. Research has shown that the genomic features of inflammatory and non-inflammatory HCC correlate with responsiveness to ICIs, but clinical biomarkers to guide treatment remain limited [4]. To enhance therapeutic efficacy, further investigation is needed to identify these biomarkers. Immunotherapies primarily target T cells and include a range of treatments such as ICIs, adoptive cell transfer therapies, and chimeric antigen receptor (CAR) T-cell therapy [5]. CAR-T cell therapy involves genetically engineering immune cells to recognize and attack tumor cells, thereby enhancing the functionality and persistence of the immune response. However, CAR-T therapy faces challenges in solid tumors, including limited tumor infiltration, restricted permeability, and an immunosuppressive tumor microenvironment, all of which hinder its effectiveness [6]. In contrast, natural killer (NK) cells offer distinct advantages in immunotherapy, including superior safety profiles, diverse tumor-killing mechanisms, broad availability, and ease of acquisition [7]. While T-cell-based therapies have garnered significant attention, NK cells have been relatively overlooked. However, the growing recognition of NK cells’ potential in immunotherapy underscores the need for further exploration of their unique role in HCC and the identification of novel molecular targets for NK cell-based treatments.
Long non-coding RNA (lncRNA) is categorized as a strain of molecules exceeding 200 nucleotides in length that cannot be translated into functional proteins and exerts a pivotal role in regulating various biological processes [8]. Immune-related lncRNAs substantially contribute to the progression of cancer. Research has demonstrated that models involving immune-related lncRNAs reliably forecast the survival rates of individuals with breast cancer [9]. An innovative prognostic risk model, predicated on lncRNAs associated with immunity, is capable of forecasting the prognosis in pancreatic cancer individuals [10]. Models grounded in immune-related lncRNAs can be employed to gauge the survival rates of colon cancer patients [11]. Nevertheless, a paucity of studies exists regarding the function of NK cell-related lncRNAs in HCC. Through further investigation of lncRNAs associated with NK cells, we may potentially identify novel targets for HCC immunotherapy.
Through analysis of intercellular interactions and enrichment analysis, we revealed a potential collaborative interaction between NK cells and CD8 + T cells in HCC. The Tumor Immune Single-Cell Hub 2 (TISCH2) platform was utilized to establish NK cell-related genes (NKGs) that exhibit differences from other immune cells in HCC. We obtained the gene levels of NKGs from The Cancer Genome Atlas Liver Hepatocellular Carcinoma (TCGA-LIHC) dataset, and through correlation analysis with a correlation coefficient threshold set at 0.6, we identified lncRNAs linked to NK cells in HCC. Subsequently, Cox regression and the Least Absolute Shrinkage and Selection Operator (LASSO) were employed to determine four lncRNAs for constructing the predictive model, which received validation through both a validation set and the comprehensive dataset. Singular and multivariable independent prognostic analyses divulged that the risk score, as well as stage, are two distinct indicators. A predictive tool in the form of a nomogram was constituted to forecast the existence outcomes of individuals with HCC.
ROC curve analysis of our clinical data validated the independent factors. We performed an analysis to separate high- and low-risk subgroups identified by the model, followed by functional enrichment of differentially expressed genes. We predicted the pathways through which NK cell-related lncRNAs influence HCC. Tumor mutational burden (TMB) analysis revealed differentially mutated genes between the subgroups. Research on immune cells highlighted the importance of NK cells within the TIME of HCC, and drug sensitivity analysis identified personalized targeted therapies. Based on the model, we proposed a novel molecular classification of HCC and explored subtype variations in responses to immunotherapy and targeted therapies, informing personalized treatment. Finally, we validated the expression and function of previously unreported model lncRNAs, enhancing the model’s reliability.
Materials and methods
Data collection and treatment
Following the links from the Gene Expression Omnibus (GEO) website on the Transcriptomic-Immunohistochemical Single-Cell Hybridization (TISCH2) platform, we downloaded and analyzed the scRNA-seq data of HCC, GSE140228_10X, which encompasses 62,530 cells [12]. Principal Component Analysis (PCA) was recruited to uncover patterns and structures within the data, achieving dimensionality reduction. The K-nearest neighbors (KNN) algorithm and the Louvain algorithm were effectively used to distinguish and classify cell populations from the single-cell data, with precise cell type annotation. Wilcoxon testing was applied to identify genes specifically expressed in the NK cell population, using a logarithmic fold change of |fold change| ≥ 1.5 and a false discovery rate (FDR) of less than 0.05 as the criteria [13].
Cell-Cell interaction analysis
We utilized the built-in Cell Chat tool from the TISCH2 platform [14], we leveraged the Cell Chat tool within the TISCH2 platform to compute the expression levels of specific ligand-receptor pairs, thereby deducing the interactions between cells. Utilizing the R packages pheatmap and Cell Chat, we visualized the cell-cell interplay networks derived from scRNA. This approach enabled us to pinpoint significant ligand-receptor pairs within each cluster, classifying them of either originating or target cells, with a significance threshold of P < 0.05 [12].
Enrichment analysis and pathway identification of immune cells
A sequence of Gene Set Enrichment Analyses (GSEA) was executed through the TISCH2 platform to evaluate alterations in gene expression. This involved sorting the gene expression data by fold change to pinpoint gene sets that show significant expression variations within different cell populations [12]. By employing various enrichment analyses, we were able to uncover significant changes in biological pathways within specific cell populations. The FDR threshold was established at 0.05 to rigorously ensure the statistical reliability of our findings [12].
Collection and analysis of gene expression data
In our exploration, we collected genetic expression profiles, comprehensive clinical profiles, and somatic mutation documentation from HCC samples in the TCGA repository. Subsequently, a meticulous differential expression analysis of mRNA levels between tumor and non-tumor tissues was carried out. We established a correlation modulus of 0.6 to pinpoint lncRNAs closely associated with NK cells. We then utilized the R.limma package [15], specifically designed for gene expression data analysis, to discern lncRNAs related to NK cells that exhibited significant expression differences among tumor and non-tumor samples.
Construction and verification of the HCC prognostic model
To develop a predictive model for HCC, participants were evenly distributed into training and validation cohorts at a 1:1 ratio. We identified 54 lncRNAs tied to HCC prognosis within the training set through the univariate COX method. To avert overfitting of the model, we further employed the LASSO strategy to refine the selection of key variables. Based on these variables, we utilized the training dataset to establish a prognostic model via multivariable Cox regression analysis. A risk score for each NK cell-related lncRNA was determined by applying the following procedure: risk score= LINC01139 * 0.27799973053339 + AL354743.2 * 0.382637825857671 + NRAV* 0.599605592290264 + LINC02870 * 0.33791658169503. Individuals were broken down into high- and low-risk categories built on their calculated risk scores, with Kaplan-Meier existence analysis then applied to evaluate the existence disparities among the distinct risk subgroups. Additionally, We gauged the model’s predictive capabilities with ROC curves and substantiated its accuracy across the validation set and the full dataset [16].
Independent prognostic analysis and nomogram development
We initially identified candidate clinical factors associated with HCC prognosis using a univariate Cox strategy. These elements were then incorporated into a multivariate Cox strategy to determine the independent prognostic variables. To translate these prognostic factors into a clinically actionable tool, we developed a nomogram. The construction of the nomogram was facilitated by the nomogram function in the rms package, which transforms the modulus from the Cox strategy into an intuitive graphical format. The predictive precision of the nomogram was assessed by employing the concordance index [17].
Differential gene expression analysis and functional enrichment study among contrary risk subgroups
To identify genes that are differentially expressed (DEGs) among contrary risk subgroups, we employed the R.limma package with stringent selection criteria (absolute |logFC| exceeding 1, with an FDR below 0.05). We successfully identified genes with significant expression differences among contrary risk subgroups. Volcano plots and heatmaps were utilized to visually display the analysis results of these DEGs. To earn further foresight into the functions and biological significance of DEGs, we conducted a Gene Ontology (GO) analysis applying the R. clusterProfiler package, covering three dimensions: biological processes (BP), cellular components (CC), and molecular functions (MF), and explored KEGG pathways. We visually presented the significantly enriched biological pathways using bubble charts. Additionally, we employed GSEA to assess differences in biological functions between different risk groups, with selection criteria of absolute |NES| exceeding 1, and FDR below 0.05.
Tumor mutational Burden, immune Infiltration, and drug sensitivity analysis
We carried out an in-depth visualization evaluation of tumor mutational burden (TMB) using the maftools package in R. This tool enabled us to quantify the frequency and number of mutations in tumor samples and to create intuitive charts to display this data. A variety of computational tools, including TIMER, XCELL, CIBERSORT-ABS, CIBERSORT, EPIC, MCPCOUNTER, and QUANTISEQ, were employed to judge the plenitude of tumor-infiltrating immune cells in TCGA-LIHC specimens [18]. The R. limma package was utilized for differential expression exploration to pinpoint markedly divergent immune checkpoint genes among risk subgroups. Subsequently, the R. reshape function was employed to restructure the differential gene data for subsequent visualization. The ssGSEA method from the R. GSVA package was applied to estimate the activity of immune cell penetration and related immune functions from gene expression data. We further explored drug sensitivity analysis, using the calcPhenotype serve from the R.oncoPredict package to estimate drug response phenotypes, and the reshape2 package to transform the data for intuitive visualization [19].
Novel molecular subtyping of HCC
We utilized the Consensus ClusterPlus package in R to conduct consensus clustering analysis within unsupervised learning. Parameters were set as follows: max K for the range of potential cluster numbers, reps for the number of resampling iterations, pItem and pFeature for the proportions of sample and feature resampling, respectively, and seed to ensure the reproducibility of the results. We identified the likely cluster members and their numbers by conducting iterative sampling of gene expression data. The consensus matrix (CM), consensus cumulative distribution function (CDF) plots, and consensus heatmaps generated from this analysis are essential for determining the optimal number of subtypes in HCC classification [20].
Cell culture
HCC cell lines HuH-7, Hep3B, HepG2, and the normal human liver cell line MIHA were purchased from Shanghai Yaji Bio. All cell lines were cultured in RPMI-1640 medium (Sigma-Aldrich, St. Louis, MO, USA) supplemented with 10% fetal bovine serum (Gibco, USA) and 1% penicillin-streptomycin-glutamine (PSG; Thermo Fisher Scientific, Dreieich, Germany) at 37 °C with 5% CO2, under standard culture conditions.
Quantitative Real-Time PCR, plasmid Construction, and transfection
Total RNA was extracted using the Trizol method (T9108, Takara, Dalian, China), and reverse transcription was performed using an enzyme kit. Quantitative real-time PCR (qRT-PCR) was carried out using ChamQ Universal SYBR qPCR Master Mix (Q711-02, Vazyme, Nanjing, China). The primer sequences used for AL354743.2 were as follows: forward: AAGCTGCTCGAGACGTCAAA; reverse: CCGAAGCTGTCAGTGTCCTT.
For overexpression experiments, the full-length AL354743.2 sequence was amplified and cloned into the pcDNA3.1(+) expression vector using standard molecular cloning techniques. The primer sequences used for plasmid construction were as follows: forward: CGTTTGAAAAGGTGACAATAGCGTATGGC; reverse: ACAGGCACAGGCACCGTG. For gene knockdown, three specific short hairpin RNA (shRNA) sequences targeting AL354743.2 were designed and inserted into the pGPU6-GFP-Neo vector for transient transfection experiments. The shRNA sequences were as follows: shRNA-1: AGAGTGACCAGAGCGTGATTA; shRNA-2: ACACTAGCGACAGGCACAAAG; shRNA-3: CCACATCCAGAAGGTAGAAAG. Among these, the shRNA exhibiting the highest knockdown efficiency was selected for subsequent functional experiments.
Cells were seeded into 6-well plates at the appropriate density, and transfection was performed 24 h later at 37 °C using Lipofectamine 3000 (Sigma, USA) according to the manufacturer’s instructions. Cells were then used for subsequent experiments.
Cell proliferation, migration and invasion assays
As previously described [24–27], cell proliferation was evaluated using clonogenic and EdU incorporation assays. For colony formation assays, transfected and control HCC cells were seeded into 6-well plates at low density and cultured for 2–3 weeks. Colonies were fixed, stained with crystal violet, and counted to assess clonogenic capacity. For EdU incorporation assays, cells were incubated with EdU (10 µM) for 2 h according to the manufacturer’s instructions, followed by fixation, nuclear counterstaining, and fluorescence imaging to evaluate proliferative activity.
Cell migration was assessed using a wound-healing assay, and cell invasion was evaluated using Matrigel-coated Transwell chambers, as previously described [24–27]. Wound closure and invading cells were imaged and quantified under a microscope.
Data statistics
Within the scope of this research, all statistical reckon were conducted employing R utensil 4.4.1. We employed the Kruskal-Wallis assay to assess differences in immune checkpoint genes, immune scores, and drug sensitivities between risk groups. The log-rank assay with R.survival package was exerted to assess survival data to perform Kaplan-Meier survival analysis. Additionally, the Cox proportional hazards model was used to determine the effect of various variables on survival. All statistical assays were conducted two-tailed, with a P-value of beneath 0.05 indicating statistical significance. The significance strata are denoted as follows: *** indicates P < 0.001, ** indicates P < 0.01, and * indicates P < 0.05.
Results
Cell interaction and transcription factor analysis of NK cell subsets in HCC
The flowchart illustrates this study’s analysis and derivation process (Fig. 1). We mined the GSE140228_10X scRNA data in the TISCH2 platform and conducted visualization analysis. The UMAP plot revealed that NK cells are composed of clusters 5 and 12 (Fig. 2A). Pie charts and bar graphs of cellular composition from five patients showed the proportion and number of NK cells, indicating that NK cells are at a normal level (Figs. 2B and C). In TIME, the mutual regulation of immune cells plays a decisive character in the mission of immune cells and the alterations in the tumor’s immune status. The internal Cell Chat algorithm of the TISCH2 platform can calculate the interactions between different cell populations, revealing a strong correlation between NK cells and primarily CD8 + T cells, as well as monocytes, macrophages, and CD4 + T cells (Figs. 2D and Supplementary Fig. 1). We also analyzed the gene interactions between NK cells and other cell types when receiving signals or emitting signals (Supplementary Fig. 2A and 2B). It was found that NK cells interact with CD8 + T cells, CD4 + T cells, dendritic cells (DCs), and monocytes in multiple genes, revealing the synergistic effects between NK cells and T cell-centric immune cells in TIME. This powerful platform can also calculate enriched transcription factors in various cell clusters according to the spatial association algorithm. The expression patterns of NK cells in various transcription factors are almost identical to those of CD8 + T cells, CD4 + T cells, DCs, and monocytes, which echoes the previous results of immune cell interactions, further revealing a certain synergistic effect among NK cells and these cells (Supplementary Fig. 3A). At the same time, we identified the top ten enriched transcription factors in NK cells, which are BRD4, BCL3, E2F1, ETS1, TAF1, NOTCH1, TBP, MAF, ZMYND8, and KMT2A (Supplementary Fig. 3B and 3 C). The expression of these transcription factors may sway the mission of NK cells, thereby governing the tumor microenvironment of HCC. This provides new clues for our grasp of the immunological regimes of HCC. To study the impact of NK cells on HCC, we used the built-in Wilcoxon test in TISCH2 (|fold change| > 1.5; FDR < 0.05) to obtain NK cell-related genes (NKGs), and a volcano plot was served to manifest the expression of NKGs (Fig. 2E).
Fig. 1. Flowchart illustrating the construction process and subsequent analysis of the NK cell-associated lncRNA model. LIHC: Liver Hepatocellular Carcinoma; TCGA: The Cancer Genome Atlas; TISCH2: Tumor Immune Single-Cell Hub 2; TF: Transcription factor; LASSO: Least absolute shrinkage and selection operator; PCA: Principal component analysis; tSNE: T-distributed stochastic neighbor embedding; NK: Natural Killer; scRNA: Single-Cell RNA Sequencing; lncRNA: Long noncoding RNA; ROC: Receiver Operating Characteristic Curve
Fig. 2. Cellular interactions of NK cell subpopulations in LIHC. A: Distribution and abundance of different cell subpopulations in LIHC; B: NK cells as a percentage of all immune cells; C: Bar graph showing the distribution of NK cells in five patients; D: Cell chat algorithm calculates the probability of NK cells interacting with other cells and visualizes it; E: Volcano plot showing NK cell-related differential genes. Green indicates fold change < 1.5, FDR < 0.05; red indicates fold change > 1.5, FDR < 0.05
Enrichment characteristics of NK cells in HCC
The TISCH2 platform, with its robust enrichment analysis capabilities, has provided insights into how NK cells may manipulate HCC tumors. Initially, through KEGG enrichment analysis, we found that NK cells primarily upregulated routes such as NK cell-mediated cytotoxicity, toll-like receptors, phagocytosis, mitogen-activated protein kinases, and spliceosomes; while downregulating pathways like the ribosome, hematopoietic cell lineage, and systemic lupus erythematosus (Supplementary Figs. 4 and 5). Subsequently, we conducted GO analysis, revealing that NK cells upregulated cytoplasmic pattern recognition receptor signaling pathways in biological processes and ATPase regulatory activity in molecular functions, and downregulated nuclear transcription of mRNA degradation processes, antigen receptor-mediated signaling pathways, pattern recognition receptor signaling pathways, cytoplasmic ribosomes, endoplasmic reticulum lumen, clathrin-coated vesicle membranes, and clathrin-coated pits in cellular components, as well as antigen binding, ribosomal RNA binding, and cytokine activity in molecular functions (Supplementary Fig. 6–8).
Building and validating a HCC model originated from NK Cell-Associated genes
To analyze NK cells, the TCGA-LIHC dataset was operated to excavate the NKGs’ expression data. Performing Pearson correlation assay, we establish a correlation modulus above 0.6 and a p-value under 0.001 to detect lncRNAs correlated with NK cells. Using the R.limma package, we identified lncRNAs exhibiting expression variance between cancerous and healthy tissues (Fig. 3A). A heatmap vividly depicted the top 50 differentially expressed lncRNAs (Fig. 3B). Next, we implemented a univariate Cox regression assay within the training dataset, yielding 54 lncRNAs that impact HCC prognosis (Fig. 3C). The heatmap’s color distinctions vividly showcased the differences between these 54 lncRNAs among HCC specimens and normal samples (Fig. 4A). The LASSO regression method was employed to prevent model overfitting (Figs. 4B and C). We then incorporated the lncRNAs obtained from the LASSO method into a multivariate Cox model for computation. Ultimately, we identified LINC01139, AL354743.2, NRAV, and LINC02870 as the lncRNAs for model construction. The risk score formula is like so: risk score = LINC01139 * 0.27799973053339 + AL354743.2 * 0.382637825857671 + NRAV * 0.599605592290264 + LINC02870 * 0.33791658169503. Next, we stratified the training cohort into high-risk and low-risk subgroups under their risk scores. The heatmap displayed the levels of the above-mentioned 4 lncRNAs in the risk subgroups (Fig. 5A), and in the training set, HCC individuals with higher risk levels exhibited progressively poorer survival status (Figs. 5B and C). The Kaplan-Meier algorithm certified that the high-risk subgroup in the training set had poorer survival outcomes (Fig. 5D). The zone below the ROC curve values about 1, 2, and 3 years, 0.827, 0.723, and 0.689, individually, indicated that the NK cell-related lncRNA model performed well in the training set (Fig. 5E). We also validated the model using the test set (Figs. 5F-J) and the entire dataset (Figs. 5K-O), concluding comparable to those from the training set, further confirming the model’s effectiveness.
Fig. 3. Enrichment characteristics of NK cells in HCC. A: Volcano plot depicting NK cell-associated lncRNAs identified in HCC (red: logFC > 1, FDR-corrected p < 0.05; green: logFC < 1, FDR-adjusted p < 0.05); B: Heatmap depicting the top 50 NK cell-associated lncRNAs with the most significant differences; C: Forest plot depicting one-way Cox regression analysis yielding 54 lncRNAs associated with prognosis in HCC patients
Fig. 4. Screening of NK cell-associated lncRNAs for model construction. A: Heatmap demonstrating the difference in expression of lncRNAs obtained in Fig. 3C in HCC samples versus normal samples (*** represents p < 0.001, ** represents p < 0.01, * represents p < 0.05.); B, C: Lasso regression analysis was used to avoid model overfitting
Fig. 5. Construction and validation of HCC prognostic model. A: Expression of four lncRNAs in high and low-risk subgroups in the training set; B: Distribution status of HCC patients in the training set; C: Survival status of HCC patients in the training set; D: Survival curves of patients in the high- and low-risk subgroups in the training set; E: Evaluation of ROC curves in the training set; F: Expression of 4 lncRNAs in the high- and low-risk subgroups in the test set; G: Distribution of HCC patients in the test set. H: Survival status of HCC patients in the test set; I: survival curves of patients in the high- and low-risk subgroups in the test set; J: ROC curve assessment of the test set; K: Expression of 4 lncRNAs in the high- and low-risk subgroups in the full dataset; L: distribution of HCC patients in the full dataset. M: Survival status of HCC patients in the full dataset; N: survival curves of patients in the high- and low-risk subgroups in the full dataset; O: ROC curve assessment in the full dataset
Independent prognostic analysis and evaluation of clinical parameters
To assess the evaluative value of the NK cell-related lncRNA model, this study included common clinical parameters and risk scores as evaluation factors. Using the univariate Cox method and criteria of Hazard ratio > 1 and p < 0.05, we identified Stage as well as risk Score as statistically significant anticipating elements for HCC (Fig. 6A). Moreover, multivariate Cox regression endorsed that the above two factors remained independent prophetic factors (Fig. 6B). Additionally, the ROC curve demonstrated that among these parameters, Stage and risk score had the best evaluative performance, with AUC numbers of 0.717 and 0.711, individually (Fig. 6C). We used these two independent factors to construct a nomogram for better disease risk prediction (Fig. 6D). Figure 6E shows that the calibration curves during 1, 2, and 3 years are well-fitted to the ideal shape, and the C-index of 0.683 also indicates the effectiveness of the nomogram. Stratified analysis of various clinical parameters revealed that the risk score has predictive effects on survival in individuals aged ≤ 65 years, > 65 years, male individuals, and early and late Grade groups, as well as in early-stage individuals. For female patients and late-stage patients, although not statistically significant, there is a clear trend that the high-risk subgroup holds a weaker existence (Fig. 6F).
Fig. 6. Independent prognostic analyses and column line graphs. A: One-way COX regression approach to screen factors for HCC prognosis; B: Multifactor COX regression approach to screen independent factors for HCC prognosis; C: ROC curves demonstrating the predictive performance of each factor; D: Column line graphs predicting the prognosis of HCC; E: Calibration curves assessing the predictive power of column-line plots; F: Survival analysis of each clinical parameter
Differential gene expression and pathway analysis among risk subgroups
We initially filtered DEGs between the contrary risk subgroups using the guideline of |log FC| > 1 and FDR adjust p < 0.05. A volcano plot visually displayed the genes distinguishingly expressed between contrary risk subgroups (Fig. 7A). The top 50 most notable differential genes were displayed using a heatmap (Fig. 7B). We comported GO analysis and observed that the DEGs were predominantly enriched in various division processes within the cell cycle, including organelle division and nuclear division, as well as biological processes related to chromosomes and cellular matrices; they were mainly associated with cellular components such as microtubule binding and glycosaminoglycan binding (Figs. 7C and 8A). Additionally, we performed GSEA analysis, which indicated that the high-risk subgroup was chiefly enriched in passageways involving cytokine-cytokine receptor interaction, extracellular matrix-receptor interaction, neuroactive ligand-receptor interaction, cell cycle, and signaling channels in cancer. The low-risk subgroup was primarily enriched in metabolic passageways including primary bile acid synthesis and glycine, serine, and threonine metabolism(Fig. 8C and D). These pathways collectively suggest that NK cell-related lncRNAs may influence HCC prognosis by regulating tumor proliferation, cell cycle progression, and tumor–microenvironment interactions.
Fig. 7. Differentially expressed genes of high and low-risk groups. A: Differentially expressed genes in different risk groups are shown by volcano plots (red: logFC > 1, FDR-corrected p < 0.05; green: logFC < 1, FDR-adjusted p < 0.05); B: Heatmap effectively visualizes and demonstrates the distribution of differential genes across various hazard groups. C: Changes of differentially expressed genes in the graphene oxide pathway
Fig. 8. Functional enrichment analysis of high and low-risk groups. A: Bubble plot demonstrating the enrichment in the graphene oxide signaling pathway; B: Bubble plots demonstrating enrichment in the KEGG signaling pathway; C: GSEA analysis demonstrating up-regulated pathways in the high-risk group; D: GSEA analysis demonstrating down-regulated pathways in the low-risk group
Discovery of the kinship among neoplasm Mutations, immune Microenvironment, and prognosis in HCC
The overall tumor mutation rates for the high-risk and low-risk subgroups were almost identical (85.47% and 85.71%), with TP53 showing a weighty distinction in mutation frequency among the high-risk subgroup (37%) and the low-risk subgroup (16%), suggesting its role in adverse prognosis in HCC (Fig. 9A). We conducted survival analysis on subgroups with high or low tumor mutational burden (TMB) and found that high TMB corresponded to poorer survival outcomes (Fig. 9B and C). Furthermore, by combining risk scores with TMB, we classified patients into four subtypes and discovered that the subtype with high TMB plus high risk had a terrible existence, while the low-risk subgroup plus low TMB was the opposite, consistent with our previous results for TMB and risk subgroups. NK cells are a pivotal ingredient of the TIME, regulating immune activities in HCC. Our TIME scoring for HCC revealed that the high-risk subgroup had inferior estimated scores (P = 0.028) and stromal scores (P = 0.0018) compared to the low-risk subgroup, with a noticeable downward trend in immune scores, although not statistically meaningful, possibly due to the tiny sample amount (Fig. 10A). This indicates that HCC individuals from the high-risk subgroup may hold a relatively suppressive TIME, which could be tied to tumor immune evasion and terrible prognosis. Using various algorithms to calculate infiltrating immune cells, we found a notable opposing correlation among risk scores and NK cell infiltration, and similarly, CD8 + T cell infiltration trends were consistent with NK cells across various algorithms (Fig. 10B). Box plots from the ssGSEA algorithm confirmed reduced NK cell infiltration in the high-risk subgroup (Fig. 10C). Moreover, we inspected significant suppression of immune processes like cytotoxic activity, pro-inflammatory, T cell co-stimulation, and type I and II interferon responses in the high-risk group. Further validation confirmed the immunosuppressive state of the high-risk group (Fig. 10D). Analysis of immune checkpoint expression showed that gene expression levels were enhanced in the high-risk subgroup for almost all outcomes, indicating more immune evasion. We conducted a drug sensitivity analysis, where higher IC50 values indicate lower sensitivity. We uncovered that the high-risk subgroup had superior sensitivity to drugs including 5-Fluorouracil, Erlotinib, Pictilisib, and Dasatinib, while the low-risk subgroup was more sensible to Entospletinib, Gemcitabine, Oxaliplatin, and Fludarabine (Fig. 10F and Supplementary Fig. 9).
Fig. 9. Tumour mutation load in high and low-risk groups. A: Waterfall diagram showing the mutated genes in the high and low-risk groups; B: Prognosis of HCC patients in different tumor mutation load groups; C: Survival analysis showing the effect of combining tumor mutation load with risk score on the prognosis of HCC patients
Fig. 10. Immune and drug sensitivity analysis in high and low-risk groups. A: ESTIMATE Score, ImmuneScore, and StromalScore in high and low-risk groups. B: Multiple algorithms for quantitative analysis of immune infiltration in HCC; C: immune cell infiltration in high and low-risk groups; D: immune function in the high and low-risk groups; E: immune checkpoint profiles in the high- and low-risk groups; F: drug sensitivity profiles of the high and low-risk groups
Analysis of novel molecular subtypes in HCC
We conducted typing of HCC tumor samples and found that setting the test value κ to 4 resulted in the flattest distribution closest to the maximum CDF, with minimal cross-contamination between subtypes (Figs. 11A-D). Based on this, we named the four subtypes of HCC samples: C1, C2, C3, and C4. Survival analysis revealed that subgroup C2 had the most favorable prognosis, whereas subgroup C4 had the least favorable (Fig. 11E). A Sankey graph matched the new subtypes of HCC with contrary risk subgroups, showing that most of C1 was in the high-risk subgroup, most of C2 was in the low-risk subgroup, and both C3 and C4 were entirely in the high-risk subgroup (Fig. 11F). PCA and tSNE analyses confirmed the validity of the tumor classification and risk typing. Therefore, we believe that this predictive tool can categorize HCC individuals into new molecular subtypes (Figs. 11G-O).
Fig. 11. New molecular typing of HCC using NK cell-associated lncRNA models. A: Consensus matrix consisting of 4 typologies; B: Consensus CDF for different numbers of typing; C: CDF curves for different numbers of typing; D: Sample distribution for different numbers of subtypes; E: Survival curves for HCC subtypes; F: Sankey diagram of the relationship between HCC subtypes and risk scores; G-J: PCA and tSNE analyses demonstrating the distribution characteristics of different risk subgroups and different molecular subtypes
The relationship between HCC molecular Subtypes, the immune Microenvironment, and drug sensitivity
We delved deeper into the TIME of HCC to further validate the position of the NK cell-linked lncRNA model in the immunotherapy of new molecular subtypes of HCC. By measuring the ESTIMATE, Immune and Stromal Score, we found that subtype C3 had the lowest ESTIMATE Score and Immune Score, followed by subtype C4, with the lowest Stromal Score in C4 (Fig. 12A). This indicates that subtypes C3 and C4 have a high tumor fineness and stronger immune suppression, which corresponds with our prior outcome on prognosis. Based on the comprehensive results of different algorithms, we found that subtype C1 had the richest immune infiltration, while C3 had the weakest (Fig. 12B). Immune checkpoint inspection was exerted to discern the most suitable ICIs for different molecular subtypes, showing that, aside from TNFSF18, TNFSF9, and CD44 being highly expressed in subtype C3, the others were highly expressed in subtypes C1 or C4 (Fig. 12C), indicating a higher level of immune suppression and poorer prognosis in subtypes C1 and C4, and their potential suitability for corresponding ICI treatments. Building on the typing, we conducted drug sensitivity scrutiny, which exhibited that: individuals in subgroup C1 were most susceptible to Gallibiscoquinazole and Carmustine; individuals in subgroup C2 were most susceptible to Axitinib and AT13148; individuals in subgroup C3 were most susceptible to Ribociclib and SB216763; and individuals in subgroup C4 were most susceptible to Alpelisib and 5-Fluorouracil (Fig. 12D and Supplementary Fig. 10). Combining the above analysis, we found that the new molecular subtypes of HCC can provide precise individual guidance in the field concerning immunotherapy and targeted therapy for individuals, helping to improve the efficacy of HCC patients.
Fig. 12. Immunological and drug sensitivity analysis of new molecular typing of HCC. A: ESTIMATE Score, ImmuneScore, and StromalScore for new molecular typing of HCC; B: Quantitative analysis of immune infiltration of different algorithms for new molecular typing of HCC; C: Immune checkpoint analysis of different HCC typing; D: Drug sensitivity analysis of different HCC typing
Expression and functional validation of model LncRNAs in HCC
This study identified four key lncRNAs (LINC01139, NRAV, AL354743.2, and LINC02870) as prognostic markers for HCC. Previous studies have reported the expression and oncogenic roles of LINC01139, NRAV, and LINC02870 in HCC progression [24–27]. However, the biological function of AL354743.2 in HCC has not been previously characterized. To address this gap, we first examined the expression of AL354743.2 in the normal human liver cell line MIHA and in HCC cell lines HuH-7, Hep3B, and HepG2. The results showed that AL354743.2 expression was significantly higher in HCC cell lines than in normal liver cells (Fig. 13A). We next established an AL354743.2 overexpression model in HuH-7 cells (Fig. 13B). Functional assays demonstrated that AL354743.2 overexpression significantly enhanced clonogenic capacity, as assessed by colony formation assays (Fig. 13C). Consistently, EdU incorporation assays further confirmed that AL354743.2 overexpression markedly promoted HCC cell proliferation (Fig. 13D). In addition, Transwell and wound-healing assays revealed that AL354743.2 overexpression significantly increased the invasive and migratory abilities of HCC cells (Fig. 13E and F). Conversely, shRNA-mediated knockdown of AL354743.2 was successfully established in HCC cells (Fig. 13G). AL354743.2 knockdown significantly suppressed clonogenic capacity, as shown by colony formation assays (Fig. 13H). EdU incorporation assays further demonstrated a pronounced reduction in cell proliferation following AL354743.2 knockdown (Fig. 13I). Moreover, Transwell and wound-healing assays showed that silencing AL354743.2 markedly impaired the invasive and migratory capabilities of HCC cells (Fig. 13J and K).
Fig. 13. Expression and functional validation of AL354743.2 in HCC. A: The expression levels of AL354743.2 in the normal human liver cell line and different HCC cell lines; B: Overexpression of AL354743.2 in HCC cells achieved by transfection with an AL354743.2 overexpression plasmid; (C) Colony formation assay assessing clonogenic capacity following AL354743.2 overexpression; (D) EdU incorporation assay evaluating cell proliferation after AL354743.2 overexpression (scale bar: 50 μm); (E) Transwell assay assessing invasive capacity after AL354743.2 overexpression (scale bar: 100 μm); (F) Wound-healing assay evaluating cell migratory ability after AL354743.2 overexpression; (G) AL354743.2 knockdown in HCC cells achieved by transfection with shAL354743.2 plasmids; (H) Colony formation assay assessing clonogenic capacity following AL354743.2 knockdown; (I) EdU incorporation assay evaluating cell proliferation after AL354743.2 knockdown (scale bar: 50 μm); (J) Transwell assay assessing invasive capacity after AL354743.2 knockdown (scale bar: 100 μm); (K) Wound-healing assay evaluating cell migratory ability after AL354743.2 knockdown. *p < 0.05, **p < 0.01, ***p < 0.001
Discussion
Given the significant therapeutic potential that NK cells have demonstrated in metastatic and hematological malignancies, their role in HCC immunotherapy is increasingly gaining attention. The rapid and potent immune capabilities of NK cells position them as a promising tool in clinical treatment [21]. As a crucial constituent of the innate immune, NK cells play a substantial task in the surveillance against immune evasion. They possess unique functions that are irreplaceable by other immune cells, which include the ability to directly kill multiple neighboring cells with oncogenic markers. This capability encompasses several mechanisms: NK cells have natural cytotoxic activity that is independent of the major histocompatibility complex and antibodies; unlike cytotoxic T lymphocytes, NK cells can non-specifically recognize targets; and through the synergy of surface inhibitory and activating receptors, they directly target tumor cells [22]. Considering the pivotal task of NK cells in tumor immune regulation, we carried out an enrichment exploration of NK cells. In addition to the classic throughways like NK cell-mediated cytotoxicity, Toll-like receptors, MAPK, and spliceosomes identified in KEGG pathway enrichment, we also uncovered understudied pathways like phagocytosis, ribosomes, hematopoietic cell lineage, and systemic lupus erythematosus. In GO analysis, we found that NK cells activate the cytoplasmic pattern recognition receptor signaling pathway and ATPase regulatory activity while inhibiting post-transcriptional mRNA degradation, antigen receptor signaling pathways, pattern recognition receptor signaling pathways, and cellular components such as cytoplasmic ribosomes, endoplasmic reticulum lumen, clathrin-coated vesicle membranes, and coated vesicles. Furthermore, NK cells showed enrichment in molecular functions related to cytoplasmic translation, free ribosomes, RNA polymerase II, ribosome structure, antigen binding, hormone binding, and membrane receptor tyrosine kinase activity. These new findings provide significant clues for exploring the position of NK cells in tumor immunity. Likewise, we discovered a highly similar expression pattern among NK cells and CD8 + T cells, providing a molecular basis for studying their synergistic interactions.
LncRNAs, as pivotal gene regulatory factors in the immune system, highly specifically modulate various aspects of both innate and adaptive immunity [23]. We employed scientifically selected four prognostic lncRNAs to construct a risk model. Among them, NRAV encourages liver cancer cell proliferation and invasion across modulation of the Wnt/β-catenin pathway [24]. LINC01139 is a crucial oncogene-associated lncRNA that is elevated in various solid and liquid tumors [25]. LINC01139 can also promote the progression of HCC by scrambling for the miR-30 family to upregulate MYBL2 [26]. LINC02870 facilitates the translation of SNAIL, thereby advancing the progression of hepatocellular carcinoma [27]. No prior studies have investigated AL354743.2, making this the first study to experimentally verify its association with HCC. Our results show that AL354743.2 is significantly overexpressed in HCC cell lines. Furthermore, its overexpression or knockdown notably enhances or inhibits HCC cell proliferation, invasion, and migration. Moreover, the interactions between these lncRNAs and NK cells have not been documented. Studying these molecules may reveal the link between NK cells and HCC. Rested on the risk score formula of this model, we isolated the population into contrary risk groups in a 1:1 ratio. We conducted enrichment exploration on distinguishingly expressed genes among subgroups. In our study, we observed that these genes are primarily engaged in key procedures of cell cycle adjustment and tumor progression, like the cell cycle, PI3K-Akt、Hippo pathway, and other signaling pathways revealed by KEGG analysis, as well as pathways connected to cytokine-cytokine receptor interaction and neuroactive ligand-receptor interaction identified by GSEA analysis. These findings uncover the mechanisms by which NK cell-associated lncRNAs regulate cell cycle control and signal transduction in tumor cells. Additionally, we identified pathways that have been less studied in the field of cellular immunity, including organelle division, nuclear division, and other cellular biological processes, together with molecular functions including microtubule binding and glycosaminoglycan binding. The discovery of these pathways provides new directions for our future research in HCC immunology and lncRNA function, particularly in areas including cell cycle regulation, cell-matrix interactions, and metabolic pathways.
This study indicates that immune cell infiltration is less pronounced in the high-risk subgroup, suggesting the presence of immune suppression in the high-risk group. The infiltration of NK cells is oppositely interrelated with the risk score and is compatible with the infiltration of CD8 + T cells. This proposes that NK cell infiltration is tied to better prognosis and that there is a synergistic effect between NK cells and CD8 + T cells. Mechanistically, NK cells and CD8 + T cells exhibit complementarity in tumor immunity. Their synergy can enhance the immune response. Most importantly, NK cells significantly participate in the recruitment of dendritic cells to tumors, which in turn boosts the activation of CD8 + T cell responses, and IL-2, released by T cells, further activates NK cells [28]. The findings from this study highlight the synergistic effects between the two cell types. This is consistent with previous research, which manifests that the TIPE2 molecule as a checkpoint factor for the auxiliary function of NK cells. Selective elimination of TIPE2 in NK cells can concurrently boost the antitumor activities of NK cells and the CD8 + T cells they support [29]. The genes predicted in this study to be linked to the synergy between NK cells and CD8 + T cells may potentially serve as targets for their mutual assistance and require further validation through multicenter experiments. As a rising contender in cell therapy following T cell therapy, despite its various advantages, our focus should also be on enhancing the cytotoxicity of NK cells and prolonging their activity [30]. By studying lncRNAs associated with NK cells, we may be able to identify effective targets to overcome the current challenges faced by NK cells.
Tumor checkpoint gene expression correlates with the tumor’s immune-suppressive state, and the use of ICIs can ameliorate this condition. The most well-known are the blockers of PD-1/PD-L1 and CTLA-4, which have shown positive clinical outcomes. Some of these drugs have been approved for specific patient populations, while others are still in clinical trial stages [31]. In this study, the high-risk subgroup exposed superior expression of almost all distinctively expressed immune checkpoint genes compared to the low-risk subgroup, suggesting a greater rank of immune evasion. However, this also indicates that ICI therapy might be particularly important for high-risk individuals, who are liable to own a higher reaction to treatment. Among these genes is the well-known CTLA-4, which is chiefly expressed during T-cell activation. It competitively binds to CD80 and CD86, similar to CD28. Extensive research has confirmed its role in various cancers, including HCC [32]. However, the roles of the majority of these checkpoint genes in HCC still require further validation.
The NK cell-related lncRNA signature identified in this study demonstrates robust predictive power. Individuals were disassociated into two subgroups through the model’s formula. Results indicate a tendency towards poorer prognosis in the high-risk subgroup. The model’s outcomes were visualized through risk heatmaps and Kaplan-Meier survival assay, and the precision of the findings was substantiated with ROC analysis. Further validation was carried out with a test set and the entire dataset to guarantee the authenticity and reproducibility of the findings. The model, in conjunction with the stage, was identified as an independent prognostic factor through independent prognostic analysis, leading to the development of a nomogram. Moreover, the ROC curves for clinical parameters indicated that these two independent prognostic factors significantly outperformed other clinical parameters.
Personalized medicine is a hallmark of modern medical advancement, and classifying tumors based on their characteristics can achieve the goal of precise treatment. We categorized HCC individuals into four subgroups based on the model characteristics. It is observable that subgroups C3 and C4 are entirely within the high-risk group, and correspondingly, these two subgroups have poorer prognoses, with C4 being the worst, presumably due to the elevated levels of immune checkpoint gene expression in subgroup C4. Conversely, subgroup C2 is largely situated within the low-risk group and consequently has the best prognosis. Immune-related scoring indicates that subgroups C3 and C4 have the lowest scores, suggesting a higher malignancy degree in these tumor subgroups. The infiltration of immune cells also shows that subgroup C3 has the weakest immune infiltration. What’s more, the exploration of immune checkpoint results indicates that subgroup C2 has the lowest expression across all immune checkpoint genes. Synthesizing all the above immune-related analyses, their prognostic outcomes are compatible with the previous typing consequences, validating the productiveness of this classification for NK cell immune-related analyses. Moreover, gleaned from the typing outcomes, we explored the most sensitive drugs for each subgroup. Considering the outcomes of immune-related analyses and drug sensitivity, we have reason to believe that this classification is beneficial for the customized treatment of HCC individuals and can refine the potency of immunotherapy and targeted treatments.
Despite the promising performance of the NK cell-related lncRNA risk model constructed in this study for predicting the prognosis of HCC, several limitations should be acknowledged: Firstly, the model was primarily built based on retrospective data from the TCGA database. Future work should involve prospective studies and larger clinical samples to further validate the robustness of the model. Secondly, this study focused mainly on bioinformatics analyses and in vitro experimental validation. We have not yet explored the molecular mechanisms of lncRNAs through in vivo experiments (e.g., tumor-bearing mouse models) or functional studies (e.g., co-culture experiments between NK cells and HCC cells). For instance, the roles of key lncRNAs such as LINC01139 and NRAV in the synergy between NK cells and CD8 + T cells need to be clarified through gene knockout/overexpression experiments to elucidate their regulatory networks. Moreover, the high heterogeneity of HCC tumors may lead to differences in the predictive efficacy of the model across various molecular subtypes or clinical characteristics (e.g., different etiologies, presence of cirrhosis). Although we preliminarily explored the differences between subtypes through molecular classification (C1-C4), further investigation is needed to combine single-cell sequencing or spatial transcriptomics to dissect the dynamic changes in the tumor microenvironment and their impact on the model.
Conclusion
The risk assessment model expanded in this research, associated with NK cell-lncRNAs, offers robust predictive capabilities for the prognosis of HCC individuals. And has unveiled the synergistic interactions among NK cells and immune cells such as CD8 T cells within the tumor microenvironment. In addition, the study has identified new molecular subtypes of HCC and explored their relationship with the TIME and drug sensitivity, offering avenues for the personalized treatment of HCC.
Supplementary Information
Below is the link to the electronic supplementary material.
Supplementary Material 1
The reference list from the paper itself. Each links out to its DOI / PubMed record.