A necroptosis-related prognostic model for predicting prognosis, immune landscape, and drug sensitivity in hepatocellular carcinoma based on single-cell sequencing analysis and weighted co-expression network

Background: Hepatocellular carcinoma (HCC) is a highly lethal cancer and is the second leading cause of cancer-related deaths worldwide. Unlike apoptosis, necroptosis (NCPS) triggers an immune response by releasing damage-related molecular factors. However, the clinical prognostic features of necroptosis-associated genes in HCC are still not fully explored.

Methods: We analyzed the single-cell datasets GSE125449 and GSE151530 from the GEO database and performed weighted co-expression network analysis on the TCGA data to identify the necroptosis genes. A prognostic model was built using COX and Lasso regression. In addition, we performed an analysis of survival, immunity microenvironment, and mutation. Furthermore, the hub genes and pathways associated with HCC were localized within the single-cell atlas.

Results: Patients with HCC in the TCGA and ICGC cohorts were classified using a necroptosis-related model with significant differences in survival times between high- and low-NCPS groups (p < 0.05). High-NCPS patients expressed more immune checkpoint-related genes, suggesting immunotherapy and some chemotherapies might prove beneficial to them. In addition, a single-cell sequencing approach was conducted to investigate the expression of hub genes and associated signaling pathways in different cell types.

Conclusion: Through the analysis of single-cell and bulk multi-omics sequencing data, we constructed a prognostic model related to necroptosis and explored the relationship between high- and low-NCPS groups and immune cell infiltration in HCC. This provides a new reference for further understanding the role of necroptosis in HCC.

Introduction

Primary liver cancer is the sixth most common cancer in the world and the second leading cause of cancer-related death (Yang et al., 2019). Hepatocellular carcinoma (HCC) is the most common type of primary liver cancer (Chaudhary et al., 2019). Most HCC patients are diagnosed at an advanced stage. The gold standard treatments, including tumor resection, local ablation with radiofrequency, and sometimes liver transplantation, have low success rates with high relapse rates and short survival times (Dhanasekaran et al., 2016). Additionally, patients with HCC who present with similar tumor, lymph node, and metastasis (TNM) stage have different clinical outcomes, and there are few current effective prognostic indicators.

Recent research has demonstrated the importance of the tumor microenvironment (TME) in promoting tumor aggressiveness (Altorki et al., 2019). The survival of patients with various malignancies can be prolonged by immune checkpoint inhibitors. However, many patients with HCC currently often respond poorly to immune checkpoint inhibitors, which may be due to low mutational loads, acquiring new immune checkpoints, and producing immunosuppressive factors (Riley et al., 2019). Therefore, there is a need to identify new biomarkers for HCC as well as to comprehend their significance in TME.

Programmed cell death has a strong impact on the characterization of the TME ecosystem (Chevrier et al., 2017). Resistance to apoptosis, a problem affecting cancer development, is one of the hallmarks of cancer (Reyna et al., 2017). In the process of cancer cell resistance to death, growth signals are overactivated, the metabolism is reprogrammed, and a change in the immune microenvironment occurs (Sahin et al., 2017). Inducing cancer cell death is becoming increasingly popular as a potential cancer treatment method (Bersuker et al., 2019). HCC cells can die by several different mechanisms, including apoptosis and necroptosis (NCPS) (Yuan et al., 2019). Both mechanisms play a significant role in homeostasis, inflammation, anti-infection, and tumorigenesis (Karki et al., 2021; Koren and Fuchs, 2021).

Necroptosis was once believed to be the “accidental death” of cells. However, current research indicates that necroptosis is distinct from conventional apoptosis (González-Juarbe et al., 2017). Necroptosis leads to membrane destabilization, which subsequently precedes swelling and lysis of cells, resulting in the release of intracellular constituents (González-Juarbe et al., 2017). Inhibited caspase 8 and receptor-interacting serine/threonine protein kinase 1 (RIPK1) are both involved in necroptosis pathway activation via recruitment and activation of receptor-interacting serine/threonine protein kinase 3 (RIPK3) (Alvarez-Diaz et al., 2016). Necroptosis occurs when caspase 8 is inactivated or absent, resulting in the activation and autophosphorylation of RIPK1 and RIPK3(Tanzer et al., 2017). During this process, the cell membrane ruptures, and the contents are released, stimulating an immune response (Kalliolias and Ivashkiv, 2016). Necroptosis becomes attractive as an alternative to apoptosis for killing tumor cells if apoptosis fails to kill them (Kalliolias and Ivashkiv, 2016). As well, the immune microenvironment is positively impacted by necroptosis (Gong et al., 2019).

Interestingly, the role of necroptosis in cancer is complex. In general, high levels of necroptosis result in strong adaptive immune responses that inhibit the progression of tumors. The recruitment of strong immune responses may also contribute to tumor progression (Koo et al., 2015; Najafov et al., 2017). Moreover, the inflammatory response may contribute to tumorigenesis and metastasis, as well as generate an immunosuppressive tumor microenvironment. Guo et al. (2022) has shown that loss of key necroptosis gene significantly reduces clinical symptoms of liver injury and fibrosis. Necroptosis has completely opposite effects on different types of cancer, the mechanism of which is still unclear. With the emergence of immune checkpoint therapy, changes in the immune microenvironment resulting from necroptosis are also important to consider. There is therefore a need to investigate the relationship between necroptosis and HCC.

Here, we downloaded the data of HCC patients from the TCGA and ICGC databases, as well as two single-cell datasets, GSE125449 (Ma et al., 2019) and GSE151530 (Ma et al., 2021), and one microarray dataset, GSE76427 (Grinchuk et al., 2018) from the GEO database. The TCGA cohort was used for model building. The ICGC cohort and GSE76427 were used to validate the results of our analysis. Two single-cell sequencing datasets, GSE125449 and GSE151530, were chosen for single-cell analysis because of their relatively large sample size and inclusion of clinical data. Through comprehensive bioinformatic analysis, we developed a prognostic model based on necroptosis and classified HCC patients into low- and high-risk groups, the results of which were significantly different. Furthermore, we explored the potential value of the signature in guiding the tumor mutational load, immune microenvironment, and drug sensitivity.

