
==== Front
Heliyon
Heliyon
Heliyon
2405-8440
Elsevier

S2405-8440(24)12592-3
10.1016/j.heliyon.2024.e36561
e36561
Research Article
Prognostic model for hepatocellular carcinoma based on necroptosis-related genes and analysis of drug treatment responses
Wu Ronghuo a
Deng Xiaoxia b
Wang Xiaomin c
Li Shanshan lishanshan33@ylu.edu.cn
c⁎⁎
Su Jing amandasu@sr.gxmu.edu.cn
de⁎
Sun Xiaoyan missharon@qq.com
f⁎⁎⁎
a Department of Economics, Jinan University, Guangzhou, 510632, China
b School of Mathematics and Statistics, Yulin Normal University, Yulin, 537000, China
c Department of Biology and Pharmacy, Yulin Normal University, Yulin, 537000, China
d Schoole of Information and Management, Guangxi Medical University, Nanning, 530002, China
e Faculty of Data Science, City University of Macau, Macao, Macao SAR, China
f Human Resources Office, Guangxi Medical University, Nanning, 530002, China
⁎ Corresponding author. Schoole of Information and Management, Guangxi Medical University, Nanning, 530002, China. amandasu@sr.gxmu.edu.cn
⁎⁎ Corresponding author. lishanshan33@ylu.edu.cn
⁎⁎⁎ Corresponding author. missharon@qq.com
22 8 2024
15 9 2024
22 8 2024
10 17 e365619 4 2024
16 8 2024
19 8 2024
© 2024 The Authors. Published by Elsevier Ltd.
2024

https://creativecommons.org/licenses/by-nc/4.0/ This is an open access article under the CC BY-NC license (http://creativecommons.org/licenses/by-nc/4.0/).
Objective

Recent studies reveal that necroptosis is pivotal in tumorigenesis, cancer metastasis, cancer immunity, and cancer subtypes. Apoptosis or necroptosis of hepatocytes in the liver microenvironment can determine the subtype of liver cancer. However, necroptosis-related genomes have rarely been analyzed in hepatocellular carcinoma (HCC). Therefore, this study aims to construct an HCC risk scoring model based on necroptosis-related genes and to validate its predictive performance in overall survival prediction and immunotherapy efficacy evaluation in HCC, as well as to analyze drug treatment responses.

Methods

This study analyzed clinical information and RNA-seq expression data of liver cancer patients from TCGA public data, identified necroptosis-related genes, and conducted GO and KEGG enrichment analyses. Using Cox regression analysis and LASSO analysis to identify independent prognostic factors, a predictive model was established and validated in clinical subgroups, and correlation analysis with immune cells and ssGSEA differential analysis were conducted. Finally, potential drugs for HCC were screened to explore the drug sensitivity of different subtypes.

Results

We identified 19 differentially expressed necroptosis-related genes and constructed a predictive model with 3 independent prognostic factors through stepwise Cox regression. Validation results from clinical subgroups showed that the constructed model performed well in risk prediction, and ssGSEA differential analysis results were significant. We analyzed 55 immunotherapy drugs, and clustered them by distinct IC50 values to guide drug selection for HCC patients. Notable, Bleomycin, Obatoclax. Mesylate, PF.562271, PF.02341066, QS11, X17. AAG, and Bl. D1870 exhibited significantly different sensitivities in different subtypes, providing references for clinical practice in HCC patients.

Keywords

Microenvironment
Necroptosis
Predictive model
Drug sensitivity
==== Body
pmc1 Introduction

Liver cancer is a highly prevalent malignant tumors with poor prognosis [1]. It includes HCC and intrahepatic cholangiocarcinoma (ICC). HCC predominates as the primary pathological variant of hepatic malignancies in China, comprising a significant majority of 80 % of the total incidence [2]. It has a low initial resectability rate, with high recurrence and mortality rates [3]. To improve the survival rate of HCC patients, effective biomarkers for early diagnosis are needed [4].

Cell death is generally considered an “active” rather than “passive” process. Its actively mediated cell suicide processes and passive death methods are collectively referred to as programmed cell death (PCD) [5]. PCD includes various forms such as apoptosis, autophagy, ferroptosis, and pyroptosis. Necroptosis is an alternative form of cell death to apoptosis and has been documented in various liver diseases. It exerts its effects by releasing DAMPs, cytokines, and chemokines in the tumor microenvironment, which activate dendritic cells, as well as CD4 and CD8 T cells. These cells then infiltrate tumor tissues to form an inflammatory immune microenvironment and generate anti-tumor effects [6,7]. Seehawer et al. published an intriguing study in the journal Nature, revealing that apoptosis or necroptosis of hepatocytes in the liver microenvironment can determine the subtype of liver cancer (HCC or ICC). The necroptosis-associated hepatocyte microenvironment determines the growth of ICC in oncogenic-transformed hepatocytes, whereas hepatocytes with the same oncogenic drivers will produce HCC if surrounded by apoptotic hepatocytes [8]. These studies also emphasize the liver microenvironment's pivotal role in disease progression through significant dynamic changes. For HCC or ICC, the authors elucidate that the liver microenvironment can regulate liver cancer type. The mechanism of cell death, whether apoptosis or necroptosis, regulates the liver microenvironment itself. Although these two cancer subtypes may originate from the same hepatocytes, they exhibit significant differences in histology, metastatic potential, and prognosis. So far, the mechanisms determining the decision of one cancer subtype over the other remain largely unknown.

