Exploration of cuprotosis-related genes for predicting prognosis and immunological characteristics in acute myeloid leukaemia based on genome and transcriptome

Data collection

The TCGA-LAML dataset was downloaded from the TGCA database, and GSE37642 (GSE37642-GPL96, GSE37642-GPL570) and GSE12417 (GSE12417-GPL96, and GSE12417-GPL570) were downloaded from the GEO database. Mutation data and copy number variation (CNV) data were also downloaded from TCGA database. We used the “sva” package to eliminate the batch effect in RNA-seq and microarray, Finally, we annotated 11,1917 protein-coding mRNA in the meta dataset. We excluded incomplete survival information, non-whole bone marrow sequencing, M3, and repeated sequencing samples. Finally, 116 patients were included in the TCGA cohort, 539 patients in the GSE37642 cohort, and 233 patients in the GSE12417 cohort. Cuprotosis-related genes (CRGs) were collected from previous references. The expression data were log2 transformed and normalized. The Gene MANIA database was used to construct protein-protein interactions (PPI) network.

Unsupervised consensus clustering

In the meta-cohort (combined TCGA dataset and GEO dataset), based on prognostic CRGs, we used unsupervised consistent clustering (“ConsensusClusterPlus” package) to divide all patients and principal component analysis (PCA) was used to determine whether each subtype was relatively independent of the others. Importantly, 1000 repetitions were performed, pltem = 0.8, to verify the stability of the subtype in unsupervised consistent clustering.

Functional enrichment analysis

Gene set variation analysis (GSVA) was used to evaluate differences in biological pathways between subtypes, using c2.cp.kegg.v7.0.symbol as a reference gene set, FDR < 0.05 was the threshold. Moreover, we used the “limma” package to analyze the differential expression genes (DEGs) between different subtypes, next, we used “clusterProfiler” package for functional enrichment based on DEGs, P-value < 0.05 and q-value < 0.05 was the threshold.

Tumor immune microenvironment

Overall, following the pipeline in previous studies, we used eight algorithms (TIMER, CIBERSORT, QUANTISEQ, MCP-counter, XCELL, EPIC, ssGSEA, and ESTIMATE) to estimate the abundances of immune cells for each sample in TCGA and GEO dataset.

Response to targeted drugs

Relevant references have confirmed the important role of the cell cycle, PI3K/mTOR pathway, and Wnt pathway in the progression of AML. Hence, we used the “pRRophetic” package to calculate the half maximal inhibitory concentration (IC50) of different targeted drugs, including CGP.60474 (cell cycle), JW.7.52.1 (PI3K/mTOR), and CHIR.99021 (Wnt).

Machine learning-derived risk cuprotosis-related risk score

We performed our previous workflow to construct a consensus prognosis model for AML patients. Firstly, we constructed 55 combination of machine learning algorithms based on nine algorithms, including LASSO, GBM, Survival-SVM, SuperPC, ridge regression, plsRcox, CoxBoost, StepCox, and Enet [15]. We used models that can perform variable filtering as the pre-model. Subsequently, we used the GSE37642 as the training set to construct signatures in prognostic DEGs matrix based on different cuprotosis subtypes. We selected the best consensus prognostic model based on the mean C-index of the three cohorts (GSE37642, TCGA-LAML and GSE12417). Risk score was calculated for each patient based on the expression of each candidate variable and the coefficients in the final model. ROC curve was used to evaluate the predictive performance of risk score for OS. The prognostic difference was analyzed by Kaplan-Meier analysis and log-rank test. In addition, we downloaded the IMvigor-210 cohort, and GSE78220 cohort from the previous references to generate a risk score and verify the effect of the risk score in the immunotherapy cohort.

Quantitative real-time polymerase chain reaction

Peripheral blood samples were collected from AML patients (12 cases) and healthy volunteers (11 cases) from Puyang Oilfield General Hospital. Patient RNA was extracted by Trizol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s protocol and cDNA was synthesized by reverse transcription kits (TAKARA, Japan). QRT-PCR was performed using a Step One Plus qRT-PCR machine (Applied Biosystems, Waltham MA, USA). All primers were provided by Suzhou Jinweizhi Biotechnology Co., Ltd (Jiangsu, China). The primer sequences as shown in Table 1.

Statistical analysis and software

The statistical analyses were conducted in the R software (version 4.1.2). Specific statistical methods have been mentioned in the bioinformatics methods above. ^***, ^**, ^*, ns refers to p < 0.001, < 0.01, < 0.05, and not significant, respectively. The Microsoft Office PowerPoint software is used to generate images for this article.

Genetic variation and expression of cuprotosis-related regulatory genes in AML