Methods

Download and processing of transcriptome data

This flowchart illustrates the key steps in the analysis (Figure 1). The data of HCC were downloaded from TCGA (https://portal.gdc.cancer.gov/) as a training cohort (Grossman et al., 2016). Count data and TPM data of HCC were extracted using R software (4.2.0), and a total of 363 tumor samples with complete clinical data were obtained. The HCC dataset was downloaded through ICGC (https://dcc.icgc.org/) database as a validation cohort, and the count data type and TPM data type of HCC were extracted, and a total of 240 tumor samples were obtained with complete clinical information (Zhang et al., 2019). GSE76427, measured using the Illumina HumanHT-12 V4.0 expression beadchip, contained 115 HCC samples (Grinchuk et al., 2018). The raw CEL files for GSE76427 were downloaded from the GEO database. More details of the data processing are in Supplementary Material S1, S2.

FIGURE 1

FIGURE 1. Flowchart of the analysis.

Download and processing of single-cell data

The single-cell datasets GSE125449 and GSE151530 for HCC were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) (Barrett et al., 2013). The GSE125449 dataset contains nine HCC samples and the GSE151530 dataset contains 32 HCC samples. We performed quality control on the data of all samples. We retained cells with genes expressed in at least 10 cells, less than 10% of mitochondrial genes, more than 200 genes, less than 5% hemoglobin genes, less than 50% ribosomal genes, and expression between 200 and 7000. We set a limit of 3000 highly variable genes. Next, we normalized all samples, removed batch effects, and integrated them by SCT. Then, using the tSNE method with the “DIMS” parameter set to 20, the dimensionality of the data was reduced. Cell clustering was then carried out using the “KNN” method with a resolution of 2.0. Subsequently, the cells were annotated with the Human Primary Cell Atlas (HPCA) from the “SingleR” package as a reference dataset (Mabbott et al., 2013). Finally, the proportion of NCPS-related genes in each cell can be calculated using the “PercateFeatureSet” function.

Identification of NCPS-related genes