Studies have shown that lncRNA directly participate in the progression of HCC and play a significant regulatory role [9]. It is closely associated with the biological behaviour of tumour cells [10,11]. Therefore, this research investigates the association between genes linked to necroptosis and HCC, and develops a survival prediction model for HCC utilizing necroptosis-related lncRNA. Additionally, drug sensitivity analysis for different cancer subtypes is conducted to identify suitable drugs and provide valuable references for effective HCC treatment.

2 Materials and methods

2.1 Data sources

This study obtained the transcriptome (RNA-seq) and related clinical data of liver cancer from the TCGA public database. Clinical information from 377 patients and RNA-seq expression data from 424 samples were downloaded, along with 67 necroptosis-related genes were sourced from the KEGG database. Drug response data were acquired from the CTRP and GDSC, covering 497 anticancer drugs. We have uploaded the key organized data and code for each analysis step to zenodo at https://zenodo.org/records/11576013.

2.2 Data processing methods

The data processing procedure in this study includes the following aspects: (1) Handling missing data: The survival status and age in the downloaded and organized clinical data were processed for missing values. Samples lacking survival time or with survival less than 30 days were excluded, resulting in 350 clinical samples with complete survival status, survival time, and age. (2) Data grouping: The integer mean age of 65 was used as the cutoff for age grouping, dividing the age into two groups. Survival status (stat), gender (Gender), distant metastasis (M), lymph node metastasis (N), primary tumour (T), and pathological stage (Stage) were organized and summarized as clinical classification variables. (3) Expression matrix processing: The transcriptome data were organized into an expression matrix. The expression matrix was transformed into an ID-based format utilizing the human genome annotation file. The data after ID conversion were separated into lncRNA and mRNA, resulting in 19,573 mRNAs and 14,056 lncRNAs. The expression matrix and lncRNA matrix after ID conversion were subjected to co-expression analysis. The filtering criteria comprised a correlation coefficient below 0.35 and a significance threshold of P < 0.001.

2.3 Co-expression analysis methods

The sample data which consists of 374 tumor and 50 normal tissue samples was divided into 19,573 mRNAs and 14,056 lncRNAs. The Pearson correlation coefficient (PCC) was calculated, and lncRNA-mRNA pairs with an absolute value of the correlation coefficient (|correlation|, |cor|) > 0.4 and P < 0.05 were selected to construct the co-expression network.

2.4 Methods for differential gene analysis

We conducted differential gene analysis using expression data from 14,056 lncRNAs (NRGLncExp.txt). We applied filtering criteria with logFC ≥1 and FDR ≤0.05, and utilized the “limma” and “heatmap” R packages to identify differentially expressed lncRNAs. Heatmaps and volcano plots were generated to illustrate expression patterns between tumor and non-tumor tissues, and genetic alterations were examined. Key organized data and code for each analysis step have been uploaded to Zenodo.

2.5 Enrichment analysis

Limma package was utilized for screening differentially expressed necroptosis-related genes. Subsequently, gene function enrichment analysis was conducted to pinpoint significant biological attributes, encompassing both GO and KEGG enrichment analyses. The enriched terms were visualized using the GOplot package for enhanced interpretability.

2.6 Prognostic gene characterization and model construction

We integrated gene expression and survival data (stored in expTime.txt) and applied univariate Cox regression analysis to identify prognostic genes with criteria of P < 0.001.

Patients were stratified into training and test sets in a 1:1 ratio. Lasso regression analysis was used to refine the selection of genes linked to overall survival (OS) prognosis. Subsequently, stepwise multivariate Cox regression analysis pinpointed three genes prognostically significant.

The expression coefficients for each independent risk gene were obtained through the multivariate Cox proportional hazards regression model. All the coefficients were then used to formulate a prognostic model. The risk score proposed in this study was shown as equation (1).(1) RS = c1 *e1 + c2 *e2 + c3 *e3

Where RS denotes risk score, ci reprents the coefficient and ei is the expression of the ith genes.

Patients were grouped into high- and low-risk categories by optimal cutoff values evaluated using “survival” and “survminer” packages. The “survROC” package was used to plot ROC curves and compute AUC values at 1, 3, and 5 years for HCC, and at 1 year for all clinical features. Survival comparisons between risk groups were assessed using Kaplan-Meier(K-M) curves and log-rank tests via the “survival” package. The gene signature's prognostic utility was validated in both the test and full cohorts.

2.7 Model validation for clinical grouping

Using the risk score data for all HCC samples (risk.all.txt) and the clinical data of HCC cases (clinical.xls), we first sorted and merged the stages into groups I-II and III-IV. In order to validate the clinical grouping model, a p < 0.05 level of significance was used in this study. Prognostic factors determined via multivariate Cox regression were integrated into a nomogram for predicting 1-year, 3-year, and 5-year OS probabilities in HCC patients. The nomogram was validated by assessing its discrimination and calibration. Calibration curves were plotted to compare predicted probabilities with observed rates.

