Identification of cuproptosis-related subtypes, construction of a prognosis model, and tumor microenvironment landscape in gastric cancer

We downloaded RNA-Sequence data and the corresponding clinicopathological data of STAD project for stomach adenocarcinoma from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/↗), and expression levels in TPM of genes were extracted and combined for all samples. Then, series matrix file of GSE84433↗ was downloaded from the Gene Expression Omnibus database (GEO) (https://www.ncbi.nlm.nih.gov/geo/↗). The TCGA cohort contained 375 GC samples and 32 normal samples from 375 patients, and the GEO cohort (GSE84433↗) contained 357 GC samples. The detailed clinicopathological information on these GC patients is presented in Table S1. The GEO dataset was combined with TCGA-STAD dataset. Before conducting subsequent analyses, we eliminated batch effects by using “Combat” algorithm. Additionally, somatic mutation data were downloaded from TCGA and contained 431 GC samples. Somatic copy number alteration data were downloaded from TCGA and contained 440 GC samples.

Nineteen CRGs, listed in Table S2, were achieved from previous cuproptosis-related publications (17, 28). Consensus unsupervised clustering analysis was carried out by using “ConsensusClusterPlus” package in R to classify GC patients into different molecular subtypes according to the expression of CRGs. Furthermore, DEGs, derived from different molecular subtypes, were grouped into different genomic subtypes by the same way. The criteria were as follows: First, the number of samples in each group was relatively consistent. Second, the cumulative distribution function (CDF) curve rose gradually and smoothly. Third, after clustering, the intra-group link was stronger, while the inter-group link was weaker.

To explore the clinical value of different molecular subtypes, we performed correlation analysis between molecular subtypes and clinicopathological characteristics of GC patients. The clinicopathological characteristics included age, gender, grade and tumor node metastasis (TNM) stage. Survival and survminer packages in R were used for survival analysis, the same as our previous research (29). We used Kaplan–Meier plot and log-rank to test the correlations between molecular subtypes and overall survival (OS) of GC patients.

In order to verify the characteristics of TME in distinct molecular subtypes, we conducted gene set variation analysis (GSVA) with the hallmark gene set (C2.CP.KEGG (186 gene sets) and C5.GO.Gene Ontology (10561 gene sets)), achieved from the MSigDB database (https://www.gsea-msigdb.org/gsea/msigdb↗). Adjusted P-value <0.05 was recognized to be statistically significant. Gene Set Enrichment Analysis (GSEA) (4.1.0) based on different gene sets, was applied to learn the specific functional profile of different molecular subtypes. The absolute value of normalize enrichment score (NES) >1, nominal p value<0.05, FDR<0.25 were considered to be statistically significant. In addition, the deconvolution algorithm (referred to as CIBERSORT) was used to calculate the abundance of tumor-infiltrating immune cells (TICs) in each GC sample (30). The gene expression signature matrix of TICs was downloaded from the CIBERSORT platform (https://cibersortx.stanford.edu/↗). The matrix data of gene expression levels in TCGA-STAD and GSE84433↗ cohorts were compared with those of the signature matrix of TICs to generate a proportion matrix for the TICs in GC tissues. We further applied Monte Carlo sampling algorithm to obtain a p value for the deconvolution of each sample, which providing a measure of confidence for the obtained data. A CIBERSORT p<0.05 were deemed qualified for further analysis. The levels of TICs in each GC sample were also calculated by using a single-sample gene set enrichment analysis (ssGSEA) algorithm.

Package “limma” in R was used to identify DEGs between different cuproptosis molecular subtypes. A fold-change of 1.5 and an adjusted p-value of <0.05 were set up to screen DEGs. “ClusterProfiler”, “org.Hs.eg.db”, “enrichplot”, and “ggplot2” packages in R were used for gene ontology (GO) enrichment analysis of DEGs. Adjusted p value <0.05 was recognized to be statistically significant.

CRG Risk scoring system was established to identify the cuproptosis patterns of the individual tumors. First, DEGs, screened out from distinct cuproptosis molecular subtypes, were subjected to univariate Cox regression analysis to seek those related to GC patients’ survival. Second, we divided patients into different cuproptosis gene subtype, including subtype A and B, through consensus clustering analysis based on the expression of prognostic DEGs. Third, “caret” package in R was used to randomly categorized all GC patients from TCGA-STAD and GSE84433↗ database into training (n=364) and testing (n=364) groups at a ratio of 1:1. Lastly, we constructed CRG Risk scoring system in the training group and further validated the system in the testing group and the combined group. In detail, we applied “glmnet” R package to conduct logistic least absolute shrinkage and selection operator (LASSO) Cox regression analysis to minimize the risk of over-fitting. The varied trajectory of each independent variable was analyzed and cross-validated to establish the model. Multivariate Cox analysis was applied again to seek candidate cuproptosis-related Risk genes and establish prognostic CRG Risk scoring system in the training set. The CRG Risk score was calculated as follows: CRG Risk score = Σ(Expi * coefi). Expi presented the expression of key cuproptosis-related Risk gene, and coefi presented the Risk coefficient. Correlation analysis was applied to evaluate the relationship between Risk score and different molecular or gene subtypes. A total of 364 GC patients in the training group were classified into high- (n=182) and low-risk (n=182) sets according to the median Risk score. Kaplan–Meier survival analysis was further used to identify the survival difference between high- and low-risk sets.

Similarly, both the testing and combined sets were grouped into high- and low-risk sets, each of which was subjected to survival analysis and the generation of receiver operating characteristic (ROC) curves.

Ten groups of GC and corresponding normal tissues were harvested from GC patients at Nanjing Jiangning Hospital. The study was permitted by the Ethics Committee of Nanjing Jiangning Hospital. RNA isolation and RT-qPCR were carried out as our previous description (31). The primer sequences used for qRT-PCR in this study are listed in Table S3.

We identified the expression of CRGs in high- and low- Risk score groups through boxplots, and further calculated the abundance of TICs in TME of each GC sample by applying CIBERSORT in R. Correlation analyses were also carried out to study the relationships between TICs and prognostic Risk genes.

We analyzed the associations of MSI, CSC, and TMB in two Risk groups. Mutation frequency analyses of high- and low- Risk groups were conducted by using the”maftools” R package.

In order to investigate the therapeutic effects of drugs in the two groups, we calculated the semi-inhibitory concentration (IC50) values of drugs using “pRRophetic” package in R.

The clinicopathological features and CRG Risk score were incorporated to develop a nomogram using the “rms” package, based on patients’ survival. In the nomogram scoring system, a variable, such as gender, age, TNM stage and CRG Risk score, was matched with a score, and the total score was obtained by adding the scores across all variables of each sample. The subsequent calibration graph of the nomogram scoring system was performed to examine the predictive value between the predicted 1-, 3-, and 5-year survival rates and the virtually outcomes.

All statistical analyses were performed using R version 4.2.1. Statistical significance was set at p < 0.05.

A total of 19 CRGs were analyzed in our following study, such as NFE2L2, NLRP3, ATP7B, ATP7A, SLC31A1, FDX1, LIAS and so on (Table S2). We first compared the expression of CRGs between tumor samples and normal samples from TCGA-STAD database, and found that 16 CRGs were up-regulated in tumor tissues, including NLRP3, ATP7A, ATP7B, SLC31A1, FDX1, LIAS, LIPT1, LIPT2, DLD, DLAT, PDHA1, PDHB, MTF1, GLS, CDKN2A, and GCSH (Figure 1A).

Next, we conducted general analysis of the somatic mutation frequency in these 19 CRGs, the result showed a relatively high mutation frequency in GC samples from TCGA-STAD database (Figure 1B). In detail, 99 (22.97%) GC samples had mutations in CRGs. LIPT1 had the highest mutation frequency (6%), followed by CDKN2A, NLRP3, ATP7A, ATP7B, DLAT, DLD, MTF1, GLS, DLST, LIAS, NFE2L2, PDHB and DBT. The other five CRGs, including FDX1, SLC31A1, LIPT2, PDHA1 and GCSH, did not have any mutations in tumor samples. Furthermore, we calculated somatic copy number alterations in these CRGs and found that copy number alterations were pervasive in all 19 CRGs. Among them, NLRP3, LIPT2, ATP7B, SLC31A1, GLS and MTF1 had relatively elevated copy number variation (CNV), while CDKN2A, DLAT, FDX1, DBT and PDHB showed relatively decreased CNV (Figure 1C). Detailed locations of CNV alterations on chromosomes were presented in Figure 1D.

From above, we noted that CRGs in GC tissues had prevalent genetic and transcriptional alterations, which might have their roles in GC oncogenesis.

In order to explore the expression pattern of CRGs involved in GC tumorigenesis, we collected 732 GC patients from TCGA database (TCGA-STAD) and GEO database (GSE84433↗) for further analyses. Detailed clinicopathological information of patients was presented in Table S1. Survival analysis revealed that 9 CRGs (ATP7A, DLAT, DLD, FDX1, LIAS, LIPT1, MTF1, NLRP3 and SLC31A1) were significantly associated with overall survival (OS) of GC patients (Figures 2A–I, p<0.05). The result of univariate Cox regression analysis on CRGs showed that both LIAS and FDX1 were significantly associated with GC patients’ survival (Table 1). The intersection between survival analysis and multivariate Cox regression analysis indicated that both LIAS and FDX1 were significantly associated with the prognosis of GC patients. Next, a cuproptosis network was carried out to comprehensively demonstrate the association among CRGs and their prognostic value in GC patients (Figure 2J; Table S4). The network indicated there were prevalent and complicated interactions among CRGs.

Regarding the comprehensive associations among CRGs, we categorized GC patients into two groups based on the expression profiles of CRGs by using a consensus clustering algorithm. The results indicated that k=2 might be an optimal selection for clarifying patients into 2 groups, including molecular subtype A (n=339) and B (n=393) (Figure 3A, Figures S1A–I; Table S5). PCA analysis verified that there were significant differences in the cuproptosis related transcription profiles between subtype A and B (Figure 3B). Survival analysis indicated that GC patients of subtype A had a higher survival probability than those in subtype B (log-rank test, p = 0.014; Figure 3C). In addition, the associations between CRGs expression and the clinical features, such as gender, age, grade and TNM stage, were profiled according to different molecular subtypes of GC (Figure 3D). The heat-map showed that a group of CRGs were upregulated in cluster A, such as FDX1, GLS, SLC31A1, LIAS and so on. And gender, age, and T stage were significantly correlated with GC patients’ cuproptosis subtypes (Figure 3D, p<0.05).

In order to further identify the characteristics of TME in distinct subtypes, we conducted not only GO and KEGG GSVA enrichment analysis, but also GSEA enrichment analysis. The results of GO GSVA enrichment analysis displayed that subtype A was significantly enriched in metabolic related pathways, including TCA cycle, peroxisomal organization and transportation, mitochondrial membrane organization and transportation, regulation of mitochondrial gene expression, protein transmembrane import, amino acid metabolic process, and so on (Figure 4A; Table S6). GSEA enrichment analysis suggested that for C5 collection, the gene ontology sets, genes in molecular subtype A were also enriched in the above metabolic related pathways (Figure 4C; Table S7). KEGG GSVA enrichment analysis also found that subtype A was primarily enriched in metabolic related pathways, such as TCA cycle, glycan biosynthesis, ubiquitin mediated proteolysis, peroxisome, RNA degradation, cysteine, methionine, glyoxylate, dicarboxylate, pyruvate, butanoate and selenoamino acid metabolism, valine leucine and isoleucine degradation and so on (Figure 4B; Table S8). GSEA enrichment analysis indicated that for C2 collection, the KEGG gene sets database, genes in molecular subtype A were also primarily enriched in the above metabolic related pathways (Figure 4D; Table S9). In particular, both GO and KEGG GSVA enrichment analysis found that TCA cycle, peroxisome, amino acid metabolism were obviously enriched in subtype A, and TCA cycle has been proven to be closely associated with cuproptosis (17).

We further explored whether CRGs were involved in TIME of GC through correlation analysis between two subtypes. Human immune cell subsets of each GC sample were calculated by using CIBERSORT algorithm (Table S10). We observed significant differences in the infiltration of most immune cells between subtype A and B through ssGSEA (Figure 4E). In detail, the infiltration levels of activated B cell, activated CD8 T cell, activated dendritic cell, immature B cell, immature dendritic cell, myeloid-derived suppressor cell (MDSC), macrophage, mast cell, monocyte, natural killer T cell, natural killer cell, plasmacytoid dendritic cell, regulatory T cell, T follicular helper cell and Type 1 T-helper cell were obviously lower in subtype A than those in subtype B (Figure 4E). From above, we primarily speculated that subtype A was enriched in metabolic related pathways, and subtype B was closely related with TIME.

To further investigate the underlying biological behavior of different cuproptosis molecular subtypes, we identified 84 DEGs between subtype A and B (Table S11). GO enrichment analysis was carried out to seek related biological pathways. The result showed that DEGs primarily enriched in digestive system development, extracellular matrix organization, extracellular structure organization, intermediate filament cytoskeleton and organization, arginine metabolic process and so on (Figures 5A, B; Table S12). According to the above enrichment analysis, we speculated that cuproptosis played an important role in the regulation of metabolism and intracellular/extracellular structure in GC. We further performed univariate Cox regression analysis to identify the prognostic value of 84 subtype-related DEGs and finally screened out 32 genes related to OS (p < 0.05), which were analyzed in the following section (Table S13)

Next, we performed consensus clustering analysis of 32 prognosis related DEGs to validate this regulation mechanism. Patients were divided into two genomic subtypes, namely, gene subtype A and B (Figure 5C; Figure S2; Table S14). The two cuproptosis gene subtypes presented significant differences in the expressions of CRGs, consistent with the expected results of the cuproptosis patterns (Figure 5E). In addition, cuproptosis gene subtype was significantly correlated with patients’ age, grade, and T and N stage (Figure 5D, p<0.05). The result of survival analysis showed that patients in gene subtype B had a higher survival probability than those in gene subtype A (Figure 5F, p<0.001).

We established CRG Risk scoring system according to the expression of prognostic DEGs derived from distinct molecular subtypes. As shown in Figure 6A, GC patients were grouped into two cuproptosis molecular subtypes, two gene subtypes, and two CRG Risk score groups. First, we randomly classified GC patients into training (n=364) and testing (n=364) groups at a ratio of 1:1, by using “caret” package in R (Tables S15, 16). LASSO and multivariate Cox analyses of 32 prognostic DEGs were used to seek optimum prognostic signature (Figure S3). Then multivariate Cox regression analysis was carried out to construct CRG Risk scoring system in the training sets: Risk score = (-0.145406843672574 * expression of SLC27A2) + (-0.100834783610163 * expression of NAT2) + (0.108030183332151* expression of TAGLN) + (0.0800602829915149* expression of SFRP2) + (0.0876442662604965 * expression of KRT17). We calculated Risk score of each GC samples in both molecular subtypes and gene subtypes, and found that Risk score was significantly elevated in molecular subtype B and gene subtype A, compared with that in molecular subtype A and gene subtype B, respectively (Figures 6B, C). Next, the expressions of five key cuproptosis-related risk genes in the training sets were profiled in Figure 6E, based on CRG Risk score (Table S17). The results indicated that the expression of five key Risk scoring genes showed a great difference between high- and low- risk sets (Figure 6E). The levels of five cuproptosis-related Risk genes were also measured in GC tissues and adjacent normal tissues by RT-qPCR. As shown in Figure 6D, the expression of SLC27A2, SFRP2, KRT17 were up-regulated in tumor tissues (p<0.05), the expression of TAGLN was down-regulated (p<0.05), while those of NAT2 remained unchanged (p>0.05), compared with the levels in the corresponding normal tissues.

The scattergram of CRG Risk score in the training sets revealed that patients’ survival time decreased while CRG Risk score increased (Figures 6F, G). The Kaplan–Meier survival curves showed that GC patients with low Risk scores had a better overall survival compared to that in patients with high Risk scores (Figure 6H, p <0.001). Additionally, the 1-, 3-, and 5-year survival rates of CRG Risk score were represented by area under the time-concentration curve (AUC) values of 0.621, 0.654, and 0.670, respectively, indicating both moderate sensitivity and specificity (Figure 6I). The CRG Risk score predicted 1- year survival with a 70% specificity and 50% sensitivity, 3- year survival with a 42% specificity and 83% sensitivity, and 5-year survival with a 61% specificity and 65% sensitivity.

To validate the accuracy of the CRG Risk scoring system, we calculated CRG Risk score in the testing group, and combined TCGA-STAD and GSE84433↗ group (Figures S4, S5; Tables S18, 19). GC patients were stratified into high- and low-risk sets, the same as which in the training set. The expression of five key Risk scoring genes in the testing set and the combined set, were presented in Figure S4A and S5A, respectively. The relationships between patients’ survival and CRG Risk score were shown in Figures S4B, C, S5B, C. Survival analyses presented that GC patients with low CRG Risk scores had a significantly favorable overall survival compared to those in patients with high scores, which was the same in the training group (Figure S4D, p=0.002; Figure S5D, p<0.001). We also conducted prognostic prediction classification efficiency analysis and found that CRG Risk score still had relatively high AUC values (Figures S4E, S5E), suggesting that the CRG Risk scoring system was suitable to accurately predict the survival of GC patients.

In order to explore the relationship of TME and CRG Risk score, we first analyzed the expression of CRGs in both high- and low- CRG Risk score groups, and found that ten CRGs were significantly related with CRG Risk score. To be specific, ATP7B, SLC31A1, FDX1, LIAS, DLD, DLAT, PDHA1, GCSH and DLST were down-regulated in high-Risk score group, while NLRP3 were up-regulated in high-Risk score group (Figure 7A). Next, we performed correlation analyses to clarify the relationship between the abundance of immune cells and CRG Risk score, through applying CIBERSORT algorithm. As shown in Figures 7B–J, CRG Risk score was positively related with resting Mast cells, activated natural killer (NK) cells, M2 Macrophages and monocytes, while negatively associated with follicular helper T cells, activated memory CD4 + T cells, plasma cells, resting NK cells and M0 macrophages. A high CRG Risk score was positively associated with both immune and stromal score (Figure 7K). In addition, we studied the association between the abundance of immune cells and the five key cuproptosis-related Risk genes. The results showed that the majority of immune cells were significantly related with the five genes (Figure 7L). As a result, CRG Risk score may be related with tumor TIME of GC.

Accumulating evidence has implied that MSI is a potential genomic biomarker to identify patients’ sensitivity to immunotherapy (32). We assessed the MSI status in distinct sets of CRG Risk score. In the high-score group, MSI-H accounted for 14%, MSI-L accounted for 16%, the rest 70% were microsatellite stable (MSS) (Figure 8A). However, in the low-score group, the proportion of MSI-H was significantly increased (22%), and the proportion of MSI-L was decreased (14%) (Figure 8A). Further correlation analyses indicated that a low CRG score was significantly related with MSI-H status, while a high CRG score was correlated with MSS status (Figure 8B). This might be associated with better efficacy of immunotherapy.

CSCs have been recognized as promising therapeutic targets for cancer therapy according to their self-renewal capacity and differentiation potential (33). As a result, we studied the correlation between CRG Risk score and CSC index values. Figure 8C showed the results of the negative linear correlation between CRG Risk score and CSC index (R = −0.64, p <.001), suggesting that GC cells with low CRG Risk score had more different stem cell properties and a lower degree of cell differentiation (Figure 8C).

TMB, reflecting cancer mutation quantity, is also clinically related with immune checkpoint inhibitors (ICIs) outcomes (34). Patients with a higher TMB usually benefited from ICIs. Our analysis of the mutation data from TCGA-STAD cohort demonstrated that lower TMB was observed in the sets of high CRG Risk score than that in the sets of low CRG Risk score (p<0.001; Figure 8D). Spearman correlation analysis discovered that TMB was negatively associated with CRG Risk score (R = −0.24, p = 3.4e-06; Figure 8E). The above analyses indicated that the low-risk set might benefit from ICIs. We further described the distribution variations of the somatic mutations between two CRG Risk score sets in TCGA-STAD cohort through maftools. The top ten mutated genes in the high- and low-CRG Risk sets were TTN, TP53, MUC16, LRP1B, ARID1A, SYNE1, CSMD3, FAT4, FLG, PCLO and ZFHX4 (Figures 8F, G). Patients with a low CRG Risk score had obviously higher frequencies of all these mutations, except CSMD3, compared to those in patients with a high CRG Risk score.

To investigate the therapeutic effects of drugs in patients of the two groups, we applied “pRRophetic” package in R to calculate the IC50 values of drugs. We assessed not only the current clinical used chemotherapy drugs, but also drugs under clinical trials. We classified drugs into multiple groups, including AKT inhibitors, BCL-2 inhibitors, PI3K inhibitors, MEK inhibitors, ROCK inhibitors, XIAP inhibitors, Raf inhibitors, Multikinase inhibitors, and so on. Interestingly, we found that the patients in the low CRG Risk score group had higher IC50 value for AKT inhibitors (MK.2206), ALK inhibitor (NVP.TAE684), proteasome inhibitor (Bortezomib, MG.132), ATP inhibitor (PF.562271), Bcl-2 inhibitor (ABT.263, TW.37), Bcr-Abl inhibitor (Imatinib, GNF.2), BTK inhibitor (LFM.A13), CDK inhibitor (CGP.60474, PD.0332991), Chk1 inhibitor (AZD7762), DNA synthesis inhibitor (Bleomycin), FGFR inhibitor (PD.173074), FTase inhibitor (FTI.277), GSK-3 inhibitor (CHIR.99021, SB.216763), HIF-PH inhibitor (DMOG), IGF-1RIR inhibitor (BMS.536924, BMS.754807), ITK inhibitor (BMS.509744), JNK inhibitor (AS601245, JNK.9L, JNK.Inhibitor.VIII), PKC Modulator (Bryostatin.1, Midostaurin), MEK inhibitor (RDEA119), MET inhibitor (PF.02341066, PHA.665752), mTOR inhibitor (Temsirolimus, AZD8055), mTOR and PI3K inhibitor (NVP.BEZ235), ATM inhibitor (KU.55933), Bcr-Abl inhibitor (Nilotinib), Doxorubicin, Elesclomol, Docetaxel, Lck/Src inhibitor (KIN001.135), MDM2 inhibitor (JNJ.26854165), PDK1 inhibitor (BX.795), PI3K inhibitor (AZD6482, GDC0941), PK inhibitor (NU.7441), RXR activator (Bexarotene), VEGFR inhibitor (AMG.706), ROCK inhibitor (GSK269962A), RSK inhibitor (CMK), Src inhibitor (A.770041, AZD.0530, WH.4.023), WIP1 inhibitor (CCT007093) and XIAP inhibitor (Embelin). However, patients in the low CRG Risk score group had lower IC50 value for eIF2α dephosphorylation inhibitor (Salubrinal). methotrexate, TrkA inhibitor (GW.441756), TNF inhibitor (Lenalidomide), PLK inhibitor (GW843682X), Rac1 inhibitor (EHT.1864), Aurora inhibitor (VX.680), Mitomycin.C, DHFR inhibitor (Pyrimethamine), RSK inhibitor (PF.4708671) and HDAC inhibitor (Vorinostat, MS.275). Furthermore, the same type of drugs might act different roles in different Risk score groups. For instance, patients in the low CRG Risk score group had higher IC50 value for Raf inhibitor (PLX4720, AZ628), and lower IC50 value for Raf inhibitor (SB590885). In the low CRG Risk score group, MAPK inhibitor (VX.702) displayed a higher IC50 value, while MAPK inhibitor (BIRB.0796) showed a lower IC50 value. Multikinase inhibitors (Dasatinib, Pazopanib, AP.24534, CEP.701) exhibited a better drug susceptibility in the group of low CRG Risk score, while multikinase inhibitor (Sorafenib) showed the opposite effect. The opposite drug susceptibility in low and high CRG Risk score groups were also seen between PARP inhibitor AZD.2281, AG.014699 and ABT.888, EGFR/Her-1/2 inhibitor Lapatinib and BIBW2992, HSP90 inhibitor AUY922 and CCT018159. Together, these results showed that CRGs were significantly related to drug sensitivity (Figures S6, 7; Table 2).

As the importance of CRG Risk score in GC patients’ survival, we established a nomogram incorporating the CRG Risk score and clinicopathological features to predict the 1-, 3-, and 5-year OS rates (Figure 9A). Clinicopathological features contained gender, age and TNM stage. The subsequent calibration graph indicated that the proposed nomogram had a similar performance in GC patients compared to an ideal model (Figure 9B).

The original contributions presented in the study are included in the article/. Further inquiries can be directed to the corresponding author. 1