Cuproptosis-related gene signatures for predicting prognosis of lung adenocarcinoma.

Lung cancer (LC) is a common malignancy with high mortality rate, and lung adenocarcinoma (LUAD) is one of the common pathological types. Cuproptosis is a recently discovered new type of cell death dependent on mitochondria. However, the role of cuproptosis in LUAD is unknown. We obtained LUAD transcriptome data from the Cancer Genome Atlas (TCGA). Long-stranded non-coding RNA (LncRNAs) based on cuproptosis prognosis associated with LUAD were constructed for prognostic multi-LncRNA characterization. We divided TCGA-LUAD into training set and validation set to prove feasibility, and all samples were divided into high-risk group or low risk group. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were used to evaluate potential biological functions and explore the relationship between risk models and immunity. We identified 3 differentially expressed LncRNAs associated with LUAD prognosis and constructed prognostic model. Kaplan-Meier (K-M) analysis revealed prognostic model and LUAD prognosis. Our risk assessment model has a good reliability in predicting the prognosis of LUAD and was able to improve predictive ability of tumor mutational burdern. Single sample gene enrichment analysis (ssGSEA) revealed risk subgroups were associated with immune-related functions. The prognostic model based on cuproptosis lncRNA has important value in predicting the survival of LUAD patients.

1. Introduction