2.8 Correlation analysis of immune cells

We combined the immune cell infiltration data across pan-cancer (infiltration_estimation_for_tcga.xls) with the risk score data for all HCC samples (risk.all.txt) to identify intersecting samples. Correlation analysis of immune cell was executed utilizing R packages including “scales”, “ggtext”, “tidyverse” and “ggpubr”. Finally, the sources of software types were visualized using bubble charts.

2.9 ssGSEA differential analysis

Utilizing gene expression data (symbol.txt) and immune gene set data (immune.gmt), we conducted ssGSEA analysis with the “GSVA” and “GSEABase” R packages. In addition, ssGSEA differential analysis was performed based on score data (ssgseaOut.txt) and risk score data for all HCC samples(risk.all.txt). Box plots were then generated to depict the disparity in immune cells and function profiles between high- and low-risk groups.

2.10 Tumor classification

We performed unsupervised consensus clustering to HCC samples(risk.all.txt) with “limma” and “ConsensusClusterPlus” packages. Specifically, Multiple consensus partitions were created for K values ranging from 2 to 10. For each K clustering partition, hierarchical clustering with 1000 resamplings was performed using the Ward linkage method and Euclidean distance as the distance metric. Based on the cumulative distribution function (CDF) of the consensus matrix, the optimal number of clusters was determined to be 3, resulting in the cluster classification file. Using the classification results file, we wrote a program using the “survival” and “survminer” R packages to present classified survival analysis with survival curves. The risk and patient correspondence were visualized with a Sankey diagram. Additionally, principal component analysis (PCA) was employed to differentiate between high- and low-risk groups and among classifications 1, 2, and 3.

2.11 Drug sensitivity analysis

Utilizing gene expression data (symbol.txt) and classification result data (cluster.txt), we used the “limma”, “ggpubr”, “pRRophetic”, and “ggplot2” R packages with a pFilter = 0.001 to conduct drug sensitivity analysis. Specifically, drug sensitivity was quantified through the calculation of dose-response AUC values. Furthermore, the chi-square test was employed to identify disparities in drug responses among HCC subgroups. Drugs exhibiting a P-value <0.05 were deemed to demonstrate statistically significant differences in sensitivity across subgroups.

2.12 Statistical analysis

Survival curves were plotted using the K-M method and log-rank tests were conducted to assess survival disparities among risk subgroups. A Cox proportional hazards model was performed for multivariate analysis. All statistical procedure were conducted in R, with two-sided tests and a significance threshold of P < 0.05.

3 Results

3.1 Differential expression and functional enrichment of necroptosis-related lncRNA in HCC

By constructing a gene co-expression network, we visualized the correlation between necroptosis-related lncRNAs and necroptosis genes (Fig. 1a). This demonstrated the scientific validity of using necroptosis-related lncRNAs for differential expression analysis in HCC. Volcano plots (Fig. 1b) and heatmaps (Fig. 1c) revealed the expression patterns of differentially expressed genes in tumour versus non-tumour tissues, identifying differences between the two groups. Ultimately, with an FDR threshold of <0.05 and a [log2(FC)] cutoff >1, we identified 19 differentially expressed genes(Table S1), including the downregulated genes ID1 and BACH2 and 17 upregulated genes. Additionally, analysis of genetic alterations showed that truncating and missense mutations as prevalent(Fig. 1d), nine genes had mutation rates ≥3 %, with CDKN2A and TRIM11 being the most frequently mutated genes (8 %).Fig. 1 Differential Expression Analysis. (a) Co-expression network with necroptosis-related lncRNAs in green and necroptosis-related genes in red. (b) Volcano plot of necroptosis-related lncRNA expression differences in HCC, with upregulated shown in red, downregulated in green, and no difference in black. (c) Heatmap of differentially expressed necroptosis-related lncRNAs in HCC, green represents normal samples, and red represents tumour samples. (d) Mutation rate chart of differentially expressed necroptosis-related lncRNAs in HCC.

Fig. 1

Enrichment analyses, including KEGG and Gene Ontology (GO), were carried out for the differentially expressed genes. The GO enrichment analysis (Fig. S1a) indicates that the pathway most significantly enriched within the biological process category is associated with tubulin binding. In terms of cellular components, the chromosomal region shows the highest enrichment, while organelle fission dominates in molecular functions. Figs. S1b–c further indicate significant associations with pathways related to cell necrosis in the KEGG enrichment analysis. The positive Z-scores of the most enriched pathways suggest potential pathway enhancements.

3.2 Construction of a prognostic risk model based on necroptosis-related genes

Univariate Cox regression analysis of the 19 differentially expressed genes revealed that 10 genes exhibited significantly associations with OS(Table S2). Subsequently, using LASSO analysis, we screened out 6 genes (Figs. S2a–b, Table S3). Following multivariate stepwise Cox regression analysis, we finally identified 3 genes strongly associated with prognosis. These genes were utilized to develop the prognostic model: RS = 0.5228 * eSQSTM1 + 0.4019 * ePLK1 + 0.3348 * eMYCN (Table S4) where ei reprents the expression of gene i.

