The copy number variation (CNV), simple nucleotide variation (VarScan2 Variant Aggregation and Masking) were downloaded from TCGA database and UCSC Xena website, respectively. Next, 1398 HCC patients including RNA expression and corresponding clinical data were retrieved from TCGA (FPKM value, n=365) data portal, GEO database (GSE14520/NCI (17), n=221; GSE36376 (18), n=223; GSE76427 (19), n=115; GSE20140 (20), n=162; GSE27150, n=81), and International Cancer Genome Consortium (ICGC, n=231). In addition, a immunotherapy dataset (IMvigor210) with a total of 348 bladder cancer patients and corresponding clinical data were included ( Table 1 ) (21). The FPKM value was firstly transformed to TPM to value more similar to those resulting from microarrays (22). Then, all raw data in the GEO database were microarray data processed on Affymetrix and Illumina. The raw data retrieved from the Affymetrix platform were processed using the RMA algorithm of the “affy” package in R for background adjustment and normalization (23). Finally, using “ComBat” algorithm of the “sva” package to batch effects among different datasets (24).

By reviewing the previous literature, we identified 76 necroptosis-related genes ( Supplementary Table 1 ). The genes with the cut-off criteria of |log FC|≥ 1.0 and adj. P <0.05 were used as DEGs by using limma package (25). Using univariate Cox regression analysis to analyze the correlation between clinical information and necroptosis-related genes expression. Finally, 15 different and prognostic necroptosis related genes were identified.

