Identification and Validation of Prognosis-Related Necroptosis Genes for Prognostic Prediction in Hepatocellular Carcinoma

Background The prediction of hepatocellular carcinoma (HCC) survival is challenging because of its rapid progression. In recent years, necroptosis was found to be involved in the progression of multiple cancer types. However, the role of necroptosis in HCC remains unclear. Methods Clinicopathological parameters and transcriptomic data of 370 HCC patients were obtained from TCGA-LIHC dataset. Prognosis-related necroptosis genes (PRNGs) were identified and utilized to construct a LASSO risk model. The GEO cohorts (GSE54236 and GSE14520) were used for external validation. We evaluated the distribution of HCC patients, the difference in prognosis, and the accuracy of the prognostic prediction of the LASSO risk model. The immune microenvironment and functional enrichment of different risk groups were further clarified. Finally, we performed a drug sensitivity analysis on the PRNGs that constructed the LASSO model and verified their mRNA expression levels in vitro. Results: A total of 48 differentially expressed genes were identified, 23 of which were PRNGs. We constructed the LASSO risk model using nine genes: SQSTM1, FLT3, HAT1, PLK1, MYCN, KLF9, HSP90AA1, TARDBP, and TNFRSF21. The outcomes of low-risk patients were considerably better than those of high-risk patients in both the training and validation cohorts. In addition, stronger bile acid metabolism, xenobiotic metabolism, and more active immune cells and immune functions were observed in low-risk patients, and high expressions of TARDBP, PLK1, and FLT3 were associated with greater drug sensitivity. With the exception of FLT3, the mRNA expression of the other eight genes was verified in Huh7 and 97H cells. Conclusions. The PRNG signature provides a novel and effective method for predicting the outcome of HCC as well as potential targets for further research.


Introduction
Necroptosis is a unique type of inflammatory programmed cell death that occurs when apoptotic pathways are halted or inhibited, and TNF-α, Fas ligand, and tumor necrosis factorrelated apoptosis-inducing ligand stimulated death receptors-TNFR1, TNFR2, and FAS-are the most prevalent triggers [1]. e role of necroptosis in the progression of cancer is complicated. On the one hand, tumor cells can be eliminated directly by the process of necroptosis [2]. In addition, necroptosis supplies antigens and inflammatory stimuli to dendritic cells to kick-start acquired immunity, which then activates CD8 + T cells and antitumor immune responses [3]. On the other hand, inflammatory responses caused by cytokines released by necrotic cells can promote the development of tumors [4,5]. Recently, a prognosisrelated necroptosis gene signature was used to predict the prognosis and describe the immune microenvironment in human tumors [6,7].
Hepatocellular carcinoma (HCC) is one of the most aggressive tumors, and the risk factors predominantly include HBV and HCV infection, alcoholic liver disease, and nonalcoholic fatty liver disease [8]. It has been found that hepatocytes can form different types of tumors (intrahepatic cholangiocarcinoma and HCC), which are mainly determined by the mode of cell death (apoptosis and necroptosis) in the tumor microenvironment [9]. In addition, chronic inflammation is an important driving factor for the development of hepatocellular carcinoma, since the release of damage-associated molecular patterns associated with necroptosis can promote angiogenesis and cell proliferation, thus promoting tumor growth and metastasis [10]. Furthermore, the selection of targeted medicines as well as survival prediction remains major challenges because of the lack of effective molecular prognostic indicators. erefore, we identified differentially expressed genes (DEGs) associated with necroptosis in this study, which was applied to develop a predictive model and describe immunological features in various modes. Our findings will provide a thorough overview of necroptosis-related genes that may have a role in the development of HCC as well as a novel perspective on the clinical application of immunotherapy drugs.