Lung cancer (LC) is the most common malignancy worldwide with a high mortality rate.[ Since 2004, lung adenocarcinoma (LUAD) has become the most common histologic type worldwide .[ LC is an extremely heterogeneous disease, and tumors with similar histology may have different outcomes.[ Despite the tremendous progress in the development of new therapies for oncology, the prognosis is poor. Currently, several drugs are commercially available that exert anticancer effects by inducing apoptosis.[ Therefore, exploring other forms of cell death can help overcome drug resistance in tumor cells, while identifying and characterizing prognostic biomarkers of LUAD can help guide clinical work.

Copper ions are widely present in various organisms and have pro-coagulant, hormonal maturation and other effects.[ The process of cuproptosis is mainly dependent on the accumulation of copper ions, which bind to the lipoacylated components of the tricarboxylic acid cycle, leading to the aggregation and dysregulation of these proteins, blocking tricarboxylic acid cycle, and then triggering proteotoxic stress and inducing cell death.[ Cuproptosis is a novel mode of death that differs from existing cell death mechanisms.[ Copper ions are significantly elevated in a variety of tumors such as breast, colon, pancreatic, and bladder cancers.[ Several studies have demonstrated that copper ion carriers and chelators promise to be potential drug molecules for the treatment of tumors.[ Long-stranded non-coding RNA (LncRNA) is a class of RNA that is >200 nt in length and does not encode protein. LncRNA is involved in a variety of biological functions such as chromatin regulation, gene expression, growth, and differentiation.[ In recent years, differences in LncRNA expression have been found in various cancers by high-throughput sequencing. LncRNA epigenetically induced lncRNA 1 can interact with proto-oncogene proteins through its 129 to 283 nt region and co-regulate the transcription of proto-oncogene target genes,[ thus promoting tumor cell cycle progression.[ Tumor cells are characterized by infinite replication in which telomeres at the ends of intact chromosomes play a central role. LncRNA telomeric-repeat-containing RNA can bind telomeres in cis and gene targets in trans. Telomeric-repeat-containing RNA and chromatin remodeling protein alpha thalassemia/mental retardation syndrome X share hundreds of target genes and play an antagonistic role at + some motifs.[ However, the role of cuproptosis-associated LncRNAs in cancer is unclear.

In this study, cuproptosis-related prognostic models were constructed and internally validated based on expression data and clinical data from LUAD patients in the TCGA database. In addition, we performed functional enrichment analysis and tumor microenvironment analysis to explore the potential mechanisms of action.

2. Materials and Methods

2.1. Data source

LUAD transcriptome data and corresponding clinical information were downloaded from the Cancer Genome Atlas (TCGA), including 535 cancer samples and 59 normal samples. Clinical information on 535 LUAD samples included sex, age, American Joint Committee on Cancer classification and survival information of Malignancies (tumor node metastasis) stage. Then, the cuproptosis gene was searched from Pubmed: NFE2L2, NLRP3, ATP7B, ATP7A, SLC31A1, FDX1, LIAS, LIPT1, LIPT2, DLD, DLAT, PDHA1, PDHB, MTF1, GLS, CDKN2A, DBT, GCSH, DLST. Further expression data of cuproptosis genes and LncRNA were extracted. Pearson correlation assessment was used to assess the relationship between the cuproptosis gene and LncRNA, |correlation coefficient| ≥ 0.4, P-value < .001 was considered as significant correlation. The ggplots package was used to draw mulberry plots showing the relationship between cuproptosis genes and LncRNAs. Expression data of cuproptosis-related LncRNAs were extracted for further analysis.

2.2. Construction of prognostic models

Survival data (survival status, survival time) and cuproptosis-related LncRNA expression data were combined and processed. Total TCGA-LUAD was divided into test group and validation group (Table S1, Supplemental Digital Content 1, http://links.lww.com/MD/H260). In the test group, univariate proportional hazards model (COX) analysis of cuproptosis-related LncRNA was performed to identify survival-related genes, which were screened by the standard of P value < .05 and presented in the form of forest plots. Based on the minimum absolute shrinkage and selective operator regression analysis (LASSO)-Cox method, cuproptosis-related LncRNAs were identified by “glmnet” package for core prognosis. Penalty parameter values (λ) were determined based on the lowest likelihood deviation of 10-fold cross-validation, and finally a prognostic model was constructed. The associated risk score was calculated for each patient in each test group by the following equation.

Risk score = β1 * expG1 + β2 * expG2 + ... + βn * expGn

β is the regression coefficient obtained from LASSO-Cox regression, and expG is the expression quantity of core prognostic cuproptosis LncRNA.

Patients were divided into high-risk score group (>median value) and low risk score group (

2.3. Validation of prediction models in test data

Using the above formula, we obtained risk scores for the test data. The “survival” package and “survminer” package were used to plot K-M curves to compare the differences in OS between high and low risk groups. The “pheatmap” package was used to draw risk curves, survival status maps, and risk heat maps to show the relationship between risk scores and survival.

2.4. Evaluation of risk prediction models

After confirming the reliability of the prediction model, we evaluated the risk scores of total TCGA-LUAD. K-M curves were plotted using the “survival” package and “survminer” package to compare the differences in OS and PFS between the high and low risk score groups. Further analysis of the differences in OS among risk score groups in stages of I to II、III to IV was performed. The “pheatmap” package was used to draw risk curves, survival status maps, and risk heat maps to demonstrate the relationship between risk scores and survival. Univariate COX regression analysis and multivariate COX regression analysis were used to determine independent risk factors of LUAD. To assess the accuracy of the model, the “survivalROC” package was used for ROC analysis, and the area under the ROC curve was calculated.

2.5. Construction and validation of predictive nomogram

To better evaluate the overall survival of LUAD, we used the “rms” package to construct nomogram integrating prognostic features to predict the OS of LUAD patients in 1, 3, and 5 years. the consistency of the nomogram in predicting LUAD was compared by calibration curves.

2.6. Validating the value of risk scores in total TCGA-LUAD

To further determine the predictive ability of the prognostic model for copper metabolism, risk scores were calculated in the same way in the test group and grouped using the same methods to obtain the same critical values. K-M analysis, ROC analysis and plotting of risk curves, survival status maps, and risk heat maps were performed as previously described.

2.7. Functional bioanalysis

We identified differential genes (DEGs) between high and low risk score groups using the “limma” package with the screening criteria: |log2FC| ≥ 1, FDR < 0.05. Gene ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of DEGs were performed using the “clusterProfiler” package.

2.8. Single sample gene enrichment analysis

GSVA package was used to assess the immune function of LUAD using the single sample gene enrichment analysis (ssGSEA) of immune cells and immune-related pathways. “Limma” was used to identify the different immune functions between groups and the results were presented by heat map.

2.9. Relationship between tumor mutational load and risk models

Download the tumor mutation burdern (TMB) data of LUAD from the TCGA database. Analyze the TMB between high and low risk groups using Perl software. Waterfall plots were drawn using the maftools package to show the top 15 genes in terms of mutation frequency between high and low risk score groups. The best cutoff values were selected into high and low tumor mutation load groups using the “survival” package and “survminer” package, and K-M curves were plotted to compare the differences in OS between the 2 groups. To better assess the predictive value of the risk model, the K-M curves of the combined tumor mutation load risk model were further plotted using the “survival” package and “survminer” package.

2.10. Statistical method

In this study, all statistical analyses were performed in R software (version 4.0.3). Continuous variables were compared using the Wilcoxon test. Survival analysis was performed using the K-M method to compare differences in OS and PFS between subgroups, and the rank sum test was used for validation. Univariate and multivariate Cox regression was used for the selection of characteristic genes. Model construction was performed using LASSO-Cox regression analysis with cross-validation. Unless otherwise stated, P value < .05 was considered significant.

3. Results

3.1. Prognostic modeling in the TCGA test group

First, we performed correlation analysis with the threshold value of |crrelation coefficient| ≥ 0.4 and P value < .001: CDKN2A, DBT, DLST, FDX1, GCSH, GLS, LIAS, LIPT1, MTF1, PDHA1, PDHB were strongly correlated with LncRNAs (Fig. 1A). Then the cuproptosis-related LncRNAs were extracted. In the test group, we first identified 4 survival-related LncRNAs (SP2-AS1, AC087501.4, AC090541.1, AC021660.4) using univariate COX (Fig. 1B). We used LASSO-Cox regression analysis to analyze and identify the expression profiles of the three cuproptosis-associated LncRNAs (SP2-AS1, AC087501.4, AC090541.1) and construct prognostic models (Fig. 1C and D). The risk score was calculated as follows:

Figure 1.

Construction of prognostic models in the test data set. (A) Cuproptosis gene-associated lncRNA correlation analysis. (B) Univariate COX regression analysis of cuproptosis-associated lncRNAs. (C) Lasso coefficient curves of survival curves between subgroups. (D) LASSO-COX regression to identify core cuproptosis-associated LncRNAs and construct prognostic models. (E) K-M curves of high and low risk groups. (F) Risk heat map of high and low risk groups. (G) Risk curves of high and low risk groups. (H) Survival status plots between high and low risk groups. K-M = Kaplan–Meier, LncRNA = long-stranded non-coding RNA.

Patients in the train group were divided into high risk and low risk groups according to the median value of the risk score. K-M survival analysis suggested that the high-risk group had worse OS than the low-risk group (P value < .001, Fig. 1E). High risk patients were more likely to die than low risk patients (Fig. 1F–H).

3.2. Validation of test data set

We used the above formula to obtain the risk scores of the test group, and also divided the patients in the validation group into high-risk score group and low risk score group. We also obtained the above conclusion: the high-risk group had worse OS than the low risk group (P value < .001, Fig. 2A); patients in the high risk group were more likely to die than those in the low risk group (Fig. 2B–D).

Figure 2.

Validation of the test data set. (A) K-M survival curves of high and low risk groups. (B) Risk heat map between high and low risk groups. (C) Risk curves between high and low risk groups. (D) Survival status plots between high and low risk groups. COX = proportional hazards model, K-M = Kaplan–Meier.

3.3. Evaluation of prediction models for TCGA data set

We used K-M analysis in the total TCGA-LUAD data set and found that the high-risk group had worse OS as well as PFS than the low-risk group (P value < .001, Fig. 3A and B); high risk patients were more likely to die than low risk patients (Fig. 3C–E). OS was different between the high and low risk score groups regardless of early or late stage (P value < .05, Fig. 3F and G). The value of using ROC curve to evaluate the model was higher than other clinical factors (age, sex, pathological stage) and had better predictive effect (Fig. 3H and I).

Figure 3.

Evaluation of the total TCGA-LUAD data set. (A) OS survival curves of the high and low risk groups. (B) PFS survival curves between high and low risk groups. (C) Risk heat map between high and low risk groups. (D) Risk curves between subgroups. (E) Survival status plots between subgroups. (F) K-M survival curves between subgroups in Stage I to II. (G) K-M survival curves between subgroups in Stage III to IV. (H) AUC curves of high and low risk groups. (I) Multi-indicator AUC curves. AUC = area under the ROC curve, K-M = Kaplan–Meier, LUAD = lung adenocarcinoma, OS = overall survival, PFS = progression-free-survival, TCGA = the cancer genome atlas.

3.4. Identification of independent predictors

Univariate and multivariate Cox regression analyses were performed in the correlated variables to determine whether risk score was an independent prognostic predictor of OS. In the univariate Cox regression analysis, risk score was significantly associated with OS in the TCGA-LUAD cohort (hazard ratio = 1.668, 95% confidence interval = 1.366–2.036, P value < .001, Fig. 4A). After correction for other confounding factors, risk score was still shown to be an independent predictor of OS in multivariate Cox regression analysis (hazard ratio = 1.550, 95% confidence interval = 1.260–1.907, P value < .001, Fig. 4B).

Figure 4.

Independent prognostic factor analysis. (A) Univariate COX regression to determine independent prognostic factors. (B) Multivariable COX regression to determine independent prognostic factors.

3.5. Functional bioanalysis in the TCGA cohort

To elucidate the biological functions and pathways associated with risk scores, GO enrichment and KEGG enrichment analysis were performed using DEGs between high and low risk scores. GO analysis revealed that DEGs between subgroups were mainly enriched in G protein-coupled receptor signaling pathway, positive mediation of MAPK cascade, etc (Fig. 5A). KEGG analysis revealed that DEGs between subgroups were mainly enriched between subgroups in neuroactive ligand-receptor interactions, calcium signaling pathway, and ECM-receptor interactions (Fig. 5B).

Figure 5.

Functional bioassay. (A) GO functional analysis. (B) KEGG functional analysis. GO = gene ontology, KEGG = Kyoto encyclopedia of genes and genomes.

3.6. Construction of predictive nomogram

Nomogram is a quantitative approach to patient prognosis. To accurately estimate the survival probability of LUAD patients in 1, 3, and 5 years, respectively, we integrated the prognostic characteristics of cuproptosis-related genes and other clinicopathological factors, including age, pathological stage, and T stage, to construct nomogram (Fig. 6A). Actual values and predictive values of the calibration curves in 1, 3, and 5 years showed that the nomenclature we constructed is reliable and accurate. This may help clinical practitioners to make clinical decisions for patients with LUAD and provide valuable insights for individualized patient treatment (Fig. 6B).

Figure 6.

Construction of the predictive nomogram. (A) Construction of the nomogram. (B) 1-, 3-, and 5-year calibration graphs.

3.7. Immunoassay

To further explore the relationship between LUAD prognosis and immune status, we quantified the infiltration scores of immune-related functions in both groups using ssGSEA algorithm. Significant differences were found between the high and low risk score groups in terms of type I interferon response, human leukocyte antigen, and antigen-presenting cell co-suppression response (P value < .05, Fig. 7).

Figure 7.

Single-sample gene set enrichment analysis of high and low risk score groups

3.8. Tumor mutation load analysis

The high-risk group possessed higher tumor mutation rate than that of the low-risk group. To accurately assess the relationship between TMB and prognosis, the high TMB group possessed a higher survival rate compared to the low TMB group (Fig. 8A and B). To further investigate the relationship between TMB risk score and survival, OS differed between different subgroups (P value < .05), with the longest OS in the high TMB + low risk score group and the shortest OS in the low TMB + high risk score group (Fig. 8C and D).

Figure 8.

Relationship between high and low risk score groups and tumor mutation load. (A) Tumor mutation burden in the high-risk group. (B) Tumor mutation load in the low-risk score group. (C) K-M survival analysis of high and low tumor mutation burden. (D) K-M survival analysis of tumor mutation burden combined with risk score. K-M = Kaplan–Meier.

4. Conclusion

In recent years, many studies have constructed LncRNA prognosis-related models for clinical guidance. Jing et al constructed a 6-gene regulator prognostic model in the TCGA data set and successfully validated the prognostic model in the validation set.[ Liu et al used the GEO data set to build a prognostic model of 4 iron apoptosis signature genes which was used to predict patient prognosis.[ Bian et al developed a new prognostic model of 4 cuproptosis-related genes, which showed strong predictive power.[ However, few prognostic models of cuproptosis-related Lncrnas have been reported in LUAD.

In this study, a novel 3-lncrNA prognostic model was identified based on cuproptosis-related genes. Compared with the traditional prediction model, it only contains 2 Lncrnas, so it is relatively easy to use. First, we divided TCGA-LUAD into a test set and a validation set. In the test set, we used LASSO-Cox regression to identify prognostic Lncrnas (SP2-AS1, AC087501.4, AC090541.1), constructed the prognostic model and calculated the risk score. According to the risk score and risk score threshold, all samples were divided into high-risk group and low risk group. The prognostic risk model was used to evaluate the survival of LUAD patients in the test set. There were differences in OS between the high and low risk groups. The same satisfactory results were obtained in the validation set. We established COX regression model and found that risk score was an independent prognostic factor. GO and KEGG analyses were performed to explore biological functions, and genes between groups were found to be closely related to cell signaling. Unfortunately, the mechanisms and biological functions of SP2-AS1, AC087501.4 and AC090541.1 have not been studied in cancer.

In this study, we found that CDKN2A, DBT, DLST, FDX1, GCSH, GLS, LIAS, LIPT1, MTF1, PDHA1, and PDHB were strongly correlated with lncRNA. The CDKN2A genomic locus is associated with human cancer and metabolic diseases, and this locus encodes antisense lncRNA (ANRIL). ANRIL mediated gene expression by staining in cells.[ ANRIL mediated histone modification and chromatin remodeling by interacting with PRC-1 and -2, and affected epigenetic transcriptional repression of adjacent genes CDKN2A and CDKN2B.[ In cancer cells, the increased function of ANRIL can promote the proliferation, metastasis, survival and epithelial-mesenchymal transition of tumor cells, while the loss of function of ANRIL reduces tumor size and growth, invasion and metastasis, and accelerates cell apoptosis and senescence.[ GLS is a phosphate-activated amide hydrolase that catalyzes the hydrolysis of glutamine into glutamate and ammonia to maintain metabolic stability and ensure biological energy balance.[ In one study, a novel lncRNA, GLS-AS, was significantly down-regulated in pancreatic cancer and associated with worse clinical outcomes.[ However, the interaction mechanism among DBT, DLST, FDX1, GCSH, LIAS, LIPT1, MTF10, PDHA1, PDHB and LncRNA is less in cancer. It is necessary to further study the relationship between these genes and LncRNA.

The body’s immune system has a powerful anti-tumor capacity and immunotherapy may become an effective tool for tumor treatment. A large body of literature confirms the point that immune-based treatment strategies provide a survival benefit for LUAD.[ One study found that copper supplementation enhanced PD-L1 expression at the mRNA and protein levels in cancer cells, and RNA sequencing showed that copper ions regulate key signaling pathways that mediate PD-L1-driven cancer immune evasion.[ In contrast, copper chelators inhibited STAT3 and EGFR phosphorylation and promoted ubiquitin-mediated PD-L1 degradation. Copper-chelating drugs significantly increased the number of tumor-infiltrating CD8+ T cells and natural killer cells, slowed tumor growth, and improved survival in mice.[ Therefore, investigating the mechanism and function between cuproptosis and LUAD to identify effective treatments is a future research direction. ssGSEA analysis demonstrated differences in immune function between high and low risk groups. Programmed death inhibitor-1 protein or its ligand (PD-L1) has gained significant clinical efficacy in the treatment of several tumors, especially LC. However, there are fewer biomarkers for its effective prediction. The assessment of TMB as the latest marker for programmed death inhibitor-1 antibody treatment efficacy has been verified in the treatment of colorectal cancer with mismatch repair defects.[ However, its use in LC has not been extensively extrapolated. In our study, the risk score was able to greatly improve the predictive power of TMB in LUAD.

Copper is an essential trace element in human body. When cells are growing under glycolytic conditions, most copper ionophores lose their killing activity.[ The researchers found that knockdown of FDX1, the gene encoding the Elesclomol target protein, attenuated cell death induced by copper ionophores.[ FDX1 is highly associated with a variety of malignant tumors. Copper ionophore may be a potential treatment for cancer cells with such metabolic characteristics. Notably, plant extracts involved in the process of copper-related metabolism may be a new intervention for the diagnosis and treatment of cancer.[

Our research has several advantages. First, we constructed a prognostic model based on cuproptosis-related lncRNA. At present, there are few related research reports. Secondly, cuproptosis is a new type of cell death, which is different from other cell death mechanisms, and its related research reports are relatively few. This new way of dying may bring promise for new treatments. Then, we constructed the prognostic score, which showed good predictive power in LUAD patients and showed a close relationship with the immune system, and could improve the predictive reliability of TMB. In addition, we focused on disease-free survival, with attention to symptom relief and tumor burden reduction, as compared with traditional prognostic models. Above, our study can provide new insights into the diagnosis and treatment of cuproptosis mechanism.

There are some limitations to our study. Firstly, due to the lack of lncRNA and survival data in other databases, this study lacks external verification. Second, our approach is based on integrated bioinformatics, and effective experimental validation of these findings is currently lacking. In addition, because cuproptosis is a newly discovered cell death mechanism, only a few genes related to cuproptosis have been identified. This can lead to biased results.

Our approach is based on comprehensive bioinformatics and currently lacks experimental validation of these findings. In addition, the prognostic value and mechanisms of cuproptosis still need further validation by more studies as cuproptosis has been less studied in LUAD. In conclusion, we constructed a model of prognostic cuproptosis characteristics in LUAD. Further validation of the risk score as a prognostic indicator independent of other clinical factors in LUAD. This study may also provide ideas and directions for cuproptosis-related studies in LUAD, which may be conducive to new therapeutic approaches.

Author contributions

Conceptualization: Wei Ye, Xingxing Li.

Data curation: Wei Ye, Xingxing Li.

Formal analysis: Wei Ye, Xingxing Li.

Methodology: Wei Ye, Yuenuo Huang.

Software: Yuenuo Huang, Xingxing Li.

Visualization: Yuenuo Huang.

Writing – original draft: Wei Ye.

Writing – review & editing: Wei Ye.