Cuproptosis-related risk score predicts prognosis and characterizes the tumor microenvironment in colon adenocarcinoma

Transcriptome and the corresponding clinicopathological data of COAD were downloaded from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/↗) and Gene Expression Omnibus database (GEO) (https://www.ncbi.nlm.nih.gov/geo/↗). In detail, the TCGA cohort included 480 COAD tissues and 41 normal tissues. The GEO cohort containing GSE17536, GSE29623 and GSE39582, included 827 COAD samples. The detailed clinicopathological data of these COAD patients was presented in Table S1. The TCGA and GEO cohorts were combined by using “Combat” algorithm in R to eliminate batch effects before conducting subsequent analyses. Principal component analysis (PCA) was applied to validate the effect of batch effect removal by using the R package ggplot2. In order to verified the accuracy of model, we also downloaded transcriptome and the corresponding clinicopathological data of GSE40967 from GEO database, which contained 585 COAD sampes.

Additionally, we downloaded somatic mutation data of 454 tumor samples and copy number variation (CNV) data of 506 tumor samples from TCGA.

A total of 10 CRGs (PDHA1, PDHB, FDX1, DLD, DLAT, MTF1, LIAS, LIPT1, GLS and CDKN2A) were obtained from the previous well-known publication of Tsvetkov et al. (19). Difference analyses of CRGs were conducted between tumor and normal tissues. Wilcoxon test was used to for statistical analysis. Survival and survminer R packages were used for survival analysis, the same as our previous study (26). Kaplan-Meier plot and cox regression analyze were further applied to evaluate the relationships between CRGs expression and patient overall survival (OS). Schoenfeld residuals were used to check the proportional assumption of COX model. Spearman correlation analyses were conducted to explore the interactions among CRGs

ConsensusClusterPlus R package was applied for consensus unsupervised clustering analysis. Patients were grouped into distinct molecular subtypes according to the expression of CRGs, and distinct gene subtypes according to the expression of prognostic differentially expressed genes (DEGs), derived from different molecular subtypes. The criteria included that the samples size in each set was relatively consistent and the cumulative distribution function (CDF) curve increased gradually and smoothly. After consensus clustering analysis, the intra-set association became stronger, while the inter-set association became weaker.

We applied Kaplan-Meier plot and log-rank test to evaluate the associations between different molecular subtypes and patient survival. Correlation analyses between molecular subtypes and clinicopathological features were carried out to learn the clinical values of distinct molecular subtypes by using Chi-square test. The clinicopathological features contained age, gender, grade and tumor node metastasis (TNM) stage.

We downloaded the hallmark gene sets, including C2.CP.KEGG (186 gene sets) and C5.GO.Gene Ontology (10561 gene sets), from the Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb↗). Gene set variation analysis (GSVA) with the above two gene sets was conducted to explore the TME characteristics of different molecular subtypes. The adjusted P-value< 0.05 was considered statistically different. Additionally, the proportion of tumor-infiltrating immune cells (TICs) in tumor samples was calculated by using the deconvolution algorithm, which was also known as CIBERSORT (27). The gene expression signature matrix of TICs was downloaded from CIBERSORT platform (https://cibersortx.stanford.edu/↗). P-value for the deconvolution of each sample was obtained by using Monte Carlo sampling algorithm in R. A CIBERSORT P-value< 0.05 was considered suitable for further analysis. Single sample gene set enrichment analysis (ssGSEA) was used to evaluate the infiltration of TICs in different molecular subtypes.

DEGs of distinct molecular subtypes were acquired by applying limma package in R. The fold change of 1.5 and the adjusted P-value< 0.05 were considered qualified for searching DEGs. Gene Ontology (GO) and Kyoto Encylopedia of Genes and Genomes (KEGG) enrichment analysis of DEGs were carried out by using org.Hs.eg.db, ClusterProfiler, enrichplot, and ggplot2 packages in R. The adjusted P-value< 0.05 was deemed statistically significant.

Firstly, cox regression analyses of DEGs, achieved from different molecular subtypes, were carried out to seek those associated with patients’ prognosis. Secondly, patients were separated into different gene subtypes via consensus clustering analysis of prognostic DEGs expression. Thirdly, patients were randomly divided into the training (n=603) and testing (n=603) sets at a ratio of 1:1. Lastly, CRG Risk scoring system was established in the training set and verified in the testing set, GSE29263, GSE17536, GSE39582 and the combined set. Logistic least absolute shrinkage and selection operator (LASSO) cox regression analysis was carried out by applying Glmnet R package to decrease the risk of over-fitting. Next, we analyzed and cross-validated the varied trajectory of each independent variable. Multivariate Cox analysis was carried out to screen prognostic DEGs in the training group. The Risk score was calculated as follows:

In detail, Expi indicated key prognostic DEGs expression and coefi indicated the coefficient of Risk. Correlation analysis between CRG Risk score and distinct subtypes was also carried out. Survival analysis between high- and low-risk sets was conducted by Kaplan-Meier plot and log-rank test. Receiver operating characteristic (ROC) curves were utilized to learn the sensitivity and specificity of the scoring system. Similarly, all of the testing group, GSE29263, GSE17536, GSE39582 and the combined group were classified into high- and low-risk groups, respectively, and further analyzed by Kaplan-Meier survival curves and ROC curves.

A total of 8 sets of COAD and paired normal tissues were harvested from COAD patients at Nanjing Jiangning Hospital. The study was permitted by the Ethics Committee of Nanjing Jiangning Hospital (2021-03-048-K01). Total RNA extraction and RT-qPCR were performed as our previous study (28). The primers used for RT-qPCR are shown in Table S2. Slides (4μm) of formalin-fixed paraffin-embedded tissue sections were incubated with GLS (1:200; Cell Signaling Technology), NOX1 (1:200; Proteintech), HOXC6 (1:50; Affinity Biosciences), TNNT1 antibody (1:100; Invitrogen). The expression level was scored semiquantitatively based on staining intensity and distribution using the immunoreactive score (IRS) as described (29) and as following: IRS = SI (staining intensity) x PP (percentage of positive cells). SI was determined as 0, negative; 1, weak; 2, moderate; and 3, strong. PP was defined as 0, negative; 1, 1-20% positive cells; 2, 21-50% positive cells; 3, 51-100% positive cells. Ten visual fields from different areas of each sample were selected randomly for the IRS evaluation and the average IRS was calculated as final value.

Difference analyses of CRGs expression levels were carried out between high- and low- Risk groups. Wilcoxon test was used for comparison. We further conducted correlation analyses not only between TICs and risk scores, but also TICs and key prognostic Risk genes. An ESTIMATE algorithm was used to analyze the ratio of immune/stromal components in TME. The Immune Score, Stromal Score and ESTIMATE Score presented the ratio of immune component, the stromal component and the sum of the both, respectively. Difference analyses of Immune/Stromal/ESTIMATE Score were conducted between high- and low- Risk score sets. Wilcoxon test was used for comparison.

Difference and correlation analyses of MSI, TMB and CSC in distinct CRG Risk score groups were conducted to study the underlying associations. Maftools package in R was applied for the comparison of mutation frequency in different Risk score sets.

In order to study effectiveness of drugs in different Risk groups, pRRophetic package in R was used to calculate the semi-inhibitory concentration (IC50) values of drugs.

We applied Rms package in R to establish a nomogram, which combined clinicopathological characteristics, patient survival and CRG Risk score. In the nomogram, a variable matched a score and the scores for all variables were added together to get an overall score. Calibration maps of the nomogram were developed to evaluate the consistency between predicted 1, 3, and 5-year survival rates and actual outcomes. ROC curve was drawn to understand the sensitivity and specificity of the scoring system.

All statistical analyses were conducted by using R version 4.2.1. Statistical significance was set at P-value< 0.05.

We analyzed 10 CRGs in our study, including DLD, DLAT, PDHB, MTF1, PDHA1, FDX1, LIAS, LIPT1, GLS and CDKN2A. Difference analyses showed that 7 of 10 CRGs were dys-regulated in tumor samples compared with those in normal samples, among which LIPT1, PDHA1, GLS and CDKN2A were up-regulated, and FDX1, DLD and MTF1 were down-regulated (Figure 1A).

In order to further study the genetic and transcriptional alterations of CRGs in COAD, we generally analyzed the somatic mutation frequency of CRGs and found 10.13% mutation frequency in tumor samples (Figure 1B). LIPT1, DLD, PDHA1 and LIAS shared the highest mutation frequency (2%), followed by PDHB, MTF1, DLAT and GLS (1%). Both FDX1 and CDKN2A had no mutations in tumor tissues. We further examined CNV frequency in CRGs, among which DLD, CDKN2A, FDX1, DLAT, PDHB and LIAS had elevated copy number loss (Figure 1C). The detailed locations of these CRGs on chromosomes were shown in Figure 1D. As a result, we noted that CRGs had relatively common genetic and transcriptional variations in COAD tissues, which might affect oncogenesis.

To learn the role of CRGs in oncogenesis of COAD, we combined expression patterns of CRGs and clinicopathological information of TCGA-COAD, GSE17536, GSE29623 and GSE39582 by using “Combat” algorithm to eliminate batch effects. PCA indicated that batch differences were well eliminated (Figure 2A). Kaplan-Meier plot revealed 3 of 10 CRGs were closely associated with patients’ OS, among which GLS and CDKN2A were negatively related, while LIAS was positively related (Figures 2B–D). Multivariate Cox regression analyses of CRGs also indicated that both GLS and CDKN2A were closely related with the survival of COAD patients (Table 1). Cuproptosis network generally described the complex interrelations among CRGs and the prognosis of patients with COAD (Figure 2E; Table S3).

Considering the pervasive interrelations among CRGs, we used consensus clustering algorithm to divide patients into three groups based on the expression profile of CRGs. K=3 appeared to be an optimal choice for grouping samples into 3 sets, including molecular subtype A (n=511), B (n=444) and C (n=328) (Figures 3A, S1A–I, Tables S4, 5). Survival analysis revealed that patients in subtype C had the worst prognosis than those in subtype A or B (Figure 3B). The heat-map exhibited the expression profile of 10 CRGs in distinct molecular subtypes (Figure 3C). CDKN2A was obviously up-regulated in molecular subtype C, while PDHA1, FDX1, DLAT, DLD and GLS were greatly elevated in subtype A (Figure 3C). In addition, grade, N, M and stage were found to be significantly associated with cuproptosis molecular subtypes (Figure 3C).

We further performed GSVA enrichment analyses to explore the features of TME in different cuproptosis subtypes. GO GSVA enrichment analysis revealed that molecular subtype A was primarily enriched in messenger ribonucleoprotein complex, regulation of translational initiation by eif2 alpha phosphorylation and phosphatase activity, compared with subtype B (Figure 4A; Table S6). Subtype B was enriched in acyl coa binding, fatty acid derivative binding and acylcoa dehydrogenase activity, compared with subtype C (Figure 4B; Table S6). Subtype C was significantly enriched in embryonic skeletal joint morphogenesis, gap junction and connexin complex, compared with subtype A (Figure 4C; Table S6). Several biological pathways, such as endoplasmic reticulum tubular network organization, cellular response to zinc ion and mrna methylation were recurrent in the comparisons of subtype A and B, A and C, and B and C (Table S7). KEGG GSVA enrichment analysis indicated subtype A mainly participated in TGF-β signaling pathway, riboflavin metabolism and RNA degradation, compared with subtype B (Figure 5A; Table S8). Subtype B was primarily enriched in metabolic related pathways, including fatty acid metabolism, butanoate metabolism, porphyrin and chlorophyll metabolism, compared with subtype C (Figure 5B; Table S8). Subtype C was mainly enriched in glycosphingolipid biosynthesis globo series, glycosaminoglycan biosynthesis chondroitin sulfate and glycosaminoglycan biosynthesis keratan sulfate, compared with subtype A (Figure 5C; Table S8).

Regarding the complex functions of different molecular subtypes in TME, we next conducted ssGSEA between TICs and different subtypes to further identify tumor immune microenvironment (TIME) characteristics of COAD. The ratio of 23 TICs in each tumor sample was presented in Table S9. The result of ssGSEA suggested great difference between the infiltration of 19 TICs and distinct subtypes. In detail, the infiltration levels of eosinophil and plasmacytoid dendritic cell were elevated in subtype A, activated B cell, activated CD8 T cell, activated dendritic cell, monocyte and neutrophil were up-regulated in subtype B, and another 12 TICs were obviously raised in subtype C (Figure 5D).

According to above analyses, we primarily speculated that different subtypes took a different part in TME, especially TIME of COAD.

As the potential role of different molecular subtypes in TME of COAD, we further explore the underlying biological behavior of different subtypes through seeking for DEGs. We identified 114 DEGs derived from subtype A and B, 90 DEGs from subtype A and C, 49 DEGs from subtype B and C (Table S10). Finally, a total of 186 DEGs were obtained for further analyses through combination (Figure 6A; Table S11). GO enrichment analysis demonstrated that 186 DEGs mainly participated in signaling pathways associated with digestion, such as maintenance of gastrointestinal epithelium and digestive system process (Figures 6B, C; Table S12). Univariate Cox regression analysis was performed to seek DEGs of prognostic value and finally identified 86 DEGs associated with patients’ OS, which were analyzed in the following section (Table S13). According to 86 prognostic DEGs expression, consensus clustering analysis was carried out to separate patients into 3 sets, namely gene subtype A (n=310), B (n=729) and C (n=244) (Figures 6D, S2A–I; Tables S14, 15). Distinct gene subtypes showed great differences in the expression levels of both prognostic DEGs and 8 CRGs (FDX1, LIPT1, DLD, PDHA1, PDHB, MTF1, GLS and CDKN2A) (Figures 6E, F; Tables S16, 17). In addition, cuproptosis gene subtypes were closely related with age, gender, grade, and T and N stage of COAD patients (Figure 6E). Survival analysis revealed that patients of gene subtype B had a better prognosis, compared with those of subtype A or C (Figure 6G).

To study the prognostic value of CRGs in COAD, we further constructed CRG Risk scoring system based on different molecular and gene subtypes. First, we applied “caret” package in R to randomly separate COAD patients into the training (n=603) and testing (n=603) groups at a ratio of 1:1. The clinicopathological characteristics of patients in the training and testing group were consistent (Table S18). Second, LASSO and multivariate Cox analyses were conducted to identify optimum prognostic signature based on 86 DEGs expression (Figure S3). Finally, CRG Risk scoring system was established through multivariate Cox regression analysis in the training set, the formula was as follow: Risk score = (0.30346935571892* expression of GLS) + (0.285346929484159 * expression of CAB39L) + (-0.171967289741126* expression of NOX1) + (0.149406405352724 * expression of HOXC6) + (0.128828618079011 * expression of TNNT1) + (-0.305462961248901* expression of ASRGL1) + (-0.142788274274145* expression of PLA2G12B). We classified patients into two groups, namely high- and low-Risk score sets, according to the calculation of Risk score in each tumor sample. Figure 7A presented the specific classifications of patients in the training set, including three cuproptosis molecular subtypes, three gene subtypes and two CRG Risk score sets. The detailed information of 7 key cuproptosis-related risk genes, Risk score and survival features in training and testing groups was displayed in Tables S19, 20. The results of difference analyses indicated that all of the expression of GLS, NOX1, HOXC6, TNNT1 and PLA2G12B were increased in tumor tissues, compared with those in normal tissues (Figure S4). Among these five genes, GLS and HOXC6 were negatively associated with patients’ survival, while NOX1 and PLA2G12B were positively related. RT-qPCR and IHC indicated the same result (Figures 7B–G). Difference analyses in the training set showed Risk score was extremely increased in both molecular subtype C and gene subtype C and decreased in both molecular subtype B and gene subtype B (Figures 7H, I). The heat-map presented a great difference of 7 key Risk score gene expression profile between high- and low-Risk score sets in the training group (Figure 7J). The scattergram of patients’ survival in different Risk score groups revealed that COAD patients’ survival got worse, while Risk score increased (Figure 7K), which was also proven by Kaplan-Meier survival curves (Figure 7L). In addition, area under the time-concentration curve (AUC) values of 1-, 3-, and 5-year survival rates of CRG Risk score in the training set were 0.693, 0.706, and 0.703, respectively, signifying both relative high sensitivity and specificity (Figure 6M).

To verify the accuracy of the scoring system, we further calculated Risk score according to the above Risk score formula, in the testing group, individual GSE17536, GSE29623, GSE39582, GSE40967, respectively (). Patients were respectively divided into distinct cuproptosis molecular subtypes, gene subtypes and Risk score sets, the same as which in the training set (). Risk score showed a great difference in both molecular subtypes and gene subtypes of the testing group, individual GSE17536, GSE29623, GSE39582 (). The expressions of 7 key Risk scoring genes in different Risk group were shown in;, respectively. Both scattergram and Kaplan-Meier survival curves showed that high Risk score predicted poor survival in testing group, individual GSE17536, GSE39582 and GSE40967 (). However, in GSE29623, survival analysis revealed that Risk score was not associated with patients’ survival, which might be related with the small sample size (). We further plot ROC curves to confirm the sensitivity and specificity of the scoring system and found relatively high AUC values in the cohorts of validation, indicating the system as an accurate predictor for patients’ survival (). 1 1 1 1 1 1 1 1

Difference analyses of CRGs indicated that 6 CRGs showed a great difference in distinct Risk score sets. To be specific, GLS and CDKN2A expression were increased, while DLD, DLAT, PDHA1 and PDHB expression were decreased in high-Risk score group, compared with those in low-Risk group (Figure 8A). In order to learn the relationships between CRG Risk score and TICs in TME of COAD, correlation analyses were carried out and suggested that CRG Risk score was positively associated with activated NK cells, memory B cells, eosinophils, M0 macrophages, M1 macrophages, M2 macrophages, and neutrophils, while negatively associated with CD8 T cells, regulatory T cells (Tregs), naïve B cells, resting dendritic cells, plasma cells and CD4 memory resting T cells (Figures 8B–N). Furthermore, all of immune, stromal and estimate score were higher in high-Risk score set than those in low-Risk score set (Figure 8O). Most immune cells were greatly associated with seven prognostic genes (Figure 8P). Consequently, CRG Risk score might be associated with TME of COAD.

Up to data, limited molecular markers are available to lead therapeutic decisions for patients with COAD, among which MSI, CSC, TMB and somatic mutations appeared to be the most promising. An increasing number of research revealed that patients with high microsatellite instability (MSI-H) tumor might benefit from immune checkpoint inhibitors (ICIs) in COAD (30, 31). As a result, we assessed the MSI status and found that in the low-risk group, 73% were MSS, 17% were low microsatellite instability (MSI-L), and 10% were MSI-H, while in the high-risk group, 59% were MSS, 20% were MSI-L, and 20% were MSI-H (Figure 9A). The results indicated that patients with high-risk shared a higher MSI-H frequency. Figure 9B suggested that patients bearing MSI-H tumors appeared to have a higher Risk score, compared with those with MSS. This might be related with better treatment outcomes of ICIs. Additionally, crosstalk between immune cells and CSCs, another important indicator of TIME, takes a great part in tumor progression (32). As presented in Figure 9C, CRG Risk score was negatively associated with CSC index, indicating COAD cells with high CRG Risk scores had less difference in stem cell properties and higher cell differentiation than those with low-risk scores. TMB, as an indicator of the number of tumor mutations, is known to be closely associated with patients’ immunotherapy benefits (33). Differential analysis indicated that TMB in high-risk group was significantly higher than that in low-risk group (Figure 9D). Correlation analysis also suggested that TMB was positively associated CRG Risk score (Figure 9E). Maftools of somatic mutations showed that the top 10 mutant genes in the high-risk and low-risk groups were APC, TP53, TTN, KRAS, PIK3CA, SYNE1, MUC16, FAT4, RYR2 and ZFHX4, respectively (Figures 9F, G).

To investigate the predictive value of CRG Risk score in drug sensitivity, we used pRRophetic R package to calculate the IC50 values of various drugs (Figures S10, 11; Table 2). Both drugs under clinical use and clinical trials were included in our analyses. Various drugs were divided into different groups, such as AKT inhibitor, AMPK activator, Bcr-Abl inhibitor, BTK inhibitor, EGFR inhibitor, MAPK inhibitor, mTOR inhibitor, TrkA inhibitor, Topoisomerase inhibitor, Microtubule assosiated inhibitor, XIAP inhibitor, TNF inhibitor and so on. In particular, patients of low-Risk score set showed increased IC50 value for AMPK activator (AICAR), Bcl-2 inhibitor (TW.37, Obatoclax.Mesylate and ABT.263), BRAF inhibitor (PLX4720), c-Kit inhibitor (AMG.706), DNA Synthesis inhibitor (Cytarabine, Bleomycin and Gemcitabine), HSP90 inhibitor (AUY922), ITK inhibitor (BMS.509744), MEK inhibitor (CI.1040 and RDEA119), PARP inhibitor (AG.014699 and AZD.2281) and ROCK inhibitor (GSK269962A). In addition, patients of high-Risk score set showed increased IC50 value for AKT inhibitor (AKT.inhibitor.VIII and A.443654), CDK inhibitor (Roscovitine), Raf/VEGFR/c-Kit inhibitor (Sorafenib), Her-2 inhibitor (Lapatinib) and EGFR inhibitor (Erlotinib and BIBW2992).

However, drugs that target the same site may have opposite effects in different risk groups. For example, patients with low CRG risk scores had increased IC50 values for Aurora kinase inhibitors (ZM.447439) and decreased IC50 values for aurora kinase inhibitors (VX.680). HDAC inhibitors (Vorinostat) had increased IC50 values and HDAC inhibitors (MS.275) had decreased IC50 values in the low-risk score set. mTOR inhibitors (Temsirolimus, NVP.BEZ235 and AZD8055) presented better drug sensitivity in low-risk score set, while mTOR inhibitor (rapamycin) had the opposite. RSK inhibitors (BI. D1870 and CMK) and PF.4708671. RSK inhibitor (BI.D1870 and CMK) and PF.4708671 also showed the opposite drug susceptibility between different risk score sets.

Regarding the important role of Risk score in patients’ survival, we constructed a nomogram combining CRG Risk scores and clinicopathological characteristics to predict 1, 3, and 5-year survival rates of COAD patients (Figure 10A). The calibration graph showed that the nomogram functioned well in predicting patients’ survival compared to an ideal model (Figure 10B). The AUC values of 1, 3, and 5-year survival rates of nomogram were 0.873, 0.798, and 0.804, respectively, suggesting both relatively high sensitivity and specificity (Figure 10C).