A Necroptosis-Related Gene Signature for Predicting Prognosis, Immune Landscape, and Therapy Sensitivity in Hepatocellular Carcinoma

Background: Hepatocellular carcinoma (HCC) remains a growing threat to global health. Necroptosis is a newly discovered regulated cell necrosis that plays a vital role in cancer development. Thus, we conducted this study to develop a predictive signature based on necroptosis-related genes. Methods: The tumor samples in The Cancer Genome Atlas (TCGA) Liver Hepatocellular Carcinoma (LIHC) cohort were subtyped using the consensus clustering algorithm. Univariate Cox regression and LASSO-Cox analysis were performed to construct a gene signature model from differentially expressed genes between tumor clusters. Then we integrated TNM stage and the prognostic model to build a nomogram. The gene signature and the nomogram were externally validated in the GSE14520 cohort from the gene expression omnibus (GEO) and LIRP-JP cohort from the International Cancer Genome Consortium (ICGC). Predictive performance evaluation was conducted using Kaplan-Meier plot, time-dependent receiver operating characteristic curve, principal components analysis, concordance index, and decision curve analysis. The tumor microenvironment was estimated using seven published methods. Finally, we also predicted the drug responses to immunotherapy, conventional chemotherapy and molecular-targeted therapy using two algorithms and two datasets. Results: We identied two necroptosis-related clusters and a ten-gene signature (MTMR2, CDCA8, S100A9, ANXA10, G6PD, SLC1A5, SLC2A1, SPP1, PLOD2, and MMP1). The gene signature and the nomogram had good predictive ability in TCGA, ICGC, and GEO cohorts. The risk score was positively associated with the degree of necroptosis and immune inltration (especially immunosuppressive cells). The high-risk group could benet more from immunotherapy. Chemotherapy and molecular-targeted therapy should be adapted to the molecular proles of each patient. Conclusion: The necroptosis-related gene signature provides reliable evidence for prognosis prediction, comprehensive treatment, and new therapeutic targets for HCC patients. The nomogram can further improve predictive accuracy.


Background
Hepatocellular cancer (HCC), the fourth most common cause of cancer mortality, is projected to result in the death of more than 1 million people in 2025 worldwide [1]. Although substantial progress has been achieved in the surveillance and systemic treatment of HCC, it is still the second deadliest cancer with a dismal overall ve-year survival rate of 18% [2]. Due to different etiological factors and the heterogeneous nature of HCC, prognostic prediction based on TNM strategy alone is challenging [3]. Furthermore, new therapeutic targets are urgently warranted to guide therapies and optimize medical interventions. Thus, a novel prognostic model utilizing molecular pro les of HCC is highly desirable.
Necroptosis is a newly discovered regulated cell necrosis mainly mediated by necrosome comprised of mixed lineage kinase domain-like protein (MLKL), receptor-interacting protein kinases (RIPK1) and RIPK3 [4]. Different from apoptosis in morphology, it is characterized by cell swelling and membrane perforation with subsequent release of intracellular damage-associated molecular patterns (DAMPs) that induce in ammation. In the context of caspase inhibition, necroptosis is initiated and functions as a backup system to combat pathogens or cancer cells that evade the apoptotic process, which is also one of the vital mechanisms for some antineoplastic drugs, such as 5-FU, etoposide, and cisplatin [5][6][7][8]. In addition to the defensive effect against malignancy, accumulating evidence suggests that necroptosis can promote tumorigenesis, cancer metastasis and immunosuppression. For example, RIPK1 is overexpressed in glioblastoma and associated with poor prognosis via attenuating p53 activation [9].
Knockout of RIPK1, RIPK3, and MLKL in the colon and esophageal cancer cells can signi cantly suppress tumor growth for decreased activity of NF-kappaB (NF-κB) [10]. Taken together, the role of necroptosis in cancer is highly complex and ambivalent. However, the relationship between the expressions of necroptosis-related genes and overall survival (OS) remains largely unknown in HCC patients.
In the present study, a prognostic model was established according to differentially expressed genes (DEGs) between two necroptosis-related HCC subtypes. Next, the molecular signature was integrated with the TNM staging system to build a nomogram for better predictive ability. Ultimately, we assessed the performance of risk score in forecasting tumor microenvironment (TME) and sensitivities to conventional chemotherapy, molecularly targeted therapy, and immunotherapy. Our ndings provide reliable evidence for the risk strati cation, comprehensive treatment, and new therapeutic targets for HCC patients.