We used a Sankey diagram to illustrate the positive and negative regulatory relationships between lncRNAs and genes (Fig. 2a), thereby identifying necroptosis-related lncRNAs from the co-expression results file. Based on the median cut-off point of the risk score (0.8746), we stratified the 343 HCC patients into a high-risk group comprising 165 cases and a low-risk group comprising 178 cases. The distribution of patients risk scores and survival status is shown in Fig. 2b, and the expression levels of these three genes were displayed in a heatmap (Fig. 2c). The OS of patients within the high-risk cohort was markedly reduced compared to those in the low-risk cohort, as depicted in Fig. 2d. To evaluate the performance of the model, we calculated the area under the curve (AUC) (Fig. 2e); the AUCs for 1 year, 3 years, and 5 years were 0.794, 0.667, and 0.638, respectively. Meanwhile, the 95 % confidence intervals (CI) for these AUCs are [0.7234, 0.8650], [0.5887, 0.7455] and [0.5396, 0.7359].Fig. 2 Performance of the prognostic risk model. (a) Sankey diagram of necroptosis-related genes and lncRNAs, illustrating their positive and negative regulatory relationships. (b) Risk score distribution and survival outcomes. (c) Heatmap showing expression patterns of three genes in high-risk versus low-risk groups. (d) Kaplan-Meier survival curves for high-risk and low-risk HCC patients in TCGA. (e) ROC analysis for predicting mortality risk at 1, 3 and 5 years using risk scores. (f) ROC analysis of 1-year mortality risk in TCGA patients using risk scores and clinicopathological features.

Fig. 2

To compare the prognostic risk score with the predictive ability of other clinicopathological characteristics, we also calculated the 1-year AUC values under the ROC curve for all clinical features (Fig. 2f). The risk feature demonstrated an AUC value of 0.794, with a 95 % confidence interval ranging from 0.7234 to 0.8650, which was considerably greater than the characteristics associated with age, gender, grade, and stage. The results suggest that the risk feature exhibits superior predictive capability for the survival of HCC patients in comparison to clinical parameters. Overall, the findings of this study suggest that the three-gene signature performs well in survival prediction.

3.3 Clinical differentiation and validation of high and low-risk scores

To verify the effect of risk scores on clinical data, we conducted a clinicopathological assessment using a correlation diagram of clinical parameters and risk features (Fig. 3a). The analysis demonstrated a significant association between the risk score and both Tumor (P = 0.003) and Stage (P = 0.002). To further validate the independent prognostic performance of necroptosis-related genes for overall survival (OS), we conducted a univariate Cox analysis. The results indicated both survival score and stage were prognostic indicators for the survival of HCC patients (Fig. 3b, left panel). Subsequently, a multivariate Cox analysis of these factors validated that necroptosis-related genes have prognostic predictive value for HCC patients (Fig. 3b, right panel). Thus, the selected necroptosis-related differential genes were validated as independent prognostic factors in clinical practice.Fig. 3 The clinical differences of risk scores and validation results in clinical subgroups. (a) Correlation diagram between clinical parameters T and Stage and risk features, (b) Forest plot of clinical-related univariate (green) and multivariate (red) Cox analyses, (c) KM survival curves for Stage I-II and Stage III-IV, (d) Nomograms for 1-, 3-, and 5-year OS based on risk score, age, grade, and Stage, (e) Calibration plots for predicting 1-, 3-, and 5-year OS to assess the prediction model.

Fig. 3

To ascertain the applicability of the risk scores across all stages, patients were stratified into Stage I-II and Stage III-IV cohorts for discrete validation. The results revealed significant findings with P < 0.001 for both groups, indicating that the model is effective across all stages(Fig. 3c).

Furthermore, we developed a predictive nomogram to estimate the survival probabilities of HCC patients. The nomogram, which includes risk levels and clinical risk features, predicts 1-, 2-, and 3-year OS incidence rates. Risk levels of the prognostic model in the nomogram (Fig. 3d) demonstrated outstanding predictive performance over clinical factors, and the calibration plots(Fig. 3e) showed ideal consistency between observed and predicted 1-, 2-, and 3-year OS ratios.

3.4 Correlation with immune cells and ssGSEA differential analysis

Through immune cell correlation analysis, we utilized a variety of analytical tools, including XCELL, TIMER, QUANTISEQ, MCPCOUNTER, EPIC, CIBERSORT-ABS, and CIBERSORT, to obtain predictive outcomes. Using a bubble plot (Fig. S3), we clearly demonstrated the positive and negative regulatory interactions between immune cells and risk assessment values. To further explore the associations between risk scores and immune cells as well as their functions, we utilized ssGSEA to assess the abundance of 16 distinct immune cell categories and their functional profiles. The results showed notable variations in immune cell types, such as aDCs, iDCs, Macrophages, Mast cells, NK cells, T helper cells, and Tregs, when comparing the high-risk and low-risk groups (Fig. 4a).Fig. 4 The analysis results of differences in immune cell populations and immune function between the high-risk and low-risk groups. (a) Relationship between risk scores and immune cells, (b) Association between risk scores and immune functions.

Fig. 4

Additionally, the low-risk group exhibited markedly elevated scores for immune functions including but not limited to Cytolytic activity, Type II IFN Response, and Type I IFN Response, indicating that immune functions related to necroptosis are more pronounced in individuals with low-risk group (Fig. 4b).

