Identification and characterization of a novel molecular classification based on disulfidptosis-related genes to predict prognosis and immunotherapy efficacy in hepatocellular carcinoma

Data acquisition

The mRNA expression data of HCC were further retrieved by searching the GEO database with the following keywords: “hepatocellular carcinoma”, “HCC”, and “liver cancer”. After initial screening, 3 gene expression omnibus (GEO) profiles with clinical information (GSE14520 [9], GSE36376 [10], and GSE76427 [11]) were selected and downloaded. The profiles base on GPL571, GPL10558, and GPL10558 platform. 221 samples in GSE14520, 223 samples in GSE36376, 115 samples in GSE76427. We also adopted ICGC (tumor: 231) and TCGA (tumor: 365; normal: 50) public database to download the RNAseq data. A total of 1,155 HCC samples with clinical information were selected. We also downloaded transcatheter arterial chemoembolization (TACE) dataset (GSE104580) to predict the clinical value of DRG.score, and a single-cell RNA dataset (GSE140228) was downloaded to explore the expression of 2 hub genes (Table 1) [12].

Data preprocessing

A FPKM gene expression matrix was acquired from TCGA and converted into TPM format [13]. The merged expression matrix was then eliminated from batch effects and normalized using the R package “sva” [14].

Analysis of differential DRG in TCGA-LIHC

The DRG expression was extracted from the collated RNA sequence data. DEG of DRG were analyzed in HCC tissues and normal tissues using the “limma” package [15].

Cox univariate and multivariate analysis for DRG in TCGA-LIHC

The “Survival” data package was used to analyze the differential genes by Cox single factor analysis, and the hazard ratio (HR) and P value were obtained. The intersect genes with P < 0.05 in two survival analyses were selected by R software. Gender, age, TNM staging, AFP and DRG were subjected to univariate and multivariate COX regression analyses. The relationship between clinicopathological characteristics and DRG and survival prognosis of patients with HCC was investigated. The sensitivity and specificity of DRG and different clinical characteristics in predicting survival prognosis of HCC patients were assessed by the area under the curve, and the predictive ability of different clinicopathological characteristics and DRG was compared.

Molecular subtypes base on 10 DRG in meta cohort

We inferred consensus cluster based on the DRG expression using the R package ConsensusClusterPlus. The optimal cluster number k was chosen depending on the elbow and CDF curve [16].

Immune cell infiltration analysis

