We obtained data of clinical data and sample information of LUAD patients from the TCGA database (https://portal.gdc.cancer.gov/) using the R package TCGAbiolinks (19). The LUAD transcriptome fragments per kilobase million (FPKM) data from TCGA were downloaded from the UCSC Xena browser (GDC hub: https://gdc.xenahubs.net), which contained 507 LUAD tissues and 59 normal lung tissues. For TCGA-LUAD data, in order to eliminate the error caused by the quantitative mRNA abundance of FPKM in multiple samples, we converted FPKM to transcripts per kilobase million (TPM) values for standardization (20). The GSE68465, GSE11969, GSE72094 data sets were downloaded from the Gene Expression Omnibus database (http://www.ncbi.nlm.nih.gov/geo) (21, 22). The data were uniformly pre-processed using the Robust Multichip Average algorithm for background correction, quantile normalization, and log2-transformation (23). The probes were converted into corresponding genes using the annotation information available from the gene platforms (GPL96-57554).

In order to investigate the expression pattern of CRGs in the LUAD patients’ tumor tissues and adjacent nontumor tissues, as well as different clinical subgroups, we used Wilcoxon’s method to calculate differential genes in clinical subgroups of immune subtype, epidermal growth factor receptor (EGFR) mutation, age, gender, pathological combined stage, pathological nodal (pN) stage, pathological tumor (pT) stage, pathological metastasis stage, chemotherapy, echinoderm microtubule-associated protein-like 4- anaplastic lymphoma kinase (EML4-ALK) fusion, and pathological type (24). Age was a continuous variable, and six years old could be selected as the interval to change it into a categorical variable. Kruskal test was selected for multi-category clinical variables. To investigate the differences in CRGs in biological and function processes, GSVA was performed by R packages ‘GSVA’ (version 1.34.0) (25, 26). The gene sets ‘h.all.v7.5.1.symbols.gmt hallmark’ in the Molecular Signatures Database (MSigDB) was selected as the reference gene set and a p-value of <0.01 was considered the threshold (27).

Somatic mutation and copy number variation (CNV) data were downloaded from the TCGA database. The high-frequency mutations in CRGs were visualized using OncoPlot. Based on Affymetrix SNP6.0 array copy number segmented data, GISTIC2.0 was used to analyze the change of regions of CNV in different groups. Deep deletions (GISTIC2 value: −2) were defined as homozygous losses in comparison to shallow deletions (GISTIC2 value: −1), which resulted in heterozygous loss. CNA gains were defined as GISTIC2 value +1, while amplifications were defined as GISTIC2 value +2. The threshold was Q-value<0.25, and the confidence level was 0.90. We used the ‘maftools’ package to provide visualization of regions of copy number variation across high and low group samples (28). Chromosome locations of copper death genes were represented by chromosome ring diagrams. The mutation frequency was plotted using the ‘ggplot’ package. The genomic location information of CRGs was downloaded from Gencode.v29.annotation.

We downloaded the somatic mutation file and calculated each patient’s tumor mutation burden (TMB) score. The influence of TMB on patient OS was evaluated by Kaplan–Meier analysis and compared between the high and low groups by t-test. Maftools R package was used.

Based on the CRGs obtained above, a univariable Cox regression model was used to screen the DEGs associated with good prognosis (overall survival, OS; progression free interval, PFI), which was performed using the ‘survival’ package in R with p<0.05 as the threshold. The median was used as a cutoff value to separate high and low expression groups. Log-rank test was performed, and the Kaplan-Meier survival analysis was used for further survival analysis. 18 CRGs were obtained. The R package “survival” and “survminer” were employed for survival analysis and drawing.

Based on TCGA data, the expression matrix of CRGs was obtained and the interaction between CRGs was calculated. The correlation coefficient and P value were calculated by Spearman method. Based on the survival analysis results obtained above, Cytoscape was used to map the interaction network of CRGs.

Unsupervised clustering of TCGA samples was performed using the ConsensusClusterPlus algorithm based on the expression level of CRGs to identify different copper modification patterns. A consensus clustering algorithm was applied to determine the number of clusters and the stability of the discovered clusters. Euclidean distance calculation of similarity measures between clusters and K-means of unsupervised clustering were used to estimate the number of CRGs clusters. The optimal number of clusters was determined by the cumulative distribution function (CDF) and the delta area and analyzed using the ConsensusClusterPlus R package with 1,000 repeats.

Principal component analysis (PCA) was done on the CRGs with the prcomp() R function and visualized using the Factoextra R package (version 1.0.7). Eigengene values of the first dimension of PCA (DIM1) were then visualized through a dot plot. We used the function fviz_contrib() to identify the genes that contributed more than average to the DIM1 axis. We used the ‘pheatmap’ package (version 1.0.12) to visualize the differences in CRGs and clinical characteristics of TCGA samples between cluster A and B.

Based on the above cluster classification model, the TCGA data set was divided into two groups, and the DEGs between the two groups were calculated using the ‘Limma’ (Linear Models for Microarray Data) package (24). The significance criterion of differential genes was set as p<0.05 (BH-corrected), and the absolute value of log fold change was greater than 0.18 (丨log2(1.2)丨). Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed and visualized by using the R package ‘clusterProfiler’ (version 3.12.0), with a strict cutoff of p < 0.05 and false discovery rate (FDR) of less than 0.05 (25). Among the genes obtained above, the genes with significantly good prognosis were screened by the same method as before. We used the ‘pheatmap’ package (version 1.0.12) to visualize the differences in CRGs and clinical characteristics of TCGA samples between DEG-cooper-cluster A and B.

Based on the DEGs of the above two subtypes, the univariable Cox regression analysis and Kaplan-Meier survival analysis were performed to obtain 474 DEGs with a good prognosis. We constructed a scoring system to quantify individual cuproptosis patterns in LUAD patients. The cuproptosis scoring system (CSS) was established based on the cuproptosis subtype-related DEGs. CSS was developed using principal component analysis (PCA), and an overall score was obtained by calculating the principal component sum.

CSS=∑(PC1_i+PC2_i)

With the median CSS values as the cutoff value, it was divided into high and low CSS groups. Patients with CSS values lower than the median risk score were categorized into the low-risk group, whereas those with CSS values greater than the median risk score were placed in the high-risk group.

The CSS scores were evaluated in TCGA and three GEO validation sets (GSE11969, GSE72094, GSE68465) to verify the prognostic effect. The prognostic effect of CSS scoring was analyzed by multivariate regression based on clinical data. The forest diagram is drawn using the ‘forestplot’ function package. In order to analyze and verify the similarity of CRGs expression level between the GEO validation data and TCGA data, the average value of CRGs expression level in these data was calculated respectively. Pearson correlation was used to calculate the correlation of CRGs expression between different data. Correlation heat map was drawn using ‘corrplot’ function package, and then correlation coefficients and significant p-values are also shown.

Single sample Gene Set Enrichment Analysis (ssGSEA) scores of the hallmark pathways between high and low CSS score groups were calculated by the ‘GSVA Bioconductor’ package (version 3.10), using pathway definitions from Molecular Signatures Database (MSigDB) Hallmark gene sets collection developed by Broad Institute (http://www.gsea-msigdb.org/gsea/msigdb/collections.jsp) (26). The Wilcoxon test was used for comparison of data between two groups.

We used the ESTIMATE algorithm to evaluate the TME scores of each individual. In addition, the 23 human immune cell infiltration levels of a single sample were calculated by ssGSEA. Furthermore, the levels of immune cell infiltration were also determined using the CIBERSORT algorithm (29). We also analyzed the correlations between the two subtypes of PD-1, PD-L1, and CTLA-4 expression.

The six immune categories of pan-cancer (C1-C6) and immune-related molecular features were calculated by Thorsson et al. (30). By using an immunogenomic analysis of more than 1000 tumor samples from 33 cancer types, six immune subtypes were identified, including wound healing (C1); IFN-γdominant (C2); inflammatory (C3); lymphocyte depleted (C4); immunologically quiet (C5); and TGF-βdominant (C6). The definitions and characteristics of C1-C6 are as follows:

C1 (Wound Healing) had elevated expression of angiogenic genes, a high proliferation rate, and a Th2 cell bias to the adaptive immune infiltrate. C2 (IFN-γ Dominant) had the highest M1/M2 macrophage polarization, a strong CD8 signal and, together with C6, the greatest TCR diversity. C3 (Inflammatory) was defined by elevated Th17 and Th1 genes, low to moderate tumor cell proliferation, and, along with C5, lower levels of aneuploidy and overall somatic copy number alterations than the other subtypes. C4 (Lymphocyte Depleted) displayed a more prominent macrophage signature, with Th1 suppressed and a high M2 response. C5 (Immunologically Quiet) exhibited the lowest lymphocyte, and highest macrophage responses, dominated by M2 macrophages. IDH mutations were enriched in C5 over C4. C6 (TGF-β Dominant) displayed the highest TGF-β signature and a high lymphocytic infiltrate with an even distribution of Type I and Type II T cells.

To explore differences in the therapeutic effects of chemotherapeutic and targeted drugs in LUAD cell lines between the high and low CSS score groups, the data that support the findings of this study were downloaded from Genomics of Drug Sensitivity in Cancer (GDSC) database (www.cancerrxgene.org) and Cancer Cell Line Encyclopedia (CCLE) (portals.broadinstitute.org/ccle). For the above data, the CSS scores were calculated and the differences in drug sensitivity between high and low CSS groups were analyzed. IC50 was quantified via the ‘pRRophetic’ package of R. Pearson method was used for correlation between CSS scores and IC50 of different drugs analysis.

We predicted the response to immune checkpoint blockade (ICB) through the TIDE website (http://tide.dfci.harvard.edu/login/). With the use of the Wilcoxon test, the difference in TIDE score between high and low CSS groups was compared. The predictive AUC of the CSS score was also calculated using TIDE scores.

The data analyses and visualization were conducted in R (version 4.1.1), and the following packages were used: “GEOquery,” “dplyr,” “Summarized Experiment,” “TCGA biolinks,” “readr,” “stringr,” “edgeR,” “org.Hs.eg.db,” “affy,” “limma,” “data.table,” “dplyr,” “zFPKM,” “ComplexHeatmap,” “ggplot2,” “ggalt,” “ggpubr,” “maftools,” “ggthemes,” “ggsignif,” “reshape,” “tidyverse,” “viridis,” “gridExtra,” “hrbrthemes,” “ggstatsplot,” “RColorBrewer,” “ggsci,” “RCircos,” “rtracklayer, rgl,” “pca3d,” “Rtsne,” “survival,” “survminer,” “ConsensusClusterPlus,” “RColorBrewer,” “pheatmap,” “GSVA,” “ggalluvial,” “plyr,” “forestplot,” “GSEABase,” “qusage,” “ggforce,” “survivalROC,” “ggsci,” “scales,” “maftools,” “enrichplot,” “gridExtra,” “TCGAbiolinks.” The Wilcoxon test was used for comparison of data between two groups, whereas the Kruskal-Wallis test was for comparison of data among three groups. Univariate survival analysis was performed by K–M survival analysis with the log-rank test. Pearson method was used for correlation analysis. The tumor mutation burden (TMB) score was calculated by the package “maftools” in R. A p value <0.05 was considered statistically significant (* p < 0.05; ** p < 0.01; *** p < 0.001; **** p < 0.0001).

After differential expression analysis in the TCGA dataset between 512 tumor tissues and 59 normal tissues, upregulated and downregulated DEGs were obtained. The 100 CRGs, including 63 upregulated genes and 37 down-regulated genes, listed in Supplementary Table 1 were obtained from Tsvetkov’s study (13). After the intersection of the two groups of data, 86 CRGs remained. Based on TCGA gene expression data and clinical data, we further analyzed the expression differences of CRGs in LUAD patients’ tumor tissues and adjacent nontumor tissues, as well as differential genes expression levels in clinical classification samples, including immune subtype, epidermal growth factor receptor (EGFR) mutation, age, gender, pathological combined stage, pathological nodal (pN) stage, pathological tumor (pT) stage, pathological metastasis stage, chemotherapy, echinoderm microtubule-associated protein-like 4- anaplastic lymphoma kinase (EML4-ALK) fusion, and pathological type ( Figure 1A and Supplementary Figure 1 ).

By analyzing the incidence of somatic mutation in 86 CRGs above, we found a relatively high mutation frequency in the LUAD cohort ( Figure 1D ). Of the 507 LUAD samples, 375 (74%) had mutations in the CRGs ( Figure 1D ). Among them, PRKDC, RB1, and CDKN2A were the top three genes with high mutation frequency, while 12 CRGs (TMEM191B, PDHB, IDH2, CERT1, MPC1, CEPT1, CIAO2A, CCDC137, C1QBP, SGF29, GCLM, and NME3) had no mutations ( Figure 1D ). We explored CNV in these 86 CRGs. Among them, UBAP2L, CHTOP, EGLN1, AHR, and HMGCS1 had widespread CNV increases, while CDKN2A, CNN1, FZR1, LDLR, and RPAP1 showed CNV decreases ( Figure 1B ). The locations of the CNV alteration in the CRGs on their respective chromosomes are shown in Figure 1C .

Interestingly, the expression levels of the top 10 CRGs were positively correlated with CNV alteration between primary tumor and normal tissues. CRGs with CNV loss, such as LDLR, CEPT1, and CNN1, were expressed at lower levels in LUAD tissues compared to those in normal tissues, while CRGs with CNV gains, such as CHTOP and UBAP2L, were significantly increased in LUAD samples ( Figures 1A, B ). These results showed that CNV might regulate CRG mRNA expression. However, some CRGs with CNV loss, such as IDH2, MTHFD1L, and ABCE1, showed upregulated mRNA expression, while other CRGs with CNV gain, such as TPK1, showed downregulated mRNA expression between tumor and normal samples. Therefore, CNV was not the only factor to regulate gene expression, although many of the observed changes in CRG expression might be explained by CNV alteration. Like some other factors, transcription factors and DNA methylation might also be involved in CRG regulation.

Our analysis revealed significant differences in the genetic landscape and expression levels of CRGs between LUAD and normal tissues, suggesting a potential role of CRGs in LUAD tumorigenesis.

In order to comprehensively understand the expression pattern of CRGs in the LUAD tumorigenesis, based on these 86 CRGs, a univariable Cox regression analysis was performed to screen the DEGs associated with a good prognosis. With the median DEGs expression values as the cutoff value, it was divided into high and low expression groups. The OS and PFI between high and low expression groups were compared by Kaplan-Meier analysis with the log-rank test, and p<0.05 was selected as the threshold for filtering. Finally, 18 CRGs with significant OS prognostic values were identified ( Supplementary Figure 2 ). The comprehensive landscape of CRGs interactions, connection strength, and their prognostic value in LUAD patients was demonstrated in a cuproptosis network ( Figure 2A ).

In order to further explore the expression characteristics of CRGs in LUAD, a consensus clustering algorithm was used to categorize the individuals with LUAD according to the expression profiles of the 18 CRGs. Our results suggested that k = 2 appeared to be an optimal selection for dividing the entire cohort into subtypes A (n=263) and B (n=249) ( Figure 2B ). PCA analysis showed that there were significant differences in the cuproptosis transcription profiles between subtypes A and B ( Figure 2C ). Kaplan–Meier analysis revealed a longer OS and PFI in patients with subtype B than that in patients with subtype A (log-rank test, p<0.0001, p=0.013, respectively; Figure 2D ). Furthermore, comparisons of the clinicopathological characteristics of the two subtypes of LUAD revealed significant differences between clinicopathological features and CRGs expression ( Figure 2E ).

GSVA enrichment analysis showed that the signaling pathways with significant differences between the two groups were mainly enriched in the immune microenvironment signaling pathway, including PI3K-AKT signaling, mTORC1 signaling, unfolded protein response, MYC target, Notch signaling, epithelial-mesenchymal transition, and KRAS signaling ( Figure 3A ). To investigate the role of CRGs in the TME of LUAD, we compared the difference in ssGSEA scores of 23 human immune cell subsets between the two clusters. We observed significant differences in the most immune cells infiltration levels between the two subtypes ( Figure 3D , Supplementary Figure 3 ). The infiltration levels of activated CD4 T cell, CD56dim natural killer cell, gamma delta T cell, natural killer cell, neutrophil, type 2 T helper cell were obviously higher in the cluster A than those in the cluster B, while eosinophil, immature dendritic cell, mast cell, T follicular helper cell had significantly lower infiltration in cluster A compared to those in cluster B. We further analyzed two important immune checkpoints and found that the expression of PD1 and PD-L1 increased in cluster A ( Figure 3B ). We also evaluated the TME score using the ‘ESTIMATE’ package and found that immune and stromal scores were elevated in cluster A (Figure 3C), which suggested higher relative contents of immunocytes or stromal cells in the TME. In contrast, the tumor purity decreased significantly in cluster A. The results demonstrated higher TME levels in patients with cluster A.

To explore the potential biological behavior of cuproptosis pattern, we identified 3153 DEGs between two cuproptosis consensus clusters and performed a functional enrichment analysis ( Figures 4A, B ). These CRGs were significantly enriched in biological processes that were correlated with extracellular matrix organization, extracellular structure organization, external encapsulating structure organization, organelle fission, nuclear fission, DNA replication, and chromosome segregation ( Figure 4A ). KEGG analysis indicated enrichment of PI3K-Akt signaling, focal adhesion, cell cycle, extracellular matrix receptor interaction, p53 signaling pathway ( Figure 4B ), suggesting that cuproptosis plays a vital role in tumor genesis and progression.

We then conducted univariable Cox regression and Kaplan–Meier analysis, revealing the significant good prognostic values of 474 DEGs in patients with LUAD. In order to further validate the mechanism related to cuproptosis, we used a consensus clustering algorithm to divide LUAD patients into two gene subtypes (DEGs-cooper-cluster A and B) according to prognostic genes. The two cuproptosis gene subtypes showed significant differences in CRG expression ( Figure 4C ), consistent with the expected results of the cuproptosis patterns. Kaplan-Meier analysis revealed a longer OS and PFI in patients with subtype A than in patients with subtype B (log-rank test, p<0.0001, p=0.0028, respectively; Figure 4D ). Furthermore, significant differences between clinicopathological features and CRGs expression were observed in DEGs-cooper-cluster A and B ( Figure 4E ).

To quantify individual cuproptosis patterns in LUAD patients, we established a CSS based on the cuproptosis subtype-related DEGs. The CSS value was divided into high- and low- CSS groups. Using the TCGA as a training set and GEO (GSE11969, GSE72094, GSE68465) as a testing set to verify prognostic efficacy, we found that the prognosis of the high-CSS group was significantly better than that of the low-CSS group in both the data sets. The expression level of DEGs was higher in the high CSS group. The distribution plot of the risk of CSS value revealed that mortality rates decreased while survival times increased with an increase in CSS value. Heatmap demonstrated the expression of DEGs between the high- and low-CSS value groups ( Figures 5A–D ). The Kaplan-Meier survival curves revealed that patients with high CSS values had a significantly favorable overall survival compared to that patient with low CSS values (log-rank test, p<0.0001, p=0.0093, p=0.0028, p<0.0001, respectively; Figures 5B–E , Supplementary Figure 4A, B ). In addition, the 1-, 3-, and 5-year survival rates of CSS values were represented by area under curve (AUC) values of 0.698, 0.652, and 0.634, respectively. Similar results were verified by the GEO database ( Figures 5C–F ). Correlation plots showed that the similarity of TCGA samples and GEO databases (GSE11969, GSE72094, GSE68465) in terms of gene expression for CRG genes ( Supplementary Figure 4C ).

To determine whether CSS value might independently predict OS in patients with LUAD, we combined the clinical features (age, sex, and TNM stage et al.) with CSS value to perform univariable and multivariable analyses. As shown in Figures 6A, B , the TNM stage and CSS value in the training set (TCGA) showed significant differences through univariable analyses, with consistent results observed in the testing set ( Figures 6C, D ). In multivariable analyses, CSS value was an independent factor that affects the OS in patients with LUAD in both training and testing sets.

The distribution of patients in the two cuproptosis consensus clusters, two cuproptosis-related DEGs subtypes, and two CSS score groups was shown in Figure 7A . We observed a significant difference in CSS value between cuproptosis gene clusters. The CSS value of cluster A was lower than that of cluster B ( Figure 7B ). More importantly, compared to DEGs subtype B, DEGs subtype A had a significantly higher CSS value ( Figure 7C ).

To investigate the impact of the CSS value on clinical characteristics, we explored the relationship between CSS value and different clinical features, including age (≤60 and >60 years), gender (female and male), TNM stage (different T, N, M stage), pathological tumor stage (stage I, II, III, and IV), radiation therapy (yes and no), and EML4-ALK fusion (yes and no). We observed significantly higher CSS values in patients in age>60, female, relatively early T stage, N stage, pathological tumor stage, and did not receive radiation therapy subgroup relative to those in the corresponding subgroup ( Figure 7D ).

Based on MSigDB cancer hallmarks, the correlation between CSS value and cancer hallmarks was calculated by Pearson correlation analysis. We found that the CSS value was positively correlated with myogenesis, hedgehog signaling, Kras signaling, bile acid metabolism, coagulation, heme metabolism, uv response dn ( Figure 8A ), while CSS values were negatively correlated with DNA repair, G2M checkpoint, unfolded protein response, PI3K-AKT signaling, MTORC1 signaling, E2F targets, MYC targets v1, MYC targets v2, oxidative phosphorylation, glycolysis, uv response up ( Figure 8A ; P< 0.05, FDR<0.25). GSEA was performed to determine the gene sets enriched in different CSS subgroups. The gene sets of the CSS-high samples were enriched in arrhythmogenic right ventricular cardiomyopathy, cardiac muscle contraction, focal adhesion, protein digestion and absorption, systemic lupus erythematosus. And the CSS-high samples were enriched in cell cycle, cellular senescence, DNA replication, oocyte meiosis, progesterone-mediated oocyte maturation ( Figures 8B, C ).

To investigate the role of CSS value in the TME of LUAD, we compared the difference in ssGSEA scores of human immune cell subsets between the high- and low- CSS groups. We observed significant differences in the most immune cells infiltration levels between the two groups ( Figure 8D ). The infiltration levels of Dendritic cells resting, Mast cell resting, Monocytes, T cells CD4 memory resting, Mast Cells, Dendritic Cells were obviously higher in the high CSS group than those in the low CSS group, while Macrophages M0, Macrophages M1, Mast cell activated, Neutrophils, T cells CD4 memory activated, T cells CD8, T Cells Follicular Helper, Neutrophils 1 had significantly lower infiltration in the high CSS group compared to those in low CSS group. We further analyzed TME scores and immune checkpoints and then found that the lymphocyte infiltration signature score increased in the high CSS subgroup ( Figures 8D, E ).

Next, we analyzed gene mutations to gain further biological insight in the CSS subgroups. We found significantly higher mutation counts in the low-CSS subgroup than in the high-CSS subgroup. Missense variations were the most common mutation type, followed by synonymous variations ( Figures 9A, B ). The mutation rates of TP53, MUC16, TTN, RYR2, KRAS, CSMD3, LRP1B, USH2A, and ZFHX4 were higher than 20% in both groups. Compared with the low-CSS subgroup, the mutations of the TP53, MUC16, and TTN genes were more common in the high-CSS subgroup, while the mutation of TP53, TTN, and CSMD3 genes were more common in the low-CSS subgroup than in high-CSS subgroup ( Figures 9A, B ). There was no significant difference in CNV between the two groups ( Figures 9C, D , Supplementary Figure 5 ).

In order to improve the therapeutic benefit of LUAD patients from chemotherapy and targeted therapy, we further explored whether cuproptosis signature could predict the sensitivity to several drugs widely used in LUAD between two CSS groups. According to the results calculated based on the GDSC database, IC50 and AUC_IC50 values of chemotherapy and targeted therapy drugs covering Acetalax, Cytarabine, EPZ5676, Pictilisib, GNE-317, AZD8186, Erlotinib, OSI-027, Buparlisib, Sapitinib, Ipatasertib, Dactolisib, AZD5363, Lapatinib, Osimertinib, Cediranib, VSP34 8731, Afatinib, AZD7762, and BMS-536924 were evaluated. The spearman correlation analysis showed that CSS values were positively correlated with IC50 and AUC_IC50 values of most drugs, except for OSI-027 and EPZ5676 ( Figure 10A ). Compared with the high-CSS subgroup, IC50 values of AZD5363, and AZD8186 were lower in the low-CSS subgroup, which indicated that low-CSS patients were more sensitive to these drugs ( Figures 10B, C ). The above results demonstrated that the CSS value had potential predictive value for chemotherapy and targeted therapy in LUAD.

We then used the tumor immune dysfunction and exclusion (TIDE) score to assess the rise and fall of immune response in different CSS subgroups. A higher TIDE prediction score represented a higher potential for immune evasion, which suggested that the patients were less likely to benefit from immune checkpoint inhibitor (ICI) therapy. In our results, the high-CSS subgroup had a lower TIDE score than the low-CSS subgroup, implying that high-CSS patients could benefit more from ICI therapy than low-CSS patients ( Figure 11A ). Besides, a higher TIDE prediction score was associated with a worse outcome. Accordingly, the high-CSS subgroup with a low TIDE score might have a better prognosis than the low-CSS group with a high TIDE score.

In addition, we found that the low-CSS subgroup had a higher myeloid-derived suppressor cells (MDSC) and T-cell exclusion score, but there was no difference in T-cell dysfunction between the two subgroups. Moreover, we assessed the prediction of the TIDE immune response by CSS. Receiver operating characteristic (ROC) curve showed that the AUCs for CSS was 0.726, which was lower than MDSC (AUC=0.887) and higher than cancer associated fibroblasts (CAF) (AUC=0.552) ( Figure 11B ). We suggested that the CSS has great predictive value for immunotherapy.