3.5 Selection of promising drug candidates for the disease

This study evaluated the infiltration levels of 22 immune cell types in the tumor microenvironment between high-risk and low-risk groups and investigated the correlation between prognostic risk scores and immune cell infiltration through ssGSEA differential analysis of immune checkpoints. The analysis revealed significant differences in genes related to immune checkpoints between the two risk groups, with higher expression levels of TNFSF18, HHLA2, TNFRSF25, and TNFRSF18 observed in the high-risk group. Conversely, IDO2 gene expression was relatively more activated in the low tumour mutational burden (TMB) group (Fig. 5a). Additionally, using the somatic mutation data from the TCGA-HCC cohort, this study calculated TMB to investigate its association with survival outcomes. Patients were categorized into high and low TMB groups based on a specific threshold. The findings revealed that the high TMB group had a lower survival probability than the low TMB group (Fig. 5b–c). Among the groups, the high TMB high-risk group showed the lowest survival probability, while the low TMB low-risk group had the highest. This suggests a strong association between the tumor immune microenvironment (TIME) and TMB.Fig. 5 Diagrams of immune cell analysis results. (a) Infiltration rates of immune cells in risk-stratified groups, (b) Survival analysis based on varying tumour mutational burden characteristics by Kaplan-Meier method, (c) Survival analysis based on varying tumour mutational burden characteristics and different risk scores by Kaplan-Meier method.

Fig. 5

Furthermore, to further compare the differences in responsiveness to immunotherapy between groups with different risk levels, and to evaluate which patients who are more inclined to gain from chemotherapy, this study calculated the half-maximal inhibitory concentration (IC50) of immunosuppressants in patients categorized as high- and low-risk. The results indicated that the p-values of 55 drugs in the database were less than 0.001 (Table S5), suggesting significant sensitivity differences between the two groups. Specifically, the distinct clustering of the 55 immunotherapy drugs, characterized by their varied IC50 values, provides critical insights for selecting appropriate treatments for HCC patients. This outcome also contributes valuable insights for further research on immunotherapy response and enhances precision medicine for HCC patients.

3.6 Drug sensitivity analysis of classifications

Due to the suboptimal results from sensitivity analysis solely within the two groups of different risk levels, and considering samples from different cancer subtypes often exhibit distinct molecular characteristics, we classified them into subtypes. Through the application of K-means and consensus clustering algorithms, we identified that three is the ideal and most stable number of prognosis-related molecular subtypes. Ultimately, all samples were grouped into three categories(Fig. 6a), comprising 180 cases in C1, 104 cases in C2, and 59 HCC patients in C3. Kaplan-Meier curve analysis indicated a significant difference in survival rates among HCC patients corresponding to the three subtypes (P < 0.001) (Fig. 6b). In this context, the overall survival rate was highest in C1, followed by C2, with C3 having the lowest survival rate.Fig. 6 Tumour Typing Related Diagrams. (a) Consensus matrix for k = 3 clusters. (b) Kaplan-Meier curves among three types. (c) Population distribution characteristics of clustering types in two groups of different risk levels. (d) Heatmap of immune cells in different subtypes.

Fig. 6

Meanwhile, the Sankey diagram(Fig. 6c) further illustrated that C1 had a higher proportion of low-risk patients, while C2 and C3 had more high-risk patients, explaining the observed differences in overall survival rates. The study utilized principal component analysis (Figs. S4a–b) and t-SNE plots (Figs. S4c–d) to identify distinctions among the three different subtypes of necroptosis-related genes, and to compare the two groups of different risk levels. The two analytical approaches together provide conclusive evidence that genes associated with necroptosis can effectively differentiate samples from patients with HCC. The heatmap of classified immune cells shows the upregulation and downregulation of cells in the three different subtypes (Fig. 6d).

To obtain the conclusive results of the drug sensitivity analysis, we initially conducted a differential analysis of immune checkpoints across the subtypes. Box plots (Fig. 7a) were employed to illustrate the levels of immune cell infiltration across each subtype, with a specific focus on the expression of necroptosis-related genes within the three subtypes.Fig. 7 Results of Drug Sensitivity Analysis. (a) Infiltration levels of various immune cell subtypes within the subgroups. (b–h) Comparison of IC50 values of drugs in the subtypes.

Fig. 7

Subsequently, we performed a drug sensitivity analysis on the three cancer subtypes. Drugs with p-values less than 0.05 were deemed to exhibit substantially different sensitivities in the respective subgroups. The results indicate that Bleomycin, Obatoclax Mesylate, PF-562271, PF-02341066, QS11, X17-AAG, and Bl-D1870 have significantly different sensitivities in the various subgroups (Fig. 7b–h).

4 Discussion

