PDF

J Cancer 2023; 14(9):1499-1514. doi:10.7150/jca.82370 This issue Cite

Research Paper

A five-cuproptosis-related LncRNA Signature: predicting prognosis, assessing immune function & drug sensitivity in lung squamous cell carcinoma

Hongtao Zhao¹, Lei Wu², Qinyuan Liao¹, Peiluo Huang^1,3, Ruonan Sun¹, Xiuzhen Yang⁴ , Juan Du¹

1. Department of Immunology, College of Basic Medicine, Guilin Medical University, Guilin 541199, Guangxi, China.
2. College of continuing education, Guilin Medical University, Guilin 541004, China.
3. College of pharmacy, Guilin Medical University, Guilin 541199, China.
4. Department of clinical laboratory, Zibo Central Hospital, Zibo 255036, China.

✉ Corresponding authors: yangxz770com; sunnydujuanedu.cn.

Abstract

Lung squamous cell carcinoma has so far lacked effective targets for diagnosis and treatment. In cancer research, long noncoding RNAs (LncRNAs) emerge as novel therapeutic targets and biomarkers. Cuprophosis is a new death type involving multiple biological processes in tumor cells. Here, we aimed to explore whether Cuprophosis-related lncRNAs could be used to predict prognosis, assess immune function, and test drug sensitivity in LUSC patients. The Cancer Genome Map (TCGA) was used to obtain genome and clinical data, and Cuprophosis-relevant genes were found in the literature. A cuproptosis-related lncRNA risk model was built using co-expression analysis, univariate/multivariate Cox regression, and LASSO analysis. The survival analysis was used to assess the model's prognostic value. The univariate and multivariate Cox regression analyses were performed to determine whether risk score, age, gender, or clinical stages could be used as independent prognostic factors. Gene Set Enrichment Analysis and mutation analysis were performed on differentially expressed mRNA between high-risk and low-risk groups. The (TIDE) algorithm was used to conduct immunological functional analysis and drug sensitivity testing. Five cuproptosis-related LncRNAs were identified, and the selected LncRNAs constructed a prognosis model. According to the Kaplan-Meier survival analysis, the overall survival time for patients in the high-risk group was shorter than for those in the low-risk group. For LUSC patients, the risk score serves as an independent prognostic indicator. The GO and KEGG enrichment analysis revealed that the differentially expressed mRNAs between the high- and low-risk groups were enriched in several immune-related processes. The enrichment score of differentially expressed mRNAs in the high-risk group is higher than that of the low-risk group in multiple immune function pathways, including the IFN-γ and MHC I pathways. The Tumor Immune Dysfunction and Exclusion (TIDE) test revealed that the high-risk group was more likely to experience immune escape. The drug sensitivity analysis showed that patients with low-risk ratings were likely to respond to GW441756 and Salubrinal. In contrast, patients with higher risk scores were more responsive to dasatinib and Z-LLNIe CHO. The 5-Cuprophosis-related lncRNA signature can be used to predict prognosis, assess immune function, and test drug sensitivity in LUSC patients.

Keywords: Cuproptosis, LncRNA, LUSC, Prognostic model, Drug sensitivity

Introduction

Lung cancer is the most prevalent malignancy and one of the most prominent reasons for high cancer-related mortality worldwide [1]. The most common subtypes of non-small cell lung cancers (NSCLCs) are lung squamous cell carcinoma (LUSC) and lung adenocarcinoma (LUAD). Meanwhile, LUSC represents approximately 30% of NSCLC cases [2, 3]. The primary treatment for patients with early-stage NSCLC is surgical tumor removal. However, for patients with advanced stages of LUSC, targeted therapy is frequently required [4]. Even though a handful of the LUSC's oncogenic drivers have been identified, the disease lacks more targeted treatment therapies [5]. Therefore, it is crucial to identify novel treatment targets for LUSC patients.

Copper (Cu) is one of the essential metal elements for most living organisms. It involves many biological processes, including mitochondrial respiratory, Fe uptake, anti-oxidant, and detoxification [6]. Furthermore, intracellular or extracellular Cu ion dyshomeostasis leads to oxidative stress and cytotoxicity [7]. Excess Cu concentrations in the tricarboxylic acid (TCA) cycle have recently been discovered to cause cell death, which differs from all previous modalities of cell death. This phenomenon is known as "cuproptosis," according to Tsvetkov et al. [8, 9]. Copper level changes in patients' serum and tumor tissues have been discovered [10]. Cu ion encourages tumor growth and metastasis by stimulating angiogenesis [11, 12]. Therefore, studies on copper-associated genes are revealing new pathways of tumorigenesis and progression.

Long noncoding RNAs (LncRNAs) are defined as RNA molecules with longer than 200 nucleotides in length and do not possess protein-coding properties [13]. LncRNA is a component of the transcriptome that can regulate the expression of specific genes based on cell type, developmental stage, and cell function [14]. The role of LncRNAs in gene expression regulation occurs at both the transcriptional and posttranscriptional levels [15]. For a long time, tens of thousands of LncRNA were discovered to be abnormally expressed or mutated in various cancer types [16]. The abnormal expression or mutation of LncRNA influences tumor occurrence, metastasis, and stage [17]. As previously stated, these phenomena occurred primarily because LncRNAs can regulate microRNA (miRNA) and messenger RNA (mRNA) involved in cancer progression [18]. As a result, LncRNAs are emerging as novel therapeutic targets and biomarkers in cancer research.

Our study discovered differentially expressed cuproptosis-related LncRNAs in LUSC tissues using the Cancer Genome Atlas (TCGA) database. Furthermore, we developed a prognostic prediction model for predicting overall survival (OS) in LUSC patients using cuproptosis-related LncRNAs. The cuproptosis-related LncRNAs pattern could be used as a prognostic indicator for LUSC patients and could provide vital information for cancer therapy.

Materials and methods

Screen cuproptosis-related LncRNAs

