Abstract
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.
1. 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 factor-related apoptosis-inducing ligand stimulated death receptors—TNFR1, TNFR2, and FAS—are the most prevalent triggers [1]. The 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 prognosis-related 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. Therefore, 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.
2. Materials and Methods
2.1. Datasets
We obtained the transcriptomic expression data and clinicopathological information of 370 HCC patients from The Cancer Genome Atlas-Liver Hepatocellular Carcinoma (TCGA-LIHC) database (https://portal.gdc.cancer.gov/projects/TCGA-LIHC) in October 2021. The same data of 78 and 221 HCC patients were obtained from the two Gene Expression Omnibus cohorts (GEO, ID : GSE54236 and GSE14520, https://www.ncbi.nlm.nih.gov/geo/) respectively.
2.2. 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. The 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/).
2.3. Clustering Analysis of TCGA-LIHC According to DEGs
We further screened out prognostic DEGs by univariate Cox regression (). Based on these DEGs, we performed a cluster analysis (“ConsensusClusterPlus” R package) and explored clinicopathological differences between different clusters.
2.4. 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. The 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. Then, the mRNA levels of the GSE54236 and GSE14520 cohorts were standardized as a validation cohort via the “sva” function in R. The 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. The 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.
2.5. 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 high- and low-risk groups by single-sample gene set enrichment analysis (ssGSEA) [13].
2.6. Cell Culture
The 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% CO2 incubator at 37°C. The total RNA was extracted with TRIZOL (Invitrogen) reagent. Then, cDNA was obtained by reverse transcription with cDNA Synthesis Mix (E047-01B; Novoprotein) and analyzed by quantitative PCR (E096-01B; Novoprotein). The 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 were considered to be statistically significant.
3. Results
3.1. Identification of DEGs, Hub Genes, and Mutation Patterns Associated with Necroptosis
The 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, KLF9, TNFRSF1A, and FAS) and 38 upregulated genes (HDAC9, DDX58, CYLD, MAPK8, CFLAR, MPG, RIPK1, TARDBP, SIRT2, BRAF, BCL2L11, OTULIN, ATRX, DIABLO, STUB, SPATA2, MLKL, MAP3K7, HAT1, FADD, ITPK1, CASP8, HSP90AA1, RNF31, HSPA4, USP22, TNFRSF21, SLC39A7, TSC1, SQSTM1, TRIM11, TRAF2, DNMT1, LEF1, PLK1, MYCN, CDKN2A, and TERT) in tumor tissue compared with normal tissue. The visual results of the DEGs are shown in a heatmap in Figure 1(a). The 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)). These genes also had different levels and types of copy number variation (Figure 1(e)) and showed different patterns of methylation (Figure 1(f)).
(a)
(b)
(c)
(d)
(e)
(f)
3.2. Relationship Between Hub Genes and Prognosis
We used the GEPIA platform to calculate the relationship between hub genes and prognosis. For OS, FADD, MAP3K7, and TNFRSF1A were associated with worse outcomes () and CFLAR, SQSTM1, and RIPK1 showed the same trend. In addition, the high expression of TRAF2 was associated with shortened disease-free survival (), and CASP8, FADD, and SQSTM1 showed the same trend (Supplementary Figure 1). Other genes had no significant correlation with prognosis.
3.3. HCC classification according to prognostic-related necroptosis genes (PRNGs)
We identified 23 PRNGs using univariate Cox analysis. KLF9 and FLT3 were associated with better prognosis (hazard ratio <1), whereas FADD, TRIM11, CASP8, IPMK, TRAF2, USP22, MAP3K7, SQSTM1, DNMT1, BRAF, CDKN2A, HSPA4, HAT1, PLK1, MYCN, SLC39A7, SPATA2, IDH1, HSP90AA1, TARDBP, and TNFRSF21 were associated with worse prognosis (hazard ratio >1, Figure 2(a)). We performed a cluster analysis and divided TCGA-LIHC into two categories according to CDF values (Figures 2(b)–2(d)). The survival analysis results showed that cluster 1 had a better prognosis than cluster 2 (Figure 2(e)), and the clinicopathological parameters of the two clusters were significantly different, which suggested that different clinical features represent different necroptosis patterns (Figure 2(f)).
(a)
(b)
(c)
(d)
(e)
(f)
3.4. Establishment and Validation of a PRNG LASSO Risk Model in TCGA Training Cohort and GEO Test Cohorts
We performed LASSO regression analysis using 23 PRNGs in TCGA-LIHC training cohort to establish the prognostic model. To minimize overfitting, nine genes were used to generate the final TCGA LASSO risk model (Figures 3(a) and 3(b)), and the risk genes and their coefficients are shown in Table 1. The formula used to calculate risk scores was as follows: Sum of (gene expressioncoefficient). Based on the median risk score, TCGA-LIHC patients were divided into high- and low-risk groups (Figure 3(c)). There were more deaths and poorer OS in the high-risk group (Figures 3(d) and 3(e)). The subsequent ROC analysis revealed that the risk model could accurately assess and predict the survival of HCC patients (AUC at 1, 3, and 5 years was 0.789, 0.735, and 0.703, respectively; Figure 3(f)). Finally, the PCA plot and t-SNE plot revealed that the risk models could distinguish the high- and low-risk groups to some extent (Figures 3(g) and 3(h)). Heatmaps showed worse staging and pathological grade in the high-risk group (Figure 3(i)).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
The prognostic assessment of the risk model was well reproduced in the GEO validation cohort. We first used SE14520 to validate the risk model. The risk scores of 221 HCC patients were calculated according to the previous formula, and 96 patients were assigned to the low-risk group and 125 to the high-risk group (Figure 4(a)). Patients in the low-risk group had longer OS (Figures 4(b) and 4(c)). The AUC of 1-year, 3-year, and 5-year OS was 0.688, 0.637, and 0.643, respectively (Figure 4(c)), and the PCA plot and t-SNE plot revealed that the risk genes were effective in distinguishing between the two risk groups (Figure 4(d)). The GSE54236 cohort containing 78 HCC patients showed similar results, and the AUC of 1-year and 3-year OS was 0.665 and 0.631, respectively (Supplementary Figure 2).
(a)
(b)
(c)
(d)
3.5. 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 independent predictor of prognosis. For the training cohort, TNM staging (, HR = 3.084, 95% CI: 1.955–4.866) and risk score (, HR = 4.480, 95% CI: 2.975–6.746) were possible risk factors in the univariate Cox regression analysis (Figure 5(a)). In the multivariate analysis, the risk score was an independent prognostic factor (, HR = 4.010, 95% CI: 2.645–6.081; Figure 5(b)). For the GSE14520 cohort, risk score is also an independent prognostic factor (, 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)).
(a)
(b)
(c)
(d)
(e)
3.6. Functional Enrichment Analysis
We further attained DEGs (logFC >0.585,) 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, MITOTIC_SPINDLE, G2M_CHECKPOINT, E2F_TARGETS, and MYC_TARGETS_V1, were associated with the high-risk subgroup, while three hallmarks, COAGULATION, BILE_ACID_METABOLISM, and XENOBIOTIC_METABOLISM, were more enriched in the low-risk subgroup in the training cohort. The enrichment analysis results of the two validation cohorts were similar to the training cohort (Figures 6(b) and 6(c)).
(a)
(b)
(c)
3.7. 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. Therefore, 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.
(a)
(b)
(c)
3.8. 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).
3.9. 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). The 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. The 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)).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
In addition, protein expression levels of risk model genes (HAT1, SQSTM1, TARDBP, and HSP90AA1) were obtained from the CPTAC database (Figure 10(a)). Meanwhile, we verified the protein expression levels of HAT1, SQSTM1, PLK1, HSP90AA1, TARDBP, and TNFRSF21 using the HPA database (Figure 10(b)).
(a)
(b)
4. 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]. This 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.
This 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. Then, 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 immunotherapy [23]. Therefore, 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 T cells) 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 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 multi-kinase inhibitor, targeting FLT3 among others. The 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]. Therefore, 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.
There 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.
5. Conclusion
This 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
There 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. The revised manuscript was supervised by Qian-Jia Ming, Dong-Lei Sun, and Wen-Xin Wang. The 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
The 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.
Supplementary Materials
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 test cohort. (a) Distribution according to the risk scores. (b, c) High-risk patients have higher mortality rates. (d) ROC analysis. (e) PCA analysis. Supplementary Table 1. Necroptosis gene list. Supplementary Table 2. Gene primer. Supplementary Table 3. Total GSEA analysis results. (Supplementary Materials)