Information on gene expression and the clinical data of LIHC patients were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov). If there was no available information on survival, the patient was not selected. One hundred and fifty-nine necroptotic regulators (hsa04217) were acquired from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (https://www.kegg.jp/). Ethical approval and informed consent were not applicable as the data are publicly available.

Unsupervised consensus clustering was applied to the LIHC patient data based on the expression of 159 necroptotic regulators, to identify specific distinct necroptosis-related patterns. The “ConsensuClusterPlus” package was used to assess the number and stability of the clusters, in accordance with the elbow method and Gap statistic (Wilkerson and Hayes, 2010). The R package “limma” was then used for the identification of the 1,516 DENRGs through cluster comparisons (Ritchie et al., 2015), using the criteria of false discovery rate (FDR) < 0.05 and fold change (FC) absolute value > 1.5.

The 365 LIHC patients were allocated to either a training or validation set in a 1:1 ratio using R. The training set was used for subsequent analysis. DENRGs that were related to survival were examined by univariate Cox regression, and further selection was performed using the “least absolute shrinkage and selection operator” (LASSO) with the “glmnet” package to avoid overfitting (Engebretsen and Bohlin, 2019). Multivariate regression was then conducted to identify DENRGs suitable for use in the prognostic signature. The coefficients of the DERNGs used in the final signature were verified at the same time and were used to determine the risk scores of each of the patients in the training set, with the score determined as: Risk   Score   = ∑ i = 0 n β i   ∗   G i

β i is the coefficient of gene i in the multivariate analysis; G i represents the expression of gene   i ; and n represents the total number of genes.

Using the median risk score, the patients were assigned to a high-risk or low-risk group. Kaplan–Meier (K-M) survival curves with the log-rank test were established to determine the difference in overall survival (OS) between the two groups. Receiver operating characteristic (ROC) curves were created for 1-year, 3-year, and 5-year OS using the “survival ROC” package, and the areas under the curves (AUCs) were determined. Principal component analysis (PCA) and T-distributed stochastic neighbor embedding (t-SNE) were applied for dimensionality reduction, using the “Rtsne” package in R.

Data on the expression of genes used in the signature were then taken from the validation set and used for calculating risk scores, as described above. The individuals in the validation set were allocated to groups using the same method as the training set. The predictive capability of the signature was verified using ROC, K-M survival curves, PCA, and t-SNE.

Multivariate Cox regression was conducted on the significant (p < 0.05) factors derived from the univariate analysis of the signature and clinical parameters in the training set. A prognostic nomogram was constructed with the “rms” package in R using the identified independent factors for OS prediction. The agreement between the actual OS and those predicted by the nomogram was evaluated by calibration curves.

Gene Ontology (GO) annotation and KEGG pathway analyses were undertaken using the “clusterprofiler” package in R, with adjusted p-values <0.05 representing the significance threshold.

The “Single Sample Gene Set Enrichment Analysis” (ssGSEA) algorithm was used to measure the proportions of 16 immune cell types in the two groups. The infiltration of six cell types, namely, B cells, CD4⁺ T cells, CD8⁺ T cells, neutrophils, macrophages, and dendritic cells, were assessed in the TCGA samples and their association with the levels of genes used in the signature was determined using the TIMER algorithm (https://cistrome.shinyapps.io/timer/).

Gene Expression Profiling Interactive Analysis (GEPIA) is a commonly used interactive website that compares gene expression differences between patients and normal controls (http://gepia.cancer-pku.cn/). The Human Protein Atlas (HPA) was used for qualitative comparison of IHC staining of liver cancer tissue with normal liver (http://www.proteinatlas.org/).

R (version 4.1.0) was used for all data analysis. Spearman’s and distance correlations were used to determine relationships between variables. Differences in survival curves were assessed by log-rank tests. Clinical parameters were compared using t-tests and chi-square tests, while Cox regression was used to measure survival-associated parameters. Two-tailed p-values <0.05 were considered statistically significant. The optimal cut-off values were established by the X-tile program (X-tile software version 3.6.1) (Camp et al., 2004) were used for the determination of gene expression levels.

A flow chart illustrating the investigation steps is given in Supplementary Figure S1. In total, the data of 374 LIHC cancerous and 50 non-cancerous tissues were obtained from the TCGA database and included in the study. One hundred and fifty-nine necroptotic regulators (Supplementary File S1) were identified from the KEGG database, and their expression between LIHC and normal tissue was found to differ significantly (Figure 1A). As seen in the volcano plot (Figure 1B), 56 of the regulators were up-regulated (log FC > 1), while only two were down-regulated in LIHC (log FC < −1). The results indicate that most necroptotic regulators may act as oncogenes in LIHC.

To investigate the molecular heterogeneity of LIHC and identify patterns in LIHC patients based on necroptotic regulators, unsupervised consensus analysis was conducted on tumor samples, and different subgroups were established. The optimal cluster numbers (k) were calculated by setting the value of k from 2 to 9 (Supplementary Figure S2); it was found that k = 2 was the optimum value in terms of both within-cluster homogeneity and between-cluster heterogeneity (Figure 2A). Cluster 1 (C1) was associated with better OS, shown by K-M analysis (Figure 2B). A comparison between C1 and Cluster 2 (C2) yielded 1,516 DENRGs (fold change = 1.5) (Supplementary File S2). The heatmap of the distributions of the different clinicopathological features and immune scores between the clusters is shown in the panel below the heatmap (Figure 2C). Clinical variables, such as T stage, sex, and immune status differed between the clusters.

The univariate analysis showed that 142 DENRGs were associated with OS (Figure 3A). In the subsequent LASSO analysis aimed at reducing the risk of overfitting (Figures 3B, C), 19 DENRGs were found to be associated with OS. The multivariate Cox regression narrowed this count to eight DENRGs, which were then used to construct the OS prognostic signature. The coefficients of these eight DENRGs were verified at the same time and were used to evaluate the scores for individual patients; these scores were then used to divide the patients in the training set into two groups according to risk. The scores were illustrated by risk curves and scatterplots and these, together with the patient survival data, are shown in Figure 4 A and B. K-M analysis found significantly reduced OS in high-risk individuals (Figure 4C). To assess the performance of the signature, ROC curves for 1-year, 3-year, and 5-year OS were created, and the AUCs at these time points were 0.825, 0.769, and 0.751, respectively (Figure 4D). Moreover, PCA and t-SNE analyses showed separation between the patients in the two groups (Figures 4E, F).

Individuals in the validation set were likewise assigned to two groups according to the risk score. Of the 181 individuals in the set, 92 were classified as having high risk, and 89 were classified as having low risk. As was observed in the training set, the risk curves and scatterplots showed an association between high-risk scores and reduced OS (Figures 5A, B), while low risk was related to higher OS (Figure 5C). The AUCs for the 1-year, 3-year, and 5-year time points were 0.807, 0.692, and 0.701, respectively (Figure 5D). Furthermore, PCA and t-SNE also showed significant separation in the two-dimensional plane between the groups (Figures 5E, F). These findings imply that the signature was effective for predicting LIHC prognosis.

The relationships between the DENRG-based signature and clinicopathological variables were examined by correlation analysis. As seen in the heatmap (Figure 6A), tumor grade and T stage were found to be significantly associated with the risk score. Compared with the low-risk group, all eight DENRGs were highly expressed in the high-risk group and were risk factors for liver cancer. Typical immunohistochemical images of five genes were downloaded from Human Protein Atlas; however, no information was available for PFKFB4, MSC-AS1, and CHST1, which were unavailable from the database (Supplementary Figure S3). Moreover, we analyzed the expression of the eight DENRGs in several other common tumors and found that except for SPESP1, these genes were highly expressed in tumors (Supplementary Figure S4). Highly expressed genes are often risk factors for tumors, which is consistent with our results. Using both univariate and multivariate Cox regression, we observed a significant link between tumor stage and the signature in the training set (Figures 6B, C) and validation set (Figures 6D, E), indicating that these can be used as independent prognostic factors for survival, and further verifying the reliability and effectiveness of the eight-gene signature we established.

A nomogram was constructed for predicting OS at different time points in the training set (Figure 7A). It was found that the predicted OS matched the actual OS, shown by the calibration chart (Figure 7B). Furthermore, a DCA curve indicated that the nomogram was highly reliable in terms of individual parameters (Figure 7C). More meaningfully, the nomogram permitted the estimation of the OS of specific patients by entering the score of individual parameters, which greatly improves its practical application.

GO and KEGG analyses were used to examine the underlying mechanisms for the differences seen in prognosis between the groups. The GO results indicated that these genes were associated with “extracellular matrix organization,” “phagocytosis,” and “leukocyte chemotaxis”, while KEGG demonstrated enrichment in several pathways associated with carcinogenesis, including “PI3K-AKT signaling pathway,” “HIF-1 signaling pathway,” and “TNF signaling pathway”, as well as immune-related pathways, such as “IL-17 signaling pathway” and “intestinal immune network for IgA production”. These results suggested that DENRGs are involved not only in the tumor but also in the tumor microenvironment.

Variations in immune infiltration can result in varied clinical outcomes despite the presence of the same histological cancer type (Li et al., 2016). As the KEGG analysis showed enrichment in immune pathways, we investigated possible variations in immune infiltration in the training-set groups. This found higher infiltration levels in the high-risk patients, specifically, infiltration of dendritic cells, macrophages, T follicular helper cells, and regulatory T cells (p < 0.05) (Figure 8A). TIMER was then used to examine the associations between seven DENRGs (TIMER does not contain information on non-coding RNAs; thus, MSC-AS1 was excluded), tumor purity, and immune infiltration. This indicated associations between all seven DENRGs and the levels of six cell types (Figures 8B–H). In addition, we analyzed the expression of immune checkpoints (Supplementary Figure S5A) and m6a regulators (Supplementary Figure S5B) and found that the expression of 33 immune checkpoints, including CD274 (PD-L1), and most m6a regulators differed significantly between the high- and low-risk groups. These findings indicate that these DENRGs may be closely linked to the formation of the tumor microenvironment and thus may represent immunotherapeutic targets for LIHC.

Data on the expression of PFKFB4, MSC-AS1, SPESP1, BMP6, CHST1, CTSV (CTSL2), TMCO3, and PAGE2 were obtained for LIHC patients. The optimal cutoff value in terms of survival prediction was confirmed by using X-tile software. K-M survival analysis shows that the survival prognoses in the high-expression groups for the genes were relatively worse than those for low-expression groups (Figures 9A–H). This is consistent with the identification of these genes as high-risk factors (Figure 6A), verifying the accuracy of the signature.