A novel cuproptosis-related lncRNA signature to predict prognosis and immune landscape of lung adenocarcinoma

Introduction

Lung cancer accounts for the highest cancer-related mortality rate worldwide, of which lung adenocarcinoma (LUAD) is the most common subtype (1). For all stages of lung cancer combined, the 5-year survival rate is around 22%. Despite progress made in the fields of surgery, chemotherapy, radiotherapy, and immunotherapy, there is still much room for improvement in cancer treatment. Thus, novel molecular biomarkers provide additional prognostic value and better risk assessment to guide personalized therapy.

The existence and physiological significance of programmed cell death (PCD) have garnered increasing interest in recent years, including necroptosis, apoptosis, pyroptosis, and ferroptosis (2,3). PCD has shown great potential in cancer treatment, and the induction of cell deaths combined with immune checkpoint inhibitors (ICIs) has shown synergistically enhanced antitumor efficacy (3). Copper-induced cell death, termed cuproptosis, is recently recognized as a new form of PCD (4). Unlike other types of cell death, cuproptosis occurs through the direct binding of copper to lipoylated components of the tricarboxylic acid cycle, leading to the aggregation of lipoylated protein and the subsequent loss of iron-sulfur cluster protein, finally causing proteotoxic stress and cell death. Copper, a redox-active metal ion, is essential for most animals. It serves as a catalytic and structural cofactor for enzymes and plays a vital role in oxygen transport, cellular metabolism, energy generation, and signal transduction (5). Copper concentrations are maintained at a low level under normal physiological conditions by a homeostatic mechanism (6). Abnormal accumulation of copper is frequently observed in several cancers, leading to tumor proliferation, angiogenesis, and metastasis (7,8). SLC3A1, known as CTR1, is a transmembrane copper transporter protein. Research showed that SLC3A1-dependent copper level was associated with the malignant degree of pancreatic cancer (9). In addition, the combination of copper and copper ionophore disulfide has been found to enhance the sensitivity of patients with head and neck squamous cell carcinoma to cisplatin (10). Therefore, the underlying mechanism of cuproptosis discovered by Tsvetkov et al. presents novel avenues for applying cuproptosis in cancer treatment (4). Considering the great potential of cuproptosis in cancer treatment, cuproptosis-related genes are promised to be novel therapeutic targets. The investigation of molecular characteristics of cuproptosis in cancer patients may provide new prognostic and therapeutic methods.

Long non-coding RNAs (lncRNAs) are non-coding RNA transcripts longer than 200 nucleotides (11). Numerous studies have noted that lncRNAs participate in various biological processes of LUAD (12-14). For example, the novel lncRNA lung cancer associated transcript 3 (LCAT3) promoted proliferation and metastasis of lung cancer cells via the recruitment of far upstream element binding protein 1 (FUBP1) to activate c-MYC (13). Recently, there has been a surge of interest in the relationship between lncRNAs and PCD (15). Various studies have shown that PCD-related lncRNAs could be used as prognostic markers to predict the response to immunotherapy and patient outcomes (16-18). Thus, the identification of cuproptosis-related lncRNAs (CuRLs) is significant for elucidating the molecular mechanism of LUAD and helping us realize the potential of cuproptosis in cancer treatment.

The present study identified CuRLs from The Cancer Genome Atlas (TCGA) dataset and constructed a risk signature based on 10 prognostic CuRLs. We developed and validated a nomogram for clinical use that integrated the 10-CuRLs signature and clinical factors. The applicability of this risk signature was validated in the Gene Expression Omnibus (GEO) meta-cohort. In addition, gene set variation analysis (GSVA), gene set enrichment analysis (GSEA), Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were performed to explore the underlying biological mechanisms of the 10-CuRLs signature comprehensively. Then, the differences between low- and high-risk groups were compared in immune infiltration, tumor mutation burden (TMB), and drug response, shedding light on the potential role of CuRLs in cancer treatment. We present the following article in accordance with the TRIPOD reporting checklist (available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-22-500/rc).

Methods

Data acquisition

