Signature construction and molecular subtype identification based on cuproptosis-related genes to predict the prognosis and immune activity of patients with hepatocellular carcinoma

Abstract

Background: Hepatocellular carcinoma (HCC) is one of the most common malignancies in the world, with high incidence, high malignancy, and low survival rate. Cuproptosis is a novel form of cell death mediated by lipoylated TCA cycle proteins-mediated novel cell death pathway and is highly associated with mitochondrial metabolism. However, the relationship between the expression level of cuproptosis-related genes (CRGs) and the prognosis of HCC is still unclear.

Methods: Combining the HCC transcriptomic data from The Cancer Genome Atlas(TCGA) and Gene Expression Omnibus (GEO) databases, we identified the differentially expressed cuproptosis-related genes (DECRGs) and obtained the prognosis-related DECRGs through univariate regression analysis.LASSO and multivariate COX regression analyses of these DECRGs yielded four genes that were used to construct the signature. Next, we use ROC curves to evaluate the performance of signatures. The tumor microenvironment, immune infiltration, tumor mutation load, half-maximum suppression concentration, and immunotherapy effects were also compared between the low-risk and high-risk groups. Finally, we analyzed the expression level, prognosis, and immune infiltration correlation on the four genes that constructed the model.

Results: Four DECRGs s were used to construct the signature. The ROC curves indicated that signature can better assess the prognosis of HCC patients. Patients were grouped according to the signature risk score. Patients in the low-risk group had a significantly longer survival time than those in the high-risk group. Furthermore, the tumor mutation burden (TMB) values were associated with the risk score and the higher-risk group had a higher proportion of TP53 mutations than the low-risk group.ESTIMATE analysis showed significant differences in stromal scores between the two groups.N6-methyladenosine (m6A) and multiple immune checkpoints were expressed at higher levels in the high-risk group. Then, we found that signature score correlated with chemotherapeutic drug sensitivity and immunotherapy efficacy in HCC patients. Finally, we further confirmed that the four DECRGs genes were associated with the prognosis of HCC through external validation.

Conclusions: We studied from the cuproptosis perspective and developed a new prognostic feature to predict the prognosis of HCC patients. This signature with good performance will help physicians to evaluate the overall prognosis of patients and may provide new ideas for clinical decision-making and treatment strategies.

Keywords: cuproptosis; hepatocellular carcinoma; immune infiltration; immune microenvironment; prognostic signature.

Figure 1

Flow chart of this study.

Figure 2

Expression and genetic alteration of CRGs in HCC. (A) the expression of 19 CRGs in HCC and normal tissues. (B) correlations between the expression of CRGs; (C) Protein-protein interaction (PPI) networks between CRGs; (D–F) the CNV and mutation frequency and classification of 19 CRGs in HCC. *p<0.05, wfi 2***p<0.001; ns, not statistically different.

Figure 3

Prognosis significance of CRGs of HCC patients in TCGA and GEO in HCC. (A–O) K-M survival curve displays the OS of HCC patients. (P) Prognostic network of CRGs.

Figure 4

Clustering analyses of the signature. (A, B) Concordance matrix and K-M survival curve of the three clusters. (C) Complex heat maps show clinical correlations among the three clusters. (D–F) The GSVA heat map showed the differences in pathways in the three clusters. (G) The differential analyses between immune cells and the scale of fraction for cluster A and cluster B and cluster C. *p<0.05, **p<0.01, ***p<0.001.

Figure 5

Functional enrichment analysis of DECRGs. (A) Principal component analysis of three clusters. (B) The Venn diagram shows the intersection of DECRGs (C) Analysis of BP, CC, and MF terms of GO enrichment demonstrated the possible function of the DECRGs. (D) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis revealed the possible pathways.

Figure 6

Clustering analyses of the signature. (A) The cumulative distribution function based on the sign indicated that the optimal number of subtypes was 3. (B) Concordance matrix of subtypes. (C) K-M survival curve of the three clusters. (D) A complex heat map illustrated the expression patterns. (E) Expression of CRGs between cluster A, cluster B, and cluster C (H) Ggalluvial shows the construction of the prognostic model. (F, G) LASSO regression analyses for screening LASSOSigGenes. Boxplots indicate the differences in risk scores in the CRG cluster (I) and the gene cluster (J).The differential analysis of CRGs expression (K). *p<0.05, **p<0.01, ***p<0.001.

