A Novel Cuprotosis-Related lncRNA Signature Effectively Predicts Prognosis in Glioma Patients

Cuprotosis is a novel and different cell death mechanism from the existing known ones that can be used to explore new approaches to treating cancer. Just like ferroptosis and pyroptosis, cuprotosis-related genes regulate various types of tumorigenesis, invasion, and metastasis. However, the relationship between cuprotosis-related long non-coding RNA (cuprotosis-related lncRNA) in glioma development and prognosis has not been investigated. We obtained relevant data from the Genotype-Tissue Expression (GTEx), Cancer Genome Atlas (TCGA), Chinese Glioma Genome Atlas (CGGA), and published articles. First, we identified 365 cuprotosis-related lncRNAs based on 10 cuprotosis-related differential genes (|R2|> 0.4, p < 0.001). Then using Lasso and Cox regression analysis methods, 12 prognostic cuprotosis-related lncRNAs were obtained and constructed the CuLncSigi risk score formula. Our next step was to divide the tumor gliomas into two groups (high risk and low risk) based on the median risk score, and we found that patients in the high-risk group had a significantly worse prognosis. We used internal and external validation methods to simultaneously analyze and validate that the risk score model has good predictive power for patients with glioma. Next, we also performed enrichment analyses such as GSEA and aaGSEA and evaluated the relationship between immune-related drugs and tumor treatment. In conclusion, we successfully constructed a formula of cuprotosis-related lncRNAs with a powerful predictive function. More importantly, our study paves the way for exploring cuprotosis mechanisms in glioma occurrence and development and helps to find new relevant biomarkers for glioma early identification and diagnosis and to investigate new therapeutic approaches.


Introduction
Glioma is a common and fatal primary tumor in the central nervous system. It is classified by the World Health Organization (WHO) as grade I to IV according to its histological features, and the higher the malignancy, the worse the prognosis (Davis 2018). Although the current advanced surgical techniques, radiotherapy, and other comprehensive treatments bring hope to patients, their overall prognosis is still poor, getting more severe economic pressure and burden on patients themselves, their families, and society (Bielecka and Markiewicz-Żukowska 2020). Therefore, there is an urgent need for early identification and diagnosis of glioma-related biomarkers in our clinical care and novel therapeutic approaches that can improve their prognosis in the long term.
Recently, Todd Golub's team demonstrated a copperregulated cell death-dependent mechanism distinct from known death mechanisms that occur through direct binding of copper to the lipidated components of the tricarboxylic acid (TCA) cycle, naming it cuprotosis (Tsvetkov et al. 2022). This study explains the pathology associated with hereditary copper overload disorders and suggests that cuprotosis could be a new approach to treating cancer (Kahlson and Dixon 2022). The study reported ten genes closely related to copper death. After reviewing the literature, we found previous related reports demonstrating the critical role of these genes in tumors. For example, FDX1 can closely influence the prognosis of lung adenocarcinoma through its involvement in fatty acid oxidation and amino acid metabolism (Zhang et al. 2021). Dihydrothioctanamide dehydrogenase (DLD) induces ferroptosis in head and neck cancer cells by regulating cystine deprivation (Shin et al. 2020). Upregulation of PDHA1 gene expression can inhibit the Warburg effect from inducing apoptosis in hepatocellular carcinoma (Sun et al. 2019). More importantly, with the confirmation of the specific mechanism of cuprotosis, the perfect answer to the previous question about why copper inhibits the growth of tumors was given (Ren et al. 2021). For instance, copper inhibits the activity of glioblastoma cells because cuprotosis may be involved in the development of tumor cells (Ceyhan et al. 2021). The fact that copper complexes can affect the redox state of gliomas to act as antiproliferative may also be due to the presence of cuprotosis in gliomas (Illán-Cabeza et al. 2020). Therefore, cuprotosis, which reveals a novel cell death mechanism, has a great potential to become a new cancer treatment. Cuprotosisrelated genes may be critical markers for treating glioma and improving the prognosis of glioma patients. So, the construction of prognosis-related features of cuprotosis-related genes to explore their gene functions is particularly important for improving the future treatment of glioma.
Long non-coding RNA, as regulators involved in various intracellular processes, has been shown to play an essential role in epigenetics and tumorigenesis and development. It has been shown that downregulation of lncRNA PVT1 can stop tumor growth by inhibiting glioma cell proliferation and migration (Li et al. 2022a). Knockdown of lncRNA KCNQ1OT1 enhanced tumor regression by TMZ in a mouse model of U251 glioma (Wang et al. 2022). In addition, it has also been reported that lncRNAs can be used to treat gliomas and predict their prognosis (Gao et al. 2018). Many studies have recently demonstrated that immune, iron death, and pyroptosis-associated lncRNAs can be used as tumor markers. Therefore, we aimed to investigate the relevant role of cuprotosis-related lncRNAs on glioma prognosis to explore their essential mechanisms in tumor development, search for suitable new glioma biomarkers for early identification and diagnosis, and investigate other new therapeutic approaches.