HCC ranks among the most prevalent malignant tumors globally, often diagnosed at an advanced stage due to an insidious onset. With the progress in novel clinical diagnostic and therapeutic methods, alongside advancements in early diagnostic techniques for HCC, the epidemiology of HCC has undergone significant changes in recent years. This work intends to create a predictive model for HCC by utilizing necroptosis-related genes, in order to advance clinical therapy research in this field. By doing differential gene expression research, we have discovered 19 genes that are expressed differently in relation to necroptosis. Through the application of statistical methodologies, we have determined that SQSTM1, PLK1, and MYCN are independent prognostic factors. Subsequently, these genes were utilized to build the prognostic model. In order to further confirm the clinical importance of the risk scores obtained from the model, we created experimental and validation groups. Each group was then broken into subgroups based on their level of risk, either high-risk or low-risk. Through a comparison of the overall survival rates of patients classified as high-risk and low-risk, we have established that the risk score is a significant determinant of unfavorable prognosis in patients with HCC. In the clinical correlation analysis, the risk score showed a substantial association with both the histology grade and the overall categorization. The ROC analysis demonstrated that the developed risk signature is a highly sensitive predictor for forecasting the survival rates of HCC patients at 1-year, 2-year, and 3-year intervals. In addition, we conducted a comparison of the predicted efficacy of the prognostic risk score with other clinical prognostic indicators and determined that the risk score outperformed all other potential factors. The risk score was determined to be a separate risk factor for accurately predicting the prognosis of patients with HCC. The ssGSEA differential analysis results were found to be statistically significant based on the validation results of the clinical grouping. This study utilized clustering analysis to group 55 immunotherapy medications based on their distinct IC50 values. The findings of this study offer a theoretical foundation for the selection of treatments in patients with HCC. Various drugs, including Bleomycin, Obatoclax Mesylate, PF-562271, PF-02341066, QS11, X17-AAG, and BI-D1870, showed distinct levels of sensitivity in different subgroups. These findings provide valuable clinical information for the treatment of HCC patients.

The justification for this investigation is based on numerous crucial factors. Previous studies have shown that programmed cell death is a sophisticated process regulated by cytokines, pattern recognition receptors, and various genes and signalling pathways. Apoptosis is a typical form of PCD, while necroptosis, a more recent finding, has several characteristics with apoptosis [12]. Current research indicates that the process of necroptosis has emerged as a pivotal contributor to the development of diverse inflammatory diseases, including atherosclerosis, systemic inflammatory response syndrome, and acute kidney injury [13]. Recent findings suggest that tumors associated with apoptosis may correlate with necroptosis. The transcription of oncogenes and expression of protein-coding genes often drives the occurrence of cancer. Previous research has shown that human genome is remarkably skewed, with non-coding regulatory genes making up over 90 % of the total, while protein-coding genes constitute only around 2 %. lncRNAs have a wide-ranging regulatory impact on cancer formation through their influence on the expression of genes relevant to tumors and their involvement in crucial biological processes [14,15]. Due to the developments in high-throughput sequencing technology, there is a growing focus on investigating, identifying, and confirming the involvement of new lncRNAs in the evolution of cancer [16,17]. Therefore, the research and exploration of biomarkers for HCC are highly valuable for both basic research and clinical treatment. This study aims to construct a risk-scoring model for HCC based on necroptosis-related genes and to validate its predictive performance in overall survival prediction and immunotherapy efficacy evaluation in HCC. Additionally, the study conducts a drug sensitivity analysis for different tumour subtypes to identify suitable drugs.

In addition, previous studies have shown that SQSTM1 acts as a mediator of MDRL/miR-361 ceRNA activity. Initial studies have indicated that miR-361 serves as a regulator of aerobic glycolysis and proliferation in gliomas [18]. Moreover, the systemic delivery of miR-361 hinders the development of gliomas by inhibiting UBR5 [19]. PLK1 is a central kinase that coordinates diverse mitotic phases and is indispensable to its cancer-driving role in malignant cells [20]. Suppression of PLK1 function leads to mitotic catastrophe and death in several forms of cancer, whereas normal cells are less impacted since they rely less on PLK1 for cell division [21]. MYCN, a component of the MYC oncogene, produces a transcription a gene transcriptional regulator that regulates the expression level of various objective genes involved in various cellular processes, these processes encompass proliferation, apoptosis, senescence, differentiation, metabolism, DNA damage repair, and protein synthesis [22,23]. The MYC family governs various biological processes, such as cell growth, proliferation, apoptosis, energy metabolism, and differentiation [24]. It has two functions in hepatocyte proliferation and hepatocarcinogenesis [25]. Currently, the majority of research on the predictive importance of lncRNAs is focused on developing models, with insufficient investigation into the processes of these genes in living organisms and laboratory settings. As necroptosis has an effect on the development of liver cancer, genes associated with it could potentially be used as prognostic biomarkers for patients with HCC. TMB has recently been identified as a new biomarker for cancer immunotherapy in the field of immunotherapy research [26]. TMB refers to the cumulative count of somatic gene coding mistakes, such as base substitutions, gene insertions, or deletions, identified per million bases [27]. Increased TMB is linked to the generation of novel tumor-specific antigens, which can be recognized by the immune system and leveraged to identify patients most likely to respond to immune checkpoint inhibitors therapies [[28], [29], [30], [31]]. These findings also provide strong support for the clinical grouping model validation conducted in this paper. Therefore, this project seeks to develop a risk-scoring model for HCC using genes related to necroptosis. The objective is to assess the model's ability to predict overall survival and the effectiveness of immunotherapy. The study also incorporates a medication sensitivity analysis for various tumor subtypes in order to find suitable treatments, hence facilitating future comprehensive research.