We first explored the landscape of 12 annotated CRGs in meta cohorts. The location of the 12 CRGs in the chromosomes was shown (Figure 1A). We further compared the expression of these 12 CRGs in different risk categories (ELN 2017 risk category) of AML. Among AML patients in different risk categories, we found differences in DLD, LIAS, FDX1, ATP7A, and PDHD in the ELN 2017 risk category (Figure 1B). Furthermore, we compared the expression of 12 CRGs between normal and AML samples, and the results showed the expression of 11 CRGs (except for PDHB) was significantly different in normal and AML patients (Figure 1C, 1D). In detail, the expression of genes LIPT1, DBT, and DLST were significantly higher in AML patients. We also found CNV amplification in FDX1, DLAT, ATP7B, DBT, DLST and LIPT1 among these 12 CRGs, while mainly showed loss in other CRGs (Figure 1E). We also found DLD mutations with a frequency of only 1% (Figure 1F). The above results showed the expression difference maybe not be caused by somatic mutations, but by CNV.

Analysis of the characteristics of different cuprotosis subtypes

We investigated the relationships between the 12 CRGs, and the network showed that all 12 CRGs were associated (Figure 2A). Multiple data sets were normalized to increase the sample size (Supplementary Figure 1). In the network of survival and correlation analysis, the results showed that 12 CRGs were still highly associated (Figure 2B). In addition, we investigated the prognostic value of 12 CRGs, we carried out Kaplan-Meier analysis of each CRGs, and the results showed that 11 CRGs had statistical significance on the overall survival (Supplementary Figure 2). Then, we divided AML patients into 3 subtypes, called cuprotosis subtypes A (n = 158), B (n = 426), and C (n = 304) (Figure 2C). The above various subtypes were demonstrated by PCA to be independent of one another over the whole gene expression spectrum (Figure 2D). Patients with subtype A had the best survival compared to other subtypes, according to a Kaplan-Meier analysis (Figure 2E). We further analyzed the expression of 12 CRGs in the three subtypes and all 12 CRGs were significantly differentially expressed in three cuprotosis subtypes (Figure 2F).

To investigate the immunological characteristics of the three cuprotosis subtypes, we analyzed the relationship between the aforementioned cuprotosis subtypes and immune cells. The results revealed significant differences among the three subtypes in the level of infiltration of most immune cells, with high expression of activated B cells, CD4⁺ T cells, CD8⁺ T cells, dendritic cells, natural killer cells, γδ T cells, neutrophils and type II helper T cells in subtype A (Figure 3A). In addition, the three subtypes of the immune system displayed notable differences in antigen presentation, cytokine receptors, inflammatory response, and T cell activation and suppression (Figure 3B). Additionally, we explored the differences in immune checkpoints among the subtypes and found that there were differences among most immune checkpoints (Figure 3C), with significantly increased expression of CD160, CD200, and CD27 in subtype A. Importantly, we discovered differences among the subtypes in most human leukocyte antigen class I-related gene (Figure 3D). According to the ESTIMATE algorithm-based results, there were disparities in the tumour microenvironment (TME) scores between the three subtypes, with subtype A score lower on both the immune and stromal scores (Figure 3E). Relevant references have confirmed the important role of the cell cycle, PI3K/mTOR pathway, and Wnt pathway in the progression of AML. It’s interesting to note that different subtypes respond differently to specific medications (Figure 3F). Specifically, subtype A was more sensitive to CHIR.99021, subtype B to JW.7.52.1, and subtype C to CGP.60474.

Moreover, we further analyzed the biological functions of different cuprotosis subtypes using GSVA. Enrichment analysis revealed that expression of subtype B was upregulated in the majority of metabolic pathways, including gluconeogenesis, pyruvate metabolism, amino acid metabolism, etc. (Supplementary Figure 3A). Contrarily, subtype A is mainly enriched in the porphyrin and chlorophyll metabolism and cytokinesis pathways (Supplementary Figure 3B), and differences in these pathways may be the reason why subtype A has a higher prognosis for survival when compared to other subtypes. Importantly, almost pathways were enriched in nutrient metabolism and circulation pathways, which was consistent with previous references on the intracellular mechanism of action of CRGs [13].

Characterization and identification of cuprotosis-regulatory subtypes

We identified 525 common DEGs between each other, and the above genes were enriched in cell cycle and metabolic pathways (Supplementary Figure 3C). Similarly, all samples were divided into 2 subtypes by unsupervised cluster analysis performed by DEGs, cuprotosis regulatory subtype A (n = 492) and B (n = 394) (Supplementary Figure 3D), and survival analysis of the two regulatory subtypes showed that regulatory subtype B had a better median survival time than regulatory subtype A (Supplementary Figure 3E). Subsequently, we analyzed the expression of 12 CRGs in two regulatory subtypes. In detail, the results showed that FDX1, LIAS, DLD, DLAT, PDHA1, PDHB, ATP7A, and ATP7B were significantly different (Supplementary Figure 3F).