Data Collection
The data used for the study included 1152 non-tumor brain tissue RNA sequence data from the GTEx database and 698 glioma RNA sequence data and clinical information from the TCGA database; the data used for validation were 1018 glioma patient RNA and clinical data from the CGGA database. All data were normalized to fragment per kilobase million (FPKM) values. When extracting clinical data, survival times less than 30 days and the presence of undocumented data collection were excluded. Cuprotosis-associated genes were obtained from those reported in the article (Tsvetkov et al. 2022). Since the data involved in this study were obtained from public databases, ethics committee approval was not required according to the relevant regulations of the databases.

Identification of Cuprotosis-Related lncRNAs and Relationship with Prognosis
First, we analyzed the differential expression of genes associated with cuprotosis in glioma and normal tissues using the Wilcoxon method (FDR < 0.05, |log 2 FC|> 1). Subsequently, the "limma" package in R software was used to identify lncR-NAs associated with cuprotosis (|R 2 |> 0.4, p < 0.001). Univariate COX regression analysis was used to calculate the lncRNAs significantly related to the prognosis of patients, and a total of 295 lncRNAs were selected. Since the number of lncRNAs is large and collinearity between variables is prevented, machine methods (elastic net, Lasso, and ridge regression) are used to reduce the collinearity of variables. Previous studies have demonstrated and emphasized the advantage of collinearity dimension reduction of Lasso regression (Rooney et al. 2015;Yao et al. 2021). We hope to use mechanical methods to reduce the number of lncRNAs to the lowest dimension and screen out lncRNAs with more substantial explanatory power. Finally, multivariate COX regression analysis identified cuprotosis-related lncRNAs as prognostic signatures for gliomas.

Construction and Validation of the Risk Score Formula
We constructed a risk score (RS) formula based on 12 prognostic signatures associated with cuprotosis.
The "Coef CuLncSigi" in this "risk scores" formula represents the coefficient value, which is the regression coefficient of the 12 prognostic lncRNAs derived from the multifactorial regression analysis. The "EXP CuLncSigi" in the formula represents the expression levels of the 12 cuprotosis-related lncRNAs. Using the RS formula, we can get the risk value of each patient. And after getting the risk value of all patients, we can find out the median risk value of the patients. We can determine the patient's risk level according to the median value. Then, Kaplan-Meier (K-M) analysis, time-dependent ROC, principal component analysis (PCA), independent prognostic analysis, nomogram (Davis 2018;Tsvetkov et al. 2022;Zhang et al. 2021), calibration curve, and the expression level of lncRNAs are used to determine the accuracy of the model. In the CGGA validation sample, we applied the same intermediate values to assess the validity and reliability of our RS formula using the same way as coef CuLncSigi × EXP CuLncSigi above. We also use the same approach to randomly divide the TCGA data into two groups for internal validation.

Functional Annotation of Cuprotosis-Related lncRNAs
Functional enrichment analysis was performed using GSEA (version 4.1.0, (p < 0.05, FDR < 0.25)) (Rooney et al. 2015). The infiltrating fraction of 16 immune cells and the activity of 13 immune-related pathways were measured by ssGSEA (Rooney et al. 2015). We also explored the correlation of immune checkpoints between the two risk groups (Yao et al. 2021).