Overall, this study built a prognostic prediction model based on the data resource from the TCGA public database. The model is based on the corrlation between necroptosis-related lncRNAs and necroptosis-related genes. This model can be used as a reference for predicting patient prognosis and providing clinical guidance. This study enhances current research on genetics associated to necroptosis in HCC by making predictions about prognosis and the effectiveness of immunotherapy. Nevertheless, it is important to recognize and address numerous constraints. First, the results were confirmed in the TCGA cohort, while their reliability needs additional validation in other independent cohorts. Further in vivo and in vitro investigations are needed to confirm if necroptosis-related lncRNAs belonging to necroptosis-related genes can boost the effectiveness of immunotherapy in patients with HCC. Second, the expression status of lncRNAs associated with necroptosis in biological samples needs to be confirmed. Further research is required to elucidate the mechanistic pathways by which the prognostic model impacts the progression of liver cancer. This could provide new targets and treatments for clinical therapy, improve systemic treatment effects, and ultimately work to prolong the overall survival duration for individuals with liver cancer. Third, it is inadequate to use a single feature for response prediction; necroptosis-related genes are only one factor that influences the prognosis of HCC. Fourth, PCA algorithm possess intrinsic drawbacks such as reliance on linear assumptions, susceptibility to outliers, loss of interpretability of features, and the need for meticulous component selection. PCA seeks to decrease dimensionality by preserving the most crucial information, but it is inevitable that certain features may be lost in the process.

Data availability statement