Using “maftools” package, we then calculated mutation frequency of necroptosis-related genes. Moreover, “RCircos” package was used to visualized the chromosome location of necroptosis related genes CNV alteration. For drug sensitivity analysis, we collected the drug response data from the Genomics of Drug Sensitivity in Cancer (GDSC, https://www.cancerrxgene.org/) database (26) and Cancer Therapeutics Response Portal (CTRP, https://portals.broadinstitute.org/ctrp/?page=#ctd2BodyHome) database (27). The Spearman correlation analysis was performed using the “pRRophetic” package, which contained 138 drugs for analysis.

Based on 15 necroptosis-related genes identified above, we then conducted consensus cluster analysis with “ConsensusClusterPlus” package and cycle computation of 1,000 times as the threshold (28). The was followed by survival and gene-expression patterns analysis using “survival” package.

We identified two NEC clusters (NEC.cluster.A and NEC.cluster.B) base on 15 necroptosis-related genes. The DEGs was generated with “limma” package. The differential expression DEGs were assigned a |log2FC|>1 with adjusted P <0.05 as a significance threshold. Using univariate Cox regression analysis to analyze the correlation between clinical information and necroptosis related genes expression. Finally, 4000 different and prognostic genes were identified. Then, R package clusterProfiler was used to explore the functions between different and prognostic genes (adj. p< 0.05), and GO and KEGG enrichment analysis was performed.

Next, we used principal component analysis (PCA) method to uantify the necroptosis related subtypes of individual patients (29). The model was based on the meta cohort and named the necroptosis score (NEC score). A NEC score for each patient was calculated according to the following formula:

where i is the TPM value of each screened gene.

Gene set enrichment analysis (GSEA) is a computational method verifying whether a priori presetting of genes indicates statistically significant, concordant differences between two biological states (30). In order to explore the biological pathways involved in HCC progression, we divided into a NEC score.low group (n= 556) and a NEC score.high group (n= 599) based on the median of NEC score as a cut-off point and the number of permutations was 1000.

ssGSEA method, which uses gene expression profiles to infer the number of tumor infiltrating immune cells (31). A series of analysis on the NEC score in HCC and its correlation with the abundance of immune infiltrates was performed.

8 HCC patients tissues and corresponding adjacent tissues were collected from the second affiliated hospital of Guangdong Medical University to explore the expression of 15 marker genes in the tissue samples. IHC was performed on formalinfixed, paraffin-embedded tissue sections using a two-step protocol. Briefly, paraffin section of 15 patients were first deparaffinized and heated in pressure cooker for 10 min. After cooling to room temperature, the sections were immersed for 5 min in PBS three times and followed by endogenous peroxidase activity blocking with 3% H₂O₂ and non-specific staining blocking with 10% goat serum. Then these sections were incubated with primary antibodies overnight at 4°C and treated for 30 min at 37°C with second antibody. The sections were staining for 2 min by using DAB and rinsed off in deionized water to terminate DAB reaction. Then using traditional method to evaluate under the optical microscope. The protein expression was calculated by German immunohistochemical score (GIS). Percentage of positive cells was graded as 0 (negative), 1 (up to 10%), 2 (11-50%), 3 (51-80%), or 4 (>80% positive cells) and staining intensity as 0 (no staining), 1 (weak), 2 (moderate), or 3 (strong). The final immunoreactive GIS was defined as the multiplication of both grading results (percentage of positive cells * staining intensity) ( Supplementary Table 2 ).

The statistical data is consolidated and implemented by R-4.0.2. Continuous variables between two groups were compared using the unpaired Student t-test and Mann-Whitney U test for parametric data and non-parametric data, respectively. Survival curves were constructed using the Kaplan-Meier method, and the differences between the survival curves were examined by the log-rank test. Univariate Cox proportional hazards regressions were applied to estimate the individual hazard ratio (HR) for overall survival (OS). The NEC score of independent of various clinical features was performed by the univariate and multivariate Cox regression analysis. The receiver operating characteristic curve and the area under the curve (AUC) were calculated to the prediction accuracy of NEC score. All reported P values were two bsided and P<0.05 was considered statistically significant.

The flowchart of the detailed identification are shown in Figure 1 . In order to obtain differentially expressed necroptosis-related genes in HCC, the expression data of HCC were collected from the TCGA-LIHC dataset. The results suggested that 19 necroptosis related genes (CASP8, HSP90AA1, RNF31, NR2C2, HSPA4, USP22, TNFRSF21, SLC39A7, TSC1, SQSTM1, TRIM11, TRAF2, DNMT1, EZH2, LEF1, PLK1, MYCN, CDKN2A, and TERT) were upregulated while 3 necroptosis related genes (ID1, ALDH2, and BACH2) were downregulated in HCC versus paired non-tumor tissues (|log2FC|>1 and adjusted P <0.05) ( Supplementary Table 3 ). Next, to determine the prognostic of necroptosis related genes in HCC, 365 cases of HCC with sufficient survival data were analyzed using the univariate Cox analysis. The result showed that high expression of 26 necroptosis-related genes had a shorter OS than the low expression, while 4 necroptosis related genes had a better OS than the low expression ( Supplementary Table 4 ). Finally, 15 differentially expressed necroptosis-related genes that were correlated with OS ( Figure 2A ). Figures 2B, C showed the expression and prognostic of 15 necroptosis related genes in HCC.

To explore the genetic mutation of these genes in HCC, we assessed the incidence of copy number variations (CNVs) and somatic mutations. We found that CDKN2A (8%) and TRIM11 (8%) showed the highest mutation frequency followed by TNFRSF21 ( Figure 2D ). Further analysis of CNV alteration frequency showed that TRIM11, SQSTM1, MYCN, and EZH2 had a high frequency of CNV amplification, whereas DNMT1, HSP90AA1, and CDKN2A mainly showed CNV depletion ( Figure 2E ). The location of CNV alteration of the modification regulators on the chromosome is shown in Figure 2F . We found that most of necroptosis-related genes (except DNMT1, HSP90AA1, and CDKN2A) with a high frequency of CNV gain were highly expressed in HCC patients, suggesting that CNVs may be a potential contributor to the regulation of the expression of necroptosis-related genes.

Next, we assessed the correlation between the expression of necroptosis-related genes and cancer-related pathways (apotoptosis, cell cycle, DNA damage resposone, EMT, PI3K/AKT, RTK). The results showed that the expression of necroptosis-related genes was positively correlated with activite cancer-related pathways ( Figure 2G ). Research on drug sensitivity linked the drug response to possible biological effects through analysis of gene expression profiles. To gain further insight into the effects of the necroptosis-related genes on drug sensitivity, we evaluated the correlation between necroptosis-related genes and drug sensitivity, which were downloaded from the GDSC and CTRP database. The results showed the correlation and significance of most of drugs with the highest correlation coefficients in GDSC ( Figure 2H ) and CTRP databases ( Supplementary Figure 1 ). These results implied that necroptosis-related genes are associated with drug sensitivity, and necroptosis-related genes might have potential applications in the development of novel chemotherapeutic drugs.

Next, 5 datasets (TCGA-LIHC, LIRI-JP, GSE14520, GSE36376, and GSE76427) with available OS data and clinical information were enrolled into one meta-cohort (N=1155). Based on 15 necroptosis-related genes in HCC, we conducted consensus clustering analysis to differentiate HCC patients. Two clusters (NEC.cluster.A and NEC.cluster.B) were suggested as the optimal clustering stability based on the similarity displayed by 15 necroptosis-related genes expression ( Figure 3A , Supplementary Table 5 ). Interestingly, patients in NEC.cluster.B had a significantly shorter OS than in NEC.cluster.A (P< 0.001; Figure 3B ). We then explored the difference in immune microenvironment of two clusters. Figure 3C showed the immune cell infiltration landscape in two clusters of HCC, which demonstrated that NEC.cluster.B was correlated with high abundance of immune infiltration levels ( Figure 3C , Supplementary Table 6 ). Immune cells remained in the stroma around tumor cell nests rather than penetrating their parenchyma. Stroma activation in TME is thought to be T cell inhibitory (32). The results showed that patients in NEC.cluster.B group had a significantly higher stroma activity than in NEC.cluster.A group ( Figure 3D ).

In order to explore the differential expression of genes between NEC.cluster.A and NEC.cluster.B. Based on the cutoff criterion of |log FC|≥1.0 and adj.P<0.05, there are 7987 DEGs between NEC.cluster.A and NEC.cluster.B ( Supplementary Table 4 ). Next, to determine the prognostic of DEGs in HCC, 1155 cases of HCC with sufficient survival data were analyzed using the univariate Cox analysis. Finally, 4000 DEGs that were correlated with OS ( Supplementary Table 7 ). Next, we performed GO and KEGG enrichment analysis to understand the potential functions of these 4000 genes. As a consequently, these genes showed widespread association with catabolic process, metabolic process in GO analysis ( Figure 3E ). Furthermore, these genes showed widespread association with cell cycle and metabolism in KEGG pathway analysis ( Figure 3F ).

Next, based on 4000 genes in HCC, we conducted consensus clustering analysis to differentiate HCC patients. Two clusters (NEC.gene.cluster.A and NEC.gene.cluster.B) were suggested as the optimal clustering stability based on the similarity displayed by 4000 genes expression ( Figure 4A , Supplementary Table 8 ). Moreover, patients in NEC.gene.cluster.B group had a significantly shorter OS than those in NEC.gene.cluster.A group (P< 0.001; Figure 4B ). We then explored the difference in immune microenvironment of this two clusters. Figure 4C showed the immune cell infiltration landscape in two clusters of HCC, which demonstrated that NEC.gene.cluster.B was correlated with high abundance of immune infiltration levels ( Supplementary Table 9 ). Moreover, patients in NEC.gene.cluster.B group had a significantly higher stroma activity than in NEC.gene.cluster.A group ( Figure 4D ), and the significant overexpressions of 15 necroptosis-related genes in NEC.gene.cluster.B were observed relative to NEC.gene.cluster.A ( Figure 4E ). The above results proved the effectiveness and stability of NEC patterns.

Next, to quantify the NEC index of each HCC patient, we constructed a scoring system using the PCA method based on NEC phenotype-related prognostic genes ( Supplementary Table 10 ). The results showed that compare with NEC.cluster.A ( Figure 5A ), and NEC score was significantly higher in NEC.cluster.B, and compare with NEC.gene.cluster.A, NEC score was significantly higher in NEC.gene.cluster.B ( Figure 5B ). In addition, patients with high NEC score had a significantly shorter OS than those with low NEC score (P< 0.001; Figure 5C ), and AUC of the ROC was 0.754 ( Figure 5D ). We then explored the difference in immune microenvironment of this two groups. Figure 4C showed the immune cell infiltration landscape in two groups of HCC, which demonstrated that high NEC score group was correlated with high abundance of immune infiltration levels ( Figure 5E , Supplementary Table 11 ). Moreover, patients in high NEC score group had a significantly higher stroma activity than in low NEC score group ( Figure 5F ). Next, GSEA was performed to explore HCC involved signaling pathways between two groups. The results indicated that the high NEC score group was more likely to be enriched in RNA translation and cancer-specific pathways, including translational initiation, notch signaling pathway, ribosome, and spliceosome ( Figures 5G, H ).

Considering clinical characters, we then performed univariate and multivariate analysis and identified NEC score as independent factors impacting on HCC patients’ prognosis in TCGA cohort (HR 1.35; 95% CI 1.08, 1.78; P=0.021) ( Figure 6A ) and GSE14520/NCI cohort (HR 1.32; 95% CI 1.15, 1.49; P<0.001) ( Figure 6B ). To detect NEC score in TCGA immune subtypes, we downloaded the TCGA immune subtypes of HCC patients from the TCGA database, and we found patients with high NEC score were characterized wound healing, IFN-gamma dominant, and inflammatory phenotype ( Figure 6C ). We also detect the association between NEC score and clinicopathological features of patients with HCC in TCGA cohort and GSE14520/NCI cohort. We found that the NEC score increased with an increase in the HCC stage, which translates to a worse prognosis in TCGA cohort ( Figure 6D ) and GSE14520/NCI cohort ( Figure 6E ).

To validate the stability of the NEC score model, we applied the NEC score established in the meta cohort to 7 external independent GEO or TCGA HCC datasets. The result showed that patients with a high NEC score had a poor prognosis in all GEO dataset (HR 1.89; 95% CI 1.47, 2.42; P<0.001) ( Figure 7A ), and the AUC of the ROC was 0.760 ( Figure 7B ). In addition, patients with a high NEC score had a poor prognosis in all TCGA dataset (HR 2.00; 95% CI 1.38, 2.90; P<0.001) ( Figure 7C ), LIRI-JP dataset (HR 2.54; 95% CI 1.35, 4.79; P=0.004) ( Figure 7E ), GSE14520 (HR 2.16; 95% CI 1.36, 3.44; P-0.001) ( Figure 7G ), GSE20140 (HR 1.82; 95% CI 1.01, 3.26; P=0.045) ( Figure 7I ), GSE27150 (HR 2.58; 95% CI 1.17, 5.65; P=0.018) ( Figure 7K ), GSE36376 (HR 1.79; 95% CI 1.12, 2.85; P=0.015) ( Figure 7M ), and GSE76427 (HR 2.00; 95% CI 1.03, 3.90; P=0.041) ( Figure 7O ), and the AUC of the ROC was 0.747 in TCGA dataset ( Figure 7D ), 0.688 in LIRI-JP dataset ( Figure 7F ), 0.754 in GSE14520 ( Figure 7H ), 0.760 in GSE20140 ( Figure 7J ), 0.862 in GSE27150 ( Figure 7L ), 0.724 in GSE36376 ( Figure 7N ), and 0.685 in GSE76427 ( Figure 7P ). These results demonstrated that NEC score is an effective and stable model and can predict the prognosis of patients with HCC.

Firstly, after consensus clustering analysis of these samples by ConsensusClusterPlus algorithm base on 15 necroptosis-related genes in IMvigor210 cohort, it revealed that the results were most stable at K=2. Then, 1278 DEGs that were correlated with OS were identified ( Supplementary Table 12 ). Finally, we constructed a NEC score base on 1278 genes ( Supplementary Table 13 ). The result showed that patients with a high NEC score had a poor prognosis in IMvigor210 dataset (HR 3.48; 95% CI 2.56, 4.72; P<0.001) ( Figure 8A ), and the AUC of the ROC was 0.775 ( Figure 8B ). We further detect the association between NEC score and immunotherapy. We found that the proportion of CR/PR patients in low NEC score group was significantly higher than that in high NEC score group ( Figure 8C , P < 0.05), and the NEC score in CR/PR patients was significantly lower than that of SD/PD ( Figures 8D, E , Table 2 , P<0.05). Moreover, the low NEC score group had a high checkpoint expression levels. Above results can partly explain the poor prognosis in high NEC score group.

Additionally, to compare the prognostic analysis effect of NEC score in patients who had a high or low neoantigen burden, survival analysis was also performed on patients in the IMvigor210 cohort. We revealed that among both patients who had a high or low neoantigen burden, patients with a high NEC score showed lower overall survival ( Figure 8G , P<0.05), and the AUC of the ROC was 0.758 ( Figure 8H ). To identify the association of the immune subtype with the NEC score, we traced each sample from the immune subtype to the NEC score. As shown in Figure 8I , patients with the inflamed subtype were assigned a low NEC score, which was consistent with better prognosis in these groups. We then explored the difference in immune microenvironment of this two groups. Figure 8J showed the immune cell infiltration landscape in two groups, which demonstrated that high NEC score group was correlated with high abundance of immune infiltration levels ( Figure 8J ). Moreover, patients in high NEC score group had a significantly higher stroma activity than in low NEC score group ( Figure 8K ).

Next, we used a HCC immune cohort (GSE140901) to detect the relationship between NEC score and immunotherapy. NEC score also varied statistically in the PR, SD, and PD groups, NEC score was notably lower in PR group than PD or SD groups ( Figure 9A ). We also discovered that NEC score was notably lower in clinical benefit response group than no clinical benefit response group ( Figure 9B ). The abovere sults suggests that NEC score was sensitivity to immunotherapy. We also compared the differential level of NEC score in different subgroups stratified by gender and status. As shown in Figures 9C, D , NEC score was no significantly betweent male and female, however, NEC score was significantly higher in dead group than in alive group.

To further validate necroptosis-related genes expression in HCC, IHC was used to measure the expression level of necroptosis-related genes in HCC (N=8), and the result showed that compared with normal group, the CASP8, CDKN2A, DNMT1, EZH2, HSP90AA1, HSPA4, MYCN, PLK1, SLC39A7, SQSTM1, TNFRSF21, TRAF2, TRIM11, and USP22 level were significantly higher in HCC group, while ALDH2 was downregulated in HCC versus normal tissues ( Figures 10A, B ).