Prediction of Immune-Related Drugs
To find more drugs for the treatment of glioma, we focused on evaluating and predicting immune-related drugs. According to the online tool Cancer Drug Sensitivity Genomics, the IC50 of different drugs on glioma samples was predicted using the R package 'pRRophetic.' The primary use of 'pRRophetic' is to predict phenotypes from gene expression data (to predict clinical outcomes using Cancer Genome Project CGP cell line data), to predict drug sensitivity in external cell lines (CCLE), and also for clinical data prediction.

Statistical Analysis
This study was statistically analyzed using R software (version 4.1.2) and GSEA software (version 4.1.0). Wilcoxon test was used to determine the expression levels of cuprotosis-related genes in cancer and normal tissues. Machine methods (elastic net, Lasso, and ridge regression) and Cox regression assessed the prognostic impact of cuprotosisrelated lncRNA signatures. Survival curves were generated using the Kaplan-Meier (KM) method and compared using the log-rank test.
In addition, to further understand the intrinsic association between the ten cuprotosis-related genes, we also performed a correlation analysis. The results showed the most significant positive and negative correlations (Fig. 1B). The above results suggest some interaction between cuprotosis-related genes in gliomas. Then, 365 cuprotosis-related lncRNAs (Table 1) were screened (|R 2 |> 0.4, p < 0.001).
We calculated the risk score for each patient in the TCGA dataset using the RS formula.
T After obtaining a risk score for each patient, the patients were divided into two groups based on the median risk score: a high-risk group and a low-risk group ( Fig. 2A). We found that more and more patients died as the risk score increased (Fig. 2B). Figure 3C shows the 12 prognostic cuprotosisrelated lncRNAs involved in both groups by heat map. The ROC curve region shows the good predictive ability of the model based on 12 survival-related lncRNAs. In the TCGA data, the AUC values were 0.896, 0.93, and 0.901 at 1, 3, and 5 years, respectively (Fig. 2D). According to K-M analysis, patients with high RS had worse survival than those with low RS (Fig. 2E).

Internal and External Validation of Prognostic Features
Using the same cut-off for the CGGA validation data as for the TCGA data, it was possible to distinguish the highrisk group from the low-risk group. However, the number of patients in the low-risk group was significantly lower ( Fig. 2F). CGGA patients showed that high-risk patients were positively associated with poor prognosis (Fig. 2G). The expression of prognostic cuprotosis-related lncRNAs in CGGA was similar to that in TCGA samples (Fig. 2H). In CGGA samples, the AUC values were 0.754, 0.809, and 0.822 at 1, 3, and 5 years, respectively (Fig. 2I). K-M analysis performed on the CGGA data showed the same results as the TCGA data (p < 0.001, Fig. 2J). The similar validation results presented in the two validation datasets of TCGA (the AUC values were 0.921, 0.961, and 0.927 at 1, 3, and 5 years; the AUC values were 0.867, 0.893, and 0.71 at 1, 3, and 5 years) also demonstrated the excellent predictive ability of the model ( Figure S1).
To further explore the relationship between risk scores and clinical characteristics of patients, we performed subgroup analysis. We found that risk scores predicted prognosis in different patients (Fig. 3). These results also strengthen the predictive potential of the risk score.         We also evaluated the expression levels of cuprotosisrelated lncRNAs in the TCGA dataset. We found that all genes showed significant differences in different classes, and in the validation data CGGA data, all genes except one showed similar trends in other classes (Fig. 4A, B). In both TCGA and CGGA datasets, the expression of genes showed the same trend as the tumor grade increased. PCA analysis showed that high-risk and low-risk patients could be distributed in different quadrants according to RS of prognostic cuprotosis-related lncRNAs (Fig. 4C-E).