Analysis of DEGs Associated with Necroptosis.
We identified 67 genes associated with necroptosis from a previous study [11], which are presented in Supplementary  Table 1. We utilized the "Limma" R package to normalize the mRNA expression levels and screened DEGs associated with necroptosis in TCGA-LIHC. e associations of DEGs were obtained by the protein-protein interaction network (PPI network; https://string-db.org/cgi/input.pl), and Cytoscape was used to visualize the results and obtain hub genes. We analyzed the relationship between all DEGs and prognosis using Gene Expression Profiling Interactive Analysis (GEPIA) platform (http://gepia.cancer-pku.cn/) and explored their mutation characteristics using GSCALite platform (http://bioinfo.life.hust.edu.cn/GSCA/).

Clustering Analysis of TCGA-LIHC According to DEGs.
We further screened out prognostic DEGs by univariate Cox regression (p < 0.05). Based on these DEGs, we performed a cluster analysis ("ConsensusClusterPlus" R package) and explored clinicopathological differences between different clusters.

Construction and Validation of the Least Absolute
Shrinkage and Selection Operator Regression Model. We used the "GLMnet" R package to build the least absolute shrinkage and selection operator (LASSO) regression model. e specific method was described previously [12]. In short, nine genes were used to establish a risk model and divided TCGA-LIHC into high-risk and low-risk groups according to the median risk score. en, the mRNA levels of the GSE54236 and GSE14520 cohorts were standardized as a validation cohort via the "sva" function in R. e risk score of each HCC patient was calculated according to the same formula. Overall survival (OS) of the two subgroups was compared based on the Kaplan-Meier analysis. Principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) analysis were performed by the "t-SNE" package in R. e receiver operating curve (ROC) and area under the curve (AUC) of 1-, 3-, and 5-year OS were analyzed using the "survivalROC" package in R.

Functional Enrichment Analysis.
Gene set enrichment analysis (GSEA) with h.all.v7.2.symbols gene sets was used to investigate potential biological functions according to risk scores. In addition, we calculated the difference in the immune cell enrichment fraction and function between highand low-risk groups by single-sample gene set enrichment analysis (ssGSEA) [13].

Cell Culture.
e cells and culture conditions were as follows: human normal liver cell line LO-2 (RPMI-1640 containing 10% FBS), and HCC cell lines Huh7 (MEM containing 10% FBS) and 97H (DMEM containing 10% FBS) in a 5% CO 2 incubator at 37°C. e total RNA was extracted with TRIZOL (Invitrogen) reagent. en, cDNA was obtained by reverse transcription with cDNA Synthesis Mix (E047-01B; Novoprotein) and analyzed by quantitative PCR (E096-01B; Novoprotein). e mRNA expression levels were normalized against GAPDH expression. All primers used for the nine genes in this study are listed in Supplementary Table S2. 2.7. Statistical Analysis. All statistical analyses were completed by R software 4.0.1, and p < 0.05 were considered to be statistically significant.

Identification of DEGs, Hub Genes, and Mutation Patterns Associated with Necroptosis.
e expression levels of a total of 67 genes with necroptosis were calculated, and we identified 10 downregulated genes (ID1, BACH2, AXL,  MYC, GATA3, TLR3, FLT3,  e correlation network is shown in Figure 1(b). A PPI network was used to identify the interactions between DEGs, which was displayed in Cytoscape, and the MCC algorithm in Cytohubba plug-in was used to calculate the top 10 hub genes ( Figure 1(c)). A missense mutation was the main single nucleotide mutation type (Figure 1(d)). ese genes also had different levels and types of copy number variation (Figure 1(e)) and showed different patterns of methylation ( Figure 1(f )).

Independent Prognostic Analysis and Establishment of a
Prognostic Nomogram Based on TCGA-LIHC. We further evaluated whether the risk model could be used as an     Figure 5(a)). In the multivariate analysis, the risk score was an independent prognostic factor (p < 0.001, HR � 4.010, 95% CI: 2.645-6.081; Figure 5(b)). For the GSE14520 cohort, risk score is also an independent prognostic factor (p � 0.003, HR � 2.209, 95% CI: 1.307-3.734; Figures 5(c) and 5(d)). Finally, we created a novel prognostic nomogram based on the training cohort that incorporates risk scores and clinical data to provide a credible method for predicting HCC patient survival ( Figure 5(e)).

