Comprehensive genomic signature of pyroptosis-related genes and relevant characterization in hepatocellular carcinoma

Data sources

Gene expression data, clinical characteristics and survival information of patients in three independent HCC cohorts (TCGA-LIHC, n = 365; LIRI-JP, n = 208; GSE14520, n = 242) were downloaded from three public databases: The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/repository), the International Cancer Genomics Consortium (ICGC, Toronto, Canada, https://dcc.icgc.org/), and the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/). The criteria for data screening were to retain survival information for patients with non-zero survival time in three HCC cohorts for further analysis. Our research followed public data access policies and did not address ethical relationships. Data collection was in full accordance with the TCGA, ICGC and GEO database usage guidelines. Furthermore, an immunotherapy dataset (Imvigor210) was downloaded using the “IMvigor210CoreBiologies” R package (Mariathasan et al., 2018). Meanwhile, the 48 PRGs used in this study were extracted from the MsigDB database (http://www.gsea-msigdb.org/gsea/msigdb/) and published literature (Table S1) (Shao et al., 2021; Ye, Dai & Qi, 2021). The flow chart of our work is illustrated in Fig. S1.

Consensus clustering analysis of differentially expressed PRGs

Differentially expressed genes (DEGs) in 365 HCC samples and 50 adjacent non-tumor samples from the TCGA-LIHC cohort were determined using the “limma” R package with the criteria of |log2(FC)| > log2(1.5) and FDR < 0.05, and then differentially expressed PRGs were extracted from DEGs. The volcano plot of DEGs was drawn using the “ggplot2” R package.

Consensus unsupervised clustering analysis was implemented using the “Consensus ClusterPlus” R package to divide HCC patients into distinct subtypes according to differentially expressed PRGs. Consensus clustering was performed using the k-means method and 1,000 repetitions were run to ensure the stability of the classification. This clustering met the following requirements: the cumulative distribution function (CDF) curve slowly and gradually elevated; the number of patients without any group was small; after clustering, the intergroup and intragroup correlation decreased and increased, respectively. Meanwhile, the optimum k value for consensus clustering was determined by combining the proportion of ambiguous clustering (PAC). The heatmap plot of PRGs was drawn using the “pheatmap” R package. The difference in overall survival (OS) between the two subtypes was analyzed using Kaplan–Meier curves.

Determination and functional annotation of DEGs between the two subtypes

DEGs in the two pyroptosis-related subtypes were determined by the “limma” R package based on the following criteria: |log2(FC)| > 2 and FDR < 0.01. We performed Kyoto Encyclopedia of Genes and Genomes (KEGG) as well as Gene Ontology (GO) functional enrichment analysis of pyroptosis-related DEGs based on the “clusterprofiler” R package to investigate the latent functions of these DEGs and identify relevant enrichment pathways. Following, pathways and terms were defined as statistically significant using FDR < 0.05.

Construction of the pyroptosis-related prognostic PRGS

We utilized univariate Cox regression analysis to determine the pyroptosis-related DEGs correlated with OS in the training cohort (TCGA-LIHC, n = 365), screened by the criterion of p-value < 0.0001. Next, the LASSO Cox algorithm eliminated overfitting between prognostic genes to perform penalty parameter tuning by ten-fold cross validation according to the “glmnet” R package in an attempt to narrow the range of prognostic genes. The pyroptosis-related risk score for each patient with HCC was calculated using gene expression levels multiplied by the LASSO Cox regression coefficients. Based on this, we constructed an efficient prognostic prediction model PRGS to obtain the optimal prediction outcomes from the TCGA-LIHC cohort using overall survival.

The calculation formula of PRGS was: PRGS = Σ (Coefi * Expi), where Coefi and Expi denote the regression coefficient and corresponding expression level of each prognostic gene, respectively. A total of 365 patients from the TCGA-LIHC cohort were classified into low-risk (PRGS < median value) and high-risk (PRGS > median value) groups according to the median PRGS, followed by Kaplan–Meier survival and receiver operating characteristic (ROC) curves to assess the predictive power of PRGS model. Then, we performed principal component analysis (PCA) according to the “factoextra” R package. Likewise, HCC patients in the testing cohorts (LIRI-JP and GSE14520) were also classified into low- and high-risk groups, and each cohort underwent Kaplan–Meier survival and ROC curves.

Clinical correlation analysis of the prognostic PRGS

The Chi-square test was performed to assess the correlation between PRGS and multiple clinical characteristics (age, gender, stage, TNM, status, grade). Univariate and multivariate Cox regression analyses were conducted in the training and testing cohorts to assess whether PRGS was a prognostic factor independent of other available clinical characteristics. The p-value < 0.05 for pyroptosis-related signature PRGS and other clinical characteristics were considered as independent prognostic risk factors.

Development and validation of the nomogram

We developed a nomogram model to forecast the survival time of HCC patients according to the prognostic model PRGS and different clinical characteristics (age, gender, stage,grade, and T stage) using the “survival” and “rms” R packages. In the nomogram, each factor has a corresponding score, and the scores of all factors for every patient are added together to gain the total score to predict the corresponding overall survival time. The C-index was applied to appraise the discriminative ability of the nomogram. The predictive power of the nomogram was assessed using ROC curves for 1-, 3-, and 4-year survival. Furthermore, we used calibration plots of the nomogram to exhibit the consistency between the actual observed survival and the predicted 1-, 3-, and 4-year survival.

Evaluation of the tumor immune microenvironment

To assess the proportions of infiltrating immune cells in the tumor immune microenvironment, we utilized the CIBERSORT algorithm to calculate the abundance of 22 infiltrating immune cells in low- and high-risk groups. Next, we applied single-sample gene set enrichment analysis (ssGSEA), the enrichment scores for 13 immune functions and 16 immune cells in HCC patients were calculated based on the specific expressed gene signatures using the “GSVA” R package. In addition, the associations between PRGS and the enrichment scores for 13 immune functions and 16 immune cells were also investigated.

Immunotherapy efficacy and drug susceptibility analysis

Immunophenoscore (IPS) is a good predictor of ICIs response, and its score is determined by quantifying vital components of tumor immunity. Specifically, IPS scores are calculated based on z-scores of immune-related gene expression in representative cell types, ranging from 0 to 10. High IPS scores reflect enhanced immunogenicity. We obtained IPS scores for patients with HCC (from the TCGA-LIHC cohort) from The Cancer Immunome Atlas (TCIA, https://www.tcia.at/home), and compared the distribution of IPS scores in low- and high-risk groups. Additionally, the expression of several typical immune checkpoints in low- and high-risk groups was also analyzed. The IMvigor210 cohort was applied to assess the immunotherapy efficacy in low- and high-risk groups.

To assess differences in the efficacy of targeted therapy and chemotherapy drugs between patients in low- and high-risk groups, we applied the “pRRophetic” R package to predict the half-maximal inhibitory concentrations (IC50) values of chemotherapeutic drugs commonly utilized to treat HCC. The common chemotherapeutic drugs, such as cisplatin, gefitinib, sorafenib were selected for HCC.

Cell culture and real-time PCR analysis

Human normal liver LO2 cell line and liver cancer Huh7 and HepG2 cell lines (obtained from the cell bank of the Chinese Academy of Science) were applied to detect the expression of five prognostic genes of PRGS. Total RNA from LO2, Huh7 and HepG2 cell lines was isolated using Trizol reagent. Total RNA was pretreated with RNase‑free DNase and then reverse‑transcribed with AMV to synthesize cDNA. The RT-qPCR reactions were conducted using the Light Cycler 480 SYBR Green I kit with specific primers commercially synthesized by Invitrogen Biotechnology Ltd. (Table S2). The 20 μl of reaction mixture was composed of 1 μl of diluted cDNA, 2 μl of primer mix (reverse and forward primers), 7 μl of 10X buffer and 10 μl of Mix (2X PCR Buffer, 2X dNTP, 2X Taq DNA Polymerase). We performed amplification reactions using the LightCycler® 480 instrument with an initial step at 95 °C for 5 min, followed by 50 cycles of three step PCR at 95 °C for 15 s, 60 °C for 15 s and 72 °C for 30 s. Quantitative analysis was normalized using β-actin.

Statistical analysis

The Chi-square test and Student’s t-tests were utilized to calculate the significant differences between variables. For baseline clinical data, Kruskal–Wallis and Wilcoxon rank-sum tests were implemented to assess statistical differences between multiple and two groups, respectively. The log-rank test was utilized to assess differences in Kaplan–Meier curves. The p-value < 0.05 represented a statistically significant difference. All statistical analyses involving data in this study were conducted using R 4.0.3 and Prism 8 software.

Genetic and expression alterations of PRGs in HCC

In this study, we selected 48 PRGs for analysis. At the genetic level, the occurrence of somatic mutations in 48 PRGs was analyzed, and the results showed that 150 (41.21%) of the 364 HCC samples had mutations in the PRGs (Fig. 1A). Among all 48 PRGs, the gene with the highest mutation incidence was TP53 (28%), followed by NLRP3 and NLRP2. Figure 1B depicts the correlation network including all PRGs. We then evaluated somatic copy number alterations in 48 PRGs and found that such alterations were present in most PRGs. Among these PRGs, copy number variation (CNV) of GSDMC, GSDMD, AIM2, CHMP4C, NLRP3 and CHMP6 was generally increased, while those of CASP3, NLRP1, CHMP7 and IRF2 was decreased (Fig. 1C). Figure 1D exhibits the location of CNV alterations in PRGs on their corresponding chromosomes.

To compare the differential expression of PRGs in 365 HCC samples and 50 non-tumor samples from the training (TCGA-LIHC) cohort, we first conducted differential expression analysis on all genes, and then identified the corresponding expression differences for 48 PRGs (Fig. 1E). Volcano plot was utilized to describe the DEGs between HCC and non-tumor samples (Fig. 1F). We further analyzed the possible relationship between CNV alterations and differential expression of PRGs. Interestingly, PRGs with CNV gain, including GSDMD, GSDMC, CHMP4C and CHMP6, were significantly higher in HCC samples than in adjacent non-tumor samples (p < 0.001), suggesting that CNVs may regulate PRGs expression. However, compared with non-tumor samples, some PRGs with CNV gain, including NLRP3, were significantly downregulated in HCC (p < 0.01); some PRGs with CNV loss, including CASP3, were significantly upregulated in HCC (p < 0.001). Therefore, although CNV can be used to explain most of the alterations in PRGs expression, CNV is not a determinant of gene expression regulation (Stranger et al., 2007). Other important factors modulating gene expression levels include histone modifications and DNA methylation (Angeloni & Bogdanovic, 2019; Pope & Medzhitov, 2018). Our analysis revealed remarkable differences in genetic alterations and expression levels of PRGs in HCC and non-tumor samples, suggesting a latent role of PRGs in HCC tumorigenesis.

Identification of pyroptosis-related subtypes in HCC

Based on the levels of differentially expressed PRGs, we identified distinct pyroptosis-related subtypes in the training (TCGA-LIHC) cohort using a consensus clustering algorithm (Fig. S2). Our results suggested that the clustering variable k = 2 seemed to be the optimum choice for classifying the whole cohort into pyroptosis-related cluster 1 (n = 170) and pyroptosis-related cluster 2 (n = 195) (Fig. 2A, Table S3). The results of principal component analysis (PCA) suggested that two pyroptosis-related subtypes, namely cluster 1 and cluster 2, were clustered separately (Fig. 2B). Besides, Silhouette analysis showed results consistent with PCA (Fig. S3). Kaplan–Meier curves showed that patients in cluster 2 had significantly longer OS than those in cluster 1 (p = 0.037, Fig. 2C). In addition, comparing the clinical characteristics of cluster 1 and cluster 2 in HCC using the Chi-square test revealed remarkable differences between the two subtypes. As shown in Fig. 2D, cluster 2 had better survival status, lower grade, and lower T stage than cluster 1 (p < 0.05), which corresponds to the previously described cluster 2 having longer OS.

Next, we identified DEGs between two pyroptosis-related subtypes, and used the obtained pyroptosis-related DEGs for subsequent analysis. To evaluate the latent functions of pyroptosis-related subtypes, we identified 203 pyroptosis-related DEGs in cluster 1 and cluster 2, and then conducted biological function analysis, including GO and KEGG enrichment analysis (Figs. 2E and 2F). According to GO analysis, in biological process (BP) category, pyroptosis-related DEGs were primarily enriched in drug metabolic process, humoral immune response and regulation of complement activation; in cellular component (CC) term, the pyroptosis-related DEGs were primarily enriched in blood microparticle, endoplasmic reticulum lumen and immunoglobulin complex; regarding molecular function (MF), pyroptosis-related DEGs were primarily enriched in antigen binding and iron ion binding (Fig. 2E). As for KEGG pathway analysis, pyroptosis-related DEGs were significantly enriched in chemical carcinogenesis, TNF signaling pathway and IL-17 signaling pathway (Fig. 2F).

Construction and validation of the prognostic PRGS

A total of 365 HCC patients in the TCGA-LIHC cohort were utilized to construct the prognostic prediction model. First, univariate Cox analysis was performed to determine the association between 203 pyroptosis-related DEGs and OS outcomes in HCC patients. Among the 203 pyroptosis-related DEGs, five genes were screened. To avoid overfitting between genes, we performed LASSO Cox regression using a minimized lambda (Fig. S4A). The results of the LASSO Cox analysis suggested that all five genes had non-zero regression coefficients and were extracted for subsequent analysis. These five genes and their corresponding coefficients were utilized to develop a pyroptosis-related gene scoring model, named “PRGS”. These five prognostic genes were utilized to construct the model PRGS, and the calculation formula was as follows: Risk score = (−0.06198467 * expression level of CYP2C9) + (0.03381356 * expression level of MYBL2) + (0.08247219 * expression level of SPP1) + (0.20457157 * expression level of CTSV) + (0.15033404 * expression level of EPO) (Fig. 3A).

Based on the above calculation formula, we calculated the risk score for every patient in the TCGA-LIHC cohort. Then, based on the median PRGS, patients were classified into low- and high-risk groups. The results of PCA showed that low- and high-risk groups were clustered separately (Fig. S4B). The overall survival of patients in the low-risk group was significantly better than that in the high-risk group (p < 0.0001, Fig. 3B). The distribution of PRGS, survival time and survival status for low- and high-risk groups were shown in Fig. 3C, and patients with death status increased with PRGS. As exhibited in Fig. 3D, in the training cohort, the AUC values of PRGS for 1-, 3-, and 4-year were 0.750, 0.684, and 0.715, respectively. We performed univariate and multivariate Cox regression analyses to investigate in depth whether PRGS model was independent of other clinical characteristics, including age, gender, grading, T stage and stage (Table 1). The results indicated that PRGS was identified as an independent prognostic risk factor (HR = 2.33, 95% CI [1.55–3.50], p < 0.001).

To demonstrate the robustness and versatility of model PRGS, we also validated this risk score in the testing cohorts (LIRI-JP and GSE14520) and obtained similar results. We used the same formula to classify patients in the testing cohorts into low- and high-risk groups. In the LIRI-JP cohort, the OS of patients in the low-risk group was significantly prolonged compared to patients in the high-risk group (p < 0.0001, Fig. 4A). The distribution of PRGS, survival time and survival status were exhibited in Fig. 4B, and patients with death status increased with PRGS. The AUC values of PRGS for 1-, 3-, and 4-year were 0.751, 0.717, and 0.773, respectively (Fig. 4C). Likewise, in the GSE14520 cohort, patients in the low-risk group also had significantly longer overall survival than those in the high-risk group (p = 0.0001, Fig. 4D). The distribution of PRGS, survival time and survival status were shown in Fig. 4E, and patients with survival status decreased with increasing PRGS. The AUC values of PRGS for 1-, 3-, and 4-year were 0.667, 0.680, and 0.653, respectively (Fig. 4F). Similarly, univariate and multivariate Cox regression analyses of PRGS and clinical characteristics such as age, gender and stage were carried out in the testing cohorts (Table 1). As a result, PRGS was an independent prognostic factor with excellent predictive power for survival in HCC patients (HR = 3.48, 95% CI [1.63–7.43], p = 0.001; HR = 1.61, 95% CI [1.01–2.60], p = 0.047). Moreover, in the training and testing cohorts, we applied ROC curves to compare the prognostic predictive power of PRGS with other pyroptosis-related signatures. We found that PRGS consistently maintained superior performance compared to other published pyroptosis-related signatures (Figs. 4G and S5, Table S4) (Fu & Song, 2021; Jin et al., 2022; Liu et al., 2021; Wu et al., 2021b). These results suggested that PRGS was an extremely reliable prognostic prediction model and outperformed other pyroptosis-related signatures.

Validation of the expression levels of five genes for PRGS

The expression levels of five prognostic genes in human normal liver cells and liver cancer cells were detected by RT-qPCR. As illustrated in Figs. 5A and 5B, the expression levels of MYBL2, SPP1, CTSV and EPO were remarkably increased in liver cancer cells compared with normal liver cells, while the expression level of CYP2C9 was significantly decreased (p < 0.05). Notably, the expression of these five genes was not statistically different between Huh7 and HepG2 cell lines (Fig. S6). Furthermore, we assessed the protein levels of these five prognostic genes in normal and HCC liver tissues using the Human Protein Atlas (HPA, https://www.proteinatlas.org/) database. Immunohistochemical (IHC) staining information for five prognostic genes in normal and HCC liver tissues from the HPA database was reflected in Fig. 5C. We found that the protein levels of SPP1, MYBL2, CTSV and EPO were increased in HCC tissues, while the protein level of CYP2C9 was decreased.

Associations between PRGS and clinical characteristics of HCC patients

To evaluate the association between PRGS and clinical characteristics, we explored the clinical manifestations of low- and high-risk groups in both the training and testing cohorts, including age, gender, stage, TNM, status, grade. We created heatmaps of clinical characteristics for three cohorts (Figs. 6A–6C). In three cohorts, except for CYP2C9, the expression levels of the remaining of the four genes (MYBL2, SPP1, CTSV and EPO) were elevated in the in high-risk group. We observed that the survival status of HCC patients in three cohorts was extremely consistent with the overall survival advantage described in our previous study (Figs. 6A–6C). In the TCGA-LIHC and LIRI-JP cohorts, early stage (I–II) patients dominated the low-risk group, whereas advanced stage (III–IV) patients were mainly divided into the high-risk group (p < 0.01, Figs. 6A and 6B). Likewise, in the TCGA-LIHC and GSE14520 cohorts, patients with early T (or TNM) stage were predominantly distributed in the low-risk group, whereas patients with advanced T (or TNM) were mostly distributed in the high-risk group (p < 0.01, Figs. 6A, 6C). This explains why the low-risk group was associated with a better overall survival advantage. Importantly, in the GSE14520 cohort, we noted that AFP (<=300 ng/ml), tumor size (<=5 cm) and non-multinodular were associated with the low-risk group, whereas AFP (>300 ng/ml), tumor size (>5 cm) and multinodular were mainly associated with the high-risk group (p < 0.05, Fig. 6C). In addition, in the TCGA-LIHC cohort, we also assessed whether PRGS maintained its predictive power in various clinical characteristic subgroups, such as age (<=60 and >60 years), gender (female and male), stage (I–II and III–IV), grade (1–2 and 3–4), T stage (T1-2 and T3-4), and fibrosis (fibrosis & cirrhosis and no fibrosis), and the results showed that PRGS could be used to forecast OS in various subgroups of clinical characteristics (except for the female subgroup) (p < 0.05, Fig. S7). These findings illustrated that low- and high-risk groups of PRGS represented distinct clinical characteristics of HCC patients, and that PRGS could maintain excellent predictive power in clinical characteristic subgroups.

Establishment and evaluation of the predictive nomogram

Considering the importance of PRGS in the prognosis prediction of HCC patients, we then tried to investigate its feasibility in clinical applications. We established a nomogram containing PRGS and clinical characteristics that were readily available and recognized to have some degree of influence on the prognosis of HCC, to predict 1-, 3-, and 4-year survival in the TCGA-LIHC cohort (Fig. 7A). The C-index of 0.715 demonstrated the excellent predictive power of the nomogram. Calibration plot of the nomogram indicated strong agreement between observed and predicted 1-, 3-, and 4-year survival (Fig. 7B). The AUC values of the nomogram model at 3-, and 4-year showed higher accuracy in predicting survival (Fig. 7C). Specifically, the AUC values of the nomogram for 1-, 3-, and 4-year were 0.739, 0.722, and 0.723, respectively. As shown in Fig. 7D, the alluvial plot was utilized to describe changes in the characteristics of HCC patients.

Furthermore, we also constructed nomograms in the LIRI-JP and GSE14520 cohorts to forecast 1-, 3-, and 4-year survival (Figs. S8A and S8B). The C-index of the nomogram was 0.769 in the LIRI-JP cohort and 0.702 in the GSE14520 cohort, indicating superior predictive performance for both nomograms. Calibration plots showed that the nomograms could accurately estimate 1-, 3-, and 4-year survival (Figs. S8C and S8D). The 1-, 3-, and 4-year AUC values of the nomograms in the LIRI-JP and GSE14520 cohorts were 0.726, 0.758, and 0.744, as well as 0.728, 0.726, and 0.599, showing good accuracy (Figs. S8E and S8F). In addition, in the training and testing cohorts, we compared the prognostic predictive power of our nomogram with other published pyroptosis-related nomograms (Fu & Song, 2021; Jin et al., 2022; Liu et al., 2021). The results showed that our nomogram has better predictive power (Figs. S8G–S8I, Table S5). According to the above results, we observed that PRGS can play a critical role in the clinical prognosis prediction of HCC patients.

Evaluation of the immune microenvironment between low- and high-risk groups

We identified the abundance of 22 infiltrating immune cells in each HCC patient between low- and high-risk groups according to the CIBERSORT algorithm. In the TCGA-LIHC cohort, the infiltration abundance of macrophages M1, B cells naive, macrophages M2 and T cells gamma delta in the low-risk group were remarkably higher than those in the high-risk group (p < 0.01, Fig. 8A), whereas the infiltration abundance of B cells memory, T cells regulatory (Tregs), macrophages M0 and dendritic cells resting in the low-risk group were significantly lower than those in the high-risk group (p < 0.05, Fig. 8A). In the GSE14520 cohort, the abundance of B cells naive, Macrophages M1 as well as T cells gamma delta were significantly increased in the low-risk group compared to the high-risk group (p < 0.05, Fig. 8B); conversely, the abundance of macrophages M0, B cells memory and T cells regulatory (Tregs) were remarkably higher in the high-risk group compared to the low-risk group (p < 0.05, Fig. 8B). In the LIRI-JP cohort, the infiltration abundance of macrophages M1 was remarkably increased in the low-risk group compared with the high-risk group, while the infiltration abundance of neutrophils was significantly decreased (p < 0.01, Fig. 8C). Interestingly, we observed a significantly higher infiltration abundance of macrophages M1 in the low-risk group compared with the high-risk group in all three cohorts (p < 0.05).

Furthermore, we explored the prognostic value of 22 infiltrating immune cells, and in the TCGA-LIHC cohort, the infiltration abundance of T cells gamma delta, T cells CD8 and macrophages M0 were significantly correlated with OS (p < 0.05, Fig. 8D). Higher infiltration abundance of macrophages M0 was associated with worse overall survival, whereas higher infiltration abundance of T cells CD8 and T cells gamma delta were related to better OS. In the GSE14520 cohort, the infiltration abundance of macrophages M1, macrophages M0 and B cells naive were significantly associated with OS in HCC patients (p < 0.05, Fig. 8E). The higher infiltration abundance of macrophage M0 showed poorer overall survival, whereas the higher infiltration abundance of macrophages M1 and B cells naive showed better overall survival.

Next, we calculated the enrichment scores for 13 immune functions and 16 immune cells in three cohorts using ssGSEA to further evaluate the association between PRGS and the scores of 13 immune functions and 16 immune cells (Figs. 9A and 9B). The results showed that Tregs were significantly positively correlated with PRGS in three cohorts (p < 0.05). The lower PRGS was closely related to B cells and NK cells. In immune functions, the higher PRGS was strongly associated with chemokine receptor (CCR) and immune checkpoints. These results demonstrated that PRGS was associated with many immune cells and immune functions, implying that PRGS could potentially reflect the tumor immune microenvironment.

Immune checkpoints and drug susceptibility analysis

To further examine the association of PRGS with immune checkpoints, we evaluated differences in immune checkpoints in low- and high-risk groups. In the training cohorts, immune checkpoints were expressed at higher levels in the high-risk group (p < 0.05, Fig. 10A). Overexpressed immune checkpoints promote tumor immune escape, which also explains that elevated immune checkpoints expression is associated with poorer prognosis in HCC patients (Xu et al., 2018; Zheng et al., 2021). In addition, the correlation analysis between PRGS and immune checkpoints showed that PRGS was significantly positively correlated with PD1, CTLA4, LAG3, PVRIG and TIGIT (p < 0.01, Fig. 10B).

In addition, we utilized IPS scores to determine the patient response to ICIs treatment. Four subtypes of IPS scores were performed to forecast response to anti-PD1 and anti-CTLA4 therapy in HCC patients. We found higher IPS scores in the low-risk group, indicating that patients in the low-risk group may respond better to immunotherapy (Fig. S9). Due to the complexity of ICIs in the tumor immune microenvironment, we further examined whether ICIs have an impact on the prognostic prediction of PRGS. As illustrated in Fig. 10C, the low-risk group with high PD1 exhibited better OS compared with the high-risk group and high PD1, and the OS of the low-risk group with low PD1 was longer than that of the high-risk group and low PD1 (p < 0.0001). Similarly, CTLA4 exhibited a survival prediction consistent with PD1 (p < 0.0001, Fig. 10C). Then, we utilized ROC curves to compare the prognostic predictive power of PRGS and ICIs, and the results showed that PRGS was superior to PD1 and CTLA4 (Fig. 10D).

We next selected targeted therapy and chemotherapy drugs currently utilized to treat HCC to evaluate the sensitivity of low- and high-risk groups to these drugs. Interestingly, the results showed that the IC50 values of bleomycin were absolutely decreased in the high-risk group, whereas the IC50 values of targeted therapy and chemotherapy drugs such as gefitinib, sorafenib, docetaxel and cisplatin were remarkably lower in the low-risk group (p < 0.01, Fig. 10E). Taken together, the above results suggested that PRGS was associated with drug sensitivity.

The role of PRGS in immunotherapy

Given the relevance of PRGS to the tumor immune microenvironment, we further validated the effect of PRGS on immunotherapy efficacy in an immunotherapy cohort (IMvigor210). In order to make PRGS more convenient for the determination of low- and high-risk in patients receiving immunotherapy, we need a cut-off value to assist. Here, the cut-off value was 0.915, which is the optimum cut-off value obtained using the “survminer” R package. Patients in the IMvigor210 cohort were grouped into low-risk (PRGS < 0.915) and high-risk (PRGS > 0.915) groups based on the cut-off value, and we found that patients in the low-risk group exhibited better OS compared to those in the high-risk group (p = 0.0045, Fig. 11A). Not surprisingly, applying the cut-off value to the training and testing cohorts described above, the results also showed that the OS of patients in the low-risk group was superior to that in the high-risk group (p < 0.01, Fig. S10). Moreover, we also analyzed the relationship between PRGS and immunotherapy response, immune cell and tumor cell subgroups, and the results indicated that PRGS was not associated with immunotherapy response subgroups (Fig. 11B), but was associated with immune cell and tumor cell subgroups (p < 0.01, Figs. 11C and 11D). Interestingly, we found that PRGS was also positively correlated with immune checkpoints in the IMvigor210 cohort (p < 0.05, Fig. S11).

Growing evidence suggests that cancer patients with higher tumor mutation burden (TMB) might benefit more from immunotherapy because of their higher tumor neoantigen burden (TNB) (Hollern et al., 2019). In the IMvigor210 cohort, high TMB and TNB represented better overall survival (p < 0.01, Figs. 11E and 11F). ROC curve analysis showed that the AUC value (0.788) obtained by the combination of TMB, TNB and PRGS was higher than that of TMB (0.685), TNB (0.697) or PRGS (0.532) alone (Fig. 11G). Next, we further evaluated the impact of TMB and TNB on the prognostic prediction of PRGS. Figure 11H showed that the low-PRGS with high TMB exhibited better overall survival compared with high-PRGS and high TMB, and the overall survival of the low-PRGS with low TMB were superior to that of high-PRGS and low TMB (p < 0.0001). Similarly, TNB exhibited a survival prediction consistent with TMB (p = 0.00035, Fig. 11I), and PRGS combined with TMB or TNB predicted more significant results than either alone. These results suggested that considering PRGS together with TMB and TNB could improve the accuracy of assessing the efficacy of immunotherapy and determine patients who will respond better to immunotherapy.