RNA sequencing, survival data, and clinical phenotypes were collected from The Cancer Genome Atlas (TCGA) database in University of California Santa Cruz (UCSC) Xena platform (http://www.genome.ucsc.edu/). After data cleaning and normalization, a total of 442 tumor tissues and 49 paracancerous tissues were enrolled. Meanwhile, the normal tissues from Genotype-Tissue Expression (GTEx) (https://www.gtexportal.org) were obtained. In addition, a total of 332 patients with survival time and mRNA expression matrix were download from GSE31210 (246) and GSE30219 (86) in Gene Expression Omnibus database (GEO, https://www.ncbi.nlm.nih.gov/geo/), and 301 patients with survival time clinical data were collected from Affiliated Tumor Hospital of Nantong University.

Unsupervised clustering was performed using The “ConsensusClusterPlus” R package. A consensus clustering approach (Euclidean distance) was run by K-means clustering algorithm for 1000 times with a resampling rate of 80%. An empirical cumulative distribution function plot was used to Identified the optimal number of clusters.

First, univariate Cox regression identified PCRGs in the TCGA-LUAD cohort. Subsequently, least absolute shrinkage and selection operator (LASSO) algorithm were performed on the PCRGs in TCGA-LUAD cohort. Finally, the cuproptosis-related signature was constructed by the stepwise Cox regression algorithm. Risk score = ∑ i = 1 n = Coef (Gene) × Expr (Gene). Coef is the coefficient, Expr is the Fragments Per Kilobase of exon Model per million mapped fragments (FPKM) of each gene. LUAD patients with risk-score above the median were categorized as the high-risk subgroup, and the rest were included in the low-risk subgroup. For model, the time-dependent receiver operating characteristic (ROC) was calculated across validation datasets. Survival analysis was used “survival”, “survminer” and “timeROC” packages, and the nomogram was used “rms” package. The ROC curve, the Harrell’s concordance index (C-index), calibration curve and detrended correspondence analysis (DCA) were performed to evaluate the performance of the nomogram.

The threshold value for screening differential genes was set to |log2FC| > 1 and adjusted P < 0.05. The threshold for stepwise-Cox analysis was set to P < 0.1 to screen the final pcrg. Paired samples of cancer tissue and paraneoplastic tissue from TCGA were used to test spatial differences in gene expression in the same individual. Human Protein Atlas database was used for the purpose of examining the expression of cuproptosis-related signature at the protein level (https://www.proteinatlas.org/).

The co-expression network of TCGA-LUAD was generated using the weighted gene weighted gene co-expression network analysis (WGCNA) package. A suitable soft threshold β is calculated based on the criteria for scale-free networks. In the following step, the weighted adjacency matrix was converted into a topological overlap matrix (TOM), and the corresponding dissimilarity (1-TOM) was calculated. Module identification was conducted using the dynamic tree cutting approach. The modules most relevant to the clinical phenotype of the risk score were selected for subsequent analysis.

The gene network analysis was performed using GENEMANIA to analyze potential interactions between these genes (http://genemania.org/). In order to explore the potential mechanisms and pathways between the riskscore subgroups and two clusters, the Gene ontology (GO), Kyotoencyclopedia of genes and genomes (KEGG) functional enrichment analysis, and gene set enrichment analysis (GSEA), gene set variation analysis (GSVA) were conducted among DEGs between the riskscore subgroups and clusters using the R packages “clusterProfiler”, “enrichplot”, “limma”, “ggplot2”, and “org.Hs.eg.db”.

“CIBERSORT” algorithm was used to assess the infiltration of 22 immune cells. Four other algorithms such as xcell, MCP-counter, GSVA, and ESTIMATE, were used to verify the robustness of CIBERSOR algorithm. The immune checkpoints were based on published research (14). The expression of HLA gene were used to assess the capabilities of antigen-presentation response.

In order to estimate the likelihood of drug interactions, the DEGs was identified from the high- and low-risk subgroups. L1000FWD (https://maayanlab.cloud/L1000FWD/) database was used to detect the signaling pathways affected by small molecule drugs. PubChem database was used for visualizing the structure of drugs (https://pubchem.ncbi.nlm.nih.gov).

Written informed consent was obtained from all patients involved in this study. The ethics committee of the Affiliated Tumor Hospital of Nantong University approved this study.

The 3D structures of mitoxantrone were downloaded from PubChem database (https://pubchem.ncbi.nlm.nih.gov/). The 3D structure of FDX1 was downloaded from the PDB protein database (https://www.rcsb.org/). Further, the protein was dehydrated and ligand extracted. Then, we conduct molecular simulation docking for the mitoxantrone and FDX1.

A complete set of data processing, statistical analysis, and plotting was carried out in R 4.0.3 software. In the case of normally distributed data, the unpaired Student t-test was used, while in the case of non-normally distributed data, the Wilcoxon test was used. An assessment of the correlation between two continuous variables was conducted using Pearson’s correlation coefficients. Univariate Cox regression and multivariate Cox regression were used to detect the effect of factors on prognosis of LUAD,. P < 0.05 was regarded as statistically significant.

As described in the literature, 12 genes including FDX1, LIAS, LIPT1, DLD,DLAT, PDHA1, PDHB, MTF1, GLS, CDKN2A, SLC31A1, and ATP7B, were confirmed to be associated with cuproptosis (6). These 12 genes were confirmed in LUAD by comparing their expression patterns in normal and tumor tissues from the TCGA and GTEx databases. It turned out that the expression of these 12 genes was significantly different in both tumor tissues and normal tissues ( Figures 1A, B ). Kaplan-Meier plotter shows 9 genes associated with lung adenocarcinoma prognosis and screened for subsequent analysis ( Figure S2 ). Co-expression analysis combined with prognosis value indicated that the expression of 4 cuproptosis-related genes (PDHB, DLAT, PDHB, DLD) was positively co-expressed with each other, and 9 genes (FDX1, LIPT1, DLAT, PDHA1, PDHB, MTF1, GLS, CDKN2A and SLC31A1) had prognostic value in LUAD ( Figure 1C ). Thereafter, we found a significant correlation between the expression of 12 cuproptosis-related genes in the TCGA cohort of LUAD patients ( Figure 1D ). The protein-protein interactions network among the cuproptosis-associated genes were displayed in Figure 1E . In parallel, we mapped the association of prognosis-related genes in the TCGA cohort with the standardized mortality rates ( Figure 2 ).

In a consensus clustering analysis of 442 tumor specimens from the TCGA cohort, 9 cuproptosis-associated genes were compared to the expression of LUAD subtypes to investigate the relationship between the two. Based on the cumulative distribution function curves of the consensus score matrix (CDF) ( Figure 3A ) and the proportion of ambiguous clustering statistics (PAC) ( Figures 3B, C ), k=2 was the optimal number. When a pair of samples has a high consensus score, they are more inclined to be clustered into the same cluster in successive iterations. Consequently, the 442 patients were divided into cluster 1 (n = 288) and cluster 2 (n = 154). Subsequently, a significant difference in OS was observed between patients in two clusters ( Figure 3D ). The expression of CDKN2A, MTF1, PDHA1, PDHB, SLC31A1 was significantly different between clusters 1 and 2 ( Figure 3F ). The distribution patterns from PCA analysis showed that samples could completely be distinguished into cluster1 and cluster 2 ( Figure 3H ). The heatmap showed the relationship between the different expression of genes and clinical characters ( Figure 3E ), and the volcano plot showed the logFC and FDR value among these genes ( Figure 3G ). There were 169 upregulated genes as well as 137 downregulated genes between the two clusters. The GSEA enrichment analysis found different expressed genes of the two cluster were enriched in some pathways, including signaling of antigen processing presentation mediated by class I MHC, transcription regulation of TP53, signaling of B cell receptor BCR and so on ( Figure 3J ). In addition, the GO enrichment analysis of these genes showed that they were also enriched in certain molecular processes, including histone-serine phosphorylation, response to decreased oxygen levels, multivesicular body, CXCR chemokine receptor binding, etc. ( Figure 3I ).

A prognostic model was constructed based on cuproptosis-related genes. In Cox univariate analysis, nine genes associated with cuproptosis were found to be significantly associated with overall survival (OS) in LUAD patients. Eight cuproptosis-related genes were tested and screened out by LASSO analysis ( Figures 4A, B ). Five genes were extracted using stepwise regression to construct the model ( Figure 4C ). In order to develop a risk-score model, the following algorithm was used: Risk score = (0.659854097) * DLAT + (-0.204720564) * GLS + (-0.67817864) * MTF1 + (0.423421209) * PDHA1 + (-0.805217147) * PDHB. The coefficients of the five PCRGs were displayed in Figure 4D , and TCGA-LUAD patients was divided into high-riskscore (n = 221) and low-riskscore (n =221) subgroups according to the median to facilitate the next step of study. As a measure of the signature’s specificity and sensitivity, areas under curve (AUC) values were calculated for 1-, 3- and 5 years. They were 0.71, 0.68 and 0.63 in TCGA-training set ( Figure 4E ). In addition, we found that the signature had not only a prognostic value but also a good diagnostic value ( Figure 4F , AUC=0.899). We detected the MRNA expression of cuproptosis-related signature in TCGA paired samples and found the expression of DLAT, PDHB, PDHA1 were significantly increased in cancer tissues when compared with cancer-adjacent tissues, while the expression of GLS and MTF1 trended downward, the decrease was not significant between paracancerous tissues and cancer tissues ( Figure 4G ). In the Human Protein Atlas database, the expressions of GLS, MTF1 and PDHA1 in tumor tissues were lower compared with those in normal tissues, while the expressions of DLAT and PDHB were significantly higher in tumor tissues compared with those in normal tissues ( Figure 4H ). Notably, all antibodies in the database were used for each gene to avoid false positives.

The risk survival status charts and Kaplan–Meier (K-M) survival analysis showed that survival time of LUAD patients in the low-risk group was longer than that of LUAD patients in the high-risk group, not only in TCGA training set ( Figures 5A, E ), but also in TCGA internal testing set ( Figures 5B, F ) and external verification set (GSE31210 and GSE30219) ( Figures 5C, D, G, H ).

According to K-M curves, the OS of the high-risk set was significantly worse than that of the low-risk group in age > 65 subgroups. In clinical subgroups of age ≤ 65, prognosis was not significant after riskscore stratification ( Figure 6A ). After stratifying for the characteristic variable of gender, it was found that the OS of the high-risk set was worse than that of the low-risk set in gender subgroups of patients ( Figure 6B ). In addition, stratification of pathological staging yielded resulted that the OS of the high-set was significantly worse than that of the low-risk set in stages III and IV subgroup ( Figure 6C ). Even though there were no significant difference in the OS of low and high risk sets in age ≤ 65, female, stage I and II subgroups, its overall OS trend was consistent with that seen before stratification. We speculated that this might be due to the lower sample size after stratifying samples and resulting in less statistical power.

Ultimately, even in different clinical subgroups, LUAD patients in the high-risk group had a lower survival probability than those in the low-risk group. In the subgroups after K-means clustering, it showed a significant difference in the distribution of riskscore, which proved the consistency of K-means clustering resulted with cuproptosis-related signature results ( Figures 5E-H ). Combining the analysis results of the K-means algorithm, we found that the riskscore was higher in Cluster 2 compared to the Cluster 1 ( Figure 6D ). The distribution of riskscores differed significantly between survival outcomes, suggesting that those who occurring ending events were more likely to be from a high-risk group ( Figure 6E ). There was no significant difference in the distribution of the riskscore between the two subgroups (aged > 65 years or ≤ 65 years), implying that age was not a factor affecting the distribution of the riskscore ( Figure 6F ). However, gender affected the distribution of the riskscore, with female patients having lower riskscore values compared to male patients ( Figure 6G ). It was noteworthy that the riskscores were higher as the tumor development and pathological stage changed ( Figure 6H ). The above results suggest that after stratification of clinical features, our riskscore distribution influenced prognostic outcomes, and our signature has high specificity and sensitivity in predicting the prognosis of LUAD patients.

Setting the soft threshold β to 7 (unsigned, R = 0.86) in the WGCNA could provide a suitable power value for the co-expression network ( Figures 7A, B ). After identifying 8 different modules, the correlation between different color modules and clinical features was calculated separately ( Figures 7C, D ). Brown module and riskscore subgroups exhibited the highest correlation of module-trait relationships. For the brown module, the R value between GS and MM was reached 0.74 in riskscore ( Figure 7E ) and 0.44 in status ( Figure 7F ), which suggested that the cuproptosis-related module was well constructed and the module was significantly associated with prognosis. To identify Hub cuproptosis-related module derived from cuproptosis-related patterns within the brown module, 167 genes with GS > 0.7 and MM > 0.2 were considered hub cuproptosis-related genes. The GO analysis of these hub genes was enriched into the lymphocyte differentiation, T cell activation, lymphocyte proliferation, protein complex involved in cell adhesion, inflammasome complex, immunological synapse, guanyl−nucleotide exchange factor activity, GTPase regulator activity and so on ( Figure 7G ). KEGG analysis enriched these genes into Chemokine signaling pathway, B cell receptor signaling pathway, PD−L1 expression and PD−1 checkpoint pathway in cancer and so on ( Figure 7H ). This suggested that our modules closely related to the cuproptosis-related signature and mainly affected the immunomodulation of patients with LUAD.

The heatmap showed the relationship between the riskscore subgroups and clinical characters ( Figure 8A ). The volcano plot showed the logFC and FDR value among these different expressed genes ( Figure 8B ).

As part of this study, we performed enrichment analyses on GO ( Figure 8C ), GSEA ( Figure 8D ), and GSVA ( Figure 8E ) to elucidate the underlying functions and pathways associated with our prognostic features. The results of GO, GSEA, GSVA, and WGCNA all suggested that genes associated with our prognostic features were associated with immune infiltration. Meanwhile, Risk-scores were negatively correlated with CD274 (PD-L1) and CTLA4 expression ( Figure 9A ). We evaluated the differences in infiltration of 22 immune cell types in the riskscore subgroup through the Cibersort database ( Figures 9B, C ). The results showed that B cells naive, Plasma cells, etc. were significantly higher while T cells CD4 memory resting, Macrophages M2, mast cells resting were significantly lower in the high-risk group when compared to the low-risk group. To demonstrate that the functions of the two subgroups were not biased by the analytical algorithm, Xcell ( Figure 9D ), ESTIMATE ( Figures 9E-H ), and MCP-counter ( Figure 9I ) algorithms were used to verify the stability and robustness of the Cibersort result. We further found significant differences in common immune checkpoints between the two subgroups as well ( Figure 9J ). HLA genes are also closely related to tumor immunity (16). Additionally, we examined whether HLA-related genes were expressed differently in risk subgroups. High-risk individuals tended to have lower HLA gene expression than low-risk individuals ( Figure 9K ). Together, these results revealed a higher tumor purity with lower immune cell infiltration in high-risk subgroup, which may affect immunotherapy outcomes in LUAD patients.

Identification of 20 important small molecules through the L1000FWD database as drugs to improve poor prognosis in high-risk groups ( Figure 10A ). Mitoxantrone was predicted to be the most promising drug. We further visualized the 2D ( Figure 10B ) and 3D ( Figure 10C ) structure of Mitoxantrone. A nomogram was created by integrating clinical information and genetic features from TCGA and performing multivariate Cox regression models ( Figure 10D ). Calibration plots were applied in OS outcomes, which demonstrated favorable concordance between predicted and observed OS at 1-, 3- and 5-year survival ( Figure 10E ). In terms of prediction, the C-index was 0.70 (0.67-0.73), which reflects relatively good performance. Also, LUAD-patients with high score had worse survival than those with low score ( Figure 10F ). The AUC values of 1-, 3- and 5-year in nomogram were 0.72, 0.74 and 0.74, respectively ( Figure 10G ). Additionally, we evaluated the nomogram model in TCGA-LUAD using decision curve analysis (DCA) ( Figure 10H ). Altogether, risk score was an independent and good prognostic indicator, and LUAD patients could benefit more after combining with pathological stage, T-stage and risk scores.

To confirm the role of cuproptosis-signature in LUAD, we further verified their differential expression in normal and tumor samples by vitro experiments. qRT-PCR was performed in paired samples of cancer and paracancer to detect mRNA expression levels of the prognostic cuproptosis-signature. We found that the expression of PDHA1, GLS, DLAT, PDHB and MTF1 were differed between cancer and paracancer ( Figure 11A ). We have deciphered that the small molecule compound mitoxantrone could treat the risk of death due to immune tolerance induced by cuproptosis-related signature. FDX1 as an enzyme catalyzing cuproptosis-related signature (17), we further used molecular docking to verify the direct interaction of FDX1 with the small molecule compound mitoxantrone ( Figure 11B ).

Next, tissue microarray of 301 lung cancer patients from the Nantong Cancer Hospital were used in a cohort study. FDX1 expression was identified by immunohistochemical staining in 301 lung cancer samples. We observed that FDX1 was located in the cytoplasm of tumor cells ( Figure 11C ). Immunohistochemical score (H-score) for immunostaining of tumor tissues also differed in each sample ( Figure 11D ). 119 (39.5%) patients were classified into FDX1 low expression subgroup while the FDX1 high expression subgroup had 182 (60.5%) patients. In our cohort, the patients with high FDX1 expression had significantly better outcomes than those with low FDX1 expression ( Figure 11E ). Stratification by pathological type revealed that the high expression group of FDX1 had a better prognosis than the low expression group in both LUAD and other lung cancer patients ( Figures 11F, G ). Most importantly, we constructed Nomogram and found that H-score was one of the most significant independent predictors of OS ( Figure 11H ). The Harrell’ s c-index for the nomogram model to predict the overall survival was 0.81. In a word, our study showed that FDX1, a key enzyme for cuproptosis, affected the prognosis of lung cancer patients by influencing the expression of GLS, PDHA1, PDHB, DLAT, and MTF1.