3.1 Overview of genetic and transcriptional variation of necroptosis-related genes in HCC
We first 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 nonsignificant alteration. The comparison between 50 normal and 371 HCC samples indicated that 13 genes were significantly 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 tumor-normal 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 . Some mechanisms, including DNA methylation, transcription factor, m6A modification, long non-coding RNAs, and RNA binding proteins, also play vital roles in gene expression [45–48]. Fig. 1F depicts the correlations of these genes and their prognostic significances. Collectively, our results indicated that the expression of pyroptosis-related genes was significantly different between HCC and normal samples, which could be a potential contributing factor to tumorigenesis and heterogeneity.
3.2 Identification of necroptosis-related clusters in HCC
According to the results of k-means consensus clustering, we identified 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 significantly worse (P < 0.001), and most necroptosis-related genes were upregulated (Fig. 2B, C). As presented in Fig. 2E, there were significant differences between the two clusters in transcription profile, age, T stage, TNM stage, tumor grade, and gender. The top 100 of 1517 DEGs were plotted in the heatmap (Additional file 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 findings of GO enrichment (Fig. 2D, F; Additional file 9: Table S4).
3.3 Development and validation of a ten-gene signature based on DEGs between necroptosis related clusters.
After adjustment of batch effect, only genes shared in TCGA, ICGC and GEO cohorts were remained, so we only screened prognostic gene signature from 985 of 1517 DEGs (Additional file 3: Figure S3A). Univariate Cox regression identified 376 genes were significantly related to OS of HCC patients in TCGA-LIHC cohort (P < 0.001; Additional file 9: Table S5). Next, a 10-fold cross-validated Lasso-Cox regression was performed to determine the optimal number of genes and coefficients for the predictive model (Additional file 3: Figure S3B, C). The risk scores of TCGA, ICGC and GEO cohorts were calculated as follows: risk score = (0.0276 * MTMR2 expression) + (0.1250 * CDCA8 expression) + (0.0291 * SLC1A5 expression) + (0.0577 * G6PD expression) + (0.0469 * PLOD2 expression) + (0.0324 * MMP1 expression) + (0.0385 * SLC2A1 expression) + (-0.0136 * ANXA10 expression) + (0.0044 * S100A9 expression) + (0.0050 * SPP1 expression).
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 significantly worse than that of the low-risk group (P < 0.001; Fig. 3C). Subsequently, we performed time-dependent ROC analyses to evaluate the predictive ability of this prognostic model. The area under curves (AUCs) were 0.827/0.729/0.695 at 1/3/5 years in the TCGA cohort, 0.786/0.782/0.800 at 1/3/4 years in the ICGC cohort, and 0.645/0.668/0.668 at 1/3/5 years in the GEO cohort (Fig. 3D). Of note, the high-risk group was mainly composed of the C1 cluster and had significantly 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 stratification 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 file 4: Figure S3D
3.4 External validation of the expression and genetic alteration of the ten genes
We first investigated the transcription levels of the ten genes in the TCGA, ICGC, and GEO cohorts, finding 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 amplification 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, five 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 low-risk groups. S100A9 protein, also known as myeloid-related protein, typically resides in immune cells such as neutrophils, macrophages, and monocytes . 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 . So, we further compared the tumor microenvironment of tumor and adjacent normal tissues and confirmed that the fractions of neutrophils, macrophages M2, and monocytes were significantly higher in adjacent normal tissues (Additional File 4: Figure S4).
3.5 Construction and Validation of Predictive Nomogram
Univariate and multivariate Cox analyses identified 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 benefit of the nomogram was evident in most cases (Fig. 8).
3.6 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 reflecting TME and guiding individualized treatment. Seven published methods were used to estimate the correlation between risk score and immune infiltration 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 fibroblasts (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 significantly 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, finding that responders (CR/PR) to anti-PD1/PD-L1 treatment had significantly 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 efficacy. 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, fluorouracil, and methotrexate were more sensitive in the low-risk group (Fig. 11).