Identification of prognostic pyroptosis-related DEGs and Construction of a prognostic model
A total of 50 normal tissue samples and 374 tumor tissue samples of HCC were screened from the TCGA database. Table 1 shows the clinicopathological data of the patients with HCC included in the study. A total of 30 differential genes(DEGs) were identified between tumor and the adjacent nontumor tissues. The DEGs expression in tumor and the adjacent nontumor tissues was represented in the heatmap(Fig. 1a,green: low expression level;red:high expression level). Also, The relationship between these genes is presented in Fig. 1b,which suggesting the tight relationship between these genes.
The prognostic values of DEGs were evaluated by univariate Cox analysis, and 7 genes that were significantly related to the OS of patients were identified(Fig. 1c). To build a prognostic model for HCC, we used training dataset and applied the LASSO Cox regression analysis to identify stable markers from 7-survival-related genes mentioned above. To find an optimal λ, 10-fold cross validation with minimum criteria was employed, so the final model gave minimum cross validation error. Then, 3 DEGs(BAK1, GSDME, and NOD2)were identified(Fig. 2a-2b). Using the following formula, the risk score was calculated 0.045673* expression level of BAK1 + 0.005602* expression level of NOD2 + 0.143950* expression level of GSDME. A median cut-off value was applied to stratify patients into a high-risk group (n = 183) and a low-risk group (n = 182)(Fig. 3a). The PCA showed that patients in different risk groups were distributed in two directions by the constructed model cohort (Fig. 3b). In the high risk score group, the risk of mortality obviously increased with the risk score increased (Fig. 3c). Kaplan-Meier analysis showed that the probability of survival was significantly higher in low-risk patients than in high-risk patients(Fig. 3d,p < 0.05). The calibration plots showed that calibration appeared good between the predicted and observed risks(Fig. 3e). Also, the AUC of the time-dependent ROC curve was plotted to assess the prognostic effect. the AUC of the prognostic risk assessment model was 0.844 at 1 years, 0.889 at 3 years, 0.774 at 5 years (Fig. 3f).
Validation Of The 3-gene Prognostic Model
In order to evaluate the robustness of the prognostic model constructed using data from the training cohort, its performance was also assessed using validation cohort. The risk score for each patient in validation cohort was calculated with the same formula as training group. The AUC of the 3-gene signature was 0.787 at 1 years, 0.801 at 3 years, 0.820 at 5 years(Fig. 4a). Besides, patients in the high-risk group were significantly worse off the overall survival time compared to the low-risk group (Fig. 4b). Western Blot and qRT-PCR was used to detect GSDME expression. The results showed that GSDME expression in HCC was significantly higher than that in adjacent tissues(Fig. 5).
Independent Prognostic Value Of The Risk Model
To investigate the prognostic value of different clinical features in HCC, univariate Cox regression analysis was performed. According to the results of univariate Cox regression analysis, the risk model showed their prognostic value in predicting OS(Fig. 6a). The multivariable Cox regression analysis suggested risk score had prognostic value for HCC patients. Even after adjusting for confounding factors, the risk score is still a prognostic factor(Fig. 6b) for patients with HCC. Figure 6c shows heatmaps of the pyroptosis-related DEGs modules on the clinical data and found a different distribution of patients' age and survival status between low- and high-risk subgroups.
DEGs-based tumor classification.
In addition, We explore the connections between the expression of DEGs and HCC subtypes, consensus clustering analysis was used to classify HCC patients into subtypes. By increasing the clustering variable (k) from 2 to 10, we found that the highest intra-group correlation and low inter-group correlation was observed when k = 2 (Fig. 7a). The results showed that the HCC patients could be well divided into two clusters based on the DEGs.The gene expression profile were matched to the clinical data and was presented in a heatmap (Fig. 7b), We also compared the OS in 2 clusters.There were significant differences in 2 clusters (Fig. 7c).
Functional Analyses And Infiltrating Immune Cells
Meanwhile, GO function and KEGG pathway enrichment analysis were used to investigate the intrinsic pathway of DEGs. GO function analysis was performed using three ontologies: biological process (BP), cellular component (CC), and molecular function (MF). Based on GO analysis, the DEGs were most enriched in BP, and specifically in immune related processes(such as T cell activation, regulation of T cell activation, leukocyte cell − cell adhesion, Fig. 8a). In addition, KEGG pathway analyses indicated that these DEGs were highly enriched in Cytokine − cytokine receptor interaction, Neuroactive ligand − receptor interaction, Chemokine signaling pathway(Fig. 8b). To further explore relationship between immune cell infiltration and risk scores in the prognostic risk mode, ssGSEA score was used to quantify the level of enrichment of immune cell functions or pathways in cancer samples(Fig. 9a-b). One interesting finding is the many types of immune cells were significantly different between the two groups(such as macrophages, activated dendritic cell, Treg). Another intriguing fact was the immune cell functions were also significantly different between the two groups(such as APC_co_inhibition, CCR, Check − point, HLA, MHC_class_I, Type_II_IFN_Reponse). There are also significant differences in immune checkpoint-related genes in patients of different grades(Fig. 9c).