Functional Enrichment Analysis.
We further attained DEGs (logFC >0.585,p < 0.05) to distinguish the biological functions and networks associated with the risk group in the training and validation cohorts. Total GSEA analysis results are presented in Supplementary Table 3. Figure 6(a) shows that the top five hallmarks, SPERMATOGENESIS, MITO-TIC_SPINDLE, G2M_CHECKPOINT, E2F_TARGETS, and MYC_TARGETS_V1, were associated with the highrisk subgroup, while three hallmarks, COAGULATION, BILE_ACID_METABOLISM, and XEN-OBIOTIC_METABOLISM, were more enriched in the lowrisk subgroup in the training cohort. e enrichment analysis results of the two validation cohorts were similar to the training cohort (Figures 6(b) and 6(c)).

Differences in the Tumor Immune Microenvironment
Associated with Risk Subgroups. Despite the lack of effective immunotherapeutic biomarkers for HCC, the immune score of the tumor immune microenvironment is a promising indicator. erefore, we further investigated 16 types of immune cells and 13 types of immune functions in different risk subgroups according to ssGSEA. We found that low-risk patients had more activated immune cells and functions in both the training cohort and the validation cohorts (Figure 7), which may be the reason for the better prognosis of low-risk patients.

Relationship Between Risk Model Genes and Drug
Sensitivity. By analyzing the GDSC and CTRP databases, potential drugs were found to be associated with genes involved in the risk model. In general, we found that high expression of TNFRSF21 and SQSTM1 mostly reduced drug sensitivity, while the high expression of TARDBP, PLK1, and FLT3 enhanced drug sensitivity (Figure 8).

Validation of the Expression of Risk Model Genes.
We compared the mRNA levels of the risk model genes in two HCC cell lines (Huh7 and 97H) and a normal liver cell line (LO-2).
e expression of all the other genes was consistent with the previous results except for the increased expression of FLT3 in tumor cell lines (Figures 9(a)-9(i)). It is reported that FLT3 promotes the proliferation and migration of HCC, so we further investigated the reasons for the decrease in FLT3 expression. e results showed that FLT3 copy number deletion mutation exists in 37.5% of patients in TCGA-LIHC which is related to the level of mRNA expression (Figures 9(j)-9(k)).