Construction and Validated the Nomogram
We constructed a nomogram based on age, gender, grade, and risk score from TCGA data to facilitate clinical work. First, univariate and multivariate Cox regression analyses showed that age class and risk score were independent predictors of OS in patients with glioma (Fig. 5A, B). Nomograms allow precise calculation of survival probabilities based on patient grade, gender, age, and risk score (Fig. 5C). The calibration curves showed good agreement between actual OS and predicted survival rates (Fig. 5I-K). The ROC curves also validated that the risk score based on cuprotosis-related lncRNAs had the highest predictive accuracy (Fig. 5G). The reliability of the risk score was also validated by analysis of the CGGA data ( Fig. 5H and L-N).

Functional Annotation of Cuprotosis-Related lncRNAs
We used GSEA to further investigate the biological functions of prognostic features. In the KEGG analysis, the main added  features were pathogenic-escherichia-coli-infection, systemiclupus-erythematosus, cell-cycle, n-glycan-biosynthesis, aminosugar-and-nucleotide-sugar-metabolism; the reduced functions are taste-transduction, long-term-depression, long-term-potentiation, etc. (Fig. 6A). We quantified the enrichment scores of ssGSEA by measuring immune cell subpopulations and related pathways to investigate further the correlation between risk scores and immune cells and functions. We found that most cells (B cells, CD8 + T cells, DCs, Tregs, etc.) were significantly elevated in the high-risk group (Fig. 6B). T cell co-stimulation, APC co-stimulation, CCR, T cell co-stimulation, and type I IFN response were higher in the high-risk group than in the low-risk group (Fig. 6C). The above results suggest that immune function is more active in the high-risk group. Due to the importance of checkpoint-based immunotherapy, we also compared and analyzed the differences in immune checkpoint expression between the two groups (Fig. 6D).
We also analyzed the correlation between predictive characteristics and tumor immune-related drugs. The results found lower IC50 of cisplatin, etoposide, and rapamycin in the high-risk group and higher IC50 of lenalidomide and PAC-1 in the high-risk group (Fig. 6E-I), which helps to explore individualized treatment regimens suitable for the high-risk patient.

