Cuproptosis regulator-mediated patterns associated with immune infiltration features and construction of cuproptosis-related signatures to guide immunotherapy

Background: Liver hepatocellular carcinoma (HCC) is a prevalent cancer that lacks a sufficiently efficient approach to guide immunotherapy. Additionally, cuproptosis is a recently identified regulated cell death program that is triggered by copper ionophores. However, its possible significance in tumor immune cell infiltration is still unclear.

Methods: Cuproptosis subtypes in HCC were identified using unsupervised consensus cluster analysis based on 10 cuproptosis regulators expressions, and a cuproptosis-related risk signature was generated using univariate and LASSO Cox regression and validated using the ICGC data. Moreover, the relationship between signature and tumor immune microenvironment (TME) was studied through tumor immunotherapy responsiveness, immune cell infiltration, and tumor stem cell analysis. Finally, clinical specimens were analyzed using immunohistochemistry to verify the expression of the three genes in the signature.

Results: Two subtypes of cuproptosis regulation were observed in HCC, with different immune cell infiltration features. Genes expressed differentially between the two cuproptosis clusters in the TCGA were determined and used to construct a risk signature that was validated using the ICGC cohort. Greater immune and stromal cell infiltration were observed in the high-risk group and were associated with unfavorable prognosis. Elevated risk scores were linked with higher RNA stemness scores (RNAss) and tumor mutational burden (TMB), together with a greater likelihood of benefitting from immunotherapy.

Conclusion: It was found that cuproptosis regulatory patterns may play important roles in the heterogeneity of immune cell infiltration. The risk signature associated with cuproptosis can assess each patient’s risk score, leading to more individualized and effective immunotherapy.

Introduction

Liver cancer is an aggressive tumor with poor outcomes (1). Liver hepatocellular carcinoma (HCC) accounts for more than 80% of all primary liver malignancies for which surgical resection is currently the most effective treatment (2). Unfortunately, diagnosis of HCC is frequently delayed leading to unfavorable outcomes (3). Despite recent significant progress in immunotherapy and targeted therapy, the five-year survival rate remains low, at approximately 15% (4). An important reason for this is that patients with HCC vary in their response to immunotherapy, and the factors that influence and predict the response to immunotherapy in HCC remain unclear (5).

Recent research has shown that copper toxicity-mediated cell death differs from other forms of regulated cell death; this novel mechanism is termed cuproptosis (6). Copper is both necessary and potentially toxic for cells. It is an essential cofactor required by all living organisms to function properly (7–9); however, high levels of copper accumulation or improper distribution in the cell can lead to cell death. Imbalances in copper homeostasis in cells can lead to severe disease in humans, including tumor development (10, 11). Excess copper has been linked with various types of cancer, including breast (12–14), prostate (15–17), colon (18), lung (19), brain (20), and liver (21) cancer. However, the reasons underlying elevated copper levels in tumors are unclear. Recent studies revealed the mechanism associated with copper-mediated cell death. A study by Tsvetkov et al. showed that cuproptosis is associated with copper binding to fatty acylated moieties of tricarboxylic acid, resulting in the abnormal aggregation of fatty acylated proteins and the loss of iron-sulfur cluster proteins, leading to proteotoxic stress and cell death (6). These findings suggest a starting point for investigating the application of cuproptosis in disease treatment, especially, for tumor therapy.

Copper, as a key trace element, is necessary for the functioning of the immune system. Copper deficiency adversely affects immune function and exposes the organism to microbial infection (22). The immune system requires copper for a variety of functions. The metal can modulate the activation of cells associated with innate immunity such as macrophages and neutrophils during bacterial infection and leukocyte differentiation, maturation, and migration (23, 24). Copper deficiency may also affect immune cell distribution in tissues or the maturation pattern of leukocyte populations (25). In addition, intratumoral copper can regulate PD-L1 expression and affect tumor immune escape (26). These reports suggest an association between copper and the modulation of the immune response.

Here, we used data from the TCGA to investigate cuproptosis subtypes and correlate these subtypes with immune-infiltration features. Signature genes associated with cuproptosis were identified by differential and prognostic analysis, and a risk signature using these genes was generated. This cuproptosis-related risk signature identifies specific features of immune infiltration, permitting the development of individualized immunotherapy.

Methods

Collection of data