Dataset Acquisition
The analysis process was summarized in the ow chart (Additional le 1: Figure S1). Data of 371 HCC samples in Liver Hepatocellular Carcinoma (LIHC) cohort were downloaded from the Cancer Genome Atlas (TCGA) repository, including fragments per kilobase of transcript per million mapped reads (FPKM) normalized values, corresponding clinical information, single nucleotide variation (SNV) and copy number variation (CNV) data (https://portal.gdc.cancer.gov/repository). Only 341 patients with follow-up time more than 30 days were subjected to prognostic model tting. For external validation of this model, the RMA-normalized matrix of GSE14520 (including 242 tumor samples) was downloaded from Gene Expression Omnibus (GEO), and raw count data of LIRI-JP (including 231 tumor samples) from International Cancer Genome Consortium (ICGC). (https://www.ncbi.nlm.nih.gov/geo/, https://dcc.icgc.org/projects/LIRI-JP). Next, the mRNA sequencing data of TCGA-LIHC and ICGC-LIRI-JP were converted into transcripts per kilobase million (TPM) values, and then log2(x+1) transformed, which was suggested to be the most accurate quanti cation method with minimized statistical biases [11,12].
Batch effects of the above three datasets were eliminated by the "Combat" function of the "sva" package. Due to only 50 matched tumor-normal samples in TCGA-LIHC, another four datasets were included to validate the transcription levels of prognostic genes ltered by least absolute shrinkage and selection operator (LASSO)-penalized Cox regression, including GSE14520 (213 pairs), GSE25097 (243 pairs), GSE57957 (37 pairs), GSE76297 (58 pairs).

Investigation of the Expressions of Necroptosis-Related Genes
Twenty-ve necroptosis-related genes were extracted from the MSigDB team (GOBP_NECROPTOTIC_SIGNALING_PATHWAY) and previous reviews, which are presented in the Additional le 5: Table S1 [13][14][15][16][17]. SNV and CNV analyses were performed utilizing the "maftools" package and the "Rcircos" packages, respectively [18,19]. A network was plotted by the "igraph" package to summarize the correlation of the expression levels among necroptosis-related genes and their prognostic values in the TCGA cohort.

Consensus Clustering
Consensus clustering is a robust unsupervised classi cation technique achieved through multiple resampling and clustering. Similarity within each group was measured by Euclidean distance, and the whole process was repeated 1000 times using the "ConsensuClusterPlus" package [20]. The proper cluster number (k) was determined according to the cumulative distribution function (CDF) plots. The atness of the CDF curve indicates the level of consensus and the stability of clustering, which was further veri ed in principal components analysis (PCA) analysis.

Identi cation of DEGs between two necroptosis-related clusters
Assessment of DEG between two necroptosis-related clusters was identi ed by the "limma" package [21]. The false discover rate (FDR) < 0.05 and |log2FC| > 1 were criteria for statistical signi cance. The biological process of Gene Ontology (GO) enrichment was analyzed using the Metascape website (http://metascape.org). Based on the reference gene set available from the MSigDB (c2.cp.kegg.v7.4.symbols), we converted the DEG expression data into a pathway score matrix using the "GSVA" package [22]. The pathway activities were also compared by the "limma" package.

Development and Validation of Necroptosis-Related Prognostic Model
The univariate Cox analysis was carried out to select DEGs associated with overall survival. Next, LASSOpenalized Cox regression was employed to lter features to establish a risk model by the "glmnet" package [23,24]. The following formula was used to calculated risk score for each patient: . According to the median risk score, patients in training cohorts were classi ed into high-or low-risk groups. Similarly, the risk scores of external validation cohorts (ICGC and GEO) were calculated using the same formula developed from the training group and divided based on the median value of the TCGA cohort.

Establishment and Validation of Nomogram
We combined univariate and multivariate Cox analyses to determine the independent prognostic factors with statistical signi cance, which were subsequently incorporated to t a Cox proportional hazard model [25]. For the convenience of clinical utility, a nomogram was constructed to visualize the results of the Cox model. Next, we carried out external validations in ICGC and GEO cohorts to con rm the reliability of the nomogram. The total points of each patient in the validation cohorts were calculated according to the Cox model derived from the TCGA cohort. Then Cox regression were implemented again using the total points as a factor, in turn, the concordance index (C index) and calibration curve were yielded based on the Cox results [26]. In addition, we also plotted time-dependent receiver operating characteristic (ROC) curve and performed decision curve analysis (DCA).

Determination of Immune Landscape
For deciphering the immune cell in ltration, we employed single-sample gene set enrichment analysis (ssGSEA) [27] and six deconvolution algorithms to fully leverage the transcription data, including CIBERSORT [28], MCPCounter [29], xCell [30], EPIC [31], quanTIseq [32], TIMER [33]. The gene set used by ssGSEA analysis was obtained from the previous study [34], and the whole process of TME estimation was carried out using the "IOBR" package, which integrates the above seven methodologies [35].

Prediction of drug sensitivities to chemotherapy and immunotherapy
With an aim to provide personalized treatment recommendations, we explored the predictive value of risk score in drug sensitivities. First, molecular pro les and z-scored half-maximal inhibitory concentration (IC50) data of FDA-approved drugs were downloaded from the CellMiner database (https://discover.nci.nih.gov/cellminer) [36]. The log2(FPKM+1) data of 60 human cancer cell lines (NCI60) were transformed into log2(TPM+1) values, followed by risk calculation based on the model derived from TCGA and subsequent correlation analyses. Next, we employed the "oncoPredict" package [37] to build ridge regression models to examine the differences in the e cacy of commonly used agents between high-and low-risk groups. Cancer Therapeutics Response Portal V2 (CTRP V2) database was used as the training cohort, of which the log2(TPM+1) and half-maximal effective concentration (EC50) data were collected and shared by the "oncoPredict" team (https://osf.io/c6tfx).
As for immunotherapy, we rst uploaded the normalized sequencing data of TCGA to the Tumor Immune Dysfunction and Exclusion (TIDE) website (http://tide.dfci.harvard.edu/), which is a computational prediction method based on tumor pre-treatment expression pro les [41]. Due to the lack of published data on immune checkpoint inhibitor (ICI) in HCC, we calculated the risk score for each patient with metastatic urothelial cancer in the IMvigor210 cohort [42]. The risk score distribution between complete/partial response (CR/PR) and stable/progressive disease (SD/PD) groups was further explored.

Statistical Analysis
All the data analyses and visualizations were accomplished in the R environment (Version 4.12). Wilcox statistical test was used to compare differences between two groups. Categorical data was compared using the Chi-square test. The correlation coe cient was determined by the nonparametric Spearman approach. PCA was completed by the "prcomp" function of the "stats" package. Kaplan-Meier (K-M) analysis, log-rank test, and Cox regression were performed using the "survival" and the "survminer" packages. C index was calculated for each Cox model and pooled using the random-effects mixed model with maximum likelihood [43]. Time-dependent ROC curve was created by the "timeROC" package to estimate the predictivity of each prognostic factor. The "rms" package was employed to generate the nomogram and the calibrate plot. Two-tailed P < 0.05 was considered signi cant.

Overview of genetic and transcriptional variation of necroptosis-related genes in HCC
We rst explored the SNV and CNV of 25 necroptosis-related genes in the TCGA-LIHC cohort. At the genetic level, 31 of 364 patients had mutations in necroptosis-related genes, of which TLR3, RIPK3 and TYRO3 showed the highest frequency (1%), while six genes (IFNA1, IPMK, IPPK, PELI1, TRADD and TRAF2) did not display any alterations (Fig. 1A). As illustrated in Fig. 1B and 1C, genomic instability was widespread in necroptosis-related genes, including 12 CNV gains, 12 CNV losses, and one nonsigni cant alteration. The comparison between 50 normal and 371 HCC samples indicated that 13 genes were signi cantly differentially expressed, including upregulation of 12 genes and downregulation of 1 gene in tumor (Fig. 1D). The same tendency could also be observed in 50 patients with matched tumornormal samples (Fig. 1E). Additionally, we found that the gains or losses of copy number were not always positively correlated with the expression levels of pyroptosis-related genes, such as MLKL and TLR3, suggesting that CNV is not the only regulatory factor of gene expression [44]. Some mechanisms, including DNA methylation, transcription factor, m6A modi cation, long non-coding RNAs, and RNA binding proteins, also play vital roles in gene expression [45][46][47][48]. Fig. 1F depicts the correlations of these genes and their prognostic signi cances. Collectively, our results indicated that the expression of pyroptosis-related genes was signi cantly different between HCC and normal samples, which could be a potential contributing factor to tumorigenesis and heterogeneity.

Identi cation of necroptosis-related clusters in HCC
According to the results of k-means consensus clustering, we identi ed two molecular subtypes (134 patients in cluster 1 and 207 patients in cluster 2) based on necroptosis-related genes, which was also validated by PCA analysis ( Fig. 2A; Additional File 2: Figure S2; Additional File 6: Table S2). Compared with cluster 2, the overall survival (OS) of cluster 1 was signi cantly worse (P < 0.001), and most necroptosis-related genes were upregulated (Fig. 2B, C). As presented in Fig. 2E, there were signi cant differences between the two clusters in transcription pro le, age, T stage, TNM stage, tumor grade, and gender. The top 100 of 1517 DEGs were plotted in the heatmap (Additional le 7: Table S3). Additionally, GSVA analysis revealed that the pathways enriched in the C2 cluster were mainly associated with metabolism, biosynthesis, and degradation. In contrast, cluster 1 showed enrichment in tumorigenesis, cell cycle, DNA replication, DNA repair, spliceosome, and immune activation, which was consistent with ndings of GO enrichment (Fig. 2D, F; Additional le 9: Table S4).
3.3 Development and validation of a ten-gene signature based on DEGs between necroptosis related clusters.
According to the median risk value of the TCGA cohort, patients were split into low-and high-risk groups, which were distributed in two directions in the PCA plot (Fig. 3A, E). The scatter plot indicated that risk score was negatively associated with survival time (Fig. 3B). Consistently, the K-M curves suggested that the prognosis of the high-risk group was signi cantly worse than that of the low-risk group (P < 0.001;  (Fig. 3D). Of note, the high-risk group was mainly composed of the C1 cluster and had signi cantly higher expression in most necroptosis-related genes, while the low-risk group predominantly consisted of the C2 cluster (Fig. 4A, B). The Sankey diagram also illustrated that the risk strati cation based on the 10-gene signature was not entirely congruent with the result of the TNM system, which could provide new insight into HCC management. The Protein-protein interaction network between the ten selected genes and necroptosis-related genes is depicted in Additional le 4: Figure S3D 3.4 External validation of the expression and genetic alteration of the ten genes We rst investigated the transcription levels of the ten genes in the TCGA, ICGC, and GEO cohorts, nding that risk score was positively correlated with the expressions of nine genes except for AXAN10 (Fig. 5A, B, C). As depicted in the OncoPrint generated from cBioportal, 70 out of 366 (19%) showed genetic alterations in the ten genes, of which ampli cation was the most prevalent type (Fig. 5D). Typical immunohistochemical images of nine genes were downloaded from Human Protein Atlas, except for MMP1, which was unavailable from the database (Fig. 4C). Subsequently, ve datasets were used to further validate the expression levels of the ten genes between tumors and matched adjacent normal tissues (Fig. 5E-I). Interestingly, S100A9 was paradoxically overexpressed in adjacent normal tissues, while the expressions of the other nine genes were consistent with trends found between high-and lowrisk groups. S100A9 protein, also known as myeloid-related protein, typically resides in immune cells such as neutrophils, macrophages, and monocytes [49]. Although S100A9 was overexpressed in liver cancer cells, it is also intensely upregulated in myeloid cells under pathological conditions, which may account for the higher expression in adjacent normal tissue [50]. So, we further compared the tumor microenvironment of tumor and adjacent normal tissues and con rmed that the fractions of neutrophils, macrophages M2, and monocytes were signi cantly higher in adjacent normal tissues (Additional File 4: Figure S4).

Construction and Validation of Predictive Nomogram
Univariate and multivariate Cox analyses identi ed risk score and TNM stage as independent prognostic factors in TCGA, ICGC, and GEO cohorts, which were included to establish a nomogram forecasting prognosis (Fig. 6A, B, Fig. 7A). The pooled C indexes of three cohorts demonstrated that nomogram (0.74, 95% CI: 0.69-0.79) had better predictive performance over risk score (0.71, 95% CI: 0.65-0.77) and TNM stage (0.66, 95% CI: 0.63-0.70) alone (Fig. 7B). The calibration plot illustrated nomogram achieved good consistency between the predicted and observed OS outcomes (Fig. 7C). The time-dependent ROC curves showed that the nomogram had greater AUC values, further supporting the high predictive consistency of the nomogram (Fig. 7D). According to the DCA analyses, the net bene t of the nomogram was evident in most cases (Fig. 8).

Analyses of TME and Drug Sensitivity
In view of the importance of risk score in prognosis prediction of HCC patients, we next investigated its potential value in re ecting TME and guiding individualized treatment. Seven published methods were used to estimate the correlation between risk score and immune in ltration in the TCGA cohort. Most cells in tumor microenvironment were positively correlated with the risk score, especially cells facilitating HCC development and immunosuppression, including type 2 T helper (Th2) cells, type 17 T helper (Th17) cells, T regulatory (Treg) cells, myeloid-derived suppressor cells (MDSCs), neutrophils, cancer-associated broblasts (CAFs) (Fig. 9A, C) [51,52]. Macrophages M0, B cells, dendrite cells, and T follicular helper cells exhibited the same trend. By contrast, the positive correlations between risk score and tumor-killing cells were weaker such as gamma delta T cells (Tγδ), CD8+ T cells, natural killer T (NKT) cells. In addition, we found that 37 immune checkpoint-related genes were overexpressed in the high-risk group, which further impaired the anti-tumor immunity (Fig. 9B). The risk score had positive correlations with seven genes with corresponding ICIs, including PDCD1 (PD1), CD274 (PD-L1), CTLA4, LAG3, HAVCR2, TIGIT, VSIR (Fig. 9D). Subsequently, we utilized the TIDE algorithm to evaluate the response to immunotherapy. The risk score was negatively correlated with the TIDE score, and the high-risk group had a signi cantly lower TIDE score, suggesting that the high-risk group was more sensitive to immunotherapy (Fig.9E, F). Besides, we calculated risk scores for patients in the IMvigor210 cohort, nding that responders (CR/PR) to anti-PD1/PD-L1 treatment had signi cantly higher risk scores (Fig. 9G).
Based on the CellMiner database, we used our prognostic model to calculate the risk scores for 60 cancer cell lines and analyzed its correlations with 218 FDA-approved drugs. A lower EC50 or IC50 value indicates a higher drug e cacy. The top 16 drugs with the most statistical differences are presented in Fig. 10. Except for rapamycin, the other 15 drugs were associated with better therapeutic effects, including oxaliplatin, parthenopid, everolimus, nitrogen mustard, actinomycin D, daunorubicin, nelarabine, cyclophosphamide, sulfatinib, lomustine, hydroxyurea, dexrazoxane, belinostat, decitabine, crizotinib. To better link the necroptosis-related gene signature with clinical practice, we used the "oncoPredict" package to predict commonly used drugs in the clinical treatment of HCC. The results showed that regorafenib, tivantinib, cabozantinib, cediranib, olaparib, navitoclax, mitomycin, vincristine, and paclitaxel were more sensitive in the high-risk group, while erlotinib, brivanib, dasatinib, linifanib, neratinib, uorouracil, and methotrexate were more sensitive in the low-risk group (Fig. 11).

Discussion
Liver cancer remains a growing threat to global health. HCC accounts for approximately 90% of cases of liver cancer and is the fastest rising cause of cancer mortality [1]. Currently, the TNM staging is the most widely used system of prognosis prediction of HCC that only considers histopathological factors. Besides, given the excessive heterogeneity within HCC, individual-tailored treatment and more accurate risk strati cation are critical to the improvement of therapy e cacy and overall survival. Therefore, it is highly desirable to develop a panel of molecular markers to supplement the TNM system.
Necroptosis is a pro-in ammatory form of regulated cell death and serves as an alternative defensive mechanism in the context of apoptosis inhibition caused by caspase deactivation [53]. As apoptosis escape is a well-established feature of tumors, necroptosis is consequently initiated via the formation of necrosome (RIPK1-RIPK3-MLKL complex) but may exert a double-edge sword effect on carcinogenesis [14]. For example, downregulation of MLKL implies poor prognosis in ovarian and colon cancer, while increased expression of RIP1 confers a worse prognosis in glioblastoma [54][55][56]. Paradoxically, one study found the pro-tumor effect of RIPK1 in HCC patients, but most human hepatoma cell lines suppress necroptosis by methylation-dependent loss of RIPK3 expression, including Huh-7, HepG2, and Hep3B. [14,57]. So large cohort and comprehensive analysis are warranted to investigate the role of necroptosis in HCC patients to guide treatment.
Sample classi cation based on the transcription pro les of prede ned genes is a widely applied method. Based on the TCGA-LIHC cohort, we identi ed two necroptosis-related clusters and found that overexpression of these genes was related to decreased OS. To achieve better predictive ability, our prognostic model was developed from DEGs between the two subtypes rather than only from necroptosis-related genes. Consistently, the risk score was positively associated with the expression levels of most necroptosis-related genes. It has been reported that necroptosis leads to in ammation, necrosis, brosis, and eventually oncogenesis in patients with viral hepatitis, alcoholic liver disease, and non-alcoholic fatty liver disease [14]. Additionally, our results revealed that the degrees of in ltration of most immune cells and necroptosis were positively correlated, especially pro-tumor cells. MDSCs, Th2, and Treg cells can compromise normal immune surveillance and contribute to immune escape, while CAFs secrete growth factors and cytokines that favor tumor growth [51,58]. The in ux of neutrophils mediated by Th17-produced interleukin-17 can drive tumor progressions via the formation of neutrophil extracellular trap (NET) [59,60]. Collectively, we concluded that necroptosis was more inclined to cause HCC progression instead of inhibition.
In the present study, the novel signature was externally validated in ICGC and GEO cohorts and included two genes encoding downstream components of the necroptosis-related pathways, as RIPK1-dependent NF-κB activation will enhance the expression of S100A9 and MMP1 [61-63]. Except for MTMR2, the other nine genes were previously reported to involve in HCC, which could be roughly classi ed into ve categories, including cell cycle (CDCA8), chronic in ammation (S100A9), anti-oncogene (ANXA10), energy metabolism (G6PD, SLC1A5, SLC2A1), and extracellular matrix organization (SPP1, PLOD2, and MMP1). As a member of the myotubularin family of phosphoinositide lipid phosphatases, MTMR2 promotes invasion and metastasis of gastric cancer and NK/T cell lymphoma via IFNγ/STAT1/IRF1 pathway and targeting JAK1, respectively [64,65]. CDCA8 protein is a crucial regulator of mitosis and overexpressed during hepatocarcinogenesis [66]. S100A9 is a DAMP molecule that can amplify in ammation in tumor microenvironment and lead to malignancy progression [63]. G6PD encodes glucose-6-phosphate dehydrogenase that limits the rate of the pentose phosphate pathway, which can induce EMT in HCC cells through the STAT3 pathway [67]. Increased glutamine metabolism is one of the critical metabolic modes of cancer cells [68]. Upregulation of SLC1A5 is accolated with increased uptake of glutamine and dismal prognosis in hepatic cancer [69]. Glucose transporter 1 (GLUT1) encoded by SLC2A1 is the primary regulator of glucose uptake and contributes to the metastasis and chemoresistance of liver cancer [68]. Osteopontin (OPN) was rst isolated from osteoblasts and highly associated with bone mineralization [70]. In liver tissues, OPN is encoded by SPP1 and responsible for liver regeneration and correlates with the extent of liver cirrhosis and tumor growth [71]. PLOD2 protein, also known as lysyl hydroxylase, is essential for the maturation of collagen and an independent prognostic factor of HCC [72,73]. Under hypoxia environment, increased PLOD2 expression promoted migration of pancreatic cancer cells via remolding extracellular matrix [74]. Matrix metalloproteinase 1 (MMP-1) degrades the pericellular matrix, leading to enhanced invasion and metastasis [75].
There are some strengths of our studies. First, we linked necroptosis to the prognosis of HCC. Second, our prognostic model re ected the tumor immune microenvironment. Most immune cells were highly correlated with the risk score, especially pro-tumor cells. Third, our risk model could guide clinical medication. A higher risk score implied that patients could bene t more from ICI treatment. As sensitivities of drugs varied between high-and low-risk groups, chemotherapy and molecular-targeted therapy should be adapted to the sequencing result of each patient. Additionally, with the advancement of circulating tumor cell (CTC) detection and single-cell sequencing, we can also dynamically adjust the treatment regime according to the sequencing results. Fourth, the ten-gene signature and the nomogram were developed based on the TCGA cohort and further validated externally in the GEO and ICGC cohorts, which increased the reliability and stability of our results. The performance of the nomogram integrating TNM stage and 10-gene signature was superior to the two independent prognostic predictors alone.
Although we conducted comprehensive analyses and multi-database validations, several limitations inherent in this study should be noted. First, LASSO penalization would ignore some crucial genes regulating hepatocarcinogenesis. As the drug sensitivity predictions were based on bioinformatic analyses, large-scale clinical trials are warranted to con rm our ndings. Moreover, the mechanism of MTMR2 in HCC remains to be elucidated for lack of relevant research.

Conclusion
In summary, we developed a ten-gene signature based on DEGs between two necroptosis-related subtypes of HCC and further evaluated its predictive ability in prognosis, tumor microenvironment, and drug sensitivity. After external validation, we found that this signature was an independent prognostic factor that contributed to predicting OS and guiding clinical treatment. The Nomogram that combined histopathological and molecular features further improved the accuracy of prognostic prediction.

Acknowledgments
We would like to thank reviewers for helpful suggestions.    The association between gene signature and necroptosis Page 23/30 (A) Sankey diagram depicting the association between tumor cluster, risk score, and TNM stage in TCGA-LIHC cohort. (B) The expression levels of necroptosis-related genes between high-and low-risk groups in TCGA-LIHC cohort (Wilcox test, * P < 0.05; * * P < 0.01; * * * P < 0.001; * * * * P < 0.0001). Only genes with statistical signi cance were displayed in this plot. (C) Immunohistochemical images of the nine genes in the risk model obtained from the Human Protein Atlas (https://www.proteinatlas.org/).

Figure 5
Page 24/30 The expression and genetic alteration of the ten genes included in the risk model.
(A-C) The transcription levels of the ten genes and clinical characteristics between high-and low-risk groups in the TCGA-LIHC, ICGC-LIRI-JP, and GSE1450 cohorts. (D) Genetic alterations of the ten genes in TCGA-LIHC. The OncoPrint was downloaded from the cBioportal (https://www.cbioportal.org/). (E-I) Validations of expression levels of the ten genes between tumors and matched adjacent normal tissues in TCGA-LIHC, GSE14520, GSE25097, GSE57957, GSE76297cohorts.    Exploration of the tumor microenvironment (TME) and response to immunotherapy in TCGA cohort.
(A) Spearman correlation between the risk score and TME estimated by seven methods. (B) Comparisons of the expressions of immune checkpoint-related genes between high-and low-risk groups (Wilcox test, * P < 0.05; * * P < 0.01; * * * P < 0.001; * * * * P < 0.0001). (C) CIBERORT result of TME between high-and low-risk groups (Wilcox test). (D) Spearman correlation analyses between the risk score and seven immune checkpoint-related genes. (E) Spearman correlation between the Tumor Immune Dysfunction and Exclusion (TIDE) score and risk score. (F) Comparison of the TIDE score between high-and low-risk groups (Wilcox test). (G) Comparison of the risk score between responders and non-responders to anti-PD1/PD-L1 treatment in the IMvigor210 cohort.

Figure 10
Correlation analyses between drug sensitivity and risk score based on the CellMiner database