Discussion
Glioma is a common intracranial tumor, and the factors that contribute to its development and progression are not fully understood. In-depth studies at the molecular level will help to elucidate the pathogenesis of glioma and predict the prognosis of patients. Cuprotosis is associated with genetic disorders and tumors and has the potential to be a therapeutic strategy (Tsvetkov et al. 2022;Cui et al. 2021;Tsang et al. 2020;Davis et al. 2020). Recently, several studies have shown that lncRNAs are involved in the development of multiple cancers and that different types of lncR-NAs predictive features are used to predict the prognosis of cancer patients (Gao et al. 2021;Luo and Ma 2021;Zhu et al. 2021). However, the role of cuprotosis-related lncR-NAs in glioma patients has not been reported. Our study identified multiple prognostic cuprotosis-related lncRNAs for the first time and successfully constructed a formula for cuprotosis-related lncRNAs with powerful predictive features and introduced the role of cuprotosis-related lncRNAs in predicting glioma prognosis and immune status. This study first extracted ten genes involved in cuprotosis closely and then identified 12 prognostic cuprotosisrelated lncRNAs for prognostic signatures using reliable biological analysis. After analysis, AC084824.4, AC104117.3, AC121761.2, AL391834.1, DNAJC9-AS1, LINC01503, RNF219-AS1, and SNAI3-AS1 were regarded as protective factors; AL355974.2, CRNDE, LINC02328, and TMEM220-AS1 were regarded as risk factors. And among these lncR-NAs, the biological functions of some lncRNAs with were confirmed. For example, high expression of CRNDE promotes the malignant progression of gliomas (Li et al. 2017). LINC01503 plays an essential role in HCC through the MAPK/ERK pathway . TMEM220-AS1 regulates hepatocellular carcinoma by regulating the miR-484/MAGI1 axis (Cao et al. 2021). Then, we constructed an RS model containing these 12 lncRNAs. Based on the RS model, we divided glioma patients into low-risk and highrisk groups; patients in the high-risk group were shown to have a significantly worse prognosis in both TCGA and CGGA data. To evaluate the reliability of the prognostic features, on the one hand, the risk-risk model was validated by the Kaplan-Meier curve and ROC curve; the correlation analysis of clinical variables and risk scores also increased the reliability of the model. On the other hand, the onethree-and five-year nomograms reported that the model was a good predictor of OS prognosis in glioma patients. The calibration curves further confirmed that the model was accurate. We also analyzed the expression levels of 12 lncRNAs using the TCGA and CGGA databases to verify the reliability of the prognostic features. All genes showed similar trends except for AC084824.4 in the CGGA data. Most importantly, the expression levels of genes that are risk factors increased with increasing tumor grade, while the expression levels of protective genes were reversed.
GSEA shows systemic-lupus-erythematosus, cell-cycle, n-glycan-biosynthesis, amino-sugar-and-nucleotide-sugarmetabolism, and leukocyte-transendothelial-migration ware mainly enriched in at-risk groups. It is well known that intracellular metabolic processes such as nucleotides, amino acids, and glutathione play an irreplaceable role in cancer. For example, aberrant glycosylation of cells leads to abnormal expression of membrane-localized glycans, triggering a malignant transformation of cells (Thomas et al. 2021). N-glycans play an essential role in breast and oral cancers (Hirano and Furukawa 2022;Wu et al. 2022). Glutathione affects tumor development by altering the sensitivity to oxidative stress in astrocytic tumors (Moreira Franco et al. 2021). The ssGSEA results showed a significant increase in most cells (macrophages, CD8 + T cells, neutrophils, mast cells, Tregs, etc.) in the high-risk group. Some studies found that the degree of tumor immune cell infiltration correlated with the prognosis of tumor patients. For example, CD8 + T cell infiltration was associated with a poor prognosis in BC patients . High infiltration of tumor-associated macrophages is associated with poor prognosis in glioma and thyroid cancer (Li et al. 2022b;Ryder et al. 2008). The degree of MC infiltration in mouse and human gliomas is proportional to the malignancy of the tumor (Põlajeva et al. 2011(Põlajeva et al. , 2014. A high ratio of neutrophils to lymphocytes predicts a poorer OS in BC patients (Tan et al. 2018). Pathological grading of gliomas is positively correlated with infiltrating neutrophils (Khan et al. 2020). lncRNA HOXA-AS2 promotes Tregs proliferation and immune tolerance via miR-302a/KDM2A axis, thus promoting glioma progression and poor prognosis (Zhong et al. 2022). An increase in Treg and MDSC in mouse gliomas leads to decreased overall survival (Zhai et al. 2021). In addition to increased tumor immune cell infiltration, immune-related pathways were significantly higher in the high-risk group. In other words, although the immune function was more active in the high-risk group, the decreased anti-tumor immunity of patients in the high-risk group may have contributed to their poor prognosis. This also suggests a new idea for immunotherapy of glioma. Glioma immunotherapy has been a hot topic (Xu et al. 2020). Therefore, we also performed an analysis of immune checkpoint expression and CGGA. I-K Calibration plotsof the nomogram for predicting the probability of OS at 1, 3, and 5 years inthe TCGA. L-N Calibration plots of the nomogram for predicting theprobability of OS at 1, 3, and 5 years in the CGGA between high-risk and low-risk groups, hoping to provide some direction for glioma immunotherapy. We studied the sensitivity of immune-related drugs among patients and found that high-risk patients may be sensitive to Cisplatin, Etoposide, Rapamycin, and resistant to lenalidomide, PAC-1. This implies that high-risk groups may benefit from treatment with multiple immune-related drugs. We hope that the above study provides a basis for precise, individualized treatment of glioma patients. However, our work has some limitations. First, all the data and results obtained are based on public databases, and more data are needed to validate the model and improve it. Second, despite the high predictive power of this feature found in this work, the stability of this constructed model depends on the current sample size. Finally, cuprotosis is a newly discovered mode of cell death, and studies about its associated lncRNAs are limited; more studies are needed to elucidate the specific mechanism of how our discovered cuprotosis-related lncRNAs induce cuprotosis.

Conclusion
In conclusion, we screened the lncRNAs with a predictive role in cuprotosis and successfully constructed an RS model with a powerful predictive function. More importantly, this predictive model on cuprotosis provides new ideas for the therapeutic approach to glioma and offers new directions for studying cuprotosis-related lncRNAs in more fields.