Figure 7

Validation of the prognostic value of the signatures. (A–C)The KM curve of all sets, testing set, and training set. (D) Univariate Cox regression analysis and multivariate Cox regression (H) analysis risk score and clinical stage as independent predictors of OS. (E–G) ROCs for 1-year, 2-year, and 3-year OS prediction. (I) The nomogram of the risk score and clinical parameters (age, gender, and stage) of all sets. (J) The calibration curves displayed the accuracy of the nomogram in the 1st, 2nd, and 3rd years. (K, L) Comparison of the OS between the high-risk and low-risk groups of patients who are<= 65 or > 65 years, male (M) or female (N) with a stage of stage I-II (O) or stage III-IV (P).

Figure 8

Analysis of gene expression and mutation correlation. (A–C) The risk curve consists of genes expression heat map, risk score curves, and survival status point plot. (D–F) Comparison of the proportion of the TP53 mutation status in the high- and low-risk groups in the training, testing, and all sets.

Figure 9

To assess the tumor microenvironment, immune checkpoint genes, and tumor mutation burden (TMB) in different groups. (A) Comparison of ESTIMATE scores, stromal scores, and immune scores between the high-risk and low-risk groups. (B) Correlation between the model-constructed genes and immune cells. (C) Differential expression of immune checkpoint genes between the high-risk group and the low-risk groups. (D) Correlation between the stem cell content and the risk score. (E) Differential expression of m6A-related genes (F, G) and the frequency of mutations in the high-risk and low-risk groups. (H, I) The KM curve of the tumor mutation burden versus the OS. *p<0.05, **p<0.01, ***p<0.001; ns, not statistically different.

Figure 10

Signature predicts chemotherapy and immunotherapy response. (A–E) The signature showed high-risk scores were associated with a lower IC50 for chemotherapeutics such as (A) Axitinib, (B) Imatinib, (C) Lapatinib, (D) Gefitinib, and (E) Bicalutamide, whereas they were related to a higher IC50 for (F) Cisplatin, (G) Gemcitabine, (H) Doxorubicin, (I) Etoposide, (J) Rapamycin and (K) Nilotinib treatment. (L) Differences in TIDE score between high- and low-risk groups.

Figure 11

Differential expression of the signature genes. (A–D) Differential expression of CFHR4, DNASE1L3, SPP2, and TAF6 in adjacent non-tumor tissues and HCC tissues in the TCGA database. (E–H) Differential expression of CFHR4, DNASE1L3, SPP2, and TAF6 in adjacent non-tumor tissues and HCC tissues in the CPTAC dataset. (I–L) Differential expression of CFHR4, DNASE1L3, SPP2, and TAF6 of HCC patients in the wild TP53 group and the mutant TP53 groups. ***p<0.001.

Figure 12

Verify the mRNA expression levels of the four signature genes in the tissues. The qRT-PCR results showed thatDNASE1L3 (A), CFHR4 (B), and SPP2 (C) were expressed higher in non-tumor tissues (Non-Tumor) than in tumor tissues (Tumor), TAF6 (D) was expressed higher in tumor tissues (Tumor)than in non-tumor tissues (Non-Tumor). TAF6 representative IHC (E, F) stained images in HCC tissue and adjacent tissues(n = 16; magnification: left, 100×; right, 200×). Expression levels of TAF6 (G) in HCC cell lines. The bar graphs (H) indicates better silencing for si-TAF6#1 and si-TAF6#2. *p< 0.05, **p< 0.01, and ***p< 0.001.

Figure 13

Adverse effects of TAF6 on HCC in vitro. (A, D) Compared with the control group, the proliferation rate of HCCLM3 cells was significantly inhibited after TAF6 silencing by EdU staining. (B, E) Transwell experiments showed that the migratory ability of HCCLM3 was inhibited after TAF6 silencing. (C) After TAF6 silencing, the cell viability of HCCLM3 was significantly inhibited by the CCK- 8 assay. (F, G) The wound healing array showed that LTAF6-downregulated HCCLM3 cells exhibited significantly delayed wound healing compared with controls. (H–J) Effects of with or without inhibition of TAF6 expression on PD-1 protein expression levels by western blotting. Scale bar: EdU,50μm; Transwell experiments and Wound healing array,200μm. **p<0.01, ***p<0.001.