Discussion
HCC is the second most lethal tumor after lung cancer, and there were approximately 830,180 new deaths worldwide in 2020 [14]. In recent years, the traditional prognostic evaluation system based on clinicopathological parameters and staging has not been able to meet the requirements of precision medicine [15]. With the development of sequencing technology, researchers have paid more attention to the molecular typing of diseases and the seeking of new biomarkers to guide clinical diagnosis and treatment [16]. is strategy not only complements the traditional prognosis evaluation, but also reveals a new pathogenesis. Necroptosis is a unique way of cell death. Recently, its characteristics have been described in a variety of human tumors. Generally, according to different subtypes of necroptosis, the prognosis of patients can be accurately predicted [6,7]. In hepatocellular carcinoma, it has been recognized that necroptosis is a double-edged sword that provides the inflammatory environment required for carcinogenesis, while the immune response is launched to fight against  tumors [17]. However, the signature of necroptosis genes has not been fully described in HCC.
is study systematically identified DEGs related to necroptosis in patients with HCC. First, we found that 10 genes were downregulated and 38 genes were upregulated in tumor samples. en, we performed consistent clustering according to 23 PRNGs and found that patients with HCC could be divided into two subtypes between which the survival time differed substantially. Interestingly, the clusters were associated with clinicopathological parameters, which means that the strong inflammation caused by necroptosis compels tumor cells face severe natural selection, which leads to stronger invasiveness of some subclones and poor prognosis. Subsequently, nine genes were used to build a LASSO risk model according to the training cohort, which was well verified in the test cohorts. We found that risk score was an independent prognostic factor and mapped a nomogram to predict the OS of HCC patients. We further illuminated the functional enrichment characteristics of different risk subgroups. Mitotic spindle disruption [18], E2F [19], MYC pathways [20], mTORC1 signaling [21], and G2/M cell cycle [22] enriched in the high-risk group are all associated with the stronger invasiveness of HCC.
HCC is characterized by low tumor mutational burden and microsatellite stability, and the expression of immune checkpoints fails to predict the response of patients to    Journal of Oncology 9 1.00 TCGA * *** *** *** *** *** ** * * *  10 Journal of Oncology immunotherapy [23]. erefore, evaluation of the tumor immune microenvironment may be the critical index of immunotherapy in the future [24]. Recent studies revealed that immune checkpoint inhibitors cannot activate exhausted T cells, and the characteristic of "hot" tumor (innate attraction of Tcells) is the key to the effectiveness of immunotherapy [25,26]. Conceivably, the low-risk patients may benefit from immunotherapy due to more active immune cells and a stronger immune function. In addition, we also conducted drug sensitivity analysis and found that high expression of TNFRSF21 and SQSTM1 mostly reduced drug sensitivity, while high expression of TARDBP, PLK1, and FLT3 enhanced drug sensitivity.
SQSTM1, TNFRSF21, TARDBP, HAT1, PLK1, MYCN, KLF9, HSP90AA1, and FLT3 were applied in the construction of the LASSO risk model. SQSTM1 is reported as an autophagy receptor that can serve as a bridge between polyubiquitinated cargo and autophagosomes as well as mediate necroptosis by recruiting RIPK1 [27]. SQSTM1, TNFRSF21, and TARDBP have rarely been reported in HCC, but they are associated with poor outcomes (Figure 2). HAT1 is an acetyltransferase that can form a complex with RIP1/3 to reduce programmed cell death [28]. PLK1 can activate the NF-κB signaling pathway to promote HCC development; thereby, harnessing necroptosis through inhibiting PLK1 may be a promising treatment    Figure 8: Correlation between drug sensitivity and risk gene expression according to the GDSC and CTRP databases. strategy [29,30]. MYCN [31] and HSP90AA1 [32] can promote HCC by activating the MYC pathway.
FLT3 was the only gene in our validation experiment that was inconsistent with the results of the bioinformatic analysis. FLT3 belongs to the receptor tyrosine kinase family, which is more widely reported in hematological diseases [33]. Sorafenib is the first-line drug for HCC and is a multikinase inhibitor, targeting FLT3 among others. e previous results showed that the mRNA expression of FLT3 was decreased in 64% of patients with HCC because of the loss of gene copy number; however, patients with high FLT3 expression benefit from sorafenib, which improves prognosis [34]. erefore, the high expression of FLT3 promotes the proliferation and migration of HCC [35], which may be the reason for the high expression of FLT3 in vitro.
ere were some limitations that need to be clarified in this study. First, patients with HCC were not stratified according to different primary factors (such as nonalcoholic fatty liver disease-associated hepatocellular carcinoma and hepatitis virus-associated hepatocellular carcinoma), and there might be necroptosis signature heterogeneity in different population. Second, bioinformatics provides an initial   strategy for screening genes, but the function of these genes needs to be further explored via protein analysis and in vitro/ vivo experiments.

Conclusion
is study established a valuable risk model based on necroptosis genes, which can effectively predict the prognosis of patients with HCC. Our results provided some potential biomarkers and targets, and further research will assist in elucidating the role of necroptosis genes in the progression of HCC.

Data Availability
Codes and other data used to support the findings of the study are available from the corresponding author upon reasonable request.

Conflicts of Interest
ere are no potential conflicts of interest for the authors to disclose.

Authors' Contributions
Xiao-Lan Zhang designed the study. Xin Gao and Di Huang performed the bioinformatic analysis and drafted the manuscript. Di Huang and Shu-Guang Li completed the in vitro experiments. e revised manuscript was supervised by Qian-Jia Ming, Dong-Lei Sun, and Wen-Xin Wang. e final version was approved by all authors. Xin Gao and Di Huang contributed to this work equally and should be regarded as the co-first authors. Acknowledgments e National Natural Science Foundation of China (Grant no. 81870381) provided funding for this work. Specifically, the recipient was Shuguang Li, and their team was responsible for cell experiments. Figure S1. Relationship between hub genes and prognosis. (a) OS and (b) DFS outcomes using Kaplan-Meier curves. Figure S2. Validation of the LASSO risk model using the GSE54236