Data on gene expression (fragments per kilobase, FPKM), mutations, and clinical data were obtained from the Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). Patients lacking survival information were not enrolled. Xena (http://xena.ucsc.edu/) was utilized to obtain copy number variation (CNV) data and RNA stemness scores (RNAss) for TCGA samples. Validation was performed using data available from the International Cancer Genome Consortium (ICGC) (https://dcc.ICGC.org/projects/LIRI-JP) for 229 HCC patients with complete information. The 10 cuproptosis regulators identified by Tsvetkov et al. (6) were used (Table S1). Since TCGA and ICGC data are publicly available, permission from an ethics committee was not necessary for this work. Despite this, the investigation was performed in conformity with the procedure specifications established by the TCGA and the ICGC.

Profiling of cuproptosis regulators

Principal component analysis (PCA) using the expression of cuproptosis regulators was performed on both HCC and normal samples with the package “scatterplot3d” in R, with relationships between the regulators examined using “corrplot”, and the CNVs of 10 regulators on human chromosomes using “Rcircos”.

Consensus cluster analysis of cuproptosis regulators

The “ConsensusClusterPlus” package in R was employed for unsupervised consensus cluster analysis (27). The criteria for clustering were an initial gradual increase in the cumulative distribution function (CDF) curve, strong intraclass relationships and weak interclass relationships, and no unacceptably small sample sizes in any of the groups. The R packages “ggplot2” and “Rtsne” were used for PCA and t-distributed random neighbor embedding (t-SNE) analysis, respectively. Overall survival (OS) in the different clusters was analyzed with Kaplan-Meier (K-M) curves. Enrichment analysis of genes was conducted with gene set enrichment analysis (GSEA, v4.2.1; https://www.gsea-msigdb.org/gsea/downloads.jsp) (28) and Gene Ontology (GO) using the R package “clusterProfiler” (29).

Immune infiltration profiles in the different groups

The relative abundance levels of immune cell types in the HCC microenvironment were examined by single-sample GSEA (ssGSEA) (30, 31). Immune cell genomes obtained from Charoentong et al. (32) were used as markers for the various cell types. Tumor immune dysfunction and rejection (TIDE; http://tide.dfci.harvard.edu/) was also used as a marker for examining immune escape mechanisms and predicting the response to immunotherapy (33). Higher TIDE scores indicate a greater likelihood of immune escape, leading to the probability of reduced treatment efficacy. The degree of infiltration of six immune cell types associated with the expression of three differentially expressed genes (DEGs) was determined with the Tumor Immune Estimation Resource (TIMER) database (https://cistrome.shinyapps.io/timer/) (34, 35). The “SCNA” module was utilized to compare infiltration levels between tumors with different CNVs for specific genes.

Comparison of responsiveness to immunotherapy

The IMvigor 210 cohort comprising cancer patients treated with programmed death ligand 1 (PD-L1) inhibitors was used for the prediction of immunotherapy response (36). These samples were categorized based on the patient’s response, namely partial response (PR), complete response (CR), progressive disease (PD), and stable disease (SD), with PR and CR representing immunotherapy response and PD and SD representing a lack of response. The R package “pRRophetic” was utilized to compute the half-maximal inhibitory concentration (IC50) of the drugs with low IC50 values indicating greater drug sensitivity.

Construction and validation of the HCC risk signature

Two patient clusters were formed based on the expression of the 10 cuproptosis-related regulatory genes. Genes expressed differentially between the clusters were determined with the empirical Bayesian method in the “limma” package with significance set as an adjusted P-value < 0.0001. Further DEGs between HCC samples and adjacent normal tissue were analyzed to narrow down candidate signature genes (P<0.0001). Finally, univariate Cox regression as well as least absolute shrinkage and selection operator (LASSO) analysis (37) were employed to detect DEGs in the final signature. Risk scores were calculated for all HCC samples and groups with low-risk and high-risk were identified in the TCGA cohort based on median risk scores. The risk score was calculated as:

$$\text{Risk Score }=\sum _{i=0}^n\beta i*Gi$$

where βi denotes the gene coefficient; i and n are the number of genes in the signature; Gi denotes the gene expression value.

The risk scores for samples in the ICGC set were then determined using the risk coefficients of the signature DEGs in the TCGA, and the ICGC samples were allocated two groups (low risk and high risk) using the risk thresholds observed in the TCGA. K-M survival curves were utilized to evaluate OS between the groups. Risk signature sensitivity and specificity were calculated using receiver operating characteristic (ROC) curves.

Nomogram construction and verification

In both the TCGA and ICGC cohorts, hazard ratio models were developed using univariate and multivariate Cox regression to detect independent prognostic variables. A nomogram including risk scores and other clinicopathological features was then constructed in the TCGA cohort and one-, three-, and five-year calibration curves for constructed for accuracy validation. The area under the ROC curve (AUC) and decision curve analysis (DCA) were used to evaluate the nomogram’s discriminative power.

Determination of the relationships between tumor mutation and risk signature

The R package “maftools” was used to generate the waterfall charts for the low- and high-risk groups respectively. Somatic mutation analysis was used to determine the tumor mutational burden (TMB) score of each TCGA-HCC patient. Subsequently, the Pearson correlation analysis was studied to determine the link between risk score and TMB. The data was further expressed in scatterplots and boxplots. In addition, boxplots were generated to demonstrate the differences in immune cell infiltration between the low- and high-TMB groups, respectively. The overall survival (OS) rate among these groups were determined through K-M survival curves.

Comparison of immunohistochemistry in normal and HCC samples

Immunohistochemical staining was performed with antibodies against human G6PD (25413-1-AP; proteintech), CDCA8 (bs-7834R; Boaosen), and Cyclin B1 (bs-0572R; Boaosen), followed by a pathologist based on the percentage of positive cells and staining intensity to score sections. Staining intensity was scored as 0 (negative), 1 (weak), 2 (moderate), or 3 (strong), and the proportion of positive cells expressed was scored as 1 (0-25%), 2 (26-50%), 3 (51– 75%) or 4 (76–100%). The final score was obtained by multiplying the expression ratio and the signal intensity.

Expression and prognostic analysis of three DEGs

The expression of three genes differentially expressed between normal tissues and HCC was analyzed using the “limma” R package, and the “surv cutpoint” function in “survminer” was employed to examine the optimal cut-off expression values for survival. Groups with high and low expression were established and OS was compared between groups using K-M analysis.

Statistical analysis

All statistical analysis was conducted using R version 4.1.2 (http://www.R-project.org) and its accompanying packages. The Kruskal-Wallis test and one-way ANOVA were used to examine differences between three or more groups. OS comparisons were conducted using K-M analysis and the log-rank test. Univariate and multivariate Cox regression analyses were used to calculate hazard ratios and identify independent risk factors. The diagnostic value of risk scores and nomogram models was assessed using ROC curves. P<0.05 was considered to be the criterion for significant statistical difference.

Results

The landscape of cuproptosis regulators in HCC

Based on the study by Tsvetkov et al. (6),10 cuproptosis regulators (7 positive regulators and 3 negative regulators) were included in this study. The expression of these regulators differed significantly between normal and HCC tissues. Apart from the critical regulator FDX1, which was down-regulated in HCC, all other regulators were up-regulated in tumor tissue in comparison with controls (Figure 1A). PCA analysis indicated that the expression of all 10 regulatory genes clearly distinguished between HCC and normal tissue (Figure 1B). After the division of the genes into high and low expression groups using “survminer”, it was observed that low levels of FDX1 were linked with poor HCC prognosis, while the remaining genes exhibited an opposite trend (P>0.05 for DLD due to insufficient sample) (Figure S1). Analysis of co-expression of the regulators indicated that all had positive regulatory relationships apart from FDX1 where a negative relationship was observed (Figure 1C). The highest correlation (0.47) was seen between DLD and PHAR1. The regulatory network demonstrated both the interactions between the genes and their potential prognostic significance for HCC (Figure 1D). The inconsistent findings for FDX1 suggest that it may be more critical in the regulation of cuproptosis (6). Further analysis of the genetic basis of cuproptosis using the TMBs and CNVs of the 10 regulators showed alterations in 5.85% of the 364 samples (16 mutations). The highest number of mutations was seen in CDKN2 followed by DLD and MTF1 (Figure 1E). CNVs were common in the 10 genes, with LOSS occurring more frequently than GAIN. LIAS and GLS showed higher CNV gains, while CDKN2 exhibited higher loss (Figure 1F). Figure 1G shows the chromosomal location of the regulators’ CNVs. These findings marked different levels of these regulatory factors in tumor and normal tissues, implicating disturbances in cuproptosis regulator expression in HCC.

FIGURE 1

Figure 1 Expression of cuproptosis-regulatory genes in HCC. (A) Comparison of expression between normal and tumor tissues (B) PCA of the expression of the 10 genes. (C) Spearman correlations between genes. The color blue represents negative regulation; the color red represents positive regulation. (D) Comprehensive network map combining cuproptosis regulator interactions and prognosis. (E) Mutation frequencies of the 10 regulatory genes in 364 HCC s specimens. (F) CNV values for cuproptosis regulators in HCC specimens. (G) Chromosomal locations of CNV alterations in cuproptosis regulators. *P < 0.05; **P < 0.01; ***P < 0.001. HCC, hepatocellular carcinoma; PCA, principal component analysis; CNV, copy number variation.

Correlations between cuproptosis clusters and TME features

Patterns of cuproptosis were determined by unsupervised consensus cluster analysis of cuproptosis regulator expression in HCC and the results of the consensus clustering (k2–9) were visualized with a CDF plot (Figures 2A; S2). Examination of the consensus matrix showed that k=2 was the best option, and that each sample in the cluster exhibited a strong correlation (Figure 2A). The HCC patients were, therefore, assigned to one of two clusters (A or B). The clustering results are shown in Table S2. PCA (Figure 2B) and t-SNE (Figure 2C) analysis showed that there were two discrete directions within the clusters. Furthermore, we found that the prognosis of these patients was significantly different, with patients in cluster A showing better OS (Figure 2D). Multi-GSEA enrichment analysis indicated significant differences in biological processes between the clusters, specifically, in metabolism, tumor signaling pathways, and immune- and matrix-related pathways (Figure 2E). The enriched metabolic processes included those associated with fatty acid, linoleic acid, retinol, and drug metabolism cytochrome P450 showing greater enrichment in cluster A. Cluster B showed greater enrichment in the tumor-associated ERBB, MAPK, TGFβ, VEGF, and Wnt signaling pathways. Pathways related to immune activation were associated mostly with cluster B; these included leukocyte transendothelial migration, Toll-like receptor signaling, and NOD-like receptor signaling. We then investigated infiltrating immune cells in the clusters (Figure 2F), finding that the majority of immune cells were strongly associated with cluster B. At the same time, we also found that matrix-activated pathways such as focal adhesions, gap junctions, and tight junctions were enriched in cluster B (Figure 2E). It has been reported that tumor stromal cells are involved in immune regulation (38). Stromal cells are able to block immune cell entry into the tumor parenchyma and can also block T-cell killing of tumor cells (39). Tumor cells and stromal cells can induce angiogenesis, thereby promoting tumor metastasis (40).. In addition, regulatory T cells and CD 56bright NK cells highly expressed in cluster B can act as immunosuppressive cells to promote tumorigenic immune escape (41, 42). The results of the immune microenvironment analysis described above validate the poor prognosis of cluster B patients. Taken together, we found that the expression of cuproptosis regulatory proteins in HCC defines two patient clusters with significantly different immune infiltration features. Cluster A is an immune-desert phenotype with reduced immune activity and immune cell numbers while cluster B is an immune-exclusion phenotype with infiltration restricted to the peripheral matrix of tumor cells.

FIGURE 2

Figure 2 Identification of cuproptosis clusters. (A) Consensus matrix based on cuproptosis regulator expression in the HCC cohort at k = 2. (B) PCA of consensus matrix when k = 2. (C) T-SNE analysis of consensus matrix when k = 2. (D) K-M survival analysis between the two cuproptosis clusters. (E) Multiple GSEA analysis between the two cuproptosis clusters. (F) Infiltration of immune cells between the two cuproptosis clusters. *P < 0.05; **P < 0.01; ***P < 0.001. HCC, hepatocellular carcinoma; PCA, principal component analysis; T-SNE, t-distributed random neighbor embedding; K-M, Kaplan-Meier; GSEA, gene set enrichment analysis.

Characterization of cuproptosis-related phenotypes

We next identified 629 cuproptosis-associated DEGs using the R package “limma”. GO analysis showed enrichment of the DEGs in immune- and stromal-associated pathways (including neutrophil mediated immunity, neutrophil activation involved in immune response, focal adhesion, and cell−substrate junctions), confirming that cuproptosis is closely associated with the regulation of the immune microenvironment of the tumor (Figure 3A). We then further narrowed the range by comparing gene expression between the normal and HCC groups to obtain 30 DEGs. Univariate Cox regression yielded 27 prognosis-associated DEGs with HCC patients again divided into two clusters according to the unsupervised consensus cluster analysis of the expression of these genes (Figures 3B; S3). The gene clustering results are presented in Table S3. Differential expression of most of the cuproptosis regulators was visible between the clusters (Figure 3C). Additionally, RNAss values were significantly increased in gene cluster B compared with gene cluster A (Figure 3D), demonstrating the adverse effects of high RNAss values on HCC prognosis (Figure 3E). Oncogenic and immune- and matrix-activated pathways were enriched in cluster B, while cluster A showed enrichment in metabolic pathways (Figure 3F). These findings suggest the presence of two cuproptosis-associated regulatory patterns in HCC.

FIGURE 3

Figure 3 Gene cluster determination. (A) GO enrichment analysis of cuproptosis-related DEGs. (B) Consensus matrix based on cuproptosis-related gene expression in the HCC cohort at k = 2. (C) Differential expression of cuproptosis regulators in the two clusters. (D) RNAss values of cuproptosis regulators in the two clusters. (E) K-M survival analysis of cuproptosis regulators in the two clusters. (F) Multiple GSEA analysis between the two clusters. **P < 0.01; ***P < 0.001. GO, Gene Ontology; DEGs, differentially expressed genes; HCC, liver hepatocellular carcinoma; RNAss, RNA stemness scores; K-M, Kaplan-Meier; GSEA, gene set enrichment analysis.

Development of the cuproptosis-associated risk signature

The above analysis was performed on a patient population. We then investigated the precise quantification of cuproptosis patterns in individual patients. Using the 27 genes identified in the univariate analysis (Figure 4A), we performed LASSO to prevent signature overfitting (Figures 4B, C). The risk signature was finally constructed using three signature genes (CDCA8, CCNB1, and G6PD). The risk score assigned to each sample was determined as: Risk Score = (0.066092028117448) × the CDCA8 expression + (0.00250484775117196) × the CCNB1 expression + (0.00581254249797367) × the G6PD expression. This enabled the HCC patients to be entirely separated into groups with low and high risk according to the median risk score. The risk score results are provided in Table S4. PCA showed the clear separation of the 343 HCC patients and 279 DEGs in two independent clusters, trending in two different directions (Figures 4D, E). The scatter plots of risk scores and patient survival statistics are presented in Figures 4F, G. As depicted in the figure, a higher risk score was related to both decreased survival and greater mortality. The heatmap (Figure 4H) shows the expression of the three DEGs in the signature in the TCGA cohort (Figure 4H). K-M curves indicated higher OS rates in the low-risk group (Figure 4I). To examine the performance of signature, ROC curves for one-, two-, and three-year OS were generated, with AUCs of 0.783, 0.725, and 0.686, respectively (Figure 4J). The risk signature’s superiority was further established by comparing the one-year ROC curves with other clinicopathological features (Figure 4K). Univariate Cox regression was employed to examine associations, finding that stage (HR = 1.804, 95% CI = 1.456–2.234, p < 0.001) and risk score (HR = 3.935, 95% CI =2.740–5.649, p < 0.001) were related to OS in the TCGA set (Figure 4L). Multivariate regression verified that stage (HR = 1.609, 95% CI = 1.274–2.031, p < 0.001) and risk score (HR = 3.236, 95% CI = 2.176–4.811, p < 0.001) were independent prognostic factors for OS (Figure 4M).

FIGURE 4

Figure 4 Construction of the risk signature. (A) Forest plot of univariate Cox regression results for the 27 DEGs. (B) Cross-validation for selection of tuning parameters in LASSO regression. (C) LASSO coefficient profiles of candidate DEGs. (D) PCA of samples from high- and low-risk groups. (E) PCA of DEGs in the high- and low-risk groups. (F) Risk score distribution of HCC patients in the TCGA cohort. (G) Scatter plot of survival status of HCC patients in the TCGA cohort. (H) Heatmap of the expression of the three signature genes in high- and low-risk populations in the TCGA cohort. (I) K-M curves of OS in HCC patients in the risk score-based TCGA cohort. (J) ROC curves of prognostic signatures for one, two, and three years in the TCGA cohort. (K) ROC curves of prognostic signatures and other clinicopathological features for one year in the TCGA cohort. (L, M) Forest plots of univariate and multivariate Cox regression analyses of prognostic signatures and clinical features in the TCGA cohort. DEGs, differentially expressed genes; LASSO, least absolute shrinkage and selection operator; PCA, principal component analysis; HCC, hepatocellular carcinoma; TCGA, the Cancer Genome Atlas; K-M, Kaplan-Meier; OS, overall survival; ROC, receiver operating characteristic.

Validation of the risk signature

The test group samples from the ICGC were divided into low- (n = 104) and high- (n = 125) risk score groups using the cutoff values determined in the TCGA cohort. Figures 5A, B show the risk curves and scatterplots while the heatmap in Figure 5C illustrates the expression patterns of the three signature DEGs in the ICGC set. The high-risk group, as determined in the TCGA cohort, had a worse overall prognosis (Figure 5D). The AUCs of the ROC curves for one, two, and three years were 0.757, 0.760, and 0.785, respectively (Figure 5E), and the AUC of the risk score at one year was larger than those for other clinical parameters (Figure 5F). Both univariate and multivariate regression analyses found that the P-values for both risk score and sex were below 0.05, demonstrating that signature functioned as an independent predictor of HCC prognosis in the ICGC set, and further confirming its reliability in the TCGA set (Figures 5G, H).

FIGURE 5

Figure 5 Independent verification of risk signatures and construction of nomograms. (A) Distribution of HCC patient risk scores in the ICGC cohort. (B) Scatterplot showing HCC patient OS in the ICGC cohort. (C) Heatmap of the expression of the three signature genes in high- and low-risk groups in the ICGC cohort. (D) K-M curves for HCC patient OS in the risk score-based ICGC cohort. (E) ROC curves of one-, two, and three-year prognostic signatures in the ICGC cohort. (F) ROC curves of one-year prognostic signatures with other clinicopathological characteristics in the ICGC cohort. (G-H) Forest plots of univariate and multivariate Cox regression analyses of prognostic signatures and other clinicopathological features in the ICGC cohort. (I) Nomogram combining risk signatures and clinical stages for OS prediction in the TCGA cohort. (J) Calibration plots for nomograms at one, three, and five years. (K) DCA of the nomogram, risk signature, and other clinicopathological features. (L) ROC curves of the nomogram, risk signature, and clinicopathological characteristics at one year. HCC, hepatocellular carcinoma; ICGC, International Cancer Genome Consortium; K-M, Kaplan-Meier; OS, overall survival; ROC, receiver operating characteristic; DCA, decision curve analysis.

Nomogram creation and validation

To further improve the clinical utility of risk signature, the identified independent prognostic risk factors, namely, risk score and stage, were utilized to develop a nomogram for predicting OS in the TCGA cohort (Figure 5I). The calibration curves at one, three, and five years showed that the nomogram accurately predicted the outcomes of HCC patients (Figure 5J). In addition, the DCA curve showed that the nomogram had specific advantages over other clinical parameters (Figure 5K). The AUCs for the ROC curves suggested that the nomogram, with an AUC of 0.804, was superior to age (AUC=0.493), sex (AUC=0.508), pathological grade (AUC=0.490), TNM stage (AUC=0.713), and risk signature (AUC=0.788) for predicting prognosis in HCC patients (Figure 5L).

Relationships between the risk signature and immune characteristics, clinical parameters

A GSEA analysis was next used to investigate physiology-related differences between the groups (Figure 6A). The results were consistent with those of the clustering analysis, showing that many metabolism-related pathways, including those associated with fatty acid, cytochrome P450, and retinol metabolism were enriched in the low-risk group. Tumor-associated pathways (including the MAPK, NOTCH, and VEGF pathways) and immune- and matrix-related pathways (including the Toll-like receptor, B cell receptor, T cell receptor, chemotaxis factor pathways, and attachment junctions) were associated with high risk. To confirm the immune characteristics of the signature, we examined associations between immune cells and risk scores. This indicated a positive relationship between risk scores and most immune cells (Figure 6B). Taken together, it appears that high-risk scores are indeed associated with increased stromal activity and immune infiltration. Interestingly, the TIDE analysis indicated that immune escape was likely to be reduced in the high-risk group after immunotherapy (Figure 6C). It was also found that immune checkpoint levels, including those of PDCD1 and CTLA4, were raised in the high-risk group, which thus had a greater chance of benefitting from immunotherapy (Figure 6D). Calculation of the risk scores of patients in the IMvigor210 set using the TCGA set risk signature indicated that patients responding (CR/PR) to anti-PD1/PD-L1 therapy showed significantly raised risk scores (Figure 6E). We then compared HCC samples for differences in risk scores based on tumor immunophenotyping established by Thorsson et al. (43). The immunophenotyping of the HCC samples is presented in Table S5. The Wound Healing (immune C1) risk score was found to be higher than others, suggesting that high-risk scores were associated with matrix activation (Figure 6F). In addition, analysis of tumor stemness indicated that RNAss was raised in the high-risk group (Figure 6G) with a positive association between the risk scores and RNAss (Figure 6H). The RNAss values of the HCC samples are presented in Table S6. Deceased HCC patients had higher RNAss values than surviving patients (Figure 6I). The OS values in the high- and low-RNAss groups were 61% and 73%, respectively (Figure 6J), indicating both the prognostic significance of RNAss and the reliability of the risk score measure. Furthermore, the associations between clinical features and risk scores are shown in Figures S4A–F. Risk scores were found to increase significantly with TNM stage from I to III, T stage from I to III, grade from G1 to G4, and AFP expression from low to high. The risk scores did not change significantly in response to age or sex. In addition, boxplots were used to show the results for 12 drug sensitivities by estimating IC50 values between the groups. Patient groups at low risk were significantly more sensitive to Gefitinib, Sorafenib, Nilotinib, Dasatinib, Erlotinib, and Metformin (Figure S5A–F), in contrast to high-risk patients who responded to Bleomycin, Doxorubicin, Gemcitabine, Tipifarnib, Imatinib, and Mitomycin.C (Figure S5G–K). This has reference significance for guiding the clinical medication of HCC treatment.

FIGURE 6

Figure 6 The correlation between the immunity and risk signature. (A) Multiplex GSEA analysis between groups with high and low risk (B) Correlations between infiltration of immune cell levels and risk scores. The color blue implies a negative association; the color red denotes a positive association. (C) TIDE scores in the groups with high and low risk. (D) Expression of immune checkpoints in the groups with high and low risk. (E) Risk scores in anti-PD1/PD-L1-treated CR/PR samples and SD/PD samples in the IMvigor210 cohort. (F) Risk scores in relation to immune cell subtypes. (G) RNAss values in the groups with high and low risk. (H) Relationships between risk scores and RNAss values. (I) RNAss values in deceased and surviving patients. (J) Proportions of deceased and surviving patients in the groups with high and low RNAss. GSEA, gene set enrichment analysis; TIDE, tumor immune dysfunction and rejection; PR, partial response; CR, complete response; PD, progressive disease; SD, stable disease; RNAss, RNA stemness scores. The asterisks represented the statistical p value (*P < 0.05; **P < 0.01; ***P < 0.001.)

Crosstalk among cuproptosis clusters, gene clusters, risk signature, and clinicopathological features

The alluvial map shown in Figure 7A illustrated the crosstalk among cuproptosis clusters, gene clusters, and risk scores. Higher risk scores were associated with gene cluster B rather than gene cluster A (Figure 7B) and the cuproptosis cluster B also had higher risk scores than cluster A (Figure 7C). The prognosis was significantly enhanced in both cluster A and gene cluster A in comparison with their respective B clusters, again confirming the reliability and consistency of the analysis. Figure 7D shows a comprehensive heatmap of risk scores in relation to clinical features (including age, sex, grade, and TNM, T, N, and M stages), cuproptosis clustering, gene clustering, and cuproptosis regulator expression. A strong association between cuproptosis clusters, gene clusters, and risk signatures can be observed.

FIGURE 7

Figure 7 Relationship cross-links of cuproptosis clusters, gene clusters, risk signatures, and clinical features. (A) Alluvial diagram showing cuproptosis clusters, gene clusters, risk grouping, and survival status. (B) Risk scores in the cuproptosis clusters. (C) Risk scores in the gene clusters. (D) Heatmaps showing the integration of cuproptosis clusters, gene clusters, clinicopathological characteristics, and cuproptosis regulators’ expression in relation to risk groups.

An examination of the correlation between the risk signature and genetic mutations

Differences in somatic mutation distribution between low- and high-risk scores in the TCGA set were investigated with the “maftools” package. As Figures 8A, B show, the TMB was greater in the group with the high-risk scores, together with an 88.89% mutation rate (altered in 144 of 162 samples) versus 80.72% for the low-scoring group (altered in 134 of 166 samples). Quantitative analysis confirmed that tumors with high scores were correlated with higher TMB values (Figure 8C) and that risk scores and TMB values were significantly positively correlated (Figure 8D). It has been reported that high TMB is related to long-term clinical sensitivity to anti-PD1/PD-L1 therapy (44), confirmed by the findings shown in Figures 6C–E. This suggests that curoptosis-associated variations in tumors may be critical for the anti-PD-1/PD-L1 therapeutic response. High TMB values were also linked to reduced immune cell infiltration in HCC (Figure 8E) while K-M survival curves showed an association between elevated TMB and worse OS (Figure 8F), indirectly confirming the effectiveness of the risk score in the prediction of immunotherapy outcome.

FIGURE 8

Figure 8 Associations between risk signatures and gene mutations. (A, B) Mutation frequencies in the groups with high and low risk (C) TMB values in the groups with high and low risk (D) Association between risk scores and TMB. (E) Infiltration of immune cell in the groups with high and low TMB. (G) K-M curves for OS in the groups with high and low TMB. *P < 0.05; **P < 0.01; ***P < 0.001. TMB, tumor mutational burden; K-M, Kaplan-Meier.

Three DEGs in the risk signature

Finally, we analyzed the three DEGs in the risk signature. The expression of these genes was found to significantly raised in TCGA-HCC samples in comparison with normal tissue (Figures 9A–C). These results were confirmed by IHC analysis (Figures 9D–F). On the basis of gene expression, patients in the TCGA were classified into one of two groups. Low levels of expression were linked to a better prognosis, thus confirming the results in Figure 5, and all three genes represented prognostic risk factors for HCC in the risk signature (Figures 9G–I). Moreover, the DEGs expression was revealed to be positively associated with immune cell infiltration using TIMER (Figures S6A–C). Boxplots were used to compare the immune cell subset distribution with CNVs, and it can be seen that the greatest differences in arm-level gain in relation to immune infiltration were with CCNB1 and G6PD, whereas the level of immune infiltration with arm-level deletion was most significantly different in CDCA8 (Figures S6D–F).

FIGURE 9

Figure 9 Expression of the three signature DEGs and their impact on prognosis. Expression of CCNB1 (A), CDCA8 (B), and G6PD (C) in tumor and normal samples (GEPIA). IHC analysis of CCNB1 (D), CDCA8 (E), and G6PD (F) expression in tumor and normal samples. K-M survival analysis between high- and low-expression groups of CCNB1 (G), CDCA8 (H), and G6PD (I).*P < 0.05;****P < 0.0001. DEGs, differentially expressed genes.

Discussion

Recent research has shown that the homeostasis of copper is rigorously regulated and that any imbalance reduces the organism’s fitness (45). Disruptions to copper homeostasis have also been linked with tumor growth and irreversible damage (46). Thus, the role of copper in vivo has attracted much attention, specifically in the field of tumor therapy. Cuproptosis is a form of programmed cell death; the regulation of cuproptosis is complex, involving numerous regulatory factors. The role of cuproptosis in tumor development and its relationship to immunity has not been fully evaluated. Here, we systematically characterized immune cell infiltration mediated by cuproptosis, as well as the corresponding cuproptosis regulatory subtypes. In addition, a cuproptosis-related risk signature was developed to correlate individual cuproptosis subtypes with the patient’s sensitivity to immunotherapy, suggesting the potential usefulness of the signature for personalized therapy.

The study initially evaluated the expression levels, somatic mutations, and CNVs of 10 cuproptosis regulators, observing that all of the regulatory genes were shown to be expressed differently in tumor and control tissues. The highest mutation frequency was seen in CDKN2A. In addition, the copy numbers of most of the cuproptosis regulators were altered. This suggests that dysfunctional expression of these regulator genes may play an essential role in HCC. Unsupervised consensus cluster analysis of the regulators identified two distinct subtypes in HCC, termed cuproptosis clusters A and B. Surprisingly, GSEA enrichment analysis found that cluster A was mainly related to metabolic processes, while cluster B was enriched in various oncogenic and immune- and stromal-related signaling pathways. Stromal cells are documented to play immunomodulatory roles in tumors (47, 48), are able to prevent the entry of immune cells into the tumor parenchyma (49) and, even if T cells do enter the tumor microenvironment, stromal cells can surround them and prevent an effective immune response (50). Further analysis showed that cuproptosis cluster B was correlated with significant infiltration of immune cells, including regulatory cells such as regulatory T cells and CD 56bright NK cells. These findings suggest that cuproptosis cluster B contains a large number of immune cells, but activity may be suppressed by stromal and immune regulatory cells. K-M survival analysis showed that cluster B was linked with a poorer prognosis, possibly due to tumor escape resulting from immunosuppression.

To quantify the cuproptosis regulatory subtypes in individual tumors, a risk signature based on three DEGs was developed as a scoring system for individual HCC patients in the TCGA cohort. Both K-M survival and ROC curve evaluations verified the signature’s accuracy and reliability for predicting patient prognosis. In addition, we independently verified the signature using gene expression and clinical data from the ICGC database. The signature served as an independent prognostic variable in both the TCGA and ICGC cohorts, as demonstrated by multivariate Cox regression. A nomogram using a combination of the risk signature and clinical features was found to be more effective than other clinical features, greatly improving the clinical utility of the signature.

In our study, three DEGs were identified as cuproptosis-related signature genes. Glucose-6-phosphate dehydrogenase (G6PD) is a rate-limiting enzyme in the pentose phosphate pathway and is involved in energy generation through the maintenance of reduced NADPH co-enzyme levels (51). Because of its key metabolic role, G6PD is also involved in tumor pathogenesis where it has been reported to modulate proliferation (52, 53), metastasis (54), chemoresistance (55, 56), immune activation (57, 58), and tumor ferroptosis (59). Cell division cycle-associated 8 (CDCA8) forms part of the chromosomal passenger complex (CPC) and is necessary for the stabilization of the mitotic spindle. Disruption of cell cycle regulation is a hallmark of tumor development (60). CDCA8 levels are elevated in a variety of tumor types where it has been implicated in the growth and tumor progression (61–64). Cyclin B1 (CCNB1) belongs to the cyclin family and is the regulatory component of cyclin-dependent kinase 1 (65, 66). It has important functions in cell cycle regulation and dysfunctional expression promotes the development of various cancer types, including colon (67–69), cervical cancer (70), and kidney cancer (71). CCNB1 overexpression leads to unplanned entry into the cell cycle, uncontrolled cell proliferation, and tumorigenesis (69, 72–74). Studies have also shown that copper can affect the expression of G6PD and CCNB1, confirming a possible link between DEGs and cuproptosis (75, 76).

Although immunotherapy is effective for HCC (77), its efficacy has not proved consistent due to an incomplete understanding of the immune microenvironment and its variations in individual patients. It is thus necessary to determine which patients are most likely to react to immunotherapy. Here, it was found that the two risk-score subgroups had distinct immune infiltration characteristics with greater infiltration and matrix activity seen in patients with higher risk score. These patients also had lower TIDE scores, suggesting that they were more likely to be sensitive to immunotherapy and less likely to show immune escape (78). Another important result was that there was higher expression of immune checkpoints in this group. And these genes may be targeted by immunotherapy to determine the clinical response of patients (79). The accuracy of the cuproptosis-related risk signature for predicting the response to immunotherapy was confirmed by analysis in the IMvigor210 cohort. In addition, we found differences in TMB values between the groups and demonstrated an association between the risk and TMB scores. It has been reported that higher non-synonymous mutational burdens in tumors result in the formation of greater numbers of neoantigens, leading to greater immunogenicity and enhancing the immunotherapy response (80). This was an additional confirmation of the likelihood that high-risk patients will respond to immunotherapy. These findings indicate that the risk score is both reliable and effective for evaluating cuproptosis-modulating subtypes in individual patients and can also be applied effectively to determining the degree of immune cell infiltration and providing guidance for immunotherapy.

Conclusion

In conclusion, our study identified two patterns of cuproptosis regulation based on the expression of 10 cuproptosis regulators that are important contributors to the heterogeneity of immune cell infiltration. The cuproptosis-related risk signature, which is closely associated with the prognosis, immune cell infiltration characteristics, and sensitivity to immunotherapy in HCC patients, can quantify the risk of individual patients, providing new directions for individualized anti-tumor immunotherapy for HCC patients.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

Ethics statement

The affiliated hospital of Qingdao University approved the collection of clinical data for research purposes. (Approval No. QYFYWZLL26947). The patients/participants provided their written informed consent to participate in this study.

Author contributions

GW, RX, and SZ analyzed the data. GW, RX, LS, and JG wrote the manuscript. WL, YZ, and XB reviewed and comprehensively revised the manuscript. WQ and SW contributed to the design of the study. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by four grants: 1. Beijing Xisike Clinical Oncology Research Foundation (Y-XD2019-11) 2. Qingdao traditional Chinese medicine science and technology project (2021-zyym29) 3. Shandong medical and health science and technology development plan project (202103030554) 4. State Key Laboratory of Ultrasound in Medicine and Engineering (2021KFKT029).

Acknowledgments

The authors gratefully acknowledge the data provided by patients and researchers participating in TCGA and ICGC.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2022.945516/full#supplementary-material

Supplementary Figure 1 | Curves of K-M for OS in the groups with high and low expression of 10 cuproptosis regulators. (A) FDX1; (B) LITP1; (C) LIAS; (D) DLD; (E) DLAT; (F) PDHA1; (G) PDHB; (H) MTF1; (I) GLS; (J) CDKN2A. K-M Kaplan-Meier, OS overall survival

Supplementary Figure 2 | Cuproptosis clustering. (A–G) Consensus matrix based on cuproptosis regulator expression in the HCC cohort at k = 2 - 9. (H) Consensus CDF plots for cuproptosis clustering. (I) The area under the CDF curve relative changes for cuproptosis clustering

Supplementary Figure 3 | Gene clustering (A–G) Consensus matrix based on cuproptosis-related gene expression in the HCC cohort at k= 2-9. (H) Consensus CDF plots for gene clustering. (I) The area under the CDF curve relative changes for gene clustering. HCC hepatocellular carcinoma, CDF cumulative distribution function

Supplementary Figure 4 | Comparison of risk scores within different grades or between different classifications of clinicopathological features. (A) Stage; (B)T; (C) Grade; (D) AFP expression; (E) Sex; (F) Age.

Supplementary Figure 5 | Comparison of drug sensitivity between groups with high and low risk. (A) Gefitinib; (B) Sorafenib; (C) Nilotinib; (D) Dasatinib; (E) Erlotinib; (F) Metformin; (G) Bleomycin; (H) Doxorubicin; (I) Gemcitabine; (J) Tipifarnib; (K) Imatinib; (L) Mitomycin.C.

Supplementary Figure 6 | The relationship between the three DEGs in the signature and immunity and copy number. (A–C) Correlation between expression of gene and infiltration of immune cell. (D–F) Association of gene copy number alterations with infiltration of immune cells. *P < 0.05; **P < 0.01; ***P < 0.001. DEGs differentially expressed genes

References

1. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA: Cancer J Clin (2021) 71(3):209–49. doi: 10.3322/caac.21660

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Mehta N. Hepatocellular carcinoma-how to determine therapeutic options. Hepatol Commun (2020) 4(3):342–54. doi: 10.1002/hep4.1481

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Fattovich G, Stroffolini T, Zagni I, Donato F. Hepatocellular carcinoma in cirrhosis: incidence and risk factors. Gastroenterology. (2004) 127(5 Suppl 1):S35–50. doi: 10.1053/j.gastro.2004.09.014

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Yin JM, Sun LB, Zheng JS, Wang XX, Chen DX, Li N. Copper chelation by trientine dihydrochloride inhibits liver RFA-induced inflammatory responses in vivo. Inflammation Res Off J Eur Histamine Res Soc [et al] (2016) 65(12):1009–20. doi: 10.1007/s00011-016-0986-2

CrossRef Full Text | Google Scholar

5. Zheng Y, Wang T, Tu X, Huang Y, Zhang H, Tan D, et al. Gut microbiome affects the response to anti-PD-1 immunotherapy in patients with hepatocellular carcinoma. J immunother cancer. (2019) 7(1):193. doi: 10.1186/s40425-019-0650-9

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Tsvetkov P, Coy S, Petrova B, Dreishpoon M, Verma A, Abdusamad M, et al. Copper induces cell death by targeting lipoylated TCA cycle proteins. Sci (New York NY) (2022) 375(6586):1254–61. doi: 10.1126/science.abf0529

CrossRef Full Text | Google Scholar

7. Tapiero H, Townsend DM, Tew KD. Trace elements in human physiology and pathology. Copper Biomedicine pharmacotherapy = Biomedecine pharmacotherapie. (2003) 57(9):386–98. doi: 10.1016/S0753-3322(03)00012-X

CrossRef Full Text | Google Scholar

8. Vassiliev V, Harris ZL, Zatta P. Ceruloplasmin in neurodegenerative diseases. Brain Res Brain Res Rev (2005) 49(3):633–40. doi: 10.1016/j.brainresrev.2005.03.003

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Fraga CG. Relevance, essentiality and toxicity of trace elements in human health. Mol aspects Med (2005) 26(4-5):235–44. doi: 10.1016/j.mam.2005.07.013

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Squitti R, Polimanti R. Copper hypothesis in the missing hereditability of sporadic alzheimer's disease: ATP7B gene as potential harbor of rare variants. J Alzheimer's Dis JAD. (2012) 29(3):493–501. doi: 10.3233/JAD-2011-111991

CrossRef Full Text | Google Scholar

11. Squitti R, Barbati G, Rossi L, Ventriglia M, Dal Forno G, Cesaretti S, et al. Excess of nonceruloplasmin serum copper in AD correlates with MMSE, CSF [beta]-amyloid, and h-tau. Neurology. (2006) 67(1):76–82. doi: 10.1212/01.wnl.0000223343.82809.cf

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Huang YL, Sheu JY, Lin TH. Association between oxidative stress and changes of trace elements in patients with breast cancer. Clin Biochem (1999) 32(2):131–6. doi: 10.1016/S0009-9120(98)00096-4

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Kuo HW, Chen SF, Wu CC, Chen DR, Lee JH. Serum and tissue trace elements in patients with breast cancer in Taiwan. Biol Trace element Res (2002) 89(1):1–11. doi: 10.1385/BTER:89:1:1

CrossRef Full Text | Google Scholar

14. Rizk SL, Sky-Peck HH. Comparison between concentrations of trace elements in normal and neoplastic human breast tissue. Cancer Res (1984) 44(11):5390–4.

PubMed Abstract | Google Scholar

15. Habib FK, Dembinski TC, Stitch SR. The zinc and copper content of blood leucocytes and plasma from patients with benign and malignant prostates. Clinica chimica acta; Int J Clin Chem (1980) 104(3):329–35. doi: 10.1016/0009-8981(80)90390-3

CrossRef Full Text | Google Scholar

16. Nayak SB, Bhat VR, Upadhyay D, Udupa SL. Copper and ceruloplasmin status in serum of prostate and colon cancer patients. Indian J Physiol Pharmacol (2003) 47(1):108–10.

PubMed Abstract | Google Scholar

17. Chen D, Dou QP. New uses for old copper-binding drugs: Converting the pro-angiogenic copper to a specific cancer cell death inducer. Expert Opin Ther targets (2008) 12(6):739–48. doi: 10.1517/14728222.12.6.739

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Majumder S, Chatterjee S, Pal S, Biswas J, Efferth T, Choudhuri SK. The role of copper in drug-resistant murine and human tumors. Biometals an Int J role metal ions biology biochemistry Med (2009) 22(2):377–84. doi: 10.1007/s10534-008-9174-3

CrossRef Full Text | Google Scholar

19. Zhang X, Yang Q. Association between serum copper levels and lung cancer risk: A meta-analysis. J Int Med Res (2018) 46(12):4863–73. doi: 10.1177/0300060518798507

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Turecký L, Kalina P, Uhlíková E, Námerová S, Krizko J. Serum ceruloplasmin and copper levels in patients with primary brain tumors. Klinische Wochenschrift. (1984) 62(4):187–9. doi: 10.1007/BF01731643

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Moffett JR, Puthillathu N, Vengilote R, Jaworski DM, Namboodiri AM. Acetate revisited: A key biomolecule at the nexus of metabolism, epigenetics and oncogenesis-part 1: Acetyl-CoA, acetogenesis and acyl-CoA short-chain synthetases. Front Physiol (2020) 11:580167. doi: 10.3389/fphys.2020.580167

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Muñoz C, Rios E, Olivos J, Brunser O, Olivares M. Iron, copper and immunocompetence. Br J Nutr (2007) 98 Suppl 1:S24–8. doi: 10.1017/S0007114507833046

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Djoko KY, Ong CL, Walker MJ, McEwan AG. The role of copper and zinc toxicity in innate immune defense against bacterial pathogens. J Biol Chem (2015) 290(31):18954–61. doi: 10.1074/jbc.R115.647099

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Festa RA, Thiele DJ. Copper at the front line of the host-pathogen battle. Plos pathogens (2012) 8(9):e1002887. doi: 10.1371/journal.ppat.1002887

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Percival SS.Copper and immunity. Am J Clin Nutr (1998) 67 (5 Suppl):1064s–8s. doi: 10.1093/ajcn/67.5.1064S

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Voli F, Valli E, Lerra L, Kimpton K, Saletta F, Giorgi FM, et al. Intratumoral copper modulates PD-L1 expression and influences tumor immune evasion. Cancer Res (2020) 80(19):4129–44. doi: 10.1158/0008-5472.CAN-20-0471

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Wilkerson MD, Hayes DN. ConsensusClusterPlus: A class discovery tool with confidence assessments and item tracking. Bioinf (Oxford England) (2010) 26(12):1572–3. doi: 10.1093/bioinformatics/btq170

CrossRef Full Text | Google Scholar

28. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA (2005) 102(43):15545–50. doi: 10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Yu G, Wang LG, Han Y, He QY. clusterProfiler: An r package for comparing biological themes among gene clusters. Omics: J Integr Biol (2012) 16(5):284–7. doi: 10.1089/omi.2011.0118

CrossRef Full Text | Google Scholar

30. Bindea G, Mlecnik B, Tosolini M, Kirilovsky A, Waldner M, Obenauf AC, et al. Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity (2013) 39(4):782–95. doi: 10.1016/j.immuni.2013.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Finotello F, Trajanoski Z. Quantifying tumor-infiltrating immune cells from transcriptomics data. Cancer immunology immunother: CII (2018) 67(7):1031–40. doi: 10.1007/s00262-018-2150-z

CrossRef Full Text | Google Scholar

32. Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, et al. Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep (2017) 18(1):248–62. doi: 10.1016/j.celrep.2016.12.019

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med (2018) 24(10):1550–8. doi: 10.1038/s41591-018-0136-1

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Li B, Severson E, Pignon JC, Zhao H, Li T, Novak J, et al. Comprehensive analyses of tumor immunity: Implications for cancer immunotherapy. Genome Biol (2016) 17(1):174. doi: 10.1186/s13059-016-1028-7

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS, et al. TIMER: A web server for comprehensive analysis of tumor-infiltrating immune cells. Cancer Res (2017) 77(21):e108–e10. doi: 10.1158/0008-5472.CAN-17-0307

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Bell D, Berchuck A, Birrer M, Chien J, Cramer DW, Dao F. Integrated genomic analyses of ovarian carcinoma. Nature. (2011) 474(7353):609–15. doi: 10.1038/nature10166

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Crown WH. Potential application of machine learning in health outcomes research and some statistical cautions. Value Health J Int Soc Pharmacoeconomics Outcomes Res (2015) 18(2):137–40. doi: 10.1016/j.jval.2014.12.005

CrossRef Full Text | Google Scholar

38. Chen DS, Mellman I. Elements of cancer immunity and the cancer-immune set point. Nature. (2017) 541(7637):321–30. doi: 10.1038/nature21349

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Zhang B, Wu Q, Li B, Wang D, Wang L, Zhou YL. m6A regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in gastric cancer. Mol cancer. (2020) 19(1):53. doi: 10.1186/s12943-020-01170-0

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Moserle L, Casanovas O. Anti-angiogenesis and metastasis: a tumour and stromal cell alliance. J Internal Med (2013) 273(2):128–37. doi: 10.1111/joim.12018

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Groth C, Hu X, Weber R, Fleming V, Altevogt P, Utikal J, et al. Immunosuppression mediated by myeloid-derived suppressor cells (MDSCs) during tumour progression. Br J cancer. (2019) 120(1):16–25. doi: 10.1038/s41416-018-0333-1

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Tian Z, Gershwin ME, Zhang C. Regulatory NK cells in autoimmune disease. J autoimmunity (2012) 39(3):206–15. doi: 10.1016/j.jaut.2012.05.006

CrossRef Full Text | Google Scholar

43. Thorsson V, Gibbs DL, Brown SD, Wolf D, Bortone DS, Ou Yang TH, et al. The immune landscape of cancer. Immunity. (2018) 48(4):812–30.e14. doi: 10.1016/j.immuni.2018.03.023

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Barroso-Sousa R, Keenan TE, Pernas S, Exman P, Jain E, Garrido-Castro AC, et al. Tumor mutational burden and PTEN alterations as molecular correlates of response to PD-1/L1 blockade in metastatic triple-negative breast cancer. Clin Cancer research: an Off J Am Assoc Cancer Res (2020) 26(11):2565–72. doi: 10.1158/1078-0432.CCR-19-3507

CrossRef Full Text | Google Scholar

45. De la Iglesia R, Valenzuela-Heredia D, Andrade S, Correa J, González B. Composition dynamics of epilithic intertidal bacterial communities exposed to high copper levels. FEMS Microbiol ecology (2012) 79(3):720–7. doi: 10.1111/j.1574-6941.2011.01254.x

CrossRef Full Text | Google Scholar

46. Guo W, Zhang X, Lin L, Wang H, He E, Wang G, et al. The disulfiram/copper complex induces apoptosis and inhibits tumour growth in human osteosarcoma by activating the ROS/JNK signalling pathway. J Biochem (2021) 170(2):275–87. doi: 10.1093/jb/mvab045

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Cremasco V, Astarita JL, Grauel AL, Keerthivasan S, MacIsaac K, Woodruff MC, et al. FAP delineates heterogeneous and functionally divergent stromal cells in immune-excluded breast tumors. Cancer Immunol Res (2018) 6(12):1472–85. doi: 10.1158/2326-6066.CIR-18-0098

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Canzonetta C, Pelosi A, Di Matteo S, Veneziani I, Tumino N, Vacca P, et al. Identification of neuroblastoma cell lines with uncommon TAZ+/mesenchymal stromal cell phenotype with strong suppressive activity on natural killer cells. J immunother Cancer (2021) 9(1):1055. doi: 10.1136/jitc-2020-001313

CrossRef Full Text | Google Scholar

49. Turley SJ, Cremasco V, Astarita JL. Immunological hallmarks of stromal cells in the tumour microenvironment. Nat Rev Immunol (2015) 15(11):669–82. doi: 10.1038/nri3902

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Nishio N, Diaconu I, Liu H, Cerullo V, Caruana I, Hoyos V, et al. Armed oncolytic virus enhances immune functions of chimeric antigen receptor-modified T cells in solid tumors. Cancer Res (2014) 74(18):5195–205. doi: 10.1158/0008-5472.CAN-14-0697

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Pinna A, Solinas G, Giancipoli E, Porcu T, Zinellu A, D'Amico-Ricci G, et al. Glucose-6-Phosphate dehydrogenase (G6PD) deficiency and late-stage age-related macular degeneration. Int J Med Sci (2019) 16(5):623–9. doi: 10.7150/ijms.30155

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Yang HC, Wu YH, Yen WC, Liu HY, Hwang TL, Stern A, et al. The redox role of G6PD in cell growth, cell death, and cancer. Cells. (2019) 8(9):e001313. doi: 10.3390/cells8091055

CrossRef Full Text | Google Scholar

53. Ghergurovich JM, Esposito M, Chen Z, Wang JZ, Bhatt V, Lan T, et al. Glucose-6-Phosphate dehydrogenase is not essential for K-Ras-Driven tumor growth or metastasis. Cancer Res (2020) 80(18):3820–9. doi: 10.1158/0008-5472.CAN-19-2486

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Lu M, Lu L, Dong Q, Yu G, Chen J, Qin L, et al. Elevated G6PD expression contributes to migration and invasion of hepatocellular carcinoma cells by inducing epithelial-mesenchymal transition. Acta Biochim Biophys Sinica (2018) 50(4):370–80. doi: 10.1093/abbs/gmy009

CrossRef Full Text | Google Scholar

55. Yamawaki K, Mori Y, Sakai H, Kanda Y, Shiokawa D, Ueda H, et al. Integrative analyses of gene expression and chemosensitivity of patient-derived ovarian cancer spheroids link G6PD-driven redox metabolism to cisplatin chemoresistance. Cancer letters (2021) 521:29–38. doi: 10.1016/j.canlet.2021.08.018

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Ju HQ, Lu YX, Wu QN, Liu J, Zeng ZL, Mo HY, et al. Disrupting G6PD-mediated redox homeostasis enhances chemosensitivity in colorectal cancer. Oncogene. (2017) 36(45):6282–92. doi: 10.1038/onc.2017.227

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Lu C, Yang D, Klement JD, Colson YL, Oberlies NH, Pearce CJ, et al. G6PD functions as a metabolic checkpoint to regulate granzyme b expression in tumor-specific cytotoxic T lymphocytes. J immunother Cancer (2022) 10(1):e003543. doi: 10.1136/jitc-2021-003543

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Nakamura M, Nagase K, Yoshimitsu M, Magara T, Nojiri Y, Kato H, et al. Glucose-6-phosphate dehydrogenase correlates with tumor immune activity and programmed death ligand-1 expression in merkel cell carcinoma. J immunother Cancer (2020) 8(2):e001679. doi: 10.1136/jitc-2020-001679

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Cao F, Luo A, Yang C. G6PD inhibits ferroptosis in hepatocellular carcinoma by targeting cytochrome P450 oxidoreductase. Cell signalling (2021) 87:110098. doi: 10.1016/j.cellsig.2021.110098

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Hanahan D. Hallmarks of cancer: New dimensions. Cancer discovery (2022) 12(1):31–46. doi: 10.1158/2159-8290.CD-21-1059

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Gassmann R, Carvalho A, Henzing AJ, Ruchaud S, Hudson DF, Honda R, et al. Borealin: A novel chromosomal passenger required for stability of the bipolar mitotic spindle. J Cell Biol (2004) 166(2):179–91. doi: 10.1083/jcb.200404001

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Ruchaud S, Carmena M, Earnshaw WC. Chromosomal passengers: Conducting cell division. Nat Rev Mol Cell Biol (2007) 8(10):798–812. doi: 10.1038/nrm2257

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Wang Y, Zhao Z, Bao X, Fang Y, Ni P, Chen Q, et al. Borealin/Dasra b is overexpressed in colorectal cancers and contributes to proliferation of cancer cells. Med Oncol (Northwood London England) (2014) 31(11):248. doi: 10.1007/s12032-014-0248-5

CrossRef Full Text | Google Scholar

64. Dai C, Miao CX, Xu XM, Liu LJ, Gu YF, Zhou D, et al. Transcriptional activation of human CDCA8 gene regulated by transcription factor NF-y in embryonic stem cells and cancer cells. J Biol Chem (2015) 290(37):22423–34. doi: 10.1074/jbc.M115.642710

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Santos SD, Wollman R, Meyer T, Ferrell JE. Spatial positive feedback at the onset of mitosis. Cell. (2012) 149(7):1500–13. doi: 10.1016/j.cell.2012.05.028

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Nixon VL, Levasseur M, McDougall A, Jones KT. Ca(2+) oscillations promote APC/C-dependent cyclin B1 degradation during metaphase arrest and completion of meiosis in fertilizing mouse eggs. Curr biology: CB. (2002) 12(9):746–50. doi: 10.1016/S0960-9822(02)00811-4

CrossRef Full Text | Google Scholar

67. Bondi J, Husdal A, Bukholm G, Nesland JM, Bakka A, Bukholm IR. Expression and gene amplification of primary (A, B1, D1, D3, and e) and secondary (C and h) cyclins in colon adenocarcinomas and correlation with patient outcome. J Clin pathology (2005) 58(5):509–14. doi: 10.1136/jcp.2004.020347

CrossRef Full Text | Google Scholar

68. Wang A, Yoshimi N, Ino N, Tanaka T, Mori H. Overexpression of cyclin B1 in human colorectal cancers. J Cancer Res Clin Oncol (1997) 123(2):124–7. doi: 10.1007/BF01269891

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Wang Q, Moyret-Lalle C, Couzon F, Surbiguet-Clippe C, Saurin JC, Lorca T, et al. Alterations of anaphase-promoting complex genes in human colon cancer cells. Oncogene (2003) 22(10):1486–90. doi: 10.1038/sj.onc.1206224

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Zhao M, Kim YT, Yoon BS, Kim SW, Kang MH, Kim SH, et al. Expression profiling of cyclin B1 and D1 in cervical carcinoma. Exp Oncol (2006) 28(1):44–8.

PubMed Abstract | Google Scholar

71. Ikuerowo SO, Kuczyk MA, Mengel M, van der Heyde E, Shittu OB, Vaske B, et al. Alteration of subcellular and cellular expression patterns of cyclin B1 in renal cell carcinoma is significantly related to clinical progression and survival of patients. Int J cancer (2006) 119(4):867–74. doi: 10.1002/ijc.21869

PubMed Abstract | CrossRef Full Text | Google Scholar

72. Heald R, McLoughlin M, McKeon F. Human wee1 maintains mitotic timing by protecting the nucleus from cytoplasmically activated Cdc2 kinase. Cell. (1993) 74(3):463–74. doi: 10.1016/0092-8674(93)80048-J

PubMed Abstract | CrossRef Full Text | Google Scholar

73. Fotedar R, Flatt J, Gupta S, Margolis RL, Fitzgerald P, Messier H, et al. Activation-induced T-cell death is cell cycle dependent and regulated by cyclin b. Mol Cell Biol (1995) 15(2):932–42. doi: 10.1128/MCB.15.2.932

PubMed Abstract | CrossRef Full Text | Google Scholar

74. Sarafan-Vasseur N, Lamy A, Bourguignon J, Le Pessot F, Hieter P, Sesboüé R, et al. Overexpression of b-type cyclins alters chromosomal segregation. Oncogene. (2002) 21(13):2051–7. doi: 10.1038/sj.onc.1205257

PubMed Abstract | CrossRef Full Text | Google Scholar

75. Bayramoglu Akkoyun M, Temel Y, Bengü A, Akkoyun HT. Ameliorative effects of astaxanthin against copper(II) ion-induced alteration of pentose phosphate pathway and antioxidant system enzymes in rats. Environ Sci pollut Res Int (2021) 28(44):62919–26. doi: 10.1007/s11356-021-15017-8

PubMed Abstract | CrossRef Full Text | Google Scholar

76. Hanagata N, Zhuang F, Connolly S, Li J, Ogawa N, Xu M. Molecular responses of human lung epithelial cells to the toxicity of copper oxide nanoparticles inferred from whole genome expression analysis. ACS nano (2011) 5(12):9326–38. doi: 10.1021/nn202966t

PubMed Abstract | CrossRef Full Text | Google Scholar

77. Sun B, Yang D, Dai H, Liu X, Jia R, Cui X, et al. Eradication of hepatocellular carcinoma by NKG2D-based CAR-T cells. Cancer Immunol Res (2019) 7(11):1813–23. doi: 10.1158/2326-6066.CIR-19-0026

PubMed Abstract | CrossRef Full Text | Google Scholar

78. Wang S, He Z, Wang X, Li H, Liu XS. Antigen presentation and tumor immunogenicity in cancer immunotherapy response prediction. eLife (2019) 8:e49020. doi: 10.7554/eLife.49020

PubMed Abstract | CrossRef Full Text | Google Scholar

79. Hu FF, Liu CJ, Liu LL, Zhang Q, Guo AY. Expression profile of immune checkpoint genes and their roles in predicting immunotherapy response. Briefings Bioinf (2021) 22(3):bbaa176. doi: 10.1093/bib/bbaa176

CrossRef Full Text | Google Scholar

80. Rizvi NA, Hellmann MD, Snyder A, Kvistborg P, Makarov V, Havel JJ, et al. Cancer immunology. mutational landscape determines sensitivity to PD-1 blockade in non-small cell lung cancer. Sci (New York NY) (2015) 348(6230):124–8. doi: 10.1126/science.aaa1348

CrossRef Full Text | Google Scholar