Copper-Linked Gene Pattern Predicts Outlook, Tumor Environment, and Drug Response in Liver Cancer

To identify the expression differences of CRGs between the normal and tumor samples, the R limma and ggpubr package were used to generate the boxplots for the comparison.

The mutation frequency and mutation type of CRGs were visualized by the R maftools package.

The occurrence of CNV was presented in a barplot, and the distribution of CRGs in chromosomes was depicted by a circle plot with the R RCircos package. The CNV data was divided into three levels, including single deletion, normal, and single gain, and the correlations between these three categorical variables (CNV data and mRNA expression data) were tested.

To specify the prognostic value of single CRG in HCC, R limma, survival, and survminer packages were used for drawing the Kaplan-Meier (K-M) curve and statistical test. The correlations between expression levels of CRGs were calculated and visualized by the R igraph package.

The CRG set prioritization module includes T cell dysfunction, T cell exclusion score, Cox proportional-hazards (Cox-PH) regression score in the immune checkpoint inhibitor (ICB) cohort, and log-fold change (logFC) in CRISPR screens was established with the usage of the Tumor Immune Dysfunction and Exclusion (TIDE) website (http://tide.dfci.harvard.edu/↗). The T dysfunction score and T cell exclusion score indicate the interactions between gene expression level and cytotoxic T cells as well as gene expression level and immunosuppressive cells [21]. The Cox-PH score was associated with ICB therapy response, and the logFC score was employed to sift genes that can interact with lymphocyte-mediated tumor killing [21].

The nonnegative matrix factorization (NMF) approach, an unsupervised clustering analysis method, was performed to identify the presence of underlying molecular subgroups associated with cuproptosis. The R NMF package was employed in this process, and the detailed method used in NMF-clustering analysis was “brunet.” The consensus matrix heat map and cophenetic correlation coefficient were used to determine the cluster number. Principal component analysis (PCA) was adopted to verify the robustness and reliability of the molecular subtypes.

Differential gene expression analysis was initially conducted for the molecular subgroups obtained from the NMF-clustering analysis. Differentially expressed genes (DEGs) among subtypes were screened with a threshold value ∣logFC | = 0.585 and adjusted p value =0.05. Univariate Cox analysis was performed on DEGs and DEGs with a p value <0.05 were retained for subsequent analysis. The least absolute shrinkage and selection operator (LASSO) regression analysis with a repeated time of 1000 was used to improve the model stability. Multivariate Cox analysis was performed based on DEGs after LASSO regression analysis. Finally, a risk score for each patient was calculated using the following formula:

The “βi” and the “Si” in the formula represent the coefficient and expression level of genes, respectively. Patients were divided into high- and low-risk groups according to a cutoff point generated by the Akaike information criterion (AIC). R survival, survminer, and glmnet packages were used in this section.

The differences in risk scores between the molecular subtypes were analyzed with a statistical test using the R ggpubr package. A Sankey diagram was plotted to show the distribution of patients by molecular subtype, risk score, and survival status.

A K-M analysis was conducted to compare the survival curve between the molecular subtype and the risk score group. The differences in the clinical phenotype, including age, gender, histologic grade, and TNM stage, were measured as well. Receiver Operating Characteristic (ROC) curves were used to test the diagnostic performance of the signature. The ROC curves for 1-, 2-, and 3-year were plotted and further compared with the clinical ROC curve. Univariate and multivariate Cox regression analyses were adopted to determine whether the signature could be regarded as an independent prognostic factor. R packages including survival, pheatmap, pbapply, survivalROC, survminer, glmnet, and ggupbr were used for these operations.

In order to recognize the different expression of CRGs between the molecular subtypes and the risk groups, we used the R limma and the ggpubr package to generate the boxplots for the comparison.

GSVA analysis was conducted to identify the biological pathway differences between the molecular subtypes and the risk groups. The KEGG subset file with gene symbol format was downloaded from the MSigDB database (https://www.gsea-msigdb.org/gsea/msigdb/↗). The threshold values used to screen significantly differentially enriched pathways were set as follows: adjusted p value =0.05. R limma, GSEABase, GSVA, and pheatmap package were utilized.

The HCC-related immune checkpoint genes were extracted from Harkus et al.'s study, [6] and the expression levels of them were compared in accordance with diverse molecular subtypes and risk groups. The single sample gene set enrichment analysis (ssGSEA) algorithm was used to evaluate the immune cells and immune functions in tumor immune microenvironment. The infiltration abundance level of 23 types of immune cells and enrichment level of 13 types of immune functions was calculated and compared between the groups. The TIDE score, including dysfunction, exclusion, integrated TIDE score, and microsatellite instability (MSI) score of each HCC sample, were computed on the TIDE website. To evaluate the ratio of immune-stromal components in the tumor microenvironment, the R estimate package was utilized to calculate the StromalScore, ImmuneScore, and ESTIMATEScore. The immune scores were further compared between the groups. Moreover, the correlation between risk score and tumor-infiltrating immune cells, whose infiltration level was estimated by 7 algorithms including TIMER [22], quanTIseq [23], EPiC [24], CIBERSORT [25], xCell [26], MCP-counter [27], and CIBERSORT-ABS [28], was tested. The correlation between risk score and tumor mutation burden (TMB) was also tested.

Gene expression and drug sensitivity data used for drug response prediction were downloaded from the Genomics of Drug Sensitivity in Cancer website (GDSC, https://www.cancerrxgene.org↗). The R package pRRophetic was utilized to calculate the half-maximal inhibitory concentration (IC₅₀) of 4 conventional anti-HCC agents, including sorafenib, cisplatin, paclitaxel, and gemcitabine, for each TCGA-LIHC project sample. The difference in IC₅₀ between high- and low-risk groups was compared.

The LIRI-JP cohort was used to validate the constructed signature. A K-M analysis was conducted to compare the survival between the low-risk and the high-risk patients in the LIRI-JP cohort. A ROC diagram was plotted, and a corresponding area under curve (AUC) value was calculated to evaluate the accuracy of the signature. Univariate and multivariate Cox regression analyses were adopted to determine whether the signature could be regarded as an independent prognostic factor in the LIRI-JP cohort. The ability of the signature to predict ICIs' response was evaluated in the IMvigor210 cohort. The distribution of risk scores and risk status between complete response (CR)/partial response (PR) and stable disease (SD)/progressive disease (PD) was compared.

All statistical operators were conducted with R 4.0.3 software. A log-rank test was used to compare the differences between the survival curves generated by K-M analysis. The comparison and correlation analysis in continuous data were the Wilcox test and the Spearman test. Differences in categorical data were compared using χ² or the Kruskal test. A p value below 0.05 was considered a statistically significant difference.

Two clusters were identified after NMF-clustering analysis (Figure 4(a) and Supplementary Figure 2). The PCA analysis shows that HCC samples are distinctly separated into two clusters (Figure 4(b)). The K-M analysis between the C1 and the C2 suggests that the overall survival of patients in C2 is better than in C1 (Figure 4(c)). CRGs including FDX1, LIAS, DBT, DLST, SLC31A1, and ATP7B are significantly higher expressed in C2, while ATP7A, MTF1, GLS, and CDKN2A are the opposite (Figure 4(d)). According to these findings, the two molecular subtypes behave differently in cuproptosis.

The GSVA analysis shows the activation of KEGG pathways in C1/C2 HCC samples. Metabolism-related pathways such as fatty acid metabolism and drug metabolism, gene stability-related pathways such as DNA replication and mismatch repair, cell cycle, and citrate cycle tricarboxylic acid (TCA) cycle are significantly enriched (Figure 5(a) and Supplementary Table 1). The disorder of TCA cycle caused by excess copper directly result in the occurrence of cuproptosis, and the significantly enriched TCA cycle pathway indeed proves the rationale of the molecular subtype.

Immune checkpoint genes including PDCD1, CD274, CTLA4, TIGIT, LAG3, and LGALS9 are significantly expressed lower in C2 (Figure 5(b)). The results of ssGSEA analysis show that the activated CD4⁺ T cell, CD56^dim natural killer cell, eosinophil, MDSC, type 2 T helper cell, type 1 T helper cell, and type 17 T helper cell are differentially infiltrated between the molecular subtypes (Figure 5(c)). A total of 3 kinds of immune response, including APC costimulation, HLA, and type II IFN response are differentially enriched between the subtypes (Figure 5(c)). The exclusion score is significantly lower in C2, while the dysfunction score and the integrated TIDE score do not present differences between molecular subgroups (Figure 5(d)). Moreover, the StromalScore of C2 is higher than C1, while the ImmuneScore is the opposite (Figure 5(e)). These results suggest that the tumor microenvironment driven by cuproptosis of the two molecular subtypes is diverse from each other.

A 6-gene signature was constructed based on the 1479 DEGs between molecular subtypes with univariate Cox regression analysis and LASSO regression analysis (Figures 6(a)–6(d)). Multivariate Cox regression analysis of the 6 genes shows that MIR1244-2, CBX2, CLEC3B, and SLC16A11 are independent prognostic factors (Figure 6(e)). The cutoff point was set as 1.570 according to the AIC (Figure 6(f)). The AUC of 1-, 2-, and 3-year is 0.775, 0.768, and 0.757, respectively (Figure 6(g)), and the 1-year ROC was further compared with age (AUC = 0.534), gender (AUC = 0.513), histologic grade (AUC = 0.505), and TNM stage (AUC = 0.664) (Figure 6(h)).

HCC samples were divided into high- and low-risk groups in accordance with the AIC value. The overall survival of patients in low-risk group is superior to these in high-risk (Figures 7(a)–7(c)). The univariate Cox regression analysis indicates that the risk status (p < 0.001, HR = 1.478, and 95%CI [1.350–1.619]) and TNM stage (p < 0.001, HR = 1.680, and 95%CI [1.369–2.062]) are the risk factors for HCC patients, and the multivariate Cox regression analysis further proved that the risk score (p < 0.001, HR = 1.418, and 95%CI [1.279–1.572]) acts as an independent role (Figures 7(d) and 7(e)). There is no difference between age and gender in the distribution of risk scores, while increased risk scores are associated with worse histologic grade and more advanced TNM stage (Figures 7(f)–7(i)). CRGs including FDX1, DBT, GCSH, DLST, SLC31A1, and ATP7B are significantly highly expressed in the low-risk group, while LIPT1, DLAT, PDHA1, ATP7A, MTF1, GLS, and CDKN2A are the opposite (Figure 7(j)). Thus, the CRG-related signature represents an eligible AUC value as an independent prognostic factor and it is superior to other clinical phenotypes. The differential expression of CRGs between the two groups suggests that they perform differently on cuproptosis.

The results of GSVA analysis present the differentially enriched KEGG pathways in low-risk/high-risk HCC samples. Immune-related pathways such as T cell receptor signaling pathway and B cell receptor signaling pathway, gene stability-related pathways such as DNA replication and mismatch repair, cell cycle, and citrate cycle TCA cycle are significantly enriched (Figure 7(k) and Supplementary Table 2). It is apparent from these data that the signature inherits the molecular characteristics of the molecular subtype.

The results of ssGSEA analysis show that the activated CD4⁺ T cell, activated dendritic cell, eosinophil, gamma-delta T cell, MDCS, regulatory T cell, type 1 T helper cell, and type 2 T helper cell are differentially infiltrated between the risk groups (Figure 8(a)). Immune functions including APC costimulation, HLA, MHC class I, type I IFN response, and type I IFN response are differentially enriched between the risk groups (Figure 8(a)). Immune checkpoint genes including PDCD1, CTLA4, TIGIT, LAG3, PVR, and LGALS9 are significantly expressed lower in the low-risk group (Figure 8(b)).

Correlations between risk scores and tumor-infiltrating immune cells are presented in Figure 8(c) and Supplementary Table 3. The risk score is strongly and negatively related to the infiltration level of endothelial cell, hematopoietic stem cell, macrophage M2, while positively correlated to T cell CD4+ Th2, neutrophil, macrophage, myeloid dendritic cell, macrophage M1, monocyte, T cell regulatory, and NK cell. The exclusion score is significantly higher in the low-risk group, while the dysfunction score and the integrated TIDE score are the opposite (Figure 8(d)). The StromalScore of the high-risk group was lower than the low-risk group, while ImmuneScore was the opposite (Figure 8(e)). Moreover, the risk score is positively related to the TMB score (Figures 8(f) and 8(g)). The tumor microenvironment driven by cuproptosis of the two risk groups differs from each other, which indicates that the signature can be used to predict the immune status of patients with HCC.

Figures 9(a)–9(d) present the comparison of IC₅₀ for sorafenib, cisplatin, paclitaxel, and gemcitabine between the low- and high-risk groups. The results show that HCC patients in the low-risk group have a higher IC₅₀ in the four comparisons, which indicates that the patients in the high-risk group have a tendency to benefit more from the conventional agents. What stands out in this section is that the signature is able to recognize HCC patients who could benefit from conventional anti-HCC drugs.

As shown in Figure 9(e), the risk score level in C1 is significantly higher than in C2. The distribution of patients by molecular subtype, risk score, and survival status are represented in Figure 9(f).

HCC patients in the LIRI-JP cohort were divided into low- and high-risk groups based on the previously constructed signature. The prognosis of patients in the low-risk group is superior to that of those in the high-risk group (Figure 10(a)). The AUC values of the 1-, 2-, and 3-year ROC curves are 0.811, 0.741, and 0.775, respectively (Figure 10(b)). Tumor patients who positively respond to ICIs have relatively higher risk scores, and tumor patients in the high-risk group are more likely to benefit from ICIs (Figures 10(c) and 10(d)). In summary, for the informants in this section, the signature is robust in predicting the prognosis of HCC and has the potential to foresee the response to immunotherapy.