LUSC gene expression data, clinical data, and mutation data were obtained from The Cancer Genome Atlas (TCGA) data portal (https://tcga-data.nci.nih.gov/tcga/). Our analysis contains 551 samples, including 502 LUSC samples and 49 LUSC adjacent normal tissues samples. The 19 cuproptosis-related genes were identified through a literature review (Supplementary Table 1). The Pearson correlation coefficient test picked up a total of 291 cuproptosis-related LncRNAs associated with 19 genes (the 'limma' R package). We created a Sankey map for cuproptosis-related LncRNA-mRNA co-expression network using the ggalluvial package in R version 4.1.3.

Model building

To determine a cuproptosis-related LncRNA signature, we first used univariate Cox analysis to screen LncRNAs related to overall survival for LUSC. Subsequently, cuproptosis-related LncRNA with P < 0.05 (n = 23) in the univariate Cox analyses were selected for the least absolute shrinkage and selection operator (LASSO) regression analysis and multivariate Cox regression analyses. The Risk Scores were generated using a multivariate Cox regression coefficient. Risk Score = coef(LncRNA AC002467.1) × (LncRNA AC002467.1)expr + coef(LncRNA MIR3945HG) × (LncRNA MIR3945HG)expr + coef(LncRNA AC010328.1) × (LncRNA AC010328.1)expr + coef(LncRNA AC008972.2) × (LncRNA AC008972.2)expr + coef(LncRNA AC253536.6) × (LncRNA AC253536.6)expr. The correlation between 5 cuproptosis-related LncRNAs and 19 cuproptosis-related mRNAs was analyzed by the Pearson correlation coefficient.

Model evaluation

Patients lacking enough survival information (n = 6) or with a survival time of 0 days (n = 2) were excluded. Data that shared a clinical patient number (n = 1) were also consolidated into a single item. As a result, 493 samples (Supplementary Table 2) altogether were used in the clinical survival study. LUSC patients were classified into high-risk and low-risk groups according to their median risk score. Overall survival (OS) and progression-free survival (PFS) were analyzed by the Kaplan-Meier method. PFS datasets were gained from UCSC Xena (http://xena.ucsc.edu/). We also calculated the area under the ROC curve (AUC) used to assess the accuracy of the Risk Score. Afterward, the LUSC research samples were randomly divided into training (n = 198) and testing (n = 295) sets according to the ratio of 2:3. The AUC analysis was performed to verify the suitability of the Risk Score in two separate sets. The ROC curve showed with the R software 'timeROC' package. LUSC patients were randomized using the 'caret' R package. Differences in training and testing sets were compared using the chi‐square test.

Propensity Score Matching (PSM) analysis

For PSM analysis, we used R's 'wakefield,' 'MatchIt,' 'knitr,' 'captioner,' 'tidyverse,' and 'foreign' packages to reduce selection bias between high- and low-risk groups. Patients who met the following criteria were excluded: patients with incomplete age, sex, clinical stage, and survival information. The PSM study covered 484 patients in total. Gender, age, and clinical stage were used as factors in a logistic regression model with 1:1 exact matching to produce propensity scores.

Independent prognostic analysis

The univariate Cox and multivariate Cox analyses were performed based on risk score, gender, age, and clinical stage. Independent prognostic analysis was conducted using the 'survival' R package. Principal component analysis (PCA) was utilized to show the distribution of the low-risk and high-risk groups using the 'limma' and 'scatterplot3d' R packages. Survival analyses of different clinical stage LUSC patients were explored using the 'survival' and 'survminer' R package.

Nomogram construction

A nomogram for the prediction of LUSC patients was constructed based on Risk Score, age, gender, and clinical stage. Moreover, calibration curves can be used to evaluate the coherence of intended and actual results. Display the images above using the 'survival,' 'regplot,' and 'rms' R packages. The specific information of the variables for the nomogram constructions is listed in Supplementary Table 3.

Gene Set Enrichment Analysis (GSEA)

Screening criteria for the differential expression mRNAs (DEmRNAs) were |log2FC| ≥ 1 and p-value < 0.05. Analysis was done using the 'limma' R package. Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) enrichment analysis of DEmRNAs was performed by 'clusterProfiler,' 'org.Hs.eg.db', 'ggpubr' and 'RColorBrewer' R package. GO terms and KEGG pathways with p < 0.05 were considered significantly enriched. Our study chose 10 GO terms, and 18 KEGG pathways with the smallest P-value were used to display.

Tumor mutation burden (TMB) analysis, Immune-related functional analysis, and pharmaceutical screening

'Maftools' R package was used to visualize and analyze tumor mutation burden. To immunization score 493 patients with LUSC who had complete clinical information, the 'limma,' 'GSVA,' and 'GSEABase' software packages, as well as the "ssgsea" method, were used. The R software package 'pheatmap' was used to visualize the result, and P<0.05 was considered statistically significant (*, P < 0.05; **, P < 0.01; ***, P < 0.001). NSCLC dysfunction and exclusion scores were gained from TIDE (http://tide.dfci.harvard.edu/). Analysis of TIDE scores in high and low-risk groups using 'limma' and 'ggpubr' R packages. For analysis of therapeutic agents and drug sensitivity, we use the 'pRRophetic,' 'ggpubr,' and 'ggplot2' R packages with p-value <10^-12.

Cell culture and quantitative real-time PCR (qRT-PCR)

Lung squamous cell carcinoma cell lines SK-MES-1 (# FH0040) and NCI-H226 (# FH0041) cells were purchased from Fuheng Biotechnology Co., Ltd. (FuHeng Cell Center, Shanghai, China). Under rigorous sterile conditions at 37°C and 5% CO₂, NCI-H226 was grown in Gibco 1640 medium (USA) containing 10% fetal bovine serum (BBI, Shanghai, China) and 1% penicillin/streptomycin combination (Solaibao, Beijing, China), and SK-MES-1 was grown in Gibco DMEM medium (USA) containing 10% fetal bovine serum and 1% penicillin/streptomycin combination. Following a 24-hour treatment with 30 nM elesclomol (MedChemexpress Co., LTD. # HY-12040) and Cucl₂ (Sigma-Aldrich. # 10125-13-0) (1:1), the cells were detected using qRT-PCR. TRIzol (Invitrogen, USA) was used to extract total RNA from the pretreated cells. With the use of the PrimeScript® RT Reagent Kit (TaKaRa, Japan), RNA was reverse-transcribed into cDNA. Lastly, using β-actin as a reference, qRT-PCR was utilized to detect the expression of lncRNAs linked to cuproptosis-related using TB Green Premix Ex Taq II (TaKaRa, Japan). Table 1 lists the primers for lncRNAs relevant to cuproptosis. Each sample was examined three times to determine the levels of lncRNA expression. The 2^{-(△CT sample-△CT control)} was employed to investigate the differential expression of lncRNAs.

Statistical Analysis

Kaplan-Meier survival analysis was performed with the log-rank test. Differential mRNA expression analysis was performed using the Wilcoxon rank sum test. Furthermore, the analysis of therapeutic agents and drug sensitivity was also evaluated by the Wilcoxon rank sum test. On the other hand, the correlation between drug and Risk Score was analyzed using Spearman correlation. The GraphPad Prism 8 software was used to analyze the qRT-PCR data. The p-values were determined using a t-test, and the findings were presented as Mean ± SEM. The results were considered statistically significant when p < 0.05.

 Figure 1 

The diagram of the research design.

Results

Construction of a cuproptosis-related LncRNA Risk Model for LUSC

The workflow of our studies is illustrated in Figure 1. RNA-seq data for 502 LUSC and 49 normal samples were obtained from the TCGA database. First, 291 cuproptosis-related LncRNAs were identified based on 19 known cuproptosis-related mRNAs (Supplementary Table 1). A total of 291 cuproptosis-related LncRNAs associated with 19 cuproptosis-related genes were picked up by the Pearson correlation coefficient test using the 'limma' R package. The criteria for LncRNA selection is based on P < 0.001 and correlation coefficient > 0.4 ​(Figure 2A). Of the 291 cuproptosis-related LncRNAs, 23 cuproptosis-related LncRNAs (P < 0.05) as LUSC independent prognostic factors were identified by univariate Cox regression analysis​ (Figure 2B, P < 0.05). Subsequently, to precisely establish the prognostic model, LASSO analysis was performed to pick up 12 LncRNAs from 23 LncRNAs for the subsequent analyses (Figure 2C, D). Finally, multivariate Cox proportional hazards regression analysis identified five cuproptosis-related LncRNAs used to build a prognostic model (Figure 2E). Risk Score = (-0.24159 × EXP AC002467.1) + (0.523456 × EXP MIR3945HG) + (0.41216 × EXP AC010328.1) + (0.17936 × EXP AC008972.2) + (-0.16671 × EXP AC253536.6).

Survival analysis using 5-cuproptosis-related lncRNA risk model

Four hundred ninety-three patients with LUSC were grouped into low-risk or high-risk groups based on the median Risk Score. The Kaplan-Meier survival analysis revealed that the overall survival time for patients in the high-risk Score group was poorer than those in the low-Risk Score group (Figure 3A, P < 0.001). To eliminate selection bias between high-risk and low-risk groups, we recruited 484 lung squamous cell carcinoma patients with complete follow-up information for PSM analysis between high-risk and low-risk groups. After PSM analysis using 'MatchIt' R software package, 170 patients in the low-risk group and 162 patients in the high-risk group were selected. The selection bias between high-risk and low-risk groups before and after PSM analysis was shown in Table 2.

 Figure 2 

Construction of a cuproptosis-related LncRNA Risk Model for LUSC. A. Co-expression of 16 mRNAs and 291 lncRNAs relevant to cuproptosis was visualized using a Sankey diagram. B. The forest plot showed 19 cuproptosis-related lncRNAs chosen by univariate Cox regression analysis. Red nodes represented high-risk lncRNAs. Green nodes expressed low-risk lncRNAs. C. The LASSO regression analysis showed that each independent variable's trajectory varies. D. Cuproptosis-related lncRNAs were selected using LASSO regression. E. Multivariate Cox analysis was used to find five cuproptosis-related lncRNAs and create the risk model.

Reanalyzing the survival assessment of high and low risk groups after eliminating bias by PSM revealed that the overall survival rate of the high-risk group was considerably lower than that of the low-risk group (Figure 3B, P = 0.002). The results above indicated the risk score has an impact on a patient's prognosis for lung squamous cell carcinoma. Furthermore, the PFS analysis was performed based on the TCGA and UCSC Xena databases. Although patients in the high-Risk Score group presented with worse PFS compared with those in the low-Risk Score group, there was no statistical difference between the two groups (Figure 3C, P = 0.052). For clinical stage I-II patients, overall survival of the high-Risk Score group was lower than that of the low-Risk Score group (Figure 3D, P < 0.001). Similarly, for clinical stage, II-III LUSC patients, the patients in the low-Risk Score group had a significantly higher survival probability compared with those in the high-Risk Score group (Figure 3E, P = 0.023). Due to the small number of LUSC patients in clinical stages III-IV, no statistical difference was observed between the high-Risk Score group and the low-Risk Score group (Figure 3F, P = 0.689). The AUC values for 1‐year, 3‐year, and 5‐year reached 0.624, 0.609, and 0.575, respectively, which indicated that the cuproptosis-related LncRNAs' prognostic signature had a good predictive ability (Figure 3G). The Heat map displayed the expressions of 5 cuproptosis-related LncRNAs (Figure 3H). The distribution of Risk Scores, risk groups, and survival status for 493 LUSC patients are shown in Figure 3I. Figure 3J shows we analyzed the correlations between the 5-cuproptosis-related LncRNAs and 19-cuproptosis-related genes.

Validation of the risk model

We randomly divided the LUSC patients into two groups: the training group (n = 198) and the testing group (n = 295), with a 2:3 split ratio. The Chi-square test showed no difference between the two randomized groups (Table 3). The training and testing sets were divided into high-risk and low-risk groups based on the median risk score. The upper panel of figure 4A shows the distributions of LUSC patients in the training set. The lower panel of Figure 4A displayed the relationship between the patient's survival time, survival status, and Risk Score in the training set. The 5-cuproptosis-related LncRNA expressions of the training set are shown in Figure 4B. For the training set, LUSC patients in the high-risk group had significantly shorter overall survival than those in the low-risk group (Figure 4C, P = 0.003). The AUC value for the training set predicting OS at 1, 3, and 5 years were 0.600, 0.656, and 0.638, respectively (Figure 4D). Figures 5E-F showed the distributions of LUSC patients with Risk Score, survival time, survival status (Figure 4E), and 5-cuproptosis-related LncRNA expressions (Figure 4F) in the testing set. Moreover, Kaplan-Meier survival analysis showed that the training set's outcomes were similar to those of the testing set (Figure 4G, P = 0.013). The value of AUC for the testing set predicting OS at 1, 3, and 5 years were, respectively, 0.640, 0.580, and 0.544 (Figure 4H).

The risk score serves as an independent prognostic indicator for LUSC patients

We further used univariate Cox regression analysis (Figure 5A, P < 0.05) and multivariate Cox regression analysis (Figure 5B, P < 0.05) to study whether risk score, age, gender, and clinical stages can be used as prognostic factors. We discovered that the Risk Score and clinical stages of patients with LUSC were independent prognostic variables (Figure 5A, B, P < 0.05). After that, ROC curve analysis found that the AUC for Risk Score, clinical stage, age, and gender is 0.624, 0.566, 0.530, 0.502, respectively (Figure 5C); that is, our proposed risk model possessed the best ability in prediction the survival time for LUSC patients. According to the risk score, we performed PCA to observe further the distribution of patients in the high-risk and low-risk groups. As shown in Figure 5D, when the 5-cuproptosis-related LncRNAs were applied to the principal component analysis (Figure 5D), the patients were divided into high-risk and low-risk groups. However, when all genes (Figure 5E, n = 59427), 19 cuproptosis-related mRNAs (Figure 5F), and 291 cuproptosis-related LncRNAs (Figure 5G) were subjected to the principal component analysis, respectively, cases from the high-risk and low-risk groups were combined, making it impossible to classify them.

Construction and validation of a Nomogram

A nomogram for LUSC patients was constructed based on age, gender, clinical stage, and Risk Score (Figure 6A). Calibration curves showed that the nomogram-predicted OS for 1-year, 3-year, and 5-year was generally consistent with the corresponding observed OS for LUSC patients (Figure 6B).

 Figure 3 

Survival analysis using 5-cuproptosis-related lncRNA risk model. A. The overall survival (OS) curves for LUSC patients (n = 493). B. The overall survival analysis after PSM (n = 332). C. The progression-free survival (PFS) curves for LUSC patients (n = 473). D-F. For clinical stage I and II patients (D), stage II and III patients (E), and stage III and IV patients (F), the overall survival of the high-Risk group and the low-Risk group was shown. G. Overall survival ROC curves at 1, 3, and 5 years. H. The Heat map displayed the expressions of 5 cuproptosis-related lncRNAs. I. The distribution of Risk Score, risk groups, and survival status for 493 LUSC patients was conducted. J. The heat map illustrated the correlations between the 5-cuproptosis-related lncRNAs and 19-cuproptosis-related genes. Statistics were considered significant when the p-value was less than 0.05.

 Figure 4 

Validation of the risk model. A. The upper panel showed the distributions of high-risk and low-risk LUSC patients in the training set. The lower panel displayed the relationship between the patient's survival time, survival status, and Risk Score in the training set. B. The heat map displayed the 5-cuproptosis-related lncRNA expressions in the training set. C.The overall survival of the high-Risk and low-Risk groups was shown for the training set. D. The AUC value for the training set predicting OS at 1, 3, and 5 years were 0.600, 0.656, and 0.638, respectively. E. The upper panel showed the distributions of high-risk and low-risk LUSC patients in the testing set. The lower panel displayed the relationship between the patient's survival time, survival status, and Risk Score in the testing set. F. The heat map displayed the 5-cuproptosis-related lncRNA expressions in the training set. G. The overall survival of the high-Risk group and the low-Risk group for the testing set was demonstrated. H. The value of AUC for the testing set predicting OS at 1, 3, and 5 years were, respectively, 0.640, 0.580, and 0.544.

 Figure 5 

The risk score serves as an independent prognostic indicator for LUSC patients. A-B. (A) Univariate and (B) Multivariate Cox regression analysis confirmed that the Risk Score and clinical stages of patients with LUSC were independent prognostic variables. C. the AUC for Risk Score, clinical stage, age, and gender is 0.624, 0.566, 0.530, and 0.502, respectively. D-G. 5-cuproptosis-related LncRNAs (D), all genes (E), 19 cuproptosis-related mRNAs (F), and 291 cuproptosis-related LncRNAs (G) were subjected to the principal component analysis. Statistics were considered significant when the p-value was less than 0.05.

 Figure 6 

Construction and validation of a Nomogram. A. A nomogram for LUSC patients was constructed based on age, gender, clinical stage, and Risk Score. Gender: gender of patient; Age: age of patient; Stage: clinical stage of the patient; risk: risk score based on the risk model. B. Calibration curves showed that the nomogram-predicted OS for 1-year, 3-year, and 5-year was generally consistent with the corresponding observed OS for LUSC patients.

Gene Set Enrichment Analysis and mutation analysis of LUSC patients

We next evaluated the differentially expressed mRNAs between the high-and low-risk groups of LUSC patients. According to the GO enrichment analysis, the differentially expressed mRNAs were enriched in several immune-related processes, including antigen binding, immunoglobulin complex formation, and humoral immune response (Figure 7A). The KEGG enrichment analysis indicated that the differentially expressed mRNAs were enriched in asthma (hsa05310), phagosome (hsa04145), cell adhesion molecules (hsa04514), complement & coagulation cascades (hsa04610), etc. (Figure 7B). We identified the top 15 genes with the highest mutation frequency and visualized them using the 'Maftools' R package. The top 15 genes had mutations in 97% (227/234) of the high-risk category samples (Figure 7C) and 97.92% (235/240) of the low-risk category samples (Figure 7D). According to the analysis of tumor mutational burden (TMB), there was no statistically significant difference in total mRNA mutation between high-risk and low-risk patients (Figure 7E). As a result, variations in survival rates between high-risk and low-risk groups might not be caused by mutations in these genes.

Immune-related functional analysis and drug sensitivity analysis

According to the ssGSEA algorithm, the enrichment score of the high-risk group is higher than that of the low-risk group in several immune function pathways, including the IFN-γ, the MHC I pathway, and others (Figure 8A, P < 0.001). The risk of tumor immune escape was calculated using the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm. The result revealed that the high-risk group had a higher TIDE score than the low-risk group, indicating a higher probability of immune escape. (Figure 8B, P < 0.001). We further explored potentially effective therapeutic drugs for LUSC patients using the 'pRRophetic' R package (Table 4).

 Figure 7 

Gene Set Enrichment Analysis and mutation analysis of LUSC patients. A-B. The differentially expressed mRNAs between the high- and low-risk groups of LUSC patients were examined using the GO enrichment analysis (A) and the KEGG enrichment analysis (B). C. A waterfall plot was used to display the top 15 mutant genes of LUSC in the high-risk (234 samples) categories. D. A waterfall plot was used to display the top 15 mutant LUSC genes in the low-risk (240 samples) categories.

 Figure 8 

Immune-related functional analysis and drug sensitivity analysis. A. Immune-related functional analyses were displayed on heat maps. B. Tumor Immune Dysfunction and Exclusion (TIDE) analyses in high-risk and low-risk groups. C-F. Drug sensitivity analyses in (C-D) low- and (E-F) high-risk groups. G-J. Correlation analyses between drug IC50 value and risk score. Statistics were considered significant when the * p-value was less than 0.05, ** p-value less than 0.01, and *** p-value less than 0.001.

The GW441756 (Figure 8C) and Salubrinal (Figure 8D) had a lower inhibitory concentration-50 (IC50) value in the low-Risk Score group (P < 0.001) (Figure 8C, D). On the other hand, Dasatinib (Figure 8E) and Z-LLNIe-CHO (Figure 8F) drugs had a lower IC50 value in the high-Risk Score group (P < 0.001) (Figure 8E, F). As shown in figure 8G-H, the risk scores were positively correlated with the IC50 values of the two pharmaceuticals, GW441756 and Salubrinal (Figure 8G, R = 0.42; Figure 8H, R = 0.4. P < 0.001), indicating that patients with low-risk ratings were likely to respond to these medications more effectively. Meanwhile, Figure 8I-J revealed the IC50 of Dasatinib (R =-0.39) and Z-LLNIe-CHO (R = -0.36) were negatively connected with Risk Scores (Figure 8I, J, P < 0.001), suggesting that LUSC patients with higher risk scores are more responsive to dasatinib and Z-LLNIe CHO (Figure 8I, R = -0.39; Figure 8J, R = -0.36).

LUSC cuproptosis cell model qRT-PCR detection

The cells treated with 30 nM elesclomol and Cucl₂ (1:1) mixture were collected to detected the expressions of LncRNAs using qRT-PCR. AC253536.6 expression was considerably higher in NCI-H226 cells, while AC008972.2, AC010328.1, and MIR3945HG expression was significantly lower (Figure 9A-D). Similarly, the SK-MES-1 cell type had the same results (Figure 9E-H). AC002467.1 expression was not detected in either of the two cell types, which could be attributed to the cells' low expression levels.

 Figure 9 

The expression levels of 5-cuproptosis-related lncRNAs in NCI-H226 and SK-MES-1 cells. After treated with 30 nM elesclomol and Cucl₂ (1:1) mixture, the expression levels of AC008972.2 (A), AC010328.1 (B), MIR3945HG (C), and AC253536.6 (D) were detected by quantitative Real-time PCR in NCI-H226 cells. Similarly, the expression levels of AC008972.2 (E), AC010328.1 (F), MIR3945HG (G), and AC253536.6 (H) were detected by quantitative Real-time PCR in SK-MES-1 cells. Statistics were considered significant when the * p-value was less than 0.05, ** p-value less than 0.01, *** p-value less than 0.001 and **** p-value less than 0.0001.

Discussion

Lung cancer has the second-highest incidence rate yet has the most remarkable fatality rate of all malignancies [19]. Lung cancer patients typically acquire medication resistance and are vulnerable to tumor metastasis at an early stage, which makes the clinical outcome of lung cancer treatment less than ideal [20]. Lung squamous cell carcinoma, which makes up around 40% of all lung malignancies, and lung adenocarcinoma are examples of the subtype of lung cancer known as non-small cell lung cancer [21, 22]. Since lung squamous cell carcinoma lacks effective treatment [23], the predicted survival rate of patients with lung squamous cell carcinoma does not reach 5% [24]. Therefore, exploring new therapeutic targets and molecular mechanisms of lung squamous cell carcinoma is still challenging.

Copper-induced cell death, also called cuproptosis, is a novel form of death that is gaining the interest of more and more scientists. According to studies, copper ions can currently target the tricarboxylic acid cycle pathway after entering the cell, leading to the aggregation of lipoylated proteins (DBT, DLAT, GCSH, and DLST), instability of the Fe-S cluster proteins (FDX1, LIAS, ACO-2, SDHB, POLD1, and DPYD), and oligomerization of DLAT [8]. It has been discovered that malignant tumors have much higher copper content [11]. According to animal studies, chelating medicines that include copper ions can slow the growth of various cancers [25, 26]. So far, the molecular mechanism of copper death causing cancer has rarely been reported. Based on bioinformatics, many researchers have predicted the genetic markers of copper death in tumors [20, 27]. Yang et al. confirmed the increased expression of copper death-related genes (LIPT1, DLD, PDHA1, DLAT, FDX1, CDKN2A) in osteosarcoma [28]. Zhang et al. additionally foresaw the crucial function of the copper death critical gene FDX1 in hepatocellular cancer [29].

Long non-coding RNA (LncRNA) is a long-chain non-coding RNA with a length of more than 200 nucleotides, which is an essential part of the non-coding genome. They Influence cell function through different patterns [30]. Researchers found that LncRNAs play an indispensable role in various cancers, such as breast cancer, non-small cell lung cancer, etc. [31, 32]. Chen et al. discovered that the miR-877-3p/FGF2 axis controls cervical cancer migration and invasion through the LncRNA HOXD-AS1[33]. Lin et al. showed that LncRNA FGD5-AS1 could bind miR-520a-3p to target KIAA1522 and affect the proliferation of lung adenocarcinoma cells [34]. Wang et al. developed a LUAD predictive signature in the most recent study that was published by thoroughly examining the TCGA and GEO databases [35].

Our study thoroughly analyzed the TCGA database's cuproptosis-related mRNA expression profiling of LUSC patients. Pearson correlation analysis revealed 291 LncRNA that had significant relationships with 19 cuproptosis-related genes (Figure 2A, P < 0.001). Of the 291 cuproptosis-related LncRNAs, 23 cuproptosis-related LncRNAs (P < 0.05) as LUSC independent prognostic factors were identified by univariate Cox regression analysis (Figure 2B, P < 0.05). Finally, multivariate Cox proportional hazards regression analysis identified five cuproptosis-related LncRNAs (MIR3945HG, AC002467.1, AC008972.2, LncRNA AC010328.1, and AC253536.6) that were used to build a prognostic model (Figure 2E). According to Chen et al., high-expressed MIR3945HG levels were strongly connected with shorter survival of patients with LUSC [36]. Wu et al. discovered that the LncRNA AC002467.1, which was found to be abnormally expressed in ovarian cancer, had a negative relationship with the survival of ovarian cancer patients [37]. AC008972.2 has been linked to m7G RNA methylation in colon cancer, according to Liu et al. [38]. However, the role of LncRNA AC010328.1 and AC253536.6 in diseases has rarely been reported.

Four hundred ninety-three patients with LUSC were grouped into low-risk or high-risk groups based on the median Risk Score acquired from the risk model we built. The Kaplan-Meier survival analysis revealed that the survival time for patients in the high-risk Score group was poorer than those in the low-Risk Score group (Figure 3A, P < 0.001). Next, all LUSC patients were divided into three groups according to clinical stages (clinical stage I-II, II-III o low- and high-Risk Score subgroups. For clinical stage I-II patients, the overall survival of III-IV and each group was subgrouped in the high-Risk Score subgroup was lower than that of the low-Risk Score subgroup both for clinical stage I-II patients and II-III patients (Figure 3D, P < 0.001; 3E, P = 0.023). However, because only a few patients were enrolled in clinical stages III-IV, there was no statistically significant difference in overall survival between the high and low-Risk Score subgroups. (Figure 3F, P = 0.689). Additionally, an examination of five LncRNA profiles with mRNA correlations (Figure 3J) reveals LncRNA AC002467.1 was highly negatively correlated with NLRP3 (P < 0.001) and positively correlated with cuprotosis-related mRNA DLD (P < 0.001). In contrast, LncRNA AC008972.2 positively correlated with NLRP3 (P < 0.001). NLRP3 is an inflammatory body activated during mitochondrial dysfunction, resulting in alcoholic and nonalcoholic hepatitis [39]. The copper transporter ATP7B and the LncRNA AC010328.1 positively correlated (P < 0.001). ATP7B interacted with the autophagy marker LC3 to induce autophagy, and ATP7B on the autophagic membrane promoted Cu ion transfer into autophagosome [40]. MIR3945HG was found to be negatively correlated with several-copper death-related genes (NFE2L2, LIPT1, GCSH) (P < 0.001). LIPT1 has been found to promote hepatocellular carcinoma proliferation, invasion, and migration [41]. Since LncRNAs regulate gene expression through various complex pathways, the mechanism of the above five LncRNAs affecting 19-copper-death-related genes needs to be investigated further.

All of the LUSC samples were randomly divided into training (40%) and testing (60%) sets to verify the credibility of our risk model. The training and testing sets were divided into high-risk and low-risk groups based on the median risk score. For the training set, LUSC patients in the high-risk group had significantly shorter overall survival compared with those in the low-risk group (P = 0.003) (Figure 4C). Moreover, the training set's outcomes were similar to those of the testing set (P = 0.013) (Figure 4G). Following that, univariate and multivariate COX regression analysis revealed that clinical stage and risk score are independent prognostic indicators for LUSC patients (Figure 5A, B, P < 0.05). More importantly, risk scores outperform the clinical stage in terms of predictive power (Figure 5C). Following that, we created a nomogram to predict the prognosis of LUSC patients (Figure 6A, B). Each patient's survival rate of 1/3/5 year can be predicted based on age, gender, clinical stage, and Risk Score (See Figure 6A). Our nanogram will be a valuable tool for clinicians in predicting the prognosis of patients with LUSC.

The enrichment analysis of differentially expressed genes in high and low-risk groups reveals that the majority of differentially expressed genes are enriched in immune-related GO terms like leukocyte-mediated immune response, humoral immune response, B cell-mediated immune response, and other immune biological processes (Figure 7A). Recently, Zhang et al. discovered that circRNA (circHMGB2) overexpression reshaped the tumor microenvironment and regulated anti-PD-1 resistance in patients with lung squamous cell carcinoma, thereby promoting the progression of LUSC [1]. According to KEGG enrichment analysis, these mRNAs were enriched in the phagosome (hsa04145) and cell adhesion molecule (hsa04514) pathways (Figure 7B). Su et al. discovered that AIM2, a novel inflammasome, was recruited into phagosomes and then up-regulated PD-L1 to cause immunosuppression [42]. Adhesion molecule expression inhibited T cell infiltration into tumors, posing a challenge to cancer immunotherapy [43]. We discovered no discernible difference in the tumor mutational burden (TMB) between the high-risk and low-risk groups (Figure 7E), implying that mutations in the differentially expressed mRNAs had little impact on the prognosis of LUAD patients.

Notably, we used the "pRRophetic" R package to investigate potentially effective therapeutic drugs for LUSC patients (Table 4). The drug sensitivity analysis showed that patients with low-risk ratings were likely to respond to GW441756 and Salubrinal. In contrast, patients with higher risk scores are more responsive to Dasatinib and Z-LLNIe CHO. Dasatinib is currently used to treat chronic myeloid leukemia (CML), but it has been linked to pulmonary pleural effusion (PE) and pulmonary hypertension (PAH) [44]. Interestingly, researchers discovered that Dasatinib inhibited the activity of the MCM7 protein, making it a promising cancer therapy [45]. The proteasome inhibitor MG-132 has a similar structure to secretase inhibitor I (GSI-I; Z-LLNIe CHO) and can inhibit the proteasome, which is a crucial cancer treatment candidate [46].

In conclusion, we developed a predictive risk model based on five copper death-related LncRNAs for LUSC patients. To forecast the prognosis of LUSC patients, we constructed a nomogram. Additionally, the drug sensitivity analysis revealed several candidates responsive to LUSC patients at high or low risk.

Supplementary Material

Supplementary tables.

Acknowledgements

Funding

This work was supported by Guangxi Science and Technology Project for Bases and Talents (Guike AD20297024), the Innovation Project of Guangxi Graduate Education (YCSW2021258), the Graduate Research Program of Guilin Medical University (GYYK2022009), the National Natural Science Foundation of China (30901712, 81872221, 82160408, 81860525).

Availability of data and materials

The data used and analyzed during our study are available from the corresponding author on reasonable request.

Author contributions

Juan Du and Xiuzhen Yang designed the study. Hongtao Zhao was a significant contributor to writing the manuscript. Juan Du, Xiuzhen Yang, and Hongtao Zhao made substantial contributions to the design of the manuscript and revised it critically for important intellectual content. Hongtao Zhao created all the figures for the article. Lei Wu, Peiluo Huang, Qinyuan Liao and Ruonan Sun collected relevant data from TCGA. All authors have read and agreed to the published version of the manuscript.

Competing Interests

The authors have declared that no competing interest exists.

References

1. Zhang LX, Gao J, Long X, Zhang PF, Yang X, Zhu SQ. et al. The circular RNA circHMGB2 drives immunosuppression and anti-PD-1 resistance in lung adenocarcinomas and squamous cell carcinomas via the miR-181a-5p/CARM1 axis. Mol Cancer. 2022;21:110

2. Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer Statistics, 2021. CA Cancer J Clin. 2021;71:7-33

3. Langer CJ, Obasaju C, Bunn P, Bonomi P, Gandara D, Hirsch FR. et al. Incremental Innovation and Progress in Advanced Squamous Cell Lung Cancer: Current Status and Future Impact of Treatment. J Thorac Oncol. 2016;11:2066-81

4. Zarogoulidis K, Zarogoulidis P, Darwiche K, Boutsikou E, Machairiotis N, Tsakiridis K. et al. Treatment of non-small cell lung cancer (NSCLC). J Thorac Dis. 2013;5(Suppl 4):S389-96

5. Yu J, Fan Z, Zhou Z, Zhang P, Bai J, Li X. et al. TP53 and LRP1B Co-Wild Predicts Improved Survival for Patients with LUSC Receiving Anti-PD-L1 Immunotherapy. Cancers (Basel). 2022;14:3382

6. Ruiz LM, Libedinsky A, Elorza AA. Role of Copper on Mitochondrial Function and Metabolism. Front Mol Biosci. 2021;8:711227

7. Ge EJ, Bush AI, Casini A, Cobine PA, Cross JR, DeNicola GM. et al. Connecting copper and cancer: from transition metal signalling to metalloplasia. Nat Rev Cancer. 2022;22:102-13

8. Tsvetkov P, Coy S, Petrova B, Dreishpoon M, Verma A, Abdusamad M. et al. Copper induces cell death by targeting lipoylated TCA cycle proteins. Science. 2022;375:1254-61

9. Li SR, Bu LL, Cai L. Cuproptosis: lipoylated TCA cycle proteins-mediated novel cell death pathway. Signal Transduct Target Ther. 2022;7:158

10. Zhang X, Yang Q. Association between serum copper levels and lung cancer risk: A meta-analysis. The Journal of international medical research. 2018;46:4863-73

11. Oliveri V. Selective Targeting of Cancer Cells by Copper Ionophores: An Overview. Front Mol Biosci. 2022;9:841814

12. Lelievre P, Sancey L, Coll JL, Deniaud A, Busser B. The Multifaceted Roles of Copper in Cancer: A Trace Metal Element with Dysregulated Metabolism, but Also a Target or a Bullet for Therapy. Cancers (Basel). 2020;12:3594

13. Nair L, Chung H, Basu U. Regulation of long non-coding RNAs and genome dynamics by the RNA surveillance machinery. Nat Rev Mol Cell Biol. 2020;21:123-36

14. Fatica A, Bozzoni I. Long non-coding RNAs: new players in cell differentiation and development. Nat Rev Genet. 2014;15:7-21

15. Peng WX, Koirala P, Mo YY. LncRNA-mediated regulation of cell signaling in cancer. Oncogene. 2017;36:5661-7

16. Bhan A, Mandal SS. Long noncoding RNAs: emerging stars in gene regulation, epigenetics and human disease. ChemMedChem. 2014;9:1932-56

17. Vitiello M, Tuccoli A, Poliseno L. Long non-coding RNAs in cancer: implications for personalized therapy. Cell Oncol (Dordr). 2015;38:17-28

18. Zhen Q, Gao LN, Wang RF, Chu WW, Zhang YX, Zhao XJ. et al. LncRNA DANCR Promotes Lung Cancer by Sequestering miR-216a. Cancer Control. 2018;25:1073274818769849

19. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A. et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin. 2021;71:209-49

20. Liu Y, Liu Y, Ye S, Feng H, Ma L. Development and validation of cuproptosis-related gene signature in the prognostic prediction of liver cancer. Front Oncol. 2022;12:985484

21. Li Y, Gu J, Xu F, Zhu Q, Ge D, Lu C. Transcriptomic and functional network features of lung squamous cell carcinoma through integrative analysis of GEO and TCGA data. Sci Rep. 2018;8:15834

22. Bezzecchi E, Ronzio M, Dolfini D, Mantovani R. NF-YA Overexpression in Lung Cancer: LUSC. Genes (Basel). 2019;10:937

23. Hirsch FR, Scagliotti GV, Mulshine JL, Kwon R, Curran WJ Jr, Wu YL. et al. Lung cancer: current therapies and new targeted treatments. Lancet. 2017;389:299-311

24. Cheng Z, Yu C, Cui S, Wang H, Jin H, Wang C. et al. circTP63 functions as a ceRNA to promote lung squamous cell carcinoma progression by upregulating FOXM1. Nat Commun. 2019;10:3200

25. Lopez J, Ramchandani D, Vahdat L. Copper Depletion as a Therapeutic Strategy in Cancer. Metal ions in life sciences. 2019;19:9783110527872-018

26. Shanbhag VC, Gudekar N, Jasmer K, Papageorgiou C, Singh K, Petris MJ. Copper metabolism as a unique vulnerability in cancer. Biochim Biophys Acta Mol Cell Res. 2021;1868:118893

27. Zhang G, Sun J, Zhang X. A novel Cuproptosis-related LncRNA signature to predict prognosis in hepatocellular carcinoma. Sci Rep. 2022;12:11325

28. Yang M, Zheng H, Xu K, Yuan Q, Aihaiti Y, Cai Y. et al. A novel signature to guide osteosarcoma prognosis and immune microenvironment: Cuproptosis-related lncRNA. Front Immunol. 2022;13:919231

29. Zhang Z, Zeng X, Wu Y, Liu Y, Zhang X, Song Z. Cuproptosis-Related Risk Score Predicts Prognosis and Characterizes the Tumor Microenvironment in Hepatocellular Carcinoma. Front Immunol. 2022;13:925618

30. Gutschner T, Diederichs S. The hallmarks of cancer: a long non-coding RNA point of view. RNA Biol. 2012;9:703-19

31. Alkhathami AG, Hadi A, Alfaifi M, Alshahrani MY, Verma AK, Beg MMA. Serum-Based lncRNA ANRIL, TUG1, UCA1, and HIT Expressions in Breast Cancer Patients. Dis Markers. 2022;2022:9997212

32. Tao Y, Tang Y, Yang Z, Wu F, Wang L, Yang L. et al. Exploration of Serum Exosomal LncRNA TBILA and AGAP2-AS1 as Promising Biomarkers for Diagnosis of Non-Small Cell Lung Cancer. Int J Biol Sci. 2020;16:471-82

33. Chen S, Li K. HOXD-AS1 facilitates cell migration and invasion as an oncogenic lncRNA by competitively binding to miR-877-3p and upregulating FGF2 in human cervical cancer. BMC Cancer. 2020;20:924

34. Lin J, Liao S, Liu Z, Li E, Wu X, Zeng W. LncRNA FGD5-AS1 accelerates cell proliferation in pancreatic cancer by regulating miR-520a-3p/KIAA1522 axis. Cancer Biol Ther. 2021;22:257-66

35. Wang F, Lin H, Su Q, Li C. Cuproptosis-related lncRNA predict prognosis and immune response of lung adenocarcinoma. World J Surg Oncol. 2022;20:275

36. Chen WJ, Tang RX, He RQ, Li DY, Liang L, Zeng JH. et al. Clinical roles of the aberrantly expressed lncRNAs in lung squamous cell carcinoma: a study based on RNA-sequencing and microarray data mining. Oncotarget. 2017;8:61282-304

37. Wu W, Gao C, Chen L, Zhang D, Guo S. Comprehensive analysis of competitive endogenous RNAs networks reveals potential prognostic biomarkers associated with epithelial ovarian cancer. Oncol Lett. 2021;22:843

38. Liu L, Wu Y, Chen W, Li Y, Yu J, Zhang G. et al. The m7G-Related Long Noncoding RNA Signature Predicts Prognosis and Indicates Tumour Immune Infiltration in Colon Cancer. Front Genet. 2022;13:892589

39. Torres S, Segales P, Garcia-Ruiz C, Fernandez-Checa JC. Mitochondria and the NLRP3 Inflammasome in Alcoholic and Nonalcoholic Steatohepatitis. Cells. 2022;11:1475

40. Pantoom S, Pomorski A, Huth K, Hund C, Petters J, Krezel A. et al. Direct Interaction of ATP7B and LC3B Proteins Suggests a Cooperative Role of Copper Transportation and Autophagy. Cells. 2021;10:3118

41. Yan C, Niu Y, Ma L, Tian L, Ma J. System analysis based on the cuproptosis-related genes identifies LIPT1 as a novel therapy target for liver hepatocellular carcinoma. J Transl Med. 2022;20:452

42. Su S, Zhao J, Xing Y, Zhang X, Liu J, Ouyang Q. et al. Immune Checkpoint Inhibition Overcomes ADCP-Induced Immunosuppression by Macrophages. Cell. 2018;175:442-57 e23

43. Harjunpaa H, Llort Asens M, Guenther C, Fagerholm SC. Cell Adhesion Molecules and Their Roles and Regulation in the Immune and Tumor Microenvironment. Front Immunol. 2019;10:1078

44. Nekoukar Z, Moghimi M, Salehifar E. A narrative review on adverse effects of dasatinib with a focus on pharmacotherapy of dasatinib-induced pulmonary toxicities. Blood Res. 2021;56:229-42

45. Samad A, Huq MA, Rahman MS. Bioinformatics approaches identified dasatinib and bortezomib inhibit the activity of MCM7 protein as a potential treatment against human cancer. Sci Rep. 2022;12:1539

46. Meng X, Matlawska-Wasowska K, Girodon F, Mazel T, Willman CL, Atlas S. et al. GSI-I (Z-LLNle-CHO) inhibits gamma-secretase and the proteosome to trigger cell death in precursor-B acute lymphoblastic leukemia. Leukemia. 2011;25:1135-46

Author contact

Corresponding authors: yangxz770com; sunnydujuanedu.cn.

Received 2023-1-4
Accepted 2023-5-12
Published 2023-5-21

Citation styles

©2024 Ivyspring International Publisher. Terms of use