RNA sequencing data and clinical information for LUAD were obtained from TCGA (https://cancergenome.nih.gov/) database, which included 526 tumor samples and 59 normal samples. Patients with incomplete follow-up information and survival time less than 30 days were excluded to reduce statistical bias in the following analysis. A total of 472 LUAD patients were included in the further analysis. The mutation data of the included LUAD patients were also obtained from TCGA. For external validation, gene expression and clinical features of GSE31210 (N=226), GSE37745 (N=196), GSE50081 (N=181) were downloaded from GEO database (https://www.ncbi.nlm.nih.gov/geo/). They were integrated into a GEO meta-cohort for further analysis. To correct the batch effects and normalization, the R package ‘sva’ was used. Patients with other pathological types and survival time less than 30 days were excluded and finally 458 LUAD patients were included in the GEO meta-cohort. For subsequent models to be validated in the GEO meta-cohort, lncRNAs expressed in both TCGA and GEO datasets were identified. The clinical features of included LUAD patients in TCGA and GEO meta-cohort are shown in Table S1. Ethics approval was not required, as TCGA and GEO are public databases. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Identification of CuRLs

Tsvetkov et al. recently discovered several genes essential for cuproptosis (4). Those genes identified in either Cu-DDC or elesclomol-Cu group were also included. In addition, the copper importer SLC31A1, copper exporters ATP7A/ATP7B, and lipoylated proteins were considered cuproptosis-related genes. Although experiments have demonstrated that Fe-S cluster proteins were downregulated after elesclomol treatment, the particular mechanism was unclear. Thus, Fe-S cluster proteins were not included in our analysis. Following that, Pearson correlation analysis was performed between cuproptosis-related genes and lncRNAs. Those lncRNAs with |correlation coefficient| >0.3 and P<0.001 were identified as CuRLs.

Establishment of the risk model

Combined with LUAD survival data from TCGA, univariate Cox regression was utilized to screen for prognostic cuproptosis-related lncRNAs. The lncRNAs with P<0.05 were further analyzed. Least Absolute Shrinkage and Selection Operator (LASSO) Cox regression as implemented in the R package ‘glmnet’ was then used to identify lncRNAs associated with the prognosis of LUAD with lambda.min as the cut-off. Stepwise multivariate Cox regression was further applied to examine the prognostic CuRLs. The stepwise approach was based on the Akaike information criterion (AIC) to define the optimal set of variables to retain in each model minimizing the AIC value (19). Finally, a 10-CuRLs signature was established. Each patient’s risk score was calculated using the following formula:

Riskscore=∑i=1nExpression(lncRNAn)×Coefficient(lncRNAn)[1]

Bioinformatic analysis of identified CuRLs

The co-expression network was visualized by Cytoscape (version 3.9.1, http://www.cytoscape.org). Sankey diagram was obtained using the R package ‘ggalluvial’. A heatmap of genes between tumor and normal samples was constructed by R package ‘ComplexHeatmap’ and ‘circlize’ (20,21).

Evaluation of the risk model

TCGA-LUAD patients were divided into low- and high-risk groups using the median risk score as the cut-off. Kaplan-Meier survival was applied to evaluate the differences in overall survival (OS) between the low- and high-risk patients. Univariate and multivariate Cox regression analyses were used to investigate whether the risk model based on identified cuproptosis-related lncRNAs was an independent risk factor considering other clinical signatures (age, gender, and stage) in LUAD patients. A predictive nomogram was developed using the R package ‘rms’. The calibration curve and time-dependent receiver operating characteristic curve (ROC) were used to evaluate the predictive ability of the model. Principal component analysis (PCA) was performed to evaluate the grouping ability of the whole genome, cuproptosis genes, all CuRLs, and the identified 10-CuRLs signature.

External validation of the risk model

The GEO meta-cohort was considered a validation set to evaluate the model’s performance. The risk score for each patient was calculated by the same formula used for the TCGA cohort. Patients were also divided into low- and high-risk groups according to the median risk score. The risk model was analyzed to determine whether it was an independent prognostic factor in the validation cohort using similar methods in TCGA-LUAD training analysis.

Functional enrichment analysis

GSVA was applied to explore differences in enrichment scores between low-risk and high-risk groups, using pathways from the Molecular Signatures Database (MSigDB) hallmark set as a reference. GSEA was also performed for each lncRNA of the 10 CuRLs. The correlation between the expression of each lncRNA and other genes was computed. Gene lists were sorted by the value of correlation and then GSEA was performed based on KEGG pathway database with P value cutoff set to 0.05. Protein-coding genes co-expressed with the 10 CuRLs were identified using |Pearson correlation coefficient| >0.5 and P<0.001 as the cut-off. GO and KEGG functional enrichment analyses were carried out to understand the functions of these correlated mRNAs. These analyses were conducted mainly based on R packages ‘GSVA’, ‘msigdbr’, ‘clusterProfiler’ (22), ‘org.Hs.eg.db’, and ‘enrichplot’.

Estimation of tumor immune microenvironment and TMB

Widely acknowledged methods were used for the evaluation of immune infiltration levels in LUAD patients, including the tumor immune estimation resource (TIMER), (23), cell-type identification by estimating relative subsets of RNA transcripts (CIBERSORT) (24), xCELL (25), QuanTIseq (26), microenvironment cell populations-counter (MCP-counter) (27), and estimating the proportions of immune and cancer cells (EPIC) algorithms (28). The Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) algorithm was also used to calculate the immune and stromal activity between low- and high-risk groups (29). To explore the differential expression of immune checkpoint genes (ICGs) between different risk patients, 79 ICGs were obtained from an article focusing on the roles of ICGs in predicting immunotherapy response (30). The R package ‘maftools’ was utilized to calculate the TMB of LUAD patients (31).

Significance of the CuRLs signature in drug sensitivity

To evaluate the differences in drug sensitivity between low- and high-risk groups, the R package ‘pRRophetic’ was applied which could yield a drug sensitivity for each patient based on the gene expression (32). The dataset ‘cpg2016’ in ‘pRRophetic’ included 251 types of drugs and we analyzed common anti-cancer drugs according to our clinical experience and review of related literature (18). The differences in drug sensitivity between low- and high-risk groups were represented by boxplots.

Statistical analysis

Data analysis was accomplished with R (version 4.1.2, the R Foundation for Statistical Computing, Vienna, Austria). Comparisons between 2 groups were evaluated using Wilcoxon test. Kruskal-Wallis test was used to evaluate the differences among more than 2 groups. A P<0.05 was considered significant.

Results

Identification of CuRLs in patients with LUAD

The flow chart in Figure 1 depicts the primary design of our study. First, we identified 7,151 lncRNAs in LUAD samples from TCGA by excluding genes expressed in less than half of the cases. By intersecting these lncRNAs from TCGA and genes obtained in the GEO meta-cohort, 1,219 lncRNAs were identified for later investigation. Based on the excellent research by Tsvetkov et al. (4), we identified 17 cuproptosis-related genes (ATP7A, ATP7B, CDKN2A, DBT, DLAT, DLD, DLST, FDX1, GCSH, GLS, LIAS, LIPT1, MTF1, PDHA1, PDHB, SLC31A1, LIPT2) and retrieved their expression data in TCGA-LUAD. Then, 201 CuRLs were identified using Pearson correlation analysis (|correlation coefficient| >0.3 and P<0.001, Table S2). The co-expression network of cuproptosis-lncRNA was displayed as a Sankey diagram (Figure 2A).

Figure 1 Flow chart of the present study. CuRLs, cuproptosis-related lncRNAs; LASSO, least absolute shrinkage and selection operator; LUAD, lung adenocarcinoma; OS, overall survival; TCGA, The Cancer Genome Atlas; TMB, tumor mutation burden.

Figure 2 Identification of prognostic CuRLs in LUAD patients. (A) Sankey diagram for the cuproptosis genes and CuRLs. (B) Selection of tuning parameter lambda in the LASSO Cox regression model using ten-fold cross-validation. (C) LASSO coefficient profiles of 36 CuRLs. (D) Stepwise multivariate Cox regression analysis constructed a 10-CuRLs prognostic signature. (E) Heatmap for the correlation between cuproptosis genes and the 10 CuRLs. (F) A co-expression network of cuproptosis genes (orange diamond) and 10 CuRLs (blue rectangle). Red lines represent positive correlation, and blue lines represent negative correlation. (G) The correlation between the 10 CuRLs using Pearson analysis. Red squares represent positive correlation, blue squares represent negative correlation and white squares represent no significant correlation. (H) Hierarchical clustering for expression of the 10 CuRLs between normal and tumor samples. ***, P<0.001; **, P<0.01; *, P<0.05. CI, confidence interval; CuRLs, cuproptosis-related lncRNAs; HR, hazard ratio; LASSO, least absolute shrinkage and selection operator; LUAD, lung adenocarcinoma.

Identification of prognostic CuRLs

Using univariate Cox regression, we identified 11 lncRNAs that were associated with poor prognosis and 25 lncRNAs that were associated with better survival outcomes (P<0.05, Table S3). LASSO Cox regression, as a common method of multiple regression method, can prevent overfitting and enhance the generalization ability of the model by adding penalties. Using the lambda. min as the cut-off threshold, we identified 21 lncRNAs that significantly correlated to LUAD prognosis (Figure 2B,2C). Then, stepwise multivariate Cox analysis was used to analyze these lncRNAs, 10 of which were finally identified as a prognostic CuRLs signature (Figure 2D). LINC01137 and TMPO-AS1 were identified as poor prognostic factors, whereas SEPSECS-AS1, CARD8-AS1, and ZSCAN16-AS1 were considered protective. The correlation between these 10 lncRNAs and cuproptosis genes is shown in Figure 2E,2F. We noted that most lncRNAs were both positively and negatively correlated with cuproptosis genes. Furthermore, we analyzed the correlations of the 10 CuRLs and found positive correlations between APTR and ELOA-AS1, ZSCAN16-AS1 and ELOA-SA1, SEPSECS-AS1 and DLGAP1-AS1, SEPSECS-AS1 and ELOA-AS1 (|correlation coefficient| >0.3 and P<0.05, Figure 2G). These indicated a complex regulatory network among cuproptosis genes and identified 10 CuRLs. We also examined the expression of these 10 CuRLs between normal and tumor tissues of LUAD (Figure 2H, Figure S1). RUNDC3A-AS1, TMPO-AS1, ELOA-AS1, LINC01137, DLGAP1-AS1, MIR31HG, and APTR were highly expressed in tumor samples, whereas the expression of CARD8-AS1 was higher in normal tissues. The result suggested a role of these 10 CuRLs in classifying LUAD patients.

Construction and evaluation of the risk model based on prognostic CuRLs

The identified 10 prognostic CuRLs were used to construct a risk model. LUAD samples were divided into low-risk (N=236) and high-risk (N=236) groups with the median risk score as the threshold value. The distribution of the risk score and survival outcomes between the low-risk and high-risk groups are shown in Figure 3A,3B. The relative expression of the 10 prognostic CuRLs for each patient is shown in Figure 3C. The survival analysis demonstrated that the OS of the low-risk group was significantly longer than that of the high-risk group (P<0.001, Figure 3D). An up-regulated expression of ATPR, SEPSECS-AS1, ELOA-AS1, ZSCAN16-AS1, CARD8-AS1, and RUNDC3A-AS1 was observed in the low-risk group whereas the expression of MIR31HG, TMPO-AS1, DLGAP1-AS1, and LINC01137 was higher in the high-risk group (Figure S2).

Figure 3 Prognostic value of the risk model based on 10 CuRLs. (A) Risk score distribution of patients in TCGA-LUAD cohort. (B) Survival status and time in the low- and high-risk groups. (C) Hierarchical clustering analysis of 10 CuRLs between the low- and high-risk groups. (D) Kaplan-Meier curve for OS in the low- and high-risk groups. Principal component analysis of low- and high-risk groups based on the whole-genome (E), cuproptosis genes (F) and all CuRLs (G). (H) The risk model based on 10 CuRLs. CuRLs, cuproptosis-related lncRNAs; LUAD, lung adenocarcinoma; OS, overall survival; TCGA, The Cancer Genome Atlas.

We conducted PCA to compare the low-risk and high-risk groups in order to verify the grouping capability of the risk model based on the identified 10 CuRLs. Figure 3E-3H display the distributions of low-risk and high-risk samples based on the entire genome, cuproptosis genes, all CuRLs, and the risk model. The results showed that the constructed risk model could distinguish the low-risk and high-risk patients much more effectively, indicating that the 10-CuRLs signature was a significant prognostic factor for LUAD patients.

Next, we performed univariate and multivariate Cox regression to evaluate the independence of the risk model based on the 10 CuRLs. The hazard ratio (HR) of the risk score was 2.718, and the 95% confidence interval (CI) was 2.179–3.391 (P<0.001) in the univariate Cox regression analysis (Figure 4A). In the multivariate Cox regression analysis, the HR of the risk score was 2.367 (95% CI: 1.892–2.963) (P<0.001), demonstrating that the 10-CuRLs signature was an independent prognostic factor (Figure 4B). A nomogram was constructed to assess the 1-, 3-, and 5-year OS incidences (Figure 4C). Calibration plots of the nomogram for 1-, 3-, and 5-years demonstrated the great agreement between nomogram-predicted and observed results (Figure 4D). The area under the curve (AUC) values of the risk score, age, gender, and tumor stage are shown in Figure 4E, indicating the discriminative capability of the risk model with respect to LUAD survival. In addition, the AUC values at 1, 3, and 5 years reached 0.778, 0.768, and 0.731, respectively (Figure 4F).

Figure 4 Independent prognostic analysis of the risk model. Forest plots showing (A) univariate and (B) multivariate Cox regression analysis of associations between clinical characteristics (including the 10-CuRLs signature) and OS. (C) A prognostic nomogram constructed to predict the 1-, 3-, and 5-year survival. (D) Calibration plot of the nomogram between the predicted and observed outcomes. (E) ROC curves of the clinical characteristics and risk score. (F) Time-dependent ROC curves of OS at 1-, 3- and 5-year. ***, P<0.001. AUC, area under the curve; CI, confidence interval; CuRLs, cuproptosis-related lncRNAs; HR, hazard ratio; OS, overall survival; ROC, receiver operating characteristic curve; TPR, true positive rate; FPR, false positive rate.

Validation of the risk model as an independent prognostic factor for LUAD

To validate the risk model based on 10 prognostic CuRLs obtained in the TCGA training cohort, we calculated risk scores for each patient in the GEO meta-cohort. Patients (N=458) were divided into low-risk and high-risk groups according to the median risk score. The distribution of risk score, survival outcome and expression of CuRLs are displayed in Figure 5A-5C. Kaplan-Meier survival analysis also showed that the OS of patients in the high-risk group was worse than that in the low-risk group (Figure 5D). We also included the risk model and clinical signatures in the univariate and multivariate Cox regression analysis (Figure 5E,5F). The univariate Cox regression analysis indicated that the 10-CuRLs signature was an independent prognostic factor (HR: 1.468, 95% CI: 1.331–1.620, P<0.001). The associations between the risk score and OS remained significant after adjusting for potential confounding factors in the multivariate Cox regression analysis (HR: 1.398, 95% CI: 1.258–1.553, P<0.001). The risk score showed better discriminatory ability over other clinical features (Figure 5G). The predictive accuracy was also evaluated by time-dependent ROC analysis, and the AUC value of 1-, 3-, and 5-year reached 0.765, 0.751, and 0.709, respectively (Figure 5H). These results demonstrated that the 10-CuRLs signature had a good-prognostic value in both TCGA training cohort and GEO validation meta-cohort.

Figure 5 Validation of the risk model in the GEO cohort. (A) Distribution of 10-CuRLs risk model-based risk score in the GEO meta-cohort. (B) Survival status and time in the low- and high-risk groups. (C) Hierarchical clustering analysis of 10 CuRLs between the low- and high-risk groups. (D) Kaplan-Meier curve for OS in the low- and high-risk groups. Forest plots showing univariate (E) and multivariate (F) Cox regression analysis of associations between clinical characteristics (including the 10-CuRLs signature) and OS. (G) ROC curves of the clinical characteristics and risk score. (H) Time-dependent ROC curves of OS at 1-, 3- and 5-year. AUC, area under the curve; CI, confidence interval; CuRLs, cuproptosis-related lncRNAs; HR, hazard ratio; OS, overall survival; ROC, receiver operating characteristic curve; TPR, true positive rate; FPR, false positive rate.

Identification of the 10 CuRLs-associated biological pathways

We used different approaches to fully explore biological functions associated with the risk model based on the 10 CuRLs. We applied GSVA to identify biological processes that differed significantly between low- and high-risk groups (Figure 6A). Many cancer-related pathways were differentially enriched between the two risk groups, such as glycolysis, p53 pathway, MTORC1 signaling, and apoptosis. The low-risk group showed enrichment for inflammatory response and angiogenesis. We then performed GSEA analysis for each lncRNA of the 10 CuRLs based on the KEGG pathway database. These CuRLs were enriched in pathways associated with cell death and immunity, as shown in GSEA enrichment plots of APTR and CARD8-AS1 (Figure 6B,6C). Figure S3 delineated the top 5 differentially regulated pathways of DLGAP1-AS1, ELOA-AS1, LINC01137, MIR31HG, TMPO-AS1, and ZSCAN16-AS1 (all pathways and their GSEA enrichment scores were not shown). These results indicated potential links between CuRLs and other types of cell death, such as apoptosis, ferroptosis, and necroptosis. In addition, these lncRNAs tended to be involved in immune-associated pathways, such as programmed death ligand-1 (PD-L1) expression and programmed cell death protein 1 (PD-1) checkpoint pathway in cancer, T cell receptor signaling pathway, and antigen processing and presenting. Protein-coding genes co-expressed with the 10 CuRLs were identified using |Pearson correlation coefficient| >0.5 and P<0.001 as the cut-off (Table S4). GO and KEGG analyses were then performed based on these correlated messenger RNAs (mRNAs). The results showed that the 10 CuRL-related mRNAs were enriched in immune receptor activity, T cell activation, and T cell receptor signaling pathway (Figure 6D,6E). These results suggested that the 10-CuRLs signature was not only relevant to traditional cancer-related pathways but also associated with immune response.

Figure 6 Functional analysis of the risk model based on 10 CuRLs. (A) GSVA analysis between low- and high-risk groups. GSEA analysis based on KEGG pathway database of lncRNA (B) APTR and (C) CARD8-AS1. (D) GO enrichment analysis of CuRLs-related mRNAs. (E) KEGG pathway enrichment analysis of CuRLs-related mRNAs. BP, biological process; CC, cellular component; CuRLs, cuproptosis-related lncRNA; GO, Gene Ontology; GSEA, gene set enrichment analysis; GSVA, gene set variation analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function; PD-1, programmed death-1; PD-L1, programmed cell death-ligand 1; Th 1/2/17, T helper cell 1/2/17.

Evaluation of tumor-immune landscape and TMB based on the risk model

To explore the roles of the risk model based on the 10 CuRLs in TIME, we evaluated the landscape of immune cell infiltration in tumor samples. Infiltration levels of immune cells between low- and high-risk groups were calculated by TIMER, CIBERSORT, QuanTIseq, MCP-counter, xCELL, and EPIC algorithms (Figure 7 and Table S5). In low-risk samples, multiple immune cells were highly expressed, including B cells and CD4+ T cells. According to the xCELL algorithm, we observed more types of infiltrated immune cells in low-risk samples, such as activated myeloid dendritic cells, CD8+ T cells, monocytes, and so on. We found that the risk score was negatively related to the immune and stroma scores based on the xCELL algorithm, suggesting that increased risk may lead to fewer intra-tumoral immune cells. However, a significant difference in stromal activity was not observed based on the ESTIMATE algorithm (Figure 8A,8B). Considering the development of ICIs, we further examined the expression of 79 ICGs between low- and high-risk groups, and 36 of them were found to be differentially expressed (Figure 8C). Notably, there was a significant difference in CTLA4, which is dysregulated in tumors and tumor-associated T cells and is a promising immunotherapeutic target. The higher expression of CTLA4 in low-risk patients indicated that these patients may be more likely to benefit from anti-CTLA4 treatment. The relation between the risk group and CD276 also attracted our attention. CD276 has been recently discovered to enable cancer stem cells to evade immune surveillance (33). The high expression of CD276 in high-risk groups indicted that immune suppression in TIME of high-risk samples could permit tumor cells to escape, leading to drug resistance. The correlation heatmap also showed significant associations between risk score and ICGs (Figure 8D).

Figure 7 Comparisons of infiltrated immune cells between low-risk and high-risk groups based on TIMER, CIBERSORT, QUANTISEQ, MCP-counter, xCELL, and EPIC algorithms. CIBERSORT, cell-type identification by estimating relative subsets of RNA transcripts; EPIC, estimating the proportions of immune and cancer cells; MCP-counter, microenvironment cell populations-counter; TIMER, tumor immune estimation resource.

Figure 8 Differences in the TIME between low- and high-risk groups. (A,B) Comparison of the stromal score and immune score between low- and high-risk samples. (C) The differential expression of ICGs in low- and high-risk groups. (D) Correlation heatmap showing the associations between risk score and ICGs. White squares represent no significant correlation. (E) Comparison of TMB between low- and high-risk groups. (F) Kaplan-Meier analysis of the effect of TMB status on OS. (G) Kaplan-Meier analysis for OS of patients categorized by combing TMB status and risk score. ICG, immune checkpoint gene; OS, overall survival; TIME, tumor immune microenvironment; TMB, tumor mutation burden.

Emerging evidence suggests that TMB is a biomarker of response to ICI treatment (34). Therefore, we explored the differences in TMB between low- and high-risk groups. It was found that patients in the high-risk group exhibited a higher TMB compared with the low-risk group (Figure 8E). Higher TMB was associated with a favorable prognosis in LUAD patients (Figure 8F). Intriguingly, the two results appear to contradict the poor prognosis in high-risk groups. Thus, we divided patients into 4 groups based on TMB and risk score (Figure 8G). Patients with a low TMB and high risk score had the least favorable OS. Correspondingly, low-risk patients with a high TMB had the best prognosis. These results indicated that risk group combined with TMB might provide a more precise assessment of patient prognosis.

Comparison of the sensitivity to anti-cancer drugs between low- and high-risk patients

The different landscape of TIME and TMB may lead to differential susceptibility to drug response. Therefore, we utilized the pRRophetic algorithm to detect suitable targeted drugs for patients in low- and high-risk groups. The sensitivity of 15 drugs that are currently in clinical use or in preclinical trials for lung cancer showed significant differences between low- and high-risk groups (P<0.05, Figure 9A). Among drugs commonly used in the treatment of non-small cell lung cancer clinically, the sensitivity of cisplatin, docetaxel, gemcitabine, gefitinib, and paclitaxel was higher in patients with low risk (Figure 9B-9G), suggesting that these drugs may be more effective in high-risk patients. However, low-risk patients might benefit more from imatinib. Some other drugs in preclinical use also exhibited potential for lung cancer treatment in different risk groups (Figure S4).

Figure 9 Predicted drug responses in low- and high-risk groups. (A) Schematic diagram for anti-cancer drugs currently in clinical use or in preclinical trials for lung cancer (P<0.05). Differences in drug sensitivity of (B) Cisplatin, (C) Docetaxel, (D) Gemcitabine, (E) Gefitinib, (F) Imatinib, (G) Paclitaxel.

Discussion

Despite progress made in the field of surgery, chemotherapy, radiotherapy, and immunotherapy, the current status of LUAD treatment is far from satisfactory. Novel molecular biomarkers, which could provide additional prognostic value and better risk assessment, may help clinicians to develop personalized therapy.

Various types of PCD, such as necroptosis, pyroptosis, and ferroptosis, have shown vital roles in inhibiting cancer progression and clinical deterioration (35-38). Although the introduction of ICI represents a milestone in cancer treatment, only one-third of patients are responsive (39). The activation of cell death has also been discovered to boost the efficacy of anti-tumor immune therapy (3). For example, the acute activation of pyroptosis leads to the release of inflammatory cytokines and immunostimulatory alarmins, thus increasing the infiltration of various immune cells (40). Wang et al. found that pyroptosis could turn ‘cold’ tumors into ‘hot’ tumors and exhibited an anti-tumor effect synergistically with ICI (41).

Copper binding agents acting as copper ionophores have been found to lead to cell death in tumor cells and provided a rationale for the exploitation of copper toxicity as an anti-tumor tool (42). However, the molecular mechanism of the induction of copper and its role in anti-tumor activity is still unclear. A novel type of cell death, cuproptosis mediated by protein lipoylation has been recently elucidated, opening the door for cuproptosis as a novel approach in cancer treatment (4).

The role of lncRNAs in PCD cannot be ignored (15). The overexpression of lncRNA LINC0036 up-regulated cystathionine β synthase (CBS) which is a ferroptosis marker in lung cancer cells (43). It suggests that combining lncRNAs with PCD may lead to a better understanding of the underlying tumor biology. Several recent studies constructed clinical models based on PCD-related lncRNAs to predict the prognosis and immune response of patients (17,18,44,45). Therefore, it is quite worthy to explore the prognostic value of CuRLs in LUAD, which is a new hotspot in the field of cell death. In this study, using Cox regression analysis and a machine-learning algorithm based on LASSO analysis, we developed a novel 10-CuRLs signature to predict patient prognosis combined with conventional clinical risk factors. The CuRLs signature showed great prediction accuracy, which was further validated in a GEO meta-cohort.

A number of studies have focused on the crosstalk between cell death and TIME (3,46,47). In our research, we also fully explored the biological functions associated with the risk model based on the 10 CuRLs. It was found that the 10 CuRLs-related mRNAs were enriched in immune receptor activity, T cell activation, and T cell receptor signaling pathway. These results suggested that the 10-CuRLs signature was not only relevant to traditional cancer-related pathways but also associated with immune response. Therefore, we further evaluated the tumor-immune landscape based on our CuRLs risk mode and sequentially observed significant differences in the immune landscape between low-risk and high-risk groups. This may guide clinicians in making individualized treatment decisions between different risk groups.

Multiple differentially expressed ICGs between low- and high-risk groups bring about new prospects for the administration of ICIs for LUAD patients with different risk scores. In addition to CTLA4 and CD276 already mentioned above, some other ICGs are worthy of attention. Butyrophilins (BTNs), which are structurally similar to the immunosuppressive B7 family, are of great importance in restarting the anti-tumor efficacy of γδ T cells (48,49). The higher expression of BTN2A1, BTN2A2, and BTNL9 in the low-risk group indicated that these patients may have a better response to checkpoint immunotherapy. In low-risk patients, we also observed higher expression of genes in HLA class II which could present antigens to CD4+ T cells (50,51). These HLA-II genes hold promise as a predictive biomarker of response to anti-PD-1 or anti-PD-L1 therapy (52).

TMB, defined as the number of somatic mutations per megabase, is emerging as a potential biomarker for the response to immunotherapy (53). However, there is currently no consensus over the optimal cut-off of TMB for patient stratification, leading to limited further application in the clinic (54). In our present study, the combination of TMB and the CuRLs risk model brought about a more accurate analysis of patient survival, indicating that TMB combined with the CuRLs risk model may serve as a promising prognostic factor in addition to predicting the immune response. Our study further showed the ability of the CuRLs risk model to predict the sensitivity to anti-cancer drugs, which may provide a novel method for clinical medication-related decisions.

Although our findings have been validated in an independent meta-cohort, our study had some limitations. On the one hand, our study was retrospective based on publicly available TCGA and GEO databases and these results need to be verified in prospective studies for clinical use. On the other hand, further biological experiments are needed to confirm the expression of these key genes and elucidate their biological functions in lung cancer.

Conclusions

To our knowledge, the present study is the first to construct a CuRLs prognostic signature in patients with LUAD and validate it in an external meta-cohort. The risk model exhibited a high diagnostic accuracy in predicting patient survival outcomes combined with clinical features. The differences in characteristics between low- and high-risk groups, including immune infiltration, TMB, and drug sensitivity, open the way to better patient stratification and explore novel drugs in different risk groups. We are looking forward to furthering studies into the mechanisms behind how these CuRLs regulate cuproptosis and affect patient outcomes. We hope the CuRLs risk model could be validated in more clinical cohorts and further benefit patients in the future.

Acknowledgments

We thank Dr. Jianming Zeng (University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes. We also would like to acknowledge the work by TCGA, GSE31210, GSE37745, and GSE50081 for providing LUAD data.

Funding: This work was supported by the National Natural Science Foundation of China (No. 82072557 to HL, No. 81871882 to HL, No. 82001976 to HJ); National Key Research and Development Program of China (No. 2021YFC2500900 to HL); Shanghai Municipal Education Commission- Gaofeng Clinical Medicine Grant (No. 20172005, the 2nd round of disbursement to HL); Program of Shanghai Academic Research Leader from Science and Technology Commission of Shanghai Municipality (No. 20XD1402300 to HL); and Shanghai Sailing Program (No. 20YF1440000 to HJ).

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-22-500/rc

Peer Review File: Available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-22-500/prf

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-22-500/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Siegel RL, Miller KD, Fuchs HE, et al. Cancer statistics, 2022. CA Cancer J Clin 2022;72:7-33. [Crossref] [PubMed]
2. Snyder AG, Oberst A. The Antisocial Network: Cross Talk Between Cell Death Programs in Host Defense. Annu Rev Immunol 2021;39:77-101. [Crossref] [PubMed]
3. Tang R, Xu J, Zhang B, et al. Ferroptosis, necroptosis, and pyroptosis in anticancer immunity. J Hematol Oncol 2020;13:110. [Crossref] [PubMed]
4. Tsvetkov P, Coy S, Petrova B, et al. Copper induces cell death by targeting lipoylated TCA cycle proteins. Science 2022;375:1254-61. [Crossref] [PubMed]
5. Kim BE, Nevitt T, Thiele DJ. Mechanisms for copper acquisition, distribution and regulation. Nat Chem Biol 2008;4:176-85. [Crossref] [PubMed]
6. Williams DM. Copper deficiency in humans. Semin Hematol 1983;20:118-28. [PubMed]
7. Oliveri V. Selective Targeting of Cancer Cells by Copper Ionophores: An Overview. Front Mol Biosci 2022;9:841814. [Crossref] [PubMed]
8. Margalioth EJ, Schenker JG, Chevion M. Copper and zinc levels in normal and malignant tissues. Cancer 1983;52:868-72. [Crossref] [PubMed]
9. Yu Z, Zhou R, Zhao Y, et al. Blockage of SLC31A1-dependent copper absorption increases pancreatic cancer cell autophagy to resist cell death. Cell Prolif 2019;52:e12568. [Crossref] [PubMed]
10. Yao W, Qian X, Ochsenreither S, et al. Disulfiram Acts as a Potent Radio-Chemo Sensitizer in Head and Neck Squamous Cell Carcinoma Cell Lines and Transplanted Xenografts. Cells 2021;10:517. [Crossref] [PubMed]
11. Volders PJ, Anckaert J, Verheggen K, et al. LNCipedia 5: towards a reference set of human long non-coding RNAs. Nucleic Acids Res 2019;47:D135-9. [Crossref] [PubMed]
12. Chen J, Zhang K, Zhi Y, et al. Tumor-derived exosomal miR-19b-3p facilitates M2 macrophage polarization and exosomal LINC00273 secretion to promote lung adenocarcinoma metastasis via Hippo pathway. Clin Transl Med 2021;11:e478. [Crossref] [PubMed]
13. Qian X, Yang J, Qiu Q, et al. LCAT3, a novel m6A-regulated long non-coding RNA, plays an oncogenic role in lung cancer via binding with FUBP1 to activate c-MYC. J Hematol Oncol 2021;14:112. [Crossref] [PubMed]
14. Loewen G, Jayawickramarajah J, Zhuo Y, et al. Functions of lncRNA HOTAIR in lung cancer. J Hematol Oncol 2014;7:90. [Crossref] [PubMed]
15. Jiang N, Zhang X, Gu X, et al. Progress in understanding the role of lncRNA in programmed cell death. Cell Death Discov 2021;7:30. [Crossref] [PubMed]
16. Tanzhu G, Li N, Li Z, et al. Molecular Subtypes and Prognostic Signature of Pyroptosis-Related lncRNAs in Glioma Patients. Front Oncol 2022;12:779168. [Crossref] [PubMed]
17. Guo Y, Qu Z, Li D, et al. Identification of a prognostic ferroptosis-related lncRNA signature in the tumor microenvironment of lung adenocarcinoma. Cell Death Discov 2021;7:190. [Crossref] [PubMed]
18. Tang R, Wu Z, Rong Z, et al. Ferroptosis-related lncRNA pairs to predict the clinical outcome and molecular characteristics of pancreatic ductal adenocarcinoma. Brief Bioinform 2022;23:bbab388. [Crossref] [PubMed]
19. Zeidan AM, Sekeres MA, Garcia-Manero G, et al. Comparison of risk stratification tools in predicting outcomes of patients with higher-risk myelodysplastic syndromes treated with azanucleosides. Leukemia 2016;30:649-57. [Crossref] [PubMed]
20. Gu Z, Eils R, Schlesner M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics 2016;32:2847-9. [Crossref] [PubMed]
21. Gu Z, Gu L, Eils R, et al. circlize Implements and enhances circular visualization in R. Bioinformatics 2014;30:2811-2. [Crossref] [PubMed]
22. Wu T, Hu E, Xu S, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb) 2021;2:100141. [Crossref] [PubMed]
23. Li T, Fu J, Zeng Z, et al. TIMER2.0 for analysis of tumor-infiltrating immune cells. Nucleic Acids Res 2020;48:W509-14. [Crossref] [PubMed]
24. Newman AM, Liu CL, Green MR, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods 2015;12:453-7. [Crossref] [PubMed]
25. Aran D, Hu Z, Butte AJ. xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol 2017;18:220. [Crossref] [PubMed]
26. Finotello F, Mayer C, Plattner C, et al. Molecular and pharmacological modulators of the tumor immune contexture revealed by deconvolution of RNA-seq data. Genome Med 2019;11:34. [Crossref] [PubMed]
27. Becht E, Giraldo NA, Lacroix L, et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol 2016;17:218. [Crossref] [PubMed]
28. Racle J, Gfeller D. EPIC: A Tool to Estimate the Proportions of Different Cell Types from Bulk Gene Expression Data. Methods Mol Biol 2020;2120:233-48. [Crossref] [PubMed]
29. Yoshihara K, Shahmoradgoli M, Martínez E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 2013;4:2612. [Crossref] [PubMed]
30. Hu FF, Liu CJ, Liu LL, et al. Expression profile of immune checkpoint genes and their roles in predicting immunotherapy response. Brief Bioinform 2021;22:bbaa176. [Crossref] [PubMed]
31. Mayakonda A, Lin DC, Assenov Y, et al. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res 2018;28:1747-56. [Crossref] [PubMed]
32. Geeleher P, Cox N, Huang RS. pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One 2014;9:e107468. [Crossref] [PubMed]
33. Wang C, Li Y, Jia L, et al. CD276 expression enables squamous cell carcinoma stem cells to evade immune surveillance. Cell Stem Cell 2021;28:1597-1613.e7. [Crossref] [PubMed]
34. Offin M, Rizvi H, Tenet M, et al. Tumor Mutation Burden and Efficacy of EGFR-Tyrosine Kinase Inhibitors in Patients with EGFR-Mutant Lung Cancers. Clin Cancer Res 2019;25:1063-9. [Crossref] [PubMed]
35. Friedmann Angeli JP, Krysko DV, Conrad M. Ferroptosis at the crossroads of cancer-acquired drug resistance and immune evasion. Nat Rev Cancer 2019;19:405-14. [Crossref] [PubMed]
36. Jiang X, Stockwell BR, Conrad M. Ferroptosis: mechanisms, biology and role in disease. Nat Rev Mol Cell Biol 2021;22:266-82. [Crossref] [PubMed]
37. Gong Y, Fan Z, Luo G, et al. The role of necroptosis in cancer biology and therapy. Mol Cancer 2019;18:100. [Crossref] [PubMed]
38. Xia X, Wang X, Cheng Z, et al. The role of pyroptosis in cancer: pro-cancer or pro-"host"? Cell Death Dis 2019;10:650. [Crossref] [PubMed]
39. Jiang P, Gu S, Pan D, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med 2018;24:1550-8. [Crossref] [PubMed]
40. Du T, Gao J, Li P, et al. Pyroptosis, metabolism, and tumor immune microenvironment. Clin Transl Med 2021;11:e492. [Crossref] [PubMed]
41. Wang Q, Wang Y, Ding J, et al. A bioorthogonal system reveals antitumour immune function of pyroptosis. Nature 2020;579:421-6. [Crossref] [PubMed]
42. Tardito S, Bassanetti I, Bignardi C, et al. Copper binding agents acting as copper ionophores lead to caspase inhibition and paraptotic cell death in human cancer cells. J Am Chem Soc 2011;133:6235-42. [Crossref] [PubMed]
43. Wang M, Mao C, Ouyang L, et al. Long noncoding RNA LINC00336 inhibits ferroptosis in lung cancer by functioning as a competing endogenous RNA. Cell Death Differ 2019;26:2329-43. [Crossref] [PubMed]
44. Ye Y, Dai Q, Qi H. A novel defined pyroptosis-related gene signature for predicting the prognosis of ovarian cancer. Cell Death Discov 2021;7:71. [Crossref] [PubMed]
45. Yao J, Chen X, Liu X, et al. Characterization of a ferroptosis and iron-metabolism related lncRNA signature in lung adenocarcinoma. Cancer Cell Int 2021;21:340. [Crossref] [PubMed]
46. Loveless R, Bloomquist R, Teng Y. Pyroptosis at the forefront of anticancer immunity. J Exp Clin Cancer Res 2021;40:264. [Crossref] [PubMed]
47. Zeng Q, Ma X, Song Y, et al. Targeting regulated cell death in tumor nanomedicines. Theranostics 2022;12:817-41. [Crossref] [PubMed]
48. Sarter K, Leimgruber E, Gobet F, et al. Btn2a2, a T cell immunomodulatory molecule coregulated with MHC class II genes. J Exp Med 2016;213:177-87. [Crossref] [PubMed]
49. Vavassori S, Kumar A, Wan GS, et al. Butyrophilin 3A1 binds phosphorylated antigens and stimulates human γδ T cells. Nat Immunol 2013;14:908-16. [Crossref] [PubMed]
50. Neefjes J, Ovaa H. A peptide’s perspective on antigen presentation to the immune system. Nat Chem Biol 2013;9:769-75. [Crossref] [PubMed]
51. Alspach E, Lussier DM, Miceli AP, et al. MHC-II neoantigens shape tumour immunity and response to immunotherapy. Nature 2019;574:696-701. [Crossref] [PubMed]
52. Axelrod ML, Cook RS, Johnson DB, et al. Biological Consequences of MHC-II Expression by Tumor Cells in Cancer. Clin Cancer Res 2019;25:2392-402. [Crossref] [PubMed]
53. Sha D, Jin Z, Budczies J, et al. Tumor Mutational Burden as a Predictive Biomarker in Solid Tumors. Cancer Discov 2020;10:1808-25. [Crossref] [PubMed]
54. Jardim DL, Goodman A, de Melo Gagliato D, et al. The Challenges of Tumor Mutational Burden as an Immunotherapy Biomarker. Cancer Cell 2021;39:154-73. [Crossref] [PubMed]