Data merging was done using the TCGA-LUAD cohort, which included 524 LUAD samples from the TCGA database (https://tcga-data.nci.nih.gov/tcga/), and the GSE31210 cohort, which included 226 LUAD samples from the GEO database (http://www.ncbi.nlm.nih.gov/geo/). Clinical information and normalized matrix files were retrieved from the GEO database, and information on gene expression from RNA sequencing (FPKM values) and clinical information were obtained from TCGA. The FPKM values were then transformed to transcripts per kilobase million (TPM) values for additional analysis, and batch effects were removed using principal component analysis (PCA) and “ComBat” from the “SVA” R package. Samples without complete survival information were disqualified. Finally, we were able to get a complete LUAD cohort with 750 samples and 16928 genes.

LUAD tissues were obtained from patients who had undergone surgery at the Affiliated Hospital of Nantong University. In our cohort, 5 pairs of tissues were obtained between 2018 and 2020. The study was authorized by the Ethical Committee of Affiliated Hospital of Nantong University (2022-L119). RNA from tumor tissues was extracted using TRIzol Reagent (Invitrogen), and the cDNA was obtained through reverse transcription using a PrimeScript™ RT Reagent Kit (TaKaRa, Shiga, Japan). Real-time quantitative polymerase chain reaction (qPCR) was conducted in triplicate for each sample using a SYBR Premix Ex Taq II Reagent Kit (TaKaRa). The primer sequences for the target genes were as follows: CDKN2A forward 5′-CCAGGTCATGATGAT-3′, reverse5′-TGCAGCACCACCA-3′; GAPDH forward 5′-TGACTTCAACAGCGACACCCA-3′, reverse 5′-CACCCTGTTGCTGTAGCCAAA-3′.

The Human Protein Atlas (HPA) (http://www.proteinatlas.org/) is a database containing immunohistochemical staining results from a wide range of tumors and normal tissues. We have used HPA data here to explore the expression of gene.

GSCALite (http://bioinfo.life.hust.edu.cn/web/GSCALite/) provides an online cancer genomic analysis platform by integrating 33 cancers data from TCGA and normal tissue genomics data from GTEx. In this study, we analyzed the genomic level, copy number level, and methylation level of CRGs in LUAD by GSCALite.

Functional enrichment analysis of differentially expressed genes connected to CRGs in LUAD was used to investigate functional annotation and enrichment pathways. ClusterProfiler was used to evaluate Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways.

According to the positive and negative relationships between the DEGs and the cluster signature, the DEGs were divided into two groups, namely sigC1 and sigC2. Then, the “clusterProfiler” R package was used for gene annotation. We then used the Boruta algorithm combined with PCA to reduce the dimensionality of the DEGs subgroups and calculated the cuproptosis score for each sample. The LUAD comprehensive cohort was divided into the high- and low-cuproptosis score groups based on the optimal cutoff value. The cuproptosis score of each LUAD sample was calculated using the following formula: Cuproptosis score=∑PsigC1−∑PsigC2

We selected “survival status”, “T stage”, “N stage”, “M stage”, “clinical stage”, and “recurrence or metastasis” as clinical subgroup characteristics, and drew box plots to show the differences in the cuproptosis score between different clinical characteristics. A stacked histogram was drawn to show the proportion of each clinical characteristic in the high- and low-cuproptosis score groups.

We used the “survival” and “survminer” R packages to perform survival analysis to compare the differences in OS between the high- and low-cuproptosis score groups, and used the “ggalluvial” R package to draw Sankey diagrams to visualize the correspondence among cuproptosis score groups, different subtypes, and prognosis. Box plots were used to compare the differences in the cuproptosis score of different subtypes. ssGSEA was used to quantify the infiltration abundance of immune cells, and the relationship between cuproptosis score and immune cell infiltration levels was displayed using a correlation heat map.

We used the “limma” R package to compare the differences in the gene expression of several common immune targets. Next, we downloaded the immunophenoscore (IPS) data of the TCGA-LUAD cohort from The Cancer Immunome Atlas (TCIA) database to explore the differences in the efficacy of the four immune checkpoint inhibitors (ICIs) between the high- and low-cuproptosis groups (10), including ctla4_pos_pd1_pos, ctla4_neg_pd1_pos, ctla4_pos_pd1_neg, and ctla4_neg_pd1_neg. In addition, we collected data from 2 immunotherapy cohorts for immune benefit validation, including the GSE91061 (Nivolumab immunotherapy) and GSE13507 cohorts (intravenous BCG immunotherapy).

All analyses were performed using R version 4.1.1. Unless otherwise specified, Pearson’s correlation coefficient was used for correlation analysis in this study. For comparison between the two groups in the bioinformatics analysis section, the Wilcoxon test was used for difference analysis. For comparison between the two groups in the experimental section, the Students’ t-test was used for difference analysis. Kaplan-Meier survival analysis and log-rank tests were used to compare the survival of the different groups of patients. For all statistical analyses, a two-tailed p <0.05 was considered statistically significant.

First, we analyzed the transcriptomic data of patients with LUAD from TCGA database. The expression of seven cuproptosis-associated genes (CRGs) was found to be significantly higher in LUAD than in normal controls. These seven CRGs, namely, CDKN2A, DLAT, LIAS, DLD, PDHA1, MTF1, and FDX1, were deemed differentially expressed genes (DEGs) ( Figure 2A ). Of the seven genes, five (CDKN2A, DLAT, LIAS, DLD, and PDHA1) were upregulated, while two (MTF1 and FDX1) were downregulated. Gene mutations and copy number variants (CNV) are closely associated with tumorigenesis and tumor progression. Therefore, we screened LUAD patients for genetic changes in the CRGs. Figure 2B shows the variant classification, variant type, single nucleotide variant (SNV) class, and variation per sample for the CRGs. We observed that missense mutations were the most common type of genomic mutations and that CDKN2A, with a mutation frequency of 39%, harbored the most mutations of all CRGs in LUAD samples. Furthermore, CDNK2A also harbored the highest deleterious SNVs in LUAD (i.e. number of samples with at least one deleterious mutation site/number of samples with SNV mutation data) ( Figure 2C ). We then the CNVs in the CRGs and found associations between CNVs and mRNA levels of nine genes (DLD, LIAS, PDHB, DLAT, MTF1, CDKN2A, FDX1, PDHA1, and LIPT1) ( Figure 2D ). Homogeneous and heterogeneous variations, that is, homogeneous/heterogeneous deletion in CDKN2A and heterogeneous amplification in DLD/LIPT1/GLS were found ( Figure 2E ). Figure 2F shows the percentage of various types of CNVs—heterozygous amplification, heterozygous deletion, homozygous amplification, and homozygous deletion—in each CRG in LUAD. Methylation is an important epigenetic alteration that remodels genes, including those associated with cancer and may lead to uncontrolled growth. Therefore, we investigated the methylation profiles of CRGs in LUAD samples and assessed corresponding mRNA expression. We found that the methylation levels of CRGs were largely negatively correlated with mRNA levels ( Figure 2G ). In addition, six genes in LUAD (LIPT1, DLAT, DLD, PDHA1, GLS, and CDKN2A) had different methylation levels ( Figure 2H ). As CDKN2A, a pivotal contributor in the cell cycle, showed the highest ectopic expression, deleterious mutation frequency, and differentiated CNV in LUAD compared with normal controls, we tried to verify its expression in LUAD tissues. We measured CDKN2A expression in five pairs of tissues using RT-qPCR and found that CDKN2A expression was significantly higher in LUAD compared to normal lung tissue ( Figure 2I ). Immunohistochemistry results also showed that CDKN2A expression was higher in LUAD tissues compared to normal lung tissues ( Figure 2J ). In summary, our analysis revealed significant differences in the genetic profiles and expression of CRGs between LUAD and control samples, suggesting a potential role of CRGs in LUAD tumorigenesis and progression.

To investigate the potential biological functions of cuproptosis in LUAD tumorigenesis in a larger sample, we merged the TCGA-LUAD samples and the GSE31210 dataset. After batch effects removal using the R packages “limma” and “sva”, we obtained a comprehensive cohort of 750 LUAD samples and 16,928 genes ( Figure 3A ). In addition, we constructed a network map reflecting the correlation of ten CRGs with clinical outcomes, and manifesting the molecular interactions between these CRGs ( Figure 3B ). We then performed the Kaplan-Meier analysis to investigate the prognostic value of the CRGs in LUAD ( Figure 3C ). Survival curves revealed that PDHA1, DLAT, and CDKN2A were risk factors, and LIPT1, GLS, and MTF1 were protective factors in LUAD. In conclusion, we found that CDKN2A, DLAT, and PDHA1 were not only highly expressed in LUAD samples but also predicted poor clinical outcomes in the long term.

We performed unsupervised clustering and classification based on the CRGs. By increasing the clustering variable (k) from 2 to 10, we found the highest intra-group correlations and low inter-group correlations at k = 3, suggesting that LUAD patients may be well segregated into three clusters ( Figure 4A ). Survival analysis revealed that prognosis differed substantially among the three cuproptosis subtypes, and subtype C had considerable survival advantages (p < 0.001, Figure 4B ). Subtype C had a longer survival time compared to subtype A and subtype B. Moreover, the expression of CRGs in the different clusters also showed significant differences ( Figure 4C ). The clinicopathological features of the three subtypes and the expression of the CRGs were unveiled using a heat map ( Figure 4D ).

To further infer the biological characteristics of the more malignant LUAD cluster, we compared the pathways enriched in each of the LUAD subclusters. We downloaded the KEGG pathway and the HALLMARK pathway from the Msigdb database and subsequently performed GSVA scoring of these pathways and compared the differences in the pathways enriched in each of the three cuproptosis mutant subtypes ( Figures 4E, F ). We found significant differences between the three subtypes, mainly in metabolism-related pathways (propanoate metabolism, phenylalanine metabolism, fatty acid metabolism, and glutathione metabolism) and immune-related pathways (IL6/JAK/STAT3, inflammatory response, interferon alpha response, and interferon-gamma response).

The aforementioned findings suggest that 1) cuproptosis is associated with the malignancy and tumorigenesis of LUAD and that 2) LUAD subclusters grouped by cuproptosis-related mechanisms are differentiated in anti-tumor immunity. Therefore, we sought to quantify the infiltration landscape of the tumor microenvironment (TME) in patients with LUAD. We assessed TME scores (stromal score, immune score, and estimated score) of the three subtypes using the “ESTIMATE” package. For the TME score, a higher stromal score or immune score represents a higher relative abundance of stromal or immune cells in the TME, while the estimated score indicates the aggregation of stromal or immune scores in the TME. The results showed that in all three cuproptosis subtypes, patients with subtype C had the highest TME scores ( Figure 5A ). The ssGSEA algorithm was used to calculate the fraction of immune cell infiltration per LUAD sample between the three subtypes and to compare the differences in immune cell infiltration between the subtypes. We observed significant differences in the infiltration of most immune cells between the three subtypes ( Figure 5B ). The results showed that compared with clusters A and B, cluster C had more activated dendritic cells, CD56 bright natural killer cells, eosinophils, immature B cells, immature dendritic cells, macrophages, mast cells, monocytes, natural killer cells, plasmacytoid dendritic cells, T follicular helper cells, Type 1 T helper cells, and Type 17 T helper cells. This also explains why cluster C had the best survival advantage. With these results, we found that subtype C had better TME scores, immune infiltration levels, and LUAD prognosis. To further explore the potential feature of each cuproptosis subtype, we performed a differential analysis of the three cuproptosis subtypes. Discrepant analysis results of the three cuproptosis subtypes are summarized in Supplementary Table 1 and the DEGs are shown using Volcano plots (p<0.05, Figure 5C ). Notably, CDKN2A was significantly more expressed in subcluster A and subcluster B, while GLS was significantly more expressed in subcluster C, indicating a potential positive correlation between high CDKN2A and low GLS expression and worse clinical outcomes of LUAD. We then screened DEGs with a log fold change> 7 and identified the 72 most significantly differentially expressed genes among the three cuproptosis-related LUAD subtypes. GO ( Figure 5D ) and KEGG ( Figure 5E ) enrichment analyses were performed to identify the top five pathways associated with the DEGs based on adjusted p-value in KEGG analysis and their relationship networks with related genes were displayed ( Figure 5F ). The results suggested that the CRGs were mainly associated with processes related to aberrant metabolism and tumor immunity. Of the top 5 enriched pathways, immune-related processes, such as complement and coagulation cascades and IL-17 signaling pathways were further studied.

To validate our findings, we extended our sample size and covered eight single-cell transcriptomic datasets (EMTAB6149, GSE117570, GSE127465, GSE127471, GSEGSE131907, GSE139555, GSE143423, and GSE99254) to assess the correlation between CRGs and TME cellular components after cell annotation. We found that compared with CDKN2A, GLS was significantly highly expressed in several immune clusters including CD4+T, CD8+T, and mononuclear phagocytosis system (MPS) derived from the LUAD TME. ( Figure 5G ). We then investigated the prognostic relevance of CDKN2A and GLS expression in patients with LUAD. We analyzed 1925 LUAD clinical data using the Kaplan-Meier Plotter database and found that patients in the low CDKN2A/high GLS expression group had a better prognosis than those in the high CDKN2A/low GLS expression group ( Figure 5H ). These findings suggest that cuproptosis could serve as an indicator to predict tumor immunity and prognosis.

To identify CRGs with prognostic significance in LUAD, we performed a univariate Cox regression analysis on 72 DEGs identified by cuproptosis-related LUAD clusters ( Supplementary Table 2 ). Of these, 22 DEGs had prognostic value (p<0.001, Figure 6A ) and were used for further clustering. At k = 2, intra-group correlations were the lowest, indicating that LUAD patients could be well divided into two gene clusters ( Figure 6B ). Kaplan-Meier curves showed that patients from geneClusterA had worse overall survival (log-rank, p<0.001; Figure 6C ). As shown in Figure 6D , we found significant differences in the expression of the 22 DEGs in these two gene clusters. Figure 6E shows the gene expression profiles as well as clinicopathological features of the 22 DEGs in the two gene subtypes.

To further evaluate the association between cuproptosis and LUAD prognosis, we constructed a prognostic signature for cuproptosis. After dimensional reduction, we defined the Cuproptosis score that predicts the long-term prognosis as the sum of transcriptomic profiles derived from two geneClusters. LUAD patients with higher Cuproptosis scores were divided into the high-risk group and those with low Cuproptosis scores, in the low-risk group. Kaplan-Meier survival curves showed a significant difference between the two groups and verified the correlation between cuproptosis-based transcriptomic profiles and LUAD prognosis (log-rank test, p<0.001; Figure 7A ). We then explored the correlation between the Cuproptosis score and the clinicopathological features of LUAD. We observed that patients in the high-risk group were associated with higher mortality, TNM staging, and recurrent/metastatic tumors ( Figures 7B–F ). The relationship between the Cuproptosis score and cuproptosis clusters A-C and cuproptosis geneClusters A-B was further explored ( Figure 7G ). Cluster A and gene cluster A were significantly correlated with a high Cuproptosis score and a worse prognosis. It was proved that prognostic results obtained by different clustering modes were consistent. Subsequently, a Sankey diagram was drawn to illustrate the relationship between patients in the three cuproptosis subtypes, two gene subtypes, two cuproptosis score groups, and survival status ( Figure 7H ). Finally, tumor tissues from the high-risk group were enriched with immunosuppressive cells such as Treg and MDSC and lacked effector cells such as activated B cells, activated dendritic cells, eosinophils, immature dendritic cells, mast cells, monocytes, plasmacytoid dendritic cells, and T follicular helper cells ( Figure 7I ). These findings indicate that cuproptosis-based transcriptomic characteristics could contribute to identifying LUAD cases with poor immune infiltration and a high likelihood of poor clinical outcomes.

After constructing the Cuproptosis score and verifying its utility in prognosis prediction and reflecting immune infiltration, we explored whether this score could be applied to facilitate immunotherapy in clinical practice. First, we found that the expression of most immunosuppressive genes (such as PD1, PD-L1, CTLA4, LAG3, CTLA4, and TIGIT) was significantly higher in the high- than in the low-risk group ( Figure 8A ). We then used two methods to validate the ability of the Cuproptosis score to predict the benefits of immunotherapy. Recent studies have reported that IPS based on immunogenicity could effectively predict the response to immunotherapy in patients with malignancies (10). Therefore, we downloaded the IPS of the TCGA-LUAD cohort from TCGA database to explore differences in the efficacy of immunotherapy between the high- and low-cuproptosis groups. Then, we compared the correlation between the IPS and the Cuproptosis score. In the immunotherapy scoring, we found that both anti-PD-1 and anti-CTLA-4 resulted in better immunotherapy benefits for patients and that compared to patients in the high cuproptosis group, those in the low cuproptosis group benefited more from immunotherapy ( Figure 8B ). To further analyze the application of the cuproptosis score in the assessment of immunotherapy, we used the GSE91061 and GSE13507 cohorts, in which patients received nivolumab and intravesical BCG immunotherapy, to be divided into the high Cuproptosis score and low Cuproptosis score groups. Through Kaplan-Meier analysis, we found that patients with a low Cuproptosis score had a better prognosis and higher response percentage than high Cuproptosis score ( Figures 8C, D ). These results also reaffirmed that patients with low Cuproptosis scores benefit more from immunotherapy than those with high Cuproptosis scores and that the combination of the Cuproptosis score with ICI expression may help better predict sensitivity to immunotherapy. In conclusion, the Cuproptosis score has great potential in predicting prognosis and immunotherapeutic benefits.

Finally, we selected the targeted drugs currently used to treat LUAD to assess the sensitivity of these agents in patients grouped by Cuproptosis score. The “pRRophetic” package was used to predict the IC50 values for each sample for multiple targeted drugs, comparing the differences between the high- and low-score groups. In patients with a high cuproptosis score, we predicted low IC50 values and high sensitivity to targeted drugs, such as cisplatin, vinblastine, paclitaxel, docetaxel, gefitinib, and gemcitabine. In patients with a low cuproptosis score, we predicted low IC50 values and high sensitivity to targeted drugs such as imatinib, bexarotene, bicalutamide, axitinib, AKT inhibitor VIII, and erlotinib ( Figure 9 ). These results show that the Cuproptosis score can predict not only prognosis and immunotherapeutic benefits but also guide targeted drugs and comprehensive patient treatment. It is a novel tumor biomarker with outstanding potential.