Construction and validation of a risk score and prognostic model

The above construction of different subtypes provides some reference value for the prognosis of AML, but the above two kinds of subtypes only represented population. However, the absence of some clinical characteristics and individual assessments may result in clinical practice not being available. We calculated the C-index in all cohorts and selected the LASSO+Cox (stepwise) algorithm as the final model (Figure 4A). In detail, redundant genes were removed by LASSO regression, and 11 regulators were screened and the coefficients of each regulator were calculated by Cox regression (stepwise) analysis (Figure 4B). Next, the AML sample was divided into low-risk and high-risk groups based on risk scores, and patient fatalities increased as risk scores rose (Figure 4C). We further analyzed the expression of 11 regulators in high-risk and low-risk groups of patients in different datasets revealed that the expression of AKR1B1, ID1, SPINT2, and CYB5A was higher in high-risk patients (Figure 5A). Patients in the high-risk category in various data sets all had significantly shorter OS than those in the low-risk group (p < 0.05), according to Kaplan-Meier analysis, which was used to examine the impact of different subtypes on prognosis. The AUC value for the TCGA cohort at 1, 3, and 5 years were: 0.634, 0.708 and 0.683, and the AUC value for the GES37642 at 1, 3, and 5 years were: 0.723, 0.760 and 0.751. For GSE12741 cohort, the AUC value at 1 and 3 years were: 0.738 and 0.738 (Figure 5B). In addition, we again performed univariate and multivariate Cox regression analysis on important clinical features (FAB subtype; age; gender; ELN2017 risk category [16]; leukocytes; platelets; bone marrow blasts; runx1-runx1t1 fusion; runx1 mutation) and risk scores to determine whether they were independent prognostic factors for OS, and the results showed that risk score was an independent prognostic factor (Figure 5C).

To predict the cuprotosis model more accurately at the individual level, we first compared the risk scores of different cuprotosis subtypes and cuprotosis regulatory subtypes and found that cuprotosis subtype A and cuprotosis regulatory subtype B had lower risk scores (Supplementary Figure 4A, 4B), which also corresponded to the better prognostic survival described above. Moreover, we used a Sankey diagram to display the association between cuprotosis subtypes, cuprotosis regulatory subtypes, risk score groups, and prognosis. The results revealed that the majority of the subtypes with bad prognoses were significantly linked to high risk and low survival (Supplementary Figure 4C). Finally, we further explored the composition of immune cells, and the results revealed that the high-risk score group exhibited more Tregs cells, fibroblasts, endothelial cells, and other multiple cell infiltrates (Supplementary Figure 4D, 4E).

Validation of the risk score in immunotherapy cohorts

Due to the lack of information on therapeutic agents for AML patients, we selected immunotherapy cohorts (GSE78220 cohort and IMvigor cohort) as validation sets to predict the response to immunotherapy. In the GSE78220 cohort (anti-PD-L1), we discovered that the low-risk score group had a significantly longer median survival time compared to the high-risk group, while in the IMvigor cohort we found that the survival time and response to treatment was better in the high-risk score group than in the low-risk score group after treating PD-1 treatment (Figure 6A, 6B). Overall, our risk score model can, to some extent, guide the treatment of tumors in anti-PD-1 or anti-PD-L1. It’s interesting to note that prognostic markers for immunotherapy varied based on the various blocking sites. Finally, we analyzed somatic mutations between different risk score subtypes and showed that the frequency of mutations in the KIT gene was higher in the low-risk group, while the frequency of RUNX1 and TP53 mutations was higher in the high-risk group, and these two genes have been included in AML treatment guidelines as molecules with poorer prognosis, but the somatic mutation frequencies were not significantly different between the two groups (Figure 6C, 6D).

Validation of cuprotosis-related genes in peripheral blood

Considering that DBT, DLST, LIPT1 plays a key role in the cuprotosis-related pathways. Hence, we used qRT-PCR to assess three CRGs in blood samples from 12 AML patients and paired normal patients, the average age of the experimental group was 57.67 ± 14.55-years-old, with a male proportion of 50.00%. The age of the control group was 58.75 ± 10.20-years-old, with a male proportion of 50.00%. There was no statistical difference between the two groups (P < 0.05). The results showed that the expression of genes DBT, DLST, LIPT1 was significantly higher in AML patients than in normal patients (Supplementary Figure 5).

Data availability

Publicly available data sets were analyzed in this study. This data can be found here: Publicly available data sets can be obtained from the TCGA (https://portal.gdc.cancer.gov/), GEO (https://www.ncbi.nlm.nih.gov/gds/).

Code availability

The code that supports the findings of this study is available from the corresponding author upon request.