In the GeneCards database (https://www.genecards.org/), 614 genes associated with necroptosis were identified (Safran et al., 2021). A total of 92 genes were identified that had an association score of greater than 1.0 with necroptosis (Supplementary Material S3). Then, the NCPS-related genes were scored for each sample by the combined analysis of ssGSEA (Single Sample Gene Set Enrichment Analysis) and WGCNA (Weighted Co-Expression Network Analysis). The log2 processed data were used for ssGSEA analysis.

ssGSEA

Gene sets enriched in a sample are often quantified by using the ssGSEA method with “GSVA” package (version: 1.44.2) (Bindea et al., 2013; Hänzelmann et al., 2013). In this study, ssGSEA analysis was utilized to determine the NCPS-related scores of each patient with HCC.

WGCNA

WGCNA analysis is one method used in systems biology for determining patterns of genetic association among diverse samples (Langfelder and Horvath, 2008). In addition to identifying highly covariant genomes, WGCNA analysis can be used to identify potential biomarkers or therapeutic targets based on the correlation between genomes and phenotypes. In this study, gene modules associated with NCPS scores in HCC were found by “WGCNA” package (version: 1.71), and genes associated with necroptosis were obtained. Non-gray modules were identified by setting a soft threshold of eight, a minimum number of module genes of 80, and combining modules that had similarities of less than 0.3.

Construction of NCPS-related prognostic model

First, univariate COX analysis was used to identify NCPS-related genes with prognostic values by using the “survival” package (version: 3.3-1). Next, a prognostic model was developed based on the least absolute shrinkage and selection operator (LASSO) regression for NCPS-related genes by using the “glmnet” package (version: 4.1-4) (Lossos et al., 2004; Friedman et al., 2010). In this way, the NCPS score could be calculated for each HCC sample by the formula. Gene expression levels were weighted by their respective coefficients of LASSO regression to calculate the NCPS score. The formula was as follows:

$NCPSscore={\displaystyle \sum }_{i=\mathit{1}}^{n}C{oef}_i\times Ex{p}_i(1)$

where n, Exp_i, Coef_i, represented the number, the expression value, and the coefficient of each selected gene, respectively. According to the median value of the TCGA-HCC cohort, patients could be classified into low- and high-risk groups. Thereafter, we assessed the accuracy of the model by comparing prognostic differences between the two groups.

Validation of NCPS-related prognostic model

The ICGC cohort and GSE76427 were selected as the external validation cohorts. According to the formula of the prognostic model, NCPS scores for each sample were calculated, and patients were categorized based on their median NCPS scores into high-risk and low-risk groups. We then conducted a survival analysis comparing the high- and low-NCPS groups. Receiver operating characteristic (ROC) curves were utilized to evaluate the model’s accuracy by using the “timeROC” package (version: 0.4) (Li et al., 2018). To determine whether the model grouped HCC patients more effectively, principal component analysis (PCA) was performed using the “PCAtools” package (version: 2.8.0) and “scatterplot3d” package (version: 0.3-41).

Immune infiltration and mutation landscape

We performed immune infiltration analysis of HCC patients in the TCGA database using immune cell infiltration algorithms from the IOBR package (version: 0.99.9) (Zeng et al., 2021). Next, we examined the differences in the levels of immune cell infiltration between the two NCPS groups and presented the immune cells with different levels of infiltration as a heat map. Also, the expression of immune checkpoint-related genes in the various NCPS subgroups was visualized by a boxplot. We identified the top 20 genes with the highest mutation rates by comparing the mutation rates between groups with high and low NCPS scores.

Nomogram

Using clinical data and NCPS values, a nomogram was developed in this study to assess the probability of mortality in patients with HCC using the “rms” package (version: 6.3-0) and “regplot” package (version: 1.1). This nomogram was evaluated by using prognostic ROC curves and decision curve analysis (DCA) to determine its accuracy in predicting patient outcomes. The DCA analysis was performed using the “ggDCA” package (version: 1.1) (Vickers and Elkin, 2006).

Drug sensitivity, immunohistochemistry, pathways

To improve personalized treatment, we calculated half maximal inhibitory concentrations (IC50) using the “pRRophetic” package (version: 0.5) and compared these data between high-risk and low-risk groups (Geeleher et al., 2014). Low IC₅₀ values indicate greater drug effectiveness. The Human Protein Atlas (HPA) database (version: 21.1, http://www.proteinatlas.org/) is the most comprehensive database for assessing protein distribution in human tissues (Uhlén et al., 2015).

HPA database was used to obtain prognostic gene expression data. Immunohistochemical staining images of normal and HCC tissues were used to analyze the protein expression of genes. In addition, we performed an enrichment analysis of pathways associated with different cell types in the single-cell data and then mapped the significantly different pathways to tSNE plots for visualization. Pathway enrichment analysis was performed using the “irGSEA” package (version: 1.1.2). Finally, the pathways associated with HCC in the TCGA dataset were analyzed.

Statistical analysis

Statistical analysis was performed using the R software (version 4.2.0). Continuous data were analyzed using Mann-Whitney tests, and categorical data were analyzed using Fisher’s exact tests. Pearson correlation coefficient was used to estimate the correlation between continuous variables. The Kaplan-Meier method was used for survival analysis. The Log-rank test was used to determine the significance of differences. All statistical analyses were considered significant if the p-values were less than 0.05.

Results

Annotation of single-cell sequencing data and identification of differentially expressed genes associated with NCPS

We first analyzed the single-cell sequencing datasets GSE125449 and GSE151530 for HCC to integrate different samples. As shown in Figure 2A, there were no significant batch effects in the 23 samples, and further analysis of the results could be conducted. The K-Nearest Neighbor (KNN) algorithm was used to divide all cells into 43 clusters (Figure 2B). After entering 92 genes related to NCPS using the “PercateFeatureSet” function, a percentage of the genes associated with NCPS was calculated for each cell. Cells were classified based on the median percentage of NCPS genes and represented as tSNE plots (Figure 2C). We then identified eight distinct cell types based on the expression of surface markers of different cell types in different clusters. They were T cells, macrophages, endothelial cells, NK cells, monocytes, smooth muscle cells, B cells, and tumor cells (Figure 2D). The surface markers for eight types of cells are shown in Figure 2E. Furthermore, we identified 3598 genes that were differentially expressed between high- and low-NCPS groups (Supplementary Material S4). Using the WGCNA analysis of 363 samples from the TCGA cohort, we have obtained gene modules associated with necroptosis. In total, eight non-gray modules were identified by setting a soft threshold of 8 (Figure 3A). As shown in Figure 3B, MEsalmon, MEtan, and MEpurple were strongly associated with the NCPS score. Further analysis was performed on the genes in these three modules.

FIGURE 2

FIGURE 2. Single-cell analysis. (A) There were no significant batch effects in the 23 samples. (B) Dimensionality reduction and cluster analysis. (C) Percentage of genes in each cell that are involved in necroptosis. (D) Annotation of cells in accordance with their surface marker genes. (E) A list of the surface markers of the eight types of cells.

FIGURE 3

FIGURE 3. Construction of prognostic model. (A,B) WGCNA screening for modules relating to necroptosis. (C) The intersection between differential genes identified by single-cell analysis and genes identified by WGCNA. (D,E) Using Lasso regression, the final genes were selected for the prognostic model.

The NCPS-related prognostic model could be used to classify HCC patients and predict their prognosis

An intersection was drawn between differential genes derived from single-cell analysis and genes identified by WGCNA. In Figure 3C, 74 genes are shown as candidates for the next step in the analysis (Supplementary Material S5). Based on univariate COX analysis within the TCGA cohort, 45 genes have been identified as significantly associated with prognosis. The LASSO regression analysis employed a random seed of 2022, and the results indicated that gene contraction stabilized with minimal partial likelihood deviation when the number of genes included was 8 (Figures 3D,E). Table 1 summarizes the results of the Lasso regression for each of these genes. The prognostic model was constructed from eight genes, including RAD21, NBN, PRKDC, MAP2K2, RIPK2, BOP1, POLR2E, and GPX4. As follows was the prognostic model.

$\text{NCPS}=\text{RAD}21*0.07523793+\text{NBN}*0.02974632+\text{PRKDC}*0.2080255+\text{MAP}2\mathrm{K}2*0.13063826+\text{RIPK}2*0.05091183+\text{BOP}1*0.0940779+\text{POLR}2\mathrm{E}*0.18787475+\text{GPX}4*0.1444125$

TABLE 1

TABLE 1. Eight genes were identified by lasso regression to construct a prognostic model.

Based on median values, patients were divided into high- and low-risk groups. Figure 4A showed that the high-NCPS group in the training cohort had a worse prognosis (p = 0.015). Figure 4B demonstrated that patients with high-NCPS had worse outcomes than those with low-NCPS in the validation cohort (p = 0.0078). ROC curves were generated for both the training and validation cohorts to test the prognosis assessment ability. As shown in Figure 4C, the area under the curve (AUC) values were 0.722, 0.746, 0.741, and 0.763 at 1, 2, 3, and 5 years in the TCGA cohort, respectively. In the validation cohort, AUC values were 0.730, 0.653, 0.625 and 0.623 at 1, 2, 3 and 5 years, respectively (Figure 4D). In the GSE76427 cohort, the results showed that the high-NCPS group had a worse prognosis (p = 0.0037), and AUC values were 0.682, 0.696, and 0.778 at 2, 3 and 5 years, respectively (Supplementary Material S6).

FIGURE 4

FIGURE 4. Validation of prognostic model. (A) Survival analysis of the training set showed significantly poorer outcomes for NCPS high group (p = 0.015). (B) Survival analysis results in the validation set were similar to those in the training set (p = 0.0078). (C) ROC curve of the training set. (D) ROC curve of the validation set. (E,F) 3D-PCA analysis in the training set and validation set.

Based on these results, the NCPS-related prognostic model was found to be accurate in predicting the outcomes of patients in all three cohorts. Furthermore, PCA was performed on the eight genes included in all three cohorts, and the results were similar. The results showed that the model performed well in classifying HCC patients (Figures 4E,F).

The nomogram could be more reliable in predicting patient outcomes than other indicators

Combining clinical information and NCPS scores, we constructed a nomogram that allows us to assess patients’ prognoses. In Figure 5A, the estimated mortality rates for patients with the high-NCPS score “TCGA-G3-A7M9” were 0.626, 0.92, and 0.984 at 1, 3, and 5 years based on gender, age, T-stage, and total stage (Table 2). Based on the low-NCPS score, the estimated mortality rates for patients with “TCGA-DD-AADS” were 0.0389, 0.102, and 0.153 at 1, 3, and 5 years based on sex, age, T-stage, and total stage (Figure 5B). In Supplementary Material S7, NCPS scores and clinical characteristics of 363 patients from the TCGA-HCC dataset are presented. Accordingly, a clinical decision could be based on assessing a patient’s risk and guiding their subsequent treatment. Furthermore, the accuracy of the nomogram was assessed through ROC analysis, which showed AUCs of 0.75, 0.67, and 0.68 for 1, 3, and 5 years, respectively (Figure 5C). In addition, we assessed the utility of the model to support clinical decision-making by using decision curve analysis (DCA) and reported the net clinical benefit of the model. The results showed that the nomogram is better than other clinical indicators, indicating that the nomogram is effective in predicting the patient’s prognosis (Figure 5D).

FIGURE 5

FIGURE 5. Construction of the nomogram. (A) High-NCPS patient “TCGA-G3-A7M9”: Mortality rates were estimated to be 0.626, 0.92, and 0.984 at 1, 3, and 5 years, respectively. (B) Low-NCPS patient “TCGA-DD-AADS”: Mortality rates in 1, 3, and 5 years were estimated to be 0.0389, 0.102, and 0.153, respectively. (C) The ROC curve for the nomogram. (D) DCA analysis showed that the nomogram was more effective than other clinical indicators.

TABLE 2

TABLE 2. Comparison of clinical data from patients with high- and low-NCPS in the TCGA-HCC dataset.

Survival analysis and cellular localization of the eight hub genes

Survival analysis was performed for each of the eight hub genes. Compared with patients with low expression, those with high levels of RAD219 (p = 0.0078), RIPK2 (p = 0.005), BOP1 (p = 0.0038), POLR2E (p = 0.02), and MAP2K2 (p = 0.017) had significantly poorer outcomes (Figure 6A). To investigate the expression of the eight hub genes in various cell types, we conducted a single-cell sequencing analysis. As shown in Figures 6B–J, RAD21, BOP1, POLR2E, and PRKDC were mainly expressed in tumor cells, RIPK2 was mainly expressed in monocytes, MAP2K2 was mainly expressed in tumor cells and macrophages, NBN was mainly expressed in macrophages and monocytes, and GPX4 was mainly expressed in tumor cells and T cells.

FIGURE 6

FIGURE 6. Survival analysis and cellular localization of the eight hub genes. (A) The survival analysis of eight hub genes in the TCGA cohort. (B–J) The expression of eight hub genes in different types of cells.

The NCPS scores are positively correlated with the levels of immune cell infiltration and the expression of immune checkpoint genes

As shown in the above analysis, patient outcomes varied significantly within the NCPS subgroups. To explore the reasons for this and inform immunotherapy, comparisons of the levels of immune cell infiltration between the various groups were conducted.

As shown in Figure 7A, six different immune infiltration algorithms have been used to estimate the relationship between necroptosis and immune cells. Specifically, the three algorithms of MCP counter, Quanti-seq, and TIMER clearly demonstrated that there were more immune cell infiltrations in the high-NCPS group, including macrophages, NK cells, T cells, monocytes, B cells, and dendritic cells. We then investigated the expression of genes associated with immune checkpoints. Figure 7B demonstrated that many immune checkpoint genes, such as PDCD1 and CTLA4, were more highly expressed in the high NCPS group. High NCPS patients were likely to have a higher degree of immune infiltration. However, patients with high-NCPS may suffer from low response states due to high levels of immune checkpoint genes, and immune checkpoint inhibitors may be of greater benefit to patients with such conditions. In addition, we examined the immune infiltration results obtained by different algorithms. The QuantTIseq algorithm showed that patients with high-NCPS levels had more macrophages, B cells, and T cells (Figure 7C).

FIGURE 7

FIGURE 7. Immune infiltration analysis of TCGA cohort. (A) Heat map of immune cell infiltration in high and low NCPS groups with six immune infiltration algorithms. (B) Expression of immune checkpoint genes in high- and low-NCPS groups. (C) Results of the quanTIseq algorithm.

A high-NCPS score is associated with a greater incidence of gene mutations

According to the NCPS scores in the high- and low-group, 20 of the top mutated genes were identified. As shown in Figures 8A,B, the incidence of mutations in the 20 most frequently mutated genes was 89.53% (High NCPS) and 82.76 % (Low NCPS) for the two groups. In the high-NCPS group, the highest mutation rates were PT53 (40%), CTNNB1 (30%), and TTN (29%). In the low NCPS group, TTN (25%), CTNNB1 (24%), and PT53 (21%) were the mutations with the highest rates. A higher incidence of mutations was observed in the high-NCPS group as compared to the low-NCPS group. Mutations were analyzed for eight hub genes (Supplementary Material S8). The highest Variant Classification shown in Figure 8C was Missense Mutation. Single nucleotide polymorphism (SNP) was the highest Variant Type (Figure 8D). Figure 8E indicated that an average of 100 genes were mutated in each sample. Figure 8F showed that the top three base mutation types of single nucleotide variants (SNVs) were C>T, C>A, and T>C. In addition, we analyzed the correlation between pairs of mutated genes. Figure 8G showed a strong co-relation between FLG and OBSCN (p < 0.0001, OR = 8.803), FAT3 and DNAH7 (p = 0.00064, OR = 6.925). There was a strong mutually exclusive relationship between CTNNB1 and TP53 (p = 0.00811, OR = 0.459), AXIN1 and CTNNB1 (p = 0.00733, OR = 0.109) (Supplementary Material S9).

FIGURE 8

FIGURE 8. Mutation landscape of TCGA cohort. (A,B) Mutated genes in high- and low-NCPS groups. (C,D) Classification and types of variants. (E) Mutations in each sample. (F) The base mutation types of SNVs. (G) The correlation between pairs of mutated genes.

Drug sensitivity of HCC and hub gene protein expression are positively correlated with NCPS scores

Based on the “pRRophetic” package, we assessed the sensitivity of different NCPS subgroups to drugs commonly used as a treatment for HCC. The high-risk group showed higher sensitivity to cisplatin, docetaxel, paclitaxel, sunitinib, tipifarnib, bexarotene, bicalutamide, bortezomib, and bleomycin, while the low-risk group showed higher sensitivity to metformin, camptothecin, temsirolimus (Figure 9A). The immunohistochemical analysis of the HPA database showed that protein products with high NCPS-related genes were expressed at higher levels in HCC samples compared to normal tissues (Figure 9B).

FIGURE 9

FIGURE 9. Drug sensitivity and Immunohistochemical analysis. (A) Drug sensitivity analysis in. high- and low-NCPS groups. (B) Immunohistochemical analysis of eight hub genes.

Pathway enrichment and localization in single-cell sequencing data

Pathway enrichment analysis of single-cell data revealed that HALLMARK OXIDATIVE PHOSPHORYLATION was upregulated in Malignant cells but downregulated in T cells, TECs, and B cells. HALLMARK ALLOGRAFT REJECTION was downregulated in Malignant cells, upregulated in T cells, downregulated in CAFs, and upregulated in TAMS. HALLMARK-TNFA -SIGNALING-VIA-NFKB was downregulated in Malignant cells and upregulated in TAMs. HALLMARK-TGF-BETA-SIGNALING was upregulated in TECs (Figure 10). In addition, we explored the expression of these signaling pathways in different cell types by single-cell sequencing analysis (Figures 11A–D) and profiled the pathways associated with disease (Figure 11E).

FIGURE 10

FIGURE 10. Pathway enrichment analysis of single-cell data.

FIGURE 11

FIGURE 11. Pathway enrichment analysis. (A–D) Localization of different pathways in the single-cell dataset. (E) The number of pathways enriched in the TCGA cohort.

Discussion

With increasing incidence, HCC has become the second leading cause of cancer-related deaths (Bray et al., 2018). Due to lifestyle changes, HCC has become the fastest growing cancer in developed countries, but the response to antitumor therapy is relatively poor. Approximately 50% of HCC patients receive systemic therapy, traditionally with first-line sorafenib or lenvatinib. In the past 5 years, immune checkpoint inhibitors have completely altered the treatment regimen for HCC and improved the prognosis (Llovet et al., 2022). The immune microenvironment plays a significant role in the progression of HCC, and HCC with high- and low-necroptosis respond differently to immune checkpoint inhibitor therapy. However, at present, there are no validated biomarkers to aid in clinical decision-making in this regard.

Immune checkpoint inhibitors are used because immune cells can receive inhibitory signals by activating immune checkpoint molecules. By activating immune checkpoint molecules to receive inhibitory signals, their activity and proliferation are blocked (Huang and Chang, 2019). These immune checkpoints can be used by cancer cells, leading to impaired immune surveillance (Liu and Qin, 2019). PD-1, PD-L1, and cytotoxic T cell antigen 4 (CTLA-4) are the main immune checkpoints that have been targeted by monoclonal antibodies.

Utilizing comprehensive data analysis on HCC datasets from TCGA, ICGC, and GEO databases, we built a prognostic profile for NCPS-related genes associated with HCC. We calculated risk scores to identify high- and low-risk groups of patients with HCC. All three cohorts both showed that the high-risk group did significantly worse than the low-risk group in HCC. Xie et al. (2022) found similar results in triple-negative breast cancer, indicating that the higher the NCPS score, the larger the tumor and the worse the prognosis. Furthermore, the ROC curve revealed that this feature might be accurate in predicting the prognosis of patients with HCC at 1, 3, and 5 years. Based on the immune microenvironment analysis, immunotherapy was more likely to be effective in necroptosis with higher expression levels. The low response to immunotherapy of HCC could be attributed in part to the low mutational load and the generation of new immune checkpoints (Ricciuti et al., 2019; Scheiner et al., 2022). Therefore, it becomes fascinating to explore the immune microenvironment of HCC. Necroptosis may play an important role in TME by the release of inflammatory molecules during the induction of apoptosis. However, it remains unclear whether necroptosis plays a role in HCC.

Necroptosis is a necrotic programmed cell death that is powerfully immunogenic and participates in a complex interplay of autophagy and apoptosis (Gong et al., 2017). There is growing evidence that necroptosis plays an important role in prognosis, disease progression and tumor metastasis, and immune surveillance in cancer patients (Gong et al., 2019). Targeting necroptosis through immune checkpoint is also emerging as a new approach in tumor therapy.

The role of necroptosis in cancer is complex. It is still unclear exactly what role necroptosis plays in cancer. In general, high expression of necroptosis elicits strong adaptive immune responses that can inhibit tumor progression (Yatim et al., 2015). However, these recruited strong immune responses may also promote tumor progression. The inflammatory response may promote tumorigenesis and metastasis, as well as may generate an immunosuppressive tumor microenvironment (Seifert and Miller, 2017). Therefore, it is essential to investigate the molecular mechanisms and physiopathological aspects of necroptosis, as well as its interaction with immunity. In addition, it is imperative to discover the correlation between specific necroptosis markers and the prognosis of HCC. This is to unravel the confusion of necroptosis correlation in HCC and further develop targeted antitumor therapeutic drugs. In this study, combining single-cell analysis and second-generation sequence analysis, we were able to identify a significant difference between NCPS groups in terms of immune cell infiltration in HCC. Significant differences were observed between the high- and low-NCPS groups. In addition, the study findings indicated that a high level of NCPS group corresponds to a high level of immune checkpoint gene expression. Therefore, patients with HCC who have a high NCPS are more likely to respond to immunotherapy.

The datasets GSE125449 and GSE151530 have been initially explored to reveal changes in the immune microenvironment of HCC. Among the published results, GSE125449 reveals different degrees of heterogeneity of malignant cells within and between tumors and different TME landscapes by single-cell sequencing techniques. GSE151530 provides insights into the collective behavior of HCC cell communities by single-cell sequencing and potential tumor evolution in response to therapy drivers. We first classified HCC cells into two groups based on their NCPS scores by analyzing single cells of GSE125449 and GSE151530. This provided a reference for us to study the heterogeneity of necroptosis in HCC. Based on these two cell populations, we calculated the differentially expressed genes, which then served as a basis for constructing a prognostic model. For the validation of the prognostic model, survival data from the ICGC dataset was analyzed.

Our study has some limitations. First, a comprehensive analysis of HCC tissues is needed to fully validate how the eight NCPS-related genes are involved in the development of HCC. This was not examined in the current study. Second, further validation with larger patient datasets is needed better to estimate the accuracy of the model’s predictions. Finally, further experimental evidence is needed to fully understand the role and mechanisms of eight NCPS-related genes in HCC.

Conclusion

Through the analysis of single-cell and bulk multi-omics sequencing data, we constructed a prognostic model related to necroptosis and explored the relationship between high- and low-necroptosis groups and immune cell infiltration in HCC. This provides a new reference for further understanding the role of necroptosis in HCC. This may be useful in developing new therapeutic targets for the treatment of HCC. However, further molecular experiments are required to confirm the present findings.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

TY: Contributed to conception and design of the study; JL: Bioinformatics analysis and data mining; ZW and SW: Drafted the manuscript; CL: Statistical analysis; MS and JX: Literature review; YW and DF: Guidance for Bioinformatics analysis and data mining; YH and XZ: Review the manuscript; WM: Edited the manuscript. All the authors read and approved the final manuscript.

Funding

This work was supported by the Jincheng People’s Hospital (Nos. JSY-2021D007, JSY-2021D009) and the Traditional Chinese Medicine Bureau of Guangdong Province (20203002).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.984297/full#supplementary-material

References

Altorki, N. K., Markowitz, G. J., Gao, D., Port, J. L., Saxena, A., Stiles, B., et al. (2019). The lung microenvironment: An important regulator of tumour growth and metastasis. Nat. Rev. Cancer 19 (1), 9–31. doi:10.1038/s41568-018-0081-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Alvarez-Diaz, S., Dillon, C. P., Lalaoui, N., Tanzer, M. C., Rodriguez, D. A., Lin, A., et al. (2016). The pseudokinase MLKL and the kinase RIPK3 have distinct roles in autoimmune disease caused by loss of death-receptor-induced apoptosis. Immunity 45 (3), 513–526. doi:10.1016/j.immuni.2016.07.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Barrett, T., Wilhite, S. E., Ledoux, P., Evangelista, C., Kim, I. F., Tomashevsky, M., et al. (2013). NCBI GEO: Archive for functional genomics data sets-update. Nucleic Acids Res. 41, D991–D995. doi:10.1093/nar/gks1193

PubMed Abstract | CrossRef Full Text | Google Scholar

Bersuker, K., Hendricks, J. M., Li, Z., Magtanong, L., Ford, B., Tang, P. H., et al. (2019). The CoQ oxidoreductase FSP1 acts parallel to GPX4 to inhibit ferroptosis. Nature 575 (7784), 688–692. doi:10.1038/s41586-019-1705-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Bindea, G., Mlecnik, B., Tosolini, M., Kirilovsky, A., Waldner, M., Obenauf, A. C., et al. (2013). Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity 39 (4), 782–795. doi:10.1016/j.immuni.2013.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Bray, F., Ferlay, J., Soerjomataram, I., Siegel, R. L., Torre, L. A., and Jemal, A. (2018). Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. Ca. Cancer J. Clin. 68 (6), 394–424. doi:10.3322/caac.21492

PubMed Abstract | CrossRef Full Text | Google Scholar

Chaudhary, K., Poirion, O. B., Lu, L., Huang, S., Ching, T., and Garmire, L. X. (2019). Multimodal meta-analysis of 1, 494 hepatocellular carcinoma samples reveals significant impact of consensus driver genes on phenotypes. Clin. Cancer Res. 25 (2), 463–472. doi:10.1158/1078-0432.CCR-18-0088

PubMed Abstract | CrossRef Full Text | Google Scholar

Chevrier, S., Levine, J. H., Zanotelli, V. R. T., Silina, K., Schulz, D., Bacac, M., et al. (2017). An immune atlas of clear cell renal cell carcinoma. Cell 169 (4), 736–749. doi:10.1016/j.cell.2017.04.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Dhanasekaran, R., Venkatesh, S. K., Torbenson, M. S., and Roberts, L. R. (2016). Clinical implications of basic research in hepatocellular carcinoma. J. Hepatol. 64 (3), 736–745. doi:10.1016/j.jhep.2015.09.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Friedman, J., Hastie, T., and Tibshirani, R. (2010). Regularization paths for generalized linear models via coordinate descent. J. Stat. Softw. 33 (1), 1–22. doi:10.18637/jss.v033.i01

PubMed Abstract | CrossRef Full Text | Google Scholar

Geeleher, P., Cox, N., and Huang, R. S. (2014). pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PloS One 9 (9), e107468. doi:10.1371/journal.pone.0107468

PubMed Abstract | CrossRef Full Text | Google Scholar

Gong, Y.-N., Guy, C., Olauson, H., Becker, J. U., Yang, M., Fitzgerald, P., et al. (2017). ESCRT-III acts downstream of MLKL to regulate necroptotic cell death and its consequences. Cell 169 (2), 286–300. doi:10.1016/j.cell.2017.03.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Gong, Y., Fan, Z., Luo, G., Yang, C., Huang, Q., Fan, K., et al. (2019). The role of necroptosis in cancer biology and therapy. Mol. Cancer 18 (1), 100. doi:10.1186/s12943-019-1029-8

PubMed Abstract | CrossRef Full Text | Google Scholar

González-Juarbe, N., Bradley, K. M., Shenoy, A. T., Gilley, R. P., Reyes, L. F., Hinojosa, C. A., et al. (2017). Pore-forming toxin-mediated ion dysregulation leads to death receptor-independent necroptosis of lung epithelial cells during bacterial pneumonia. Cell Death Differ. 24 (5), 917–928. doi:10.1038/cdd.2017.49

PubMed Abstract | CrossRef Full Text | Google Scholar

Grinchuk, O. V., Yenamandra, S. P., Iyer, R., Singh, M., Lee, H. K., Lim, K. H., et al. (2018). Tumor-adjacent tissue co-expression profile analysis reveals pro-oncogenic ribosomal gene signature for prognosis of resectable hepatocellular carcinoma. Mol. Oncol. 12 (1), 89–113. doi:10.1002/1878-0261.12153

PubMed Abstract | CrossRef Full Text | Google Scholar

Grossman, R. L., Heath, A. P., Ferretti, V., Varmus, H. E., Lowy, D. R., Kibbe, W. A., et al. (2016). Toward a shared vision for cancer genomic data. N. Engl. J. Med. 375 (12), 1109–1112. doi:10.1056/NEJMp1607591

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, R., Jia, X., Ding, Z., Wang, G., Jiang, M., Li, B., et al. (2022). Loss of MLKL ameliorates liver fibrosis by inhibiting hepatocyte necroptosis and hepatic stellate cell activation. Theranostics 12 (11), 5220–5236. doi:10.7150/thno.71400

PubMed Abstract | CrossRef Full Text | Google Scholar

Hänzelmann, S., Castelo, R., and Guinney, J. (2013). Gsva: Gene set variation analysis for microarray and RNA-seq data. BMC Bioinforma. 14, 7. doi:10.1186/1471-2105-14-7

CrossRef Full Text | Google Scholar

Huang, P.-W., and Chang, J. W.-C. (2019). Immune checkpoint inhibitors win the 2018 Nobel Prize. Biomed. J. 42 (5), 299–306. doi:10.1016/j.bj.2019.09.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Kalliolias, G. D., and Ivashkiv, L. B. (2016). TNF biology, pathogenic mechanisms and emerging therapeutic strategies. Nat. Rev. Rheumatol. 12 (1), 49–62. doi:10.1038/nrrheum.2015.169

PubMed Abstract | CrossRef Full Text | Google Scholar

Karki, R., Sundaram, B., Sharma, B. R., Lee, S., Malireddi, R. K. S., Nguyen, L. N., et al. (2021). ADAR1 restricts ZBP1-mediated immune response and PANoptosis to promote tumorigenesis. Cell Rep. 37 (3), 109858. doi:10.1016/j.celrep.2021.109858

PubMed Abstract | CrossRef Full Text | Google Scholar

Koo, G.-B., Morgan, M. J., Lee, D.-G., Kim, W.-J., Yoon, J.-H., Koo, J. S., et al. (2015). Methylation-dependent loss of RIP3 expression in cancer represses programmed necrosis in response to chemotherapeutics. Cell Res. 25 (6), 707–725. doi:10.1038/cr.2015.56

PubMed Abstract | CrossRef Full Text | Google Scholar

Koren, E., and Fuchs, Y. (2021). Modes of regulated cell death in cancer. Cancer Discov. 11 (2), 245–265. doi:10.1158/2159-8290.CD-20-0789

PubMed Abstract | CrossRef Full Text | Google Scholar

Langfelder, P., and Horvath, S. (2008). Wgcna: an R package for weighted correlation network analysis. BMC Bioinforma. 9, 559. doi:10.1186/1471-2105-9-559

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, L., Greene, T., and Hu, B. (2018). A simple method to estimate the time-dependent receiver operating characteristic curve and the area under the curve with right censored data. Stat. Methods Med. Res. 27 (8), 2264–2278. doi:10.1177/0962280216680239

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, X., and Qin, S. (2019). Immune checkpoint inhibitors in hepatocellular carcinoma: Opportunities and challenges. Oncologist 24 (1), S3–S10. doi:10.1634/theoncologist.2019-IO-S1-s01

PubMed Abstract | CrossRef Full Text | Google Scholar

Llovet, J. M., Castet, F., Heikenwalder, M., Maini, M. K., Mazzaferro, V., Pinato, D. J., et al. (2022). Immunotherapies for hepatocellular carcinoma. Nat. Rev. Clin. Oncol. 19 (3), 151–172. doi:10.1038/s41571-021-00573-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Lossos, I. S., Czerwinski, D. K., Alizadeh, A. A., Wechser, M. A., Tibshirani, R., Botstein, D., et al. (2004). Prediction of survival in diffuse large-B-cell lymphoma based on the expression of six genes. N. Engl. J. Med. 350 (18), 1828–1837. doi:10.1056/NEJMoa032520

PubMed Abstract | CrossRef Full Text | Google Scholar

Mabbott, N. A., Baillie, J. K., Brown, H., Freeman, T. C., and Hume, D. A. (2013). An expression atlas of human primary cells: Inference of gene function from coexpression networks. BMC Genomics 14, 632. doi:10.1186/1471-2164-14-632

PubMed Abstract | CrossRef Full Text | Google Scholar

Najafov, A., Chen, H., and Yuan, J. (2017). Necroptosis and cancer. Trends Cancer 3 (4), 294–301. doi:10.1016/j.trecan.2017.03.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Reyna, D. E., Garner, T. P., Lopez, A., Kopp, F., Choudhary, G. S., Sridharan, A., et al. (2017). Direct activation of BAX by BTSA1 overcomes apoptosis resistance in acute myeloid leukemia. Cancer Cell 32 (4), 490–505. doi:10.1016/j.ccell.2017.09.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Ricciuti, B., Kravets, S., Dahlberg, S. E., Umeton, R., Albayrak, A., Subegdjo, S. J., et al. (2019). Use of targeted next generation sequencing to characterize tumor mutational burden and efficacy of immune checkpoint inhibition in small cell lung cancer. J. Immunother. Cancer 7 (1), 87. doi:10.1186/s40425-019-0572-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Riley, R. S., June, C. H., Langer, R., and Mitchell, M. J. (2019). Delivery technologies for cancer immunotherapy. Nat. Rev. Drug Discov. 18 (3), 175–196. doi:10.1038/s41573-018-0006-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Safran, M., Rosen, N., Twik, M., BarShir, R., Stein, T. I., Dahary, D., et al. (2021). “The GeneCards suite,” in Practical guide to life science databases. Editors I. Abugessaisa, and T. Kasukawa (Singapore: Springer Nature Singapore), 27–56.

CrossRef Full Text | Google Scholar

Sahin, I. H., Askan, G., Hu, Z. I., and O'Reilly, E. M. (2017). Immunotherapy in pancreatic ductal adenocarcinoma: An emerging entity? Ann. Oncol. 28 (12), 2950–2961. doi:10.1093/annonc/mdx503

PubMed Abstract | CrossRef Full Text | Google Scholar

Scheiner, B., Pomej, K., Kirstein, M. M., Hucke, F., Finkelmeier, F., Waidmann, O., et al. (2022). Prognosis of patients with hepatocellular carcinoma treated with immunotherapy - development and validation of the CRAFITY score. J. Hepatol. 76 (2), 353–363. doi:10.1016/j.jhep.2021.09.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Seifert, L., and Miller, G. (2017). Molecular pathways: The necrosome-A target for cancer therapy. Clin. Cancer Res. 23 (5), 1132–1136. doi:10.1158/1078-0432.CCR-16-0968

PubMed Abstract | CrossRef Full Text | Google Scholar

Tanzer, M. C., Khan, N., Rickard, J. A., Etemadi, N., Lalaoui, N., Spall, S. K., et al. (2017). Combination of IAP antagonist and IFNγ activates novel caspase-10- and RIPK1-dependent cell death pathways. Cell Death Differ. 24 (3), 481–491. doi:10.1038/cdd.2016.147

PubMed Abstract | CrossRef Full Text | Google Scholar

Uhlén, M., Fagerberg, L., Hallström, B. M., Lindskog, C., Oksvold, P., Mardinoglu, A., et al. (2015). Proteomics. Tissue-based map of the human proteome. Sci. (New York, N.Y.) 347 (6220), 1260419. doi:10.1126/science.1260419

CrossRef Full Text | Google Scholar

Vickers, A. J., and Elkin, E. B. (2006). Decision curve analysis: A novel method for evaluating prediction models. Med. Decis. Mak. 26 (6), 565–574. doi:10.1177/0272989X06295361

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, J., Tian, W., Tang, Y., Zou, Y., Zheng, S., Wu, L., et al. (2022). Establishment of a cell necroptosis index to predict prognosis and drug sensitivity for patients with triple-negative breast cancer. Front. Mol. Biosci. 9, 834593. doi:10.3389/fmolb.2022.834593

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, W., Ma, Y., Liu, Y., Smith-Warner, S. A., Simon, T. G., Chong, D. Q., et al. (2019). Association of intake of whole grains and dietary fiber with risk of hepatocellular carcinoma in US adults. JAMA Oncol. 5 (6), 879–886. doi:10.1001/jamaoncol.2018.7159

PubMed Abstract | CrossRef Full Text | Google Scholar

Yatim, N., Jusforgues-Saklani, H., Orozco, S., Schulz, O., Barreira da Silva, R., Reis e Sousa, C., et al. (2015). RIPK1 and NF-κB signaling in dying cells determines cross-priming of CD8⁺ T cells. Sci. (New York, N.Y.) 350 (6258), 328–334. doi:10.1126/science.aad0395

PubMed Abstract | CrossRef Full Text | Google Scholar

Yuan, J., Amin, P., and Ofengeim, D. (2019). Necroptosis and RIPK1-mediated neuroinflammation in CNS diseases. Nat. Rev. Neurosci. 20 (1), 19–33. doi:10.1038/s41583-018-0093-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeng, D., Ye, Z., Shen, R., Yu, G., Wu, J., Xiong, Y., et al. (2021). Iobr: Multi-Omics immuno-oncology biological research to decode tumor microenvironment and signatures. Front. Immunol. 12, 687975. doi:10.3389/fimmu.2021.687975

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, J., Bajari, R., Andric, D., Gerthoffert, F., Lepsa, A., Nahal-Bose, H., et al. (2019). The international cancer genome consortium data portal. Nat. Biotechnol. 37 (4), 367–369. doi:10.1038/s41587-019-0055-9

PubMed Abstract | CrossRef Full Text | Google Scholar