ESTIMATE was utilized to reveal TME in tumor tissues [17]. MCP counter, ssGSEA, cibersort (https://cibersortx.stanford.edu/), EPIC, and TIMER were utilized to reveal the immune cell infiltration [18–20].

Differential gene analysis between the 2 clusters and GO and KEGG enrichment analysis

“Limma” R package was performed to detect putative differences between 2 clusters (|log2FC|>1; adj. p<0.05) [21, 22].

WGCNA

The R language “WGCNA” package was used to construct a gene co-expression network for the normalized gene data and determine the optimal soft threshold. A scale-free network is constructed based on the optimal soft threshold, and the genes are subsequently clustered into the dynamic tree cutting algorithm is used to cluster and classify the genes into different color functional modules. Gene significance (GS) indicates the correlation between genes and traits, while module membership (MM) indicates the correlation between module feature vectors and gene expression profiles. The correlation between modules and clinical traits was analyzed by Pearson algorithm by combining GS and MM, and the module with the highest correlation with clinical traits of HCC was selected as the key module in this study [23].

DRG.score model construction

Based on the DEGs and WGCNA of the 2 clusters, we established a DRG.score by using principal component analysis (PCA) method. The DRG.score of each sample was calculated by

$$\text{DRG}\text{.score}=\sum (P\text{C}{1}_i+P\text{C}{2}_i),$$

where “i” represents the gene expression level of hub gene expression.

GSEA

GSEA judges the enrichment of gene sets based on gene expression. The difference between the two lies in that GSEA judges the enrichment of gene sets based on the contribution of genes in gene sets. To assess the signaling pathways associated with the prognostic model, we used GSEA to assess enrichment pathways in the high-risk and low-risk groups [24].

Immunotherapy efficacy analysis

Use the TIDE web server to predict each sample’s response to immunotherapy based on liver cancer transcriptomic data. Individual TIDE scores were pooled to predict the efficacy of treatment with ICIs in high- and low-DRG.score groups, with higher TIDE scores indicating that tumor cells are more prone to immune escape, implying a lower response rate to ICIs treatment [25].

Statistical analysis

Correlations between variables were explored using Spearman or Pearson coefficients. Continuous variables that conformed to the normal distribution were compared using independent t-tests for comparisons between binary groups, while continuous variables with skewed distributions were compared with the Mann–Whitney U test. Survival curves for categorical variable prognostic analyses were generated using the Kaplan–Meier method, while the log-rank test was used to estimate statistical significance. The significance level was set at P < 0.05, and all statistical tests were two-sided. All statistical data analyses were performed using R software or online analysis tools described in the relevant Materials and Methods subsections.

Data availability

All data used in the study can be downloaded from the TCGA data repository (https://portal.gdc.cancer.gov/repository), ICGC data (https://icgc.org/) and the GEO database (https://www.ncbi.nlm.nih.gov/gds/?term=).

mRNA expression levels and predictive efficiency of 10 DRG in TCGA-LIHC cohort

To confirm the biological function of 10 DRG in HCC, we first calculated the cancer major pathway score and then the correlation of 10 DRG with the cancer major pathway score. We found that 10 DRG mainly affects cancer development through the cell cycle pathway (Figure 1A). Next, the expression of 8/10 DRG were highly expressed in HCC except NUBPL and NDUFS1 (Figure 1B). Meanwhile, the expression of NCKAP1, SLC3A2, GYS1, LRPPRC were highly expressed in stage III and stage IV (Figure 1C). We also found that 10 DRG were closely positive correlation except NDUFA11 (Figure 1D). Next, to assess the accuracy of 10 DRG in predicting the prognosis of HCC, we plotted the ROC curves. The 1-year, 3-year, and 5-year AUC of NCKAP1 was 0.677, 0.605, 0.574, RPN1 was 0.652, 0.631, 0.610, SLC3A2 was 0.656, 0.590, 0.540, SLC7A11 was 0.707, 0.630, 0.574, GYS1 was 0.652, 0.583, 0.520, NDUFS1 was 0.605, 0.565, 0.508, OXSM was 0.659, 0.637, 0.608, LRPPRC was 0.711, 0.626, 0.565, NDUFA11 was 0.465, 0.499, 0.488, NUBPL was 0.548, 0.499, 0.457. The 10 DRG exhibited a low AUC area, suggesting that these gene had a worse overall performance (Figure 1E–1N). Therefore, there is an urgent need to develop a more effective model to predict the survival of HCC.

Figure 1. mRNA expression levels and predictive efficiency of 10 DRG in TCGA-LIHC cohort. (A) The relationship between 10 DRG and tumor signaling pathways. (B) The expression distribution of 10 DRG between tumor and normal. (C) The expression distribution of 10 DRG between stage III & stage IV than stage I & stage II. (D) Correlation map of 10 DRG. (E–N) ROC analysis showed the predict performance 10 DRG.

SLC7A11 and LRPPRC can be used as an independent prognosis factor in HCC in TCGA-LIHC cohort

To assess whether 10 DRG, age, gender, TNM stage, and AFP were independent prognostic factors for HCC patients, our forest plot results obtained by univariate and multivariate COX regression analysis showed that both T stage and SLC7A11 and LRPPRC were independent prognostic factors for patients with HCC (P<0.05) (Figure 2A). Correlation analysis of immune cell subpopulations in ssGSEA showed that SLC7A11 expression levels were positively correlated with T helper cells, macrophages, Th2 cells, NK CD56bright cells, Tcm, Th1 cells and negatively correlated with NK cells, B cells, cytotoxic cells, eosinophils, DC, pDC, and Th17 cells. LRPPRC expression levels were positively correlated with Tcm, T helper cells, and Th2 cells and negatively correlated with other cells (Figure 2B, 2C). Further immunohistochemistry showed that SLC7A11 and LRPPRC were higher in cancer tissues than in normal tissues (Figure 2D). Next, our single-cell dataset was validated for SLC7A11 and LRPPRC expression. Based on the GSE140228 dataset, 62,530 cells were identified. Further clustering analysis resulted in surface these cells can be classified into 12 cell types (B, CD4Tconv, CD8T, CD8Tex, DC, ILC, Mast, Mono/Macro, NK, Plasma, Tprolif, and Treg) (Figure 2E). Among them, LRPPRC is widely released in various immune cells and SLC7A11 is mainly released in DC, ILC cells (Figure 2F, 2G).

Figure 2. SLC7A11 and LRPPRC can be used as an independent prognosis factor in HCC in TCGA-LIHC cohort. (A) The univariate and multivariate Cox regression analysis between 10 DRG and overall survival. (B) The relationship between SLC7A11 and immune infiltration. (C) The relationship between LRPPRC and immune infiltration. (D) Immunohistochemistry images of SLC7A11 and LRPPRC expression in normal tissues, and HCC tissue. (E) T-distributed stochastic neighbor embedding plot of all the single cells, with each color coded for immune cell types. (F, G) The expression distribution of SLC7A11 and LRPPRC in immune cell types.

Molecular subtypes based on 10 DRG in the meta cohort

First, we merge TCGA-LIHC, LIRI-JP, GSE14520, GSE36376, and GSE76427 into a meta cohort. The consensus CDF curve and the change in area under CDF delta area curve showed that, for consensus matrix k=2, 9 DRG related expression-based classification had relatively stable clustering results (Figure 3A). We figured out that the OS time of the DRG.cluster.B had better prognosis compared with DRG.cluster.A (Figure 3B). Two clusters such as DRG.cluster.A and DRG.cluster.B were clearly identified (Figure 3C). A total of 40 mRNAs with significant differences between DRG.cluster.A and DRG.cluster.B HCC subtypes (|log2FC|>1; adj. p<0.05) (Figure 3D and Supplementary Table 1). By univariate COX regression analysis, 39 out of 40 genes were significantly associated with survival prognosis (Figure 3E). These genes were mainly enriched in metabolism-related signaling pathways (retinol metabolism, carbon metabolism, linoleic acid metabolism, tyrosine metabolism) (Figure 3F, 3G and Supplementary Tables 2, 3).

Figure 3. Molecular subtypes based on 10 DRG in the meta cohort. (A) Unsupervised consensus clustering based on 10 DRG for 1155 HCC patients in a meta cohort (GSE14520, GSE36376, GSE76427, LIRI-JP, and TCGA-LIHC). (B) Kaplan-Meier curve showed a significant difference between the 2 DRG.clusters. (C) PCA analysis between 2 DRG.clusters. (D) The different genes between the 2 DRG.clusters. (E) The univariate Cox regression analysis between 40 DEGs and overall survival. (F) GO enrichment analysis, (G) KEGG enrichment analysis for the DEGs and prognosis genes between the 2 DRG.clusters.

TME between the 2 DRG.cluster

Next, we further explored the immunological role of DRG.cluster in HCC. The results suggested that the infiltration level of most immune cells was significantly lower in DRG.cluster.B group (Figure 4A). In this analysis, we also evaluated some representative targets, and found that immune checkpoint were highly expressed in DRG.cluster.A group (Figure 4B and Supplementary Table 4). The ESTIMATE results showed that the stromal score, immune score and ESTIMATE score in DRG.cluster.A group were much higher than those in DRG.cluster.B group (Figure 4C–4F and Supplementary Table 5).

Figure 4. TME between the 2 DRG.cluster. (A) The relationship between the 2 DRG.clusters and TME. (B) The correlation of immune checkpoint condition in 2 DRG.clusters. (C–F) The ESTIMATE score, stromal score, immune score, and tumor immunity levels in 2 DRG.clusters.

Identification of hub genes between the 2 DRG.cluster

The sample clustering tree was constructed based on dynamic hybrid cuts using scale-free networks and topological overlap. Based on the scale-free topology criterion, the optimal soft threshold β=4 was determined based on the fit index and the average degree of network connectivity (Supplementary Figure 1). Based on the optimal soft threshold, the gene modules were divided, the cut height was set to 25%, and a total of 13 modules were obtained, and the module cluster tree was drawn (Figure 5A). The correlation analysis between the construct modules and clinical characteristics was performed and heat maps were drawn, and the grey module correlated most closely with HCC (r=0.31, P=4e-18) (Figure 5B). Next, the genes in the grey module were extracted, and a total of 1097 genes were obtained. A total of 6 overlapping genes were obtained by intersecting the genes in the grey module with DEGs and prognosis through the Venn diagram. These 6 genes were defined as candidate pivotal genes in this study (Figure 5C).

Figure 5. Identification of hub genes between the 2 DRG.cluster. (A) Co-expression network module clustering dendrogram, different colors represent different clusters. (B) Heat map of correlations between gene modules and clinical features of HCC, with red representing positive correlations and blue the opposite. (C) Venn diagram showing common genes between WGCNA and DEGs and prognosis.

DRG.score model was constructed based on 6 hub genes

In order to predict the prognosis of HCC patients more accurately, we built DRG.score based on 6 hub gene (Supplementary Table 6). According to the DRG.score distribution, the 2 DRG.cluster group had significant differences. The results suggested that the DRG.score in DRG.cluster.B group were much lower than those in DRG.cluster.A group (Figure 6A), and DRG.cluster.A group had much more high DRG.score patients (Figure 6B). TCGA-based molecular subtypes also showed that the DRG.scores were different for different subtypes (Figure 6C). The OS time curve showed that lower DRG.score subgroup had longer survival time (Figure 6D). The 1-year, 3-year and 5-year AUC of DRG.score was 0.847, 0.739, 0.685 (Figure 6E). The DRG.score exhibited a high AUC area, suggesting that DRG.score had a good overall performance. Meanwhile, univariate and multivariate Cox regression analysis show that DRG.score can be used as an independent factor to predict the prognosis in HCC (Figure 6F, 6G). Also, a lower mutation was discovered in the low DRG.score group (86.99% vs 83.26%) (Figure 6H, 6I). The above results show that the DRG.score can better predict the risk of prognosis.

Figure 6. DRG.score model was constructed based on 6 hub genes. (A) Differences in DRG.score among 2 DRG.clusters. (B) The number of high and low DRG.score patients in 2 DEG.clusters groups. (C) Differences in DRG.score between immune subtypes. (D) Kaplan-Meier curves for high and low DRG.score groups. (E) The predictive value of DRG.score. (F, G) The univariate and multivariate Cox regression analysis between DRG.score and overall survival. (H, I) The waterfall plot depicted the differences in frequently mutated genes of hepatocellular carcinoma among high and low DRG.score groups.

TME between the DRG.score groups

Next, the results figured out that the positive correlation between these immune cells and DRG.score (Figure 7A). The infiltration level of immune cells was significantly lower in low DRG.score group (Supplementary Figure 2). We also found the positive correlations between these immune checkpoints and DRG.score were intricate (Figure7B). In order to examine the underlying biological processes in the high and low DRG.score groups, we performed a GSEA analysis. The results showed that tumor-related pathways (cell cycle, ECM receptor interaction) were enriched in the high DRG.score group (Figure 4C, 4D). These results suggested that high DRG.score may have a significant impact in tumor development.

Figure 7. TME between the high and low DRG.score groups. (A) The correlation between DRG.score and different immune cells. (B) The correlation between DRG.score and different immune checkpoint. (C) GSEA GO identified high and low DRG.score groups related signaling pathways in HCC. (D) GSEA KEGG identified high and low DRG.score groups related signaling pathways in HCC.

DRG.score is a robust prognosis factor in HCC

To further validate the robustness of DRG.score, the prognostic implication of DRG.score was examined in multiple independent datasets. In both patient sets, patients in the high DRG.score group had a worse OS than in the low DRG.score group in TCGA-LIHC, LIRI-JP, GSE14520, GSE36376, and GSE76427 (Figure 8A–8E). Meanwhile, patients in the high DRG.score group had a worse progression free survival than in the low DRG.score group in TCGA-LIHC (Figure 8F).

Figure 8. External validation of DRG.score. (A–E) The Kaplan-Meier curves analysis between high and low DRG.score groups and overall survival in GSE14520, GSE36376, GSE76427, LIRI-JP, and TCGA-LIHC. (F) The Kaplan-Meier curves analysis between high and low DRG.score groups and progression free survival in G TCGA-LIHC.

Evaluation of the immunotherapy response between the DRG.score groups

The results showed that low DRG.score group had a higher TIDE score compared with high DRG.score group (Figure 9A). Higher TIDE scores indicate a greater likelihood of immune evasion, suggesting that patients may not benefit from immune checkpoint inhibitor therapy. Meanwhile, responder group had a higher DRG.score than non-responder group (Figure 9B), and high DRG.score group had a higher responder patients compared with low DRG.score group (Figure 9C). Moreover, a higher DRG.score was found in the TACE response group than in the TACE non-response group, and high DRG.score group had a higher TACE responder patients compared with low DRG.score group (Figure 9D, 9E). These results suggested that patients with higher DRG.score demonstrated significant immunotherapy therapeutic advantages and TACE clinical benefits.

Figure 9. DRG.score in the role of anti-PD-1/L1 immunotherapy. (A) Differences in TIDE among high and low DRG.score groups. (B) Differences in DRG.score among non-response and response groups. (C) The proportion of non-response and response patients in low or high DRG.score groups. (D) Differences in DRG.score among TACE non-response and TACE response groups. (E) The proportion of TACE non-response and TACE response patients in low or high DRG.score groups.