The transcriptome and related clinical data are openly available in TCGA public database. Clinical information was sourced from the KEGG database (https://www.kegg.jp/). Drug response data were referenced from the CTRP (https://portals.broadinstitute.org/ctrp.v2.1/) and GDSC (https://www.cancerrxgene.org/). For the convenience of readers, we have uploaded the key organized data, code, and analysis process to zenodo at https://zenodo.org/records/11576013.

CRediT authorship contribution statement

Ronghuo Wu: Writing – original draft, Data curation, Conceptualization. Xiaoxia Deng: Writing – review & editing. Xiaomin Wang: Writing – review & editing. Shanshan Li: Writing – review & editing, Conceptualization. Jing Su: Writing – review & editing, Data curation. Xiaoyan Sun: Writing – review & editing, Conceptualization.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A Supplementary data

The following is the Supplementary data to this article:Multimedia component 1

Multimedia component 1

Acknowledgements

The project was partially supported by the 10.13039/501100001809 National Natural Science Foundation of China (3260777, 12261096 ), the 10.13039/100012547 Natural Science Foundation of Guangxi (2023GXNSFAA0262292345678 , 2020GXNSFAA159155 , 2020GXNSFDA297016 ).

Appendix A Supplementary data to this article can be found online at https://doi.org/10.1016/j.heliyon.2024.e36561.
==== Refs
References

1 Díaz-González Á. Forner A. Rodríguez de Lope C. Varela M. New challenges in clinical research on hepatocellular carcinoma Rev. Esp. Enferm. Dig. 108 8 2016 Aug 485 493 26653993
2 Cheng M. Zhang J. Cao P.B. Zhou G.Q. Prognostic and predictive value of the hypoxia-associated long non-coding RNA signature in hepatocellular carcinoma Yi Chuan 44 2 2022 Feb 20 153 167 35210216
3 Chen M. Zhou X. Cai H. Li D. Song C. You H. Ma R. Dong Z. Peng Z. Feng S.T. Evaluation of hypoxia in hepatocellular carcinoma using quantitative MRI: significances, challenges, and advances J. Magn. Reson. Imag. 58 1 2023 Jul 12 25
4 Zeng F. Xu Z. Zhuang P. Integrated analysis of SKA1-related ceRNA network and SKA1 immunoassays in HCC: a study based on bioinformatic Medicine (Baltim.) 102 38 2023 Sep 22 e34826
5 Gao W. Wang X. Zhou Y. Wang X. Yu Y. Autophagy, ferroptosis, pyroptosis, and necroptosis in tumor immunotherapy Signal Transduct. Targeted Ther. 7 1 2022 Jun 20 196
6 Wang R. Li H. Wu J. Gut stem cell necroptosis by genome instability triggers bowel inflammation Nature 580 4 2020 386 390 32296174
7 Sprooten J. Wijngaert P.D. Vanmeerbeek I. Necroptosis in immunooncology and cancer immunotherapy Cells 9 8 2020 1823 32752206
8 Saeed W.K. Jun D.W. Viewpoint: necroptosis influences the type of liver cancer via changes of hepatic microenvironment Hepatobiliary Surg. Nutr. 8 5 2019 Oct 549 551 31673555
9 McFadden E.J. Hargrove A.E. Biochemical methods to investigate lncRNA and the influence of lncRNA:protein complexes on chromatin Biochemistry 55 11 2016 Mar 22 1615 1630 26859437
10 Huang J.L. Zheng L. Hu Y.W. Wang Q. Characteristics of long non-coding RNA and its relation to hepatocellular carcinoma Carcinogenesis 35 3 2014 Mar 507 514 24296588
11 Standaert L. Adriaens C. Radaelli E. Van Keymeulen A. Blanpain C. Hirose T. Nakagawa S. Marine J.C. The long noncoding RNA Neat1 is required for mammary gland development and lactation RNA 20 12 2014 Dec 1844 1849 25316907
12 Zhu P. Ke Z.R. Chen J.X. Li S.J. Ma T.L. Fan X.L. Advances in mechanism and regulation of PANoptosis: prospects in disease treatment Front. Immunol. 14 2023 Feb 9 1120034
13 Wei X. Xie F. Zhou X. Wu Y. Yan H. Liu T. Huang J. Wang F. Zhou F. Zhang L. Role of pyroptosis in inflammation and cancer Cell. Mol. Immunol. 19 9 2022 Sep 971 992 35970871
14 Statello L. Guo C.J. Chen L.L. Huarte M. Gene regulation by long non-coding RNAs and its biological functions Nat. Rev. Mol. Cell Biol. 22 2 2021 Feb 96 118 33353982
15 Zhang R. Xia L.Q. Lu W.W. Zhang J. Zhu J.S. LncRNAs and cancer Oncol. Lett. 12 2 2016 Aug 1233 1239 27446422
16 Chen B.W. Zhou Y. Wei T. Wen L. Zhang Y.B. Shen S.C. Zhang J. Ma T. Chen W. Ni L. Wang Y. Bai X.L. Liang T.B. lncRNA-POIR promotes epithelial-mesenchymal transition and suppresses sorafenib sensitivity simultaneously in hepatocellular carcinoma by sponging miR-182-5p J. Cell. Biochem. 122 1 2021 Jan 130 142 32951268
17 Sun Z. Xue S. Zhang M. Xu H. Hu X. Chen S. Liu Y. Guo M. Cui H. Aberrant NSUN2-mediated m5C modification of H19 lncRNA is associated with poor differentiation of hepatocellular carcinoma Oncogene 39 45 2020 Nov 6906 6919 32978516
18 Long N. Chu L. Jia J. Peng S. Gao Y. Yang H. Yang Y. Zhao Y. Liu J. CircPOSTN/miR-361-5p/TPX2 axis regulates cell growth, apoptosis and aerobic glycolysis in glioma cells Cancer Cell Int. 20 2020 Aug 6 374 32774168
19 Jia J. Ouyang Z. Wang M. Ma W. Liu M. Zhang M. Yu M. MicroRNA-361-5p slows down gliomas development through regulating UBR5 to elevate ATMIN protein expression Cell Death Dis. 12 8 2021 Jul 28 746 34321465
20 Lens S.M. Voest E.E. Medema R.H. Shared and separate functions of polo-like kinases and aurora kinases in cancer Nat. Rev. Cancer 10 12 2010 825 841 21102634
21 Christoph D.C. Schuler M. Polo-like kinase 1 inhibitors in mono- and combination therapies: a new strategy for treating malignancies Expert Rev. Anticancer Ther. 11 7 2011 1115 1130 21806334
22 Westermark U.K. Wilhelm M. Frenzel A. Henriksson M.A. The MYCN oncogene and differentiation in neuroblastoma Semin. Cancer Biol. 21 2011 256 266 21849159
23 Otte J. Dyberg C. Pepich A. Johnsen J.I. MYCN Function in neuroblastoma development Front. Oncol. 10 2021 624079
24 Dang C.V. MYC on the path to cancer Cell 149 2012 22 35 22464321
25 Qu A. Jiang C. Cai Y. Kim J.H. Tanaka N. Ward J.M. Shah Y.M. Gonzalez F.J. Role of Myc in hepatocellular proliferation and hepatocarcinogenesis J. Hepatol. 60 2 2014 Feb 331 338 24096051
26 Van Dijk N. Funt S.A. Blank C.U. Powles T. Rosenberg J.E. van der Heijden M.S. The cancer immunogram as a framework for personalized immunotherapy in urothelial cancer Eur. Urol. 75 3 2019 Mar 435 444 30274701
27 Yarchoan M. Hopkins A. Jaffee E.M. Tumor mutational burden and response rate to PD-1 inhibition N. Engl. J. Med. 377 25 2017 Dec 2500 2501 29262275
28 Gubin M.M. Artyomov M.N. Mardis E.R. Schreiber R.D. Tumor neoantigens: building a framework for personalized cancer immunotherapy J. Clin. Invest. 125 9 2015 Sep 3413 3421 26258412
29 Schumacher T.N. Kesmir C. van Buuren M.M. Biomarkers in cancer immunotherapy Cancer Cell 27 1 2015 Jan 12 12 14 25584891
30 Grizzi G. Caccese M. Gkountakos A. Carbognin L. Putative predictors of efficacy for immune checkpoint inhibitors in non-small-cell lung cancer: facing the complexity of the immune system Expert Rev. Mol. Diagn. 17 12 2017 Dec 1055 1069 29032709
31 Allgauer M. Budczies J. Christopoulos P. Endris V. Implementing tumor mutational burden (TMB) analysis in routine diagnostics—a primer for molecular pathologists and clinicians Transl. Lung Cancer Res. 7 6 2018 Dec 703 715 30505715
