- Department of Thoracic Surgery, The First Affiliated Hospital of Zhengzhou University, Zhengzhou, China
Background: There is a rapid increase in lung adenocarcinomas (LUAD), and studies suggest associations between cuproptosis and the occurrence of various types of tumors. However, it remains unclear whether cuproptosis plays a role in LUAD prognosis.
Methods: Dataset of the TCGA-LUAD was treated as training cohort, while validation cohort consisted of the merged datasets of the GSE29013, GSE30219, GSE31210, GSE37745, and GSE50081. Ten studied cuproptosis-related genes (CRG) were used to generated CRG clusters and CRG cluster-related differential expressed gene (CRG-DEG) clusters. The differently expressed lncRNA that with prognosis ability between the CRG-DEG clusters were put into a LASSO regression for cuproptosis-related lncRNA signature (CRLncSig). Kaplan–Meier estimator, Cox model, receiver operating characteristic (ROC), time-dependent AUC (tAUC), principal component analysis (PCA), and nomogram predictor were further deployed to confirm the model’s accuracy. We examined the model’s connections with other forms of regulated cell death, including apoptosis, necroptosis, pyroptosis, and ferroptosis. The immunotherapy ability of the signature was demonstrated by applying eight mainstream immunoinformatic algorithms, TMB, TIDE, and immune checkpoints. We evaluated the potential drugs for high risk CRLncSig LUADs. Real-time PCR in human LUAD tissues were performed to verify the CRLncSig expression pattern, and the signature’s pan-cancer’s ability was also assessed.
Results: A nine-lncRNA signature, CRLncSig, was built and demonstrated owning prognostic power by applied to the validation cohort. Each of the signature genes was confirmed differentially expressed in the real world by real-time PCR. The CRLncSig correlated with 2,469/3,681 (67.07%) apoptosis-related genes, 13/20 (65.00%) necroptosis-related genes, 35/50 (70.00%) pyroptosis-related genes, and 238/380 (62.63%) ferroptosis-related genes. Immunotherapy analysis suggested that CRLncSig correlated with immune status, and checkpoints, KIR2DL3, IL10, IL2, CD40LG, SELP, BTLA, and CD28, were linked closely to our signature and were potentially suitable for LUAD immunotherapy targets. For those high-risk patients, we found three agents, gemcitabine, daunorubicin, and nobiletin. Finally, we found some of the CRLncSig lncRNAs potentially play a vital role in some types of cancer and need more attention in further studies.
Conclusion: The results of this study suggest our cuproptosis-related CRLncSig can help to determine the outcome of LUAD and the effectiveness of immunotherapy, as well as help to better select targets and therapeutic agents.
Introduction
In China and throughout the world, lung cancer is one of the most common malignant tumors (Chen et al., 2019; Sung et al., 2021). According to WHO statistics, in 2020, 2,206,771 new lung cancer patients were diagnosed globally, accounting for 11.4% of new cancer cases; deaths were ranked the first among all cancer types, accounting for 18% of all cancer deaths (Chen et al., 2019; Sung et al., 2021). Adenocarcinoma (LUAD) accounts for approximately 40% of lung cancer cases (Siegel et al., 2021). Current treatment methods, such as surgical, radiotherapy, and chemotherapy, can hardly meet advanced LUAD patients’ survival expectations (Siegel et al., 2021; Li et al., 2022a). Therefore, the discovery of more effective prognostic models is crucial for exploring the cellular and molecular mechanism of LUAD carcinogenesis and finding treatment strategies.
Cell suicide pathways, termed regulated cell death, play a critical role in organismal development, homeostasis, and pathogenesis (Strasser and Vaux, 2020). Regulated cell death has been observed in cancer for a long time, but the initial reports were mainly in cases where necrosis was observed in hypoxic areas of growing tumors (Strasser and Vaux, 2020). During further tumor progression, cancer cells often respond to their altered state through programmed regulated cell death and are highly dependent on certain survival signals from their environment (Strasser and Vaux, 2020). Regulated cell death has long been implicated in cancer treatment, with radiation and chemotherapy killing cancer cells while destroying normal, healthy cells (Strasser and Vaux, 2020). The biology and therapeutic response of LUAD are reported shaped by various forms of regulated cell death, such as apoptosis (Wu et al., 2019), necroptosis (Lu et al., 2022), pyroptosis (Lin et al., 2021), and ferroptosis (Ma et al., 2021; Zhang et al., 2021). Excitingly, Tsvetkov and colleagues published their latest study in the journal Science, confirming the existence of copper-induced regulated cell death, which is termed cuproptosis (Tsvetkov et al., 2022). In their study, it was shown that cuproptosis is distinct from apoptosis, necroptosis, pyroptosis, and ferroptosis (Tsvetkov et al., 2022). Recently, the teams of Zhang et al. (2022a); Zhang et al. (2022b) separately developed cuproptosis-related prognostic models for hepatocellular carcinoma. Li et al. (2022b) also conducted in-depth research on cuproptosis and constructed a prediction model for oral squamous cell carcinoma. Bian et al. (2022); Ji et al. (2022) have established cuproptosis models of clear cell renal cell carcinoma, respectively, both of which claim to be effective in predicting the prognosis of the disease.
Given that some forms of regulated cell death may be more immune-targeted than others, understanding how cuproptosis is initiated, propagated, and ultimately executed in LUAD may have important implications for possible combined diagnostic and therapeutic interventions. Regarding the prognostic model about cuproptosis, no LUAD research so far has been published. Most studies of prognostic signatures target entire tumor populations without tailoring high-risk patients, making them insufficient for risk stratification or treatment (Ma et al., 2018; Ma et al., 2020a; Ma et al., 2020b; Li et al., 2020; Zheng et al., 2020; Ma et al., 2021; Zhang et al., 2021; Ma et al., 2022a; Ma et al., 2022b). To cope with the above issues, we first developed a cuproptosis-related lncRNA prognostic signature of LUAD and potential therapeutic targets and drugs for high-risk patients. We gathered proved cuproptosis-related genes to construct a lncRNA signature having the power to predict LUAD outcomes, and further validated its prognostic ability in a large independent cohort. We also tested the signature genes’ expression profile in real world and the model’s connections with other forms of regulated cell death. More importantly, we revealed that KIR2DL3, IL10, IL2, CD40LG, SELP, BTLA, and CD28, were linked closely to our signature and were potentially suitable for LUAD immunotherapy targets. We found three agents, gemcitabine, daunorubicin, and nobiletin, for those high-risk patients. Finally, we tested the ability of our signature lncRNAs in pan-cancer.
Materials and methods
Data selection and preprocessing
We downloaded the RNA-seq data and clinical phenotype of patients in the TCGA-LUAD project from the Xena Hub (https://xenabrowser.net/). Use the keyword “lung adenocarcinoma” to search in Gene Expression Omnibus (Clough and Barrett, 2016) (GEO, https://www.ncbi.nlm.nih.gov/geo/), and filter out the datasets containing total RNA and the total number of lung cancer patients with survival time greater than 80 in the results as candidate validation cohorts. GSE29013, GSE30219, GSE31210, GSE37745, and GSE50081 were selected, and their data were downloaded. The TCGA-LUAD were gathered as training cohort. For the preprocessing of GSE29013, GSE30219, GSE31210, GSE37745, and GSE50081, we used the R package “inSilicoMerging” to merge them (Taminau et al., 2012), and then we adopted the method that developed by Johnson et al. (2007). To remove the batch effect and finally obtained the data matrix treated as a validation cohort.
Identification of cuproptosis-regulated genes (CRG) subgroups using consensus clustering
In the present research, we retrospectively selected 10 CRG (FDX1, LIAS, LIPT1, DLD, DLAT, PDHA1, PDHB, MTF1, GLS and CDKN2A) from the work of Tsvetkov et al. (2022). We used the “limma” R package to determine the co-expression genes of the ten CRG with the parameters of correlation coefficient >0.4 or < −0.4 and p < 0.05. The “ConsensusClusterPlus” R package was used to cluster CRG co-expression genes in LUAD patients extracted from the training cohort into different subtypes. The counts of the subtypes were determined by the optimal k value. The Kaplan-Meier estimator (KM) curve was constructed using R packages “survival” and “survminer” to measure difference in overall survival of patients in the CRG clusters. The “scatterplot3d” R package was applied to perform principal component analysis (PCA) to see the clusters differences. Packages “reshape2,” “ggpubr,” “limma,” “GSEABase,” and “GSVA” were used to conduct the single-sample gene set enrichment analysis (ssGSEA) and visualization. We then screened the CRG cluster-related differential expressed genes (CRG-DEGs) between clusters using the R package “limma” with Log FC > 0.15 or < −0.15 and adjust p < 0.05. These CRG-DEGs were the put into an KEGG analysis for the potential pathways finding.
Constructions of CRG-DEG clusters and a cuproptosis-regulated lncRNA signature (CRLncSig)
With the help of the “ConsensusClusterPlus” R package and the CRG-DEGs, we were able to group the patients in the training cohort into different CRG-DEG clusters. The survival difference of CRG-DEG clusters was measured using the KM curve calculated using the R packages “survival” and “survminer”. The “scatterplot3d” R package was applied to perform PCA to see the clusters differences. Packages “reshape2,” “ggpubr,” “limma,” “GSEABase,” and “GSVA” were used to conduct the ssGSEA analysis and visualization. We then conducted GSVA to screen the most important KEGG pathways by comparing CRG-DEG clusters by R packages “limma,” “GSEABase,” “GSVA,” and “pheatmap”. We investigated the differently expressed lncRNA (DEL) between the CRG-DEG clusters with Log FC > 0.15 or < −0.15 and adjust p < 0.05. Then, these DELs were subjected to univariate Cox and KM analyses for choosing the ones with potential prognostic power with p < 0.05. In order to prevent overfitting, we performed the least absolute shrinkage and selection operator (LASSO) using the R package “glmnet” on the prognostic DELs. A 10-fold cross-validation was performed in the training cohort at 1 SE above the minimum partial likelihood deviation to estimate the penalty parameter (Tibshirani, 1997; Sauerbrei et al., 2007; Friedman et al., 2010; Goeman, 2010). For each patient, a risk score was calculated according to the following formula:
with βi denoting the coefficient, Expi denoting the relative expression level of each lncRNA normalized by z-score, and n denoting each lncRNA in the CRLncSig.
Validation of the CRLncSig in an indenpendent cohort
High-risk and low-risk patients were divided using median risk scores. In both the training and validation cohorts, the CRLncSig’s predictive capability was evaluated using methods such as the KM curve, univariate and multivariable Cox analysis (Cao and Lopez-de-Ullibarri, 2019), receiver operating characteristic (ROC) curve, time-dependent AUC (tAUC) analysis, and PCA. ROC curve and tAUC analyzes were implemented with the help of the “timeROC” and “survival” R packages. The “scatterplot3D” R package was used to evaluate the distribution of patients with different risk scores by PCA. Additionally, we used the R packages “survival,” “survminer,” “rms,” and “regplot,” to construct nomograms that predicted 1-, 3-, and 5-year overall survival and calibrated the model to determine whether the model’s predictions matched up with the actual consistency.
Correlations between the cuproptosis-related lncRNA signature and apoptosis, necroptosis, pyroptosis, and ferroptosis
For better knowing the interactions between our CRLncSig and other types of “cell death”, we adopted the Pearson analysis. Apoptosis, necroptosis, and pyroptosis-related genes were extracted from the GeneCard and Gene Set Enrichment Analysis (GSEA) online databases, respectively, by applying the following steps: 1) search the GeneCard using the corresponding keyword; 2) search the GSEA using the corresponding keyword: 3) merge the above results and take the unique genes. FerrDb is the first database dedicated to ferroptosis regulators and ferroptosis-related diseases (Zhou and Bao, 2020). Ferroptosis-related genes were obtained from FerrDb (http://www.datjar.com:40013/bt2104/).
GSEA
We first downloaded the GSEA program from http://www.gsea-msigdb.org/gsea/downloads.jsp, the GSEA website. By dividing LUADs by their median risk scores, we classified them as low-risk and high-risk. From the Molecular Signatures Database (Liberzon et al., 2011) (http://www.gsea-msigdb.org/gsea/downloads.jsp), we downloaded “c2.cp.kegg.v7.4.symbols.gmt” to evaluate the assess of KEGG pathways between high-risk and low-risk groups using GSEA. GSEA parameter settings were set to five minimum genes, 5,000 maximum genes, and 1,000 resamplings. A statistically significant value was a p-value of 0.05 and a false discovery rate of 0.25.
Identification of the immunological status of the CRLncSig
We calculated stromal, immune, and ESTIMATE scores for each patient based on gene expression in the training cohort using the R package “ESTIMATE” (Yoshihara et al., 2013). We calculated the association between CRLncSig and stromal, immune, and ESTIMATE scores using the Pearson coefficient and Wilcoxon rank sum. Multi-omics data can be leveraged with “IOBR” to facilitate immuno-oncology exploration, revealing tumor-immune interactions and accelerating precision immunotherapy. To calculate the immune-infiltrating cell scores for each sample of the training cohort, we used the R package “IOBR” and its methods of CIBERSORT, CIBERSORT-ABS, quantIseq, TIMER, IPS, MCPCounter, xCell, and EPIC. CRLncSig’s relationships with immune cell types of the eight algorithms were calculated using Pearson coefficients and Wilcoxon rank-sum, and lollipop plots and heatmap were applied for the visualization. Venn and cloud diagrams were used to summarize the results we obtained. As part of the “gsva” R package, the “ ssGSEA” function was used to assess 13 immune-related pathways of the CRLncSig.
Identification of the immunotherapy role and immune checkpoint target of CRLncSig
Using the “maftools,” the top 20 mutated genes were identified, and the mutations and their frequencies were visualized across all training cohort samples. Based on the median risk score of LUADs, we divided them into two groups. To compare gene mutation frequencies between low- and high-risk LUAD populations, we used the chi-square test. Tumor mutational burden (TMB) is an emerging therapeutic measure of immunotherapy sensitivity. TMB is defined as the frequency of certain mutations within a tumor’s genes (Chalmers et al., 2017). The TMB rank score of each case with LUAD was obtained as previously described (Chalmers et al., 2017). We deployed the Pearson coefficient together with the Wilcoxon rank-sum to calculate the connections between the risk score and TMB. Tumor Immune Dysfunction and Exclusion (TIDE) is a computational framework that integrates T cell dysfunction expression signatures and T cell exclusion to model tumor immune evasion. Using TIDE, tumor immune evasion can be modeled in two ways and can be used to predict immunotherapy outcomes (Jiang et al., 2018; Fu et al., 2020; Chen et al., 2021a). More importantly, we tested if our signature correlated with the TIDE. A total of 60 immune checkpoints were selected from previous studies (Thorsson et al., 2018) (Supplementary Table S1), including 24 inhibitory and 36 stimulatory checkpoints. We measured the correlation between our signature and these 60 immune checkpoints based on Pearson and Wilcoxon rank-sum analysis. To further test whether our CRLncSig could guide immunotherapy, we deployed the Kaplan–Meier estimator and Cox regression to test the prognostic ability of the 60 immune checkpoints. Venn diagram was used to summarize the above results to find good checkpoints that may have the targeting ability of the CRLncSig. We collected data from multiple immune datasets to further evaluate the potential impact of these promising checkpoint genes on the immune system and immunotherapy. We presented them in a heatmap format with the help of the regulator prioritization module of the TIDE online tool (Fu et al., 2020).
Identification of drugs for high risk score LUADs
High-risk and low-risk patients were divided using median risk scores. Data on expression profiles and somatic mutations in human cancer cell lines (CCLs) were acquired from the Broad Institute Cancer Cell Line Encyclopedia (CCLE) (https://portals.broadinstitute.org/ccle/) (Ghandi et al., 2019). The dependency map (DepMap, https://de.pmap.org/portal/) portal was used to collect the CERES scores of genome-scale CRISPR knockout screens on 18,333 genes in 739 cell lines. CERES scores indicate how dependent a gene is on certain CCLs. The lower the CERES score, the greater the likelihood that the gene plays an important role in the growth and survival of a given CCL. Based on PRISM Repurposing (https://depmap.org/portal/prism/) and the Cancer Therapeutics Response Portal (https://portals.broadinstitute.org/ctrp), drug sensitivity data of CCLs were gathered. In the CTRP, 481 compounds have been evaluated over 835 CCLs, while in the PRISM, 1,448 compounds have been evaluated over 482 CCLs. In both datasets, the area under the dose-response curve (area under the curve - AUC) represents drug sensitivity, and lower values indicate greater sensitivity. To identify drug candidates with higher drug sensitivity in patients with high risk scores, CTRP and PRISM derived drug response data were analyzed. To identify compounds with lower estimated AUC values in the high risk score group, a differential drug response analysis (log2FC > 0.9) was conducted between the groups of the top decile score and bottom decile score (Yang et al., 2021). As a next step, Spearman correlation analysis between AUC value and risk score was conducted to select compounds with a negative correlation coefficient (Spearman’s r < −0.09) (Yang et al., 2021).
Validation of drugs using connectivity map (CMap) analysis
These above results were then put into multiple perspective analyses, including clinical trial data, published experiment evidence, and CMap for further confirming their effective in LUAD (Yang et al., 2021). CMap analysis was performed as a complement to investigate the therapeutic potential of candidate agents in LUAD. We first conducted differential expression analysis between tumor and normal samples. Next, 300 genes with the most significant fold changes (150 upregulated genes and 150 downregulated genes) were submitted to the CMap website (https://clue.io/query). Gene expression signature of perturbations in this website are derived from both CMap v2 and Library of Integrated Network-Based Cellular Signatures (LINCS) database, and a total of 2,429 compounds are available for CMap analysis. CMap calculates a connectivity score for each perturbation, which ranges from −100 to 100. Specific perturbations with negative scores have opposite gene expression patterns to a particular disease, suggesting they have therapeutic potential.
Validation of the CRLncSig’s expression profile using the real-time PCR and human LUAD tissues
In order to verify the expression levels of each lncRNA of CRLncSig, nine pairs of LUAD tissues and adjacent normal tissues were examined using real-time PCR. The Ethics Review Committee of the First Affiliated Hospital of Zhengzhou University approved this study. Informed consent was obtained from all patients before surgery, no history of chemotherapy or radiotherapy was present before surgery. After extracting the tissues, we froze and stored them in liquid nitrogen. Total RNA was extracted from the sample tissues via Trizol lysate (Thermo Fisher Scientific). We reverse-transcribed total RNA from clinical samples into cDNA using the HiScript III RT SuperMix Kit (R32301, Vazyme, Nanjing, China). The real-time PCR was performed using the CFX96 system (BIO-RAD Laboratories, Inc., Hercules, CA, United States). Based on the 2−ΔΔCT method, expression levels of target RNAs were normalized to GAPDH. The mean value was used as the final experimental result for replicated wells. The manufacturer’s instructions were followed during all procedures. Student’s t-test was applied to identify genes differentially expressed between normal and tumor samples. Differentially expressed genes were defined as adjusted p-value < 0.05.
Pan-cancer ability determination of each lncRNA of the CRLncSig
First, we downloaded the TCGA TARGET GTEx dataset from the UCSC database (https://xenabrowser.net/) to determine whether CRLncSig lncRNAs are differentially expressed in tumors and normal tissues. In order to complete this step, we also excluded cancer types with fewer than three samples, giving us data on 34 cancer types. A Wilcoxon rank sum test was used in R software to determine the expression difference between normal and tumor samples in each tumor type.
For prognostic assessment, in addition to the TCGA TARGET GTEx dataset we got, we also obtained a high-quality TCGA prognostic dataset from the TCGA prognostic study previously published on Cell (Liu et al., 2018). A supplement of TARGET follow-up data was obtained from UCSC, and samples with a short than 30 days follow-up period were eliminated. Additionally, we eliminated cancer types with fewer than ten samples each and obtained data on 44 cancer types, including expression and overall survival. Based on the Cox proportional hazards regression model built with the R package “survival,” we analyzed the relationship between the lncRNAs and the prognosis of each cancer type.
In the final step, we evaluated the staging ability of each lncRNA. The TCGA TARGET GTEx dataset was obtained above, and the cancer types with fewer than three samples were eliminated, as were cancer types without tumor staging data. Finally, we obtained information on 30 cancer types. For each cancer type’s different clinical stages, we calculated expression distributions of CRLncSig lncRNAs using R software and utilized variance analysis to investigate the differences.
Results
Patient characteristics
As Figure 1B demonstrates, 500 LUAD patients of project TCGA-LUAD were treated as a training cohort for model construction. 554 LUADs from the GSE29013, GSE30219, GSE31210, GSE37745, and GSE50081 datasets were taken and merged as a validation cohort for model certifying. The effect of data merging for the validation cohort is shown in Figure 2. According to the UMAP plot (Figure 2A), each dataset’s samples are separated before the batch effect is removed. A clustering and intertwining of each data set were observed after Johnson et al. (2007) method had been applied, indicating that the batch effect is well removed by this method. The clinical parameters of the enrolled patients in each cohort are shown in Table 1.
FIGURE 1. Schematic illustration indicates the mechanism of cuproptosis induction (A) (Hu et al., 2022; Tsvetkov et al., 2022) and research design (B). DLAT, dihydrolipoamide S-acetyltransferase; FDX1, ferredoxin-1; Fe–S, iron–sulfur; TCA, tricarboxylic acid; GSH, Glutathione; BSO, buthionine sulfoximine; CRLncSig, cuproptosis-regulated lncRNA signature; TCGA, The Cancer Genome Atlas; LUAD, lung adenocarcinoma; HR, hazard ratio; ROC, receiver operating characteristic; TP, true positive rate; FP, false positive rate; AUC, area under the ROC curve; PC, principal component; KEGG, Kyoto Encyclopedia of Genes and Genomes; GSEA, Gene Set Enrichment Analysis.
FIGURE 2. Data preprocessing and the CRG clusters construction. (A) UMAP plot shows validation cohort merged before batch effect removing, displaying that the samples of each dataset are separated from each other. UMAP plot shows validation cohort merged after batch effect removing, displaying that samples of each dataset are clustered and intertwined with each other. (B) The KM survival curve demonstrated that the overall survival of patients in the CRG clusters were significantly different. The cluster C had a more favor prognosis better than that of patients in cluster B and C, while the cluster B population had a promising outcome than that in the cluster C. (C) Principal component analysis scatter plot shows three clusters were obviously separated. (D) Infiltration distribution of immune cells in the three CRG clusters. (E) Heatmap shows the relationship between CRG clusters, clinical parameters, and 1175 CRG co-expressed genes. Each row represents a gene, and each column represents a sample. (F) Venn diagram shows the process of obtaining the 344 CRG-DEGs. (G) KEGG analysis shows the top 10 enriched pathways. CRG: cuproptosis-regulated gene; UMAP: Uniform Manifold Approximation and Projection; KM: Kaplan–Meier estimator; DEGs: differentially expressed genes; KEGG: Kyoto Encyclopedia of Genes and Genomes; p-value < 0.05 was considered statistically significant; *: p-value <0.05; ***: p-value < 0.001.
Evaluation of the CRG clusters in LUADs using consensus clustering
We first explored the co-expressed genes of 10 CRGs according to our set screening criteria. After removing duplicate results, we finally identified 1,175 genes as CRG co-expressed genes (Supplementary Table S2). An algorithm of consensus clustering classified LUADs from the training cohort into three CRG clusters based on CRG co-expressed genes. A significant difference in the overall survival of patients in the CRG cluster was found with the KM survival curves. Cluster C had a better prognosis than patients in clusters A and B. Populations in cluster A had a better prognosis than those in cluster B (Figure 2B). The PCA revealed that the cluster A, B, and C were obviously separated from each other (Figure 2C). The infiltration level of the different immune cell populations in CRG clusters was determined by ssGSEA. As shown in Figure 2D, 14 immune cells, including Activated B cell, Activated CD4 T cell, Activated dendritic cell, CD56dim natural killer cell, Eosinophil, Immature B cell, Immature dendritic cell, Mast cell, Monocyte, Natural killer cell, Plasmacytoid dendritic cell, T follicular helper cell, Type 1 T helper cell, and Type 2 T helper cell were recognized significantly distributed in three CRG clusters. Furthermore, there was a statistically significant difference between the three clusters of LUAD patients in terms of age and tumor stage (Figure 2E). For exploring the mechanism of the cluster pattern, we first investigated the different expressed genes between each two clusters, finding there were 26,517 between cluster A and B, 15,995 between cluster A and C, and 10,222 between cluster B and C. We take the intersection of the above three parts and learned 344 CRG-DEGs (Figure 2F; Supplementary Table S3). Then, the KEGG method was used to assess the signaling mechanisms using the CRG-DEGs, which of the top 10 pathways were shown in Figure 2G, including Glycolysis/Gluconeogenesis, Cell cycle, Vascular smooth muscle contraction, Amino sugar and nucleotide sugar metabolism, Retinol metabolism, Central carbon metabolism in cancer, p53 signaling pathway, Apoptosis, Metabolic pathways, and Fc gamma R-mediated phagocytosis.
Two CRG-DEG clusters identification and a CRLncSig generated
Based on the CRG-DEGs, the LUADs from the training cohort were divided into two CRG-DEG clusters using consensus clustering methods. In the CRG-DEG clusters, the KM survival curve showed significant differences in overall survival. The cluster A had a more favor prognosis better than that of patients in cluster B (Figure 3A). The PCA revealed that the cluster A and B were obviously separated from each other (Figure 3B). Through the use of ssGSEA, we determined the degree of infiltration of the different immune cell populations in CRG-DEG clusters. As shown in Figure 3C, 16 immune cells, including Activated B cell, Activated CD4 T cell, Activated dendritic cell, CD56dim natural killer cell, Eosinophil, Gamma delta T cell, Immature B cell, Immature dendritic cell, Macrophage, Mast cell, Monocyte, Natural killer T cell, Natural killer cell, T follicular helper cell, Type 17 T helper cell, and Type 2 T helper cell were recognized significantly distributed in the CRG-DEG clusters. In addition, LUAD patients in the CRG-DEG clusters differed significantly in terms of age, gender, and tumor stage (Figure 3D). By comparing two CRG-DEG clusters, we carried out GSVA to identify KEGG pathways that are most important (Figure 3E; Supplementary Table S4). Surprisingly, KEGG_ALPHA_LINOLENIC_ACID_METABOLISM, KEGG_LINOLEIC_ACID_METABOLISM, KEGG_ARACHIDONIC_ACID_METABOLISM, KEGG_DNA_REPLICATION, KEGG_TAURINE_AND_HYPOTAURINE_METABOLISM, KEGG_CELL_CYCLE, KEGG_ABC_TRANSPORTERS, KEGG_PRIMARY_BILE_ACID_BIOSYNTHESIS, KEGG_FATTY_ACID_METABOLISM, and KEGG_ETHER_LIPID_METABOLISM ranked as the top 10 pathways. Notably, we looked at the distribution of ten CRGs (FDX1, LIAS, LIPT1, DLD, DLAT, PDHA1, PDHB, MTF1, GLS, and CDKN2A) across CRG-DEG clusters and found seven (FDX1, LIAS, LIPT1, PDHB, MTF1, GLS, and CDKN2A) were related to CRG-DEG clusters (Figure 3F). Interestingly, compared to CRG-DEG cluster A, only CDKN2A showed an upregulated situation, while the remaining six displayed downregulated.
FIGURE 3. Two CRG-DEG clusters identified. (A) The KM survival curve demonstrated that the overall survival of patients in the CRG-DEG clusters were significantly different. (B) Principal component analysis scatter plot shows three clusters were obviously separated. (C) Infiltration distribution of immune cells in the two CRG-DEG clusters. (D) Heatmap shows the relationship between CRG-DEG clusters, clinical parameters, and 344 CRG-DEGs. Each row represents a gene, and each column represents a sample. (E) GSVA screens the most important KEGG pathways by comparing two CRG-DEG clusters. (F) Boxplots shows the distribution of ten CRGs (FDX1, LIAS, LIPT1, DLD, DLAT, PDHA1, PDHB, MTF1, GLS, and CDKN2A) across CRG-DEG clusters. CRG: cuproptosis-regulated gene; DEGs: differentially expressed genes; KM: Kaplan–Meier estimator; GSVA: gene set variation analysis; KEGG: Kyoto Encyclopedia of Genes and Genomes; p-value < 0.05 was considered statistically significant; ns: not significant; *: p-value < 0.05; **: p-value < 0.01; ***: p-value < 0.001.
For constructing the prognosis model of the LUAD, we first investigated the DEL between the two CRG-DEG clusters, finding there were 3,711 DELs between them. We take these DELs into KM and Cox analysis to further screen and found 17 DELs meet our criteria (Supplementary Table S5). Then the LASSO analysis was carried out using the 17 DELs for in-depth shrinkage and selection. Finally, 9 lncRNAs were identified (Figures 4A, B), and each gene’s coefficient was obtained (Table 2). In addition, we built Sankey diagrams to descript the relationship of CRG cluster, CRG-DEG cluster, risk, and vital status (Figure 4C). Using boxplots, we found that the distribution of risk scores within CRG cluster was significantly different. Consistently, the distribution of risk scores also differed significantly across the CRG-DEG clusters (Figures 4D, E). We next examined the expression pattern of the ten CRGs (FDX1, LIAS, LIPT1, DLD, DLAT, PDHA1, PDHB, MTF1, GLS, and CDKN2A) in the high- and low-risk groups. In the high-risk group, seven genes (FDX1, LIAS, LIPT1, PDHA1, PDHB, MTF1, and GLS) had downregulated expression values, which statistically significant (Figure 4F). In addition, we displayed the correlations between each of the nine lncRNAs and the ten CRGs, which are shown in Supplementary Figure S1A.
FIGURE 4. Risk model CRLncSig development. (A) LASSO coefficient profiles of prognostic lncRNAs enrolled. (B) LASSO regression with ten-fold cross-validation obtained nine prognostic genes using the minimum Lambda. (C) Sankey diagrams describes the relationship of CRG cluster, CRG-DEG cluster, risk, and vital status. (D) Boxplots shows the distribution of risk score within the CRG cluster was significantly different between each two clusters. (E) Boxplots displays the distribution of risk score was significantly different within the CRG-DEG cluster. (F) Boxplots exhibits expression pattern of the ten CRGs (FDX1, LIAS, LIPT1, DLD, DLAT, PDHA1, PDHB, MTF1, GLS, and CDKN2A) in the high- and low-risk groups. Seven genes (FDX1, LIAS, LIPT1, PDHA1, PDHB, MTF1, and GLS) were found having downregulated expression values in the high-risk group than that in the low-risk group, which statistics significantly. CRLncSig: cuproptosis-regulated lncRNA signature; LASSO: least absolute shrinkage and selection operator; FDR: false discovery rate. CRG: cuproptosis-regulated gene; DEGs: differentially expressed genes; KM: Kaplan–Meier estimator; p-value < 0.05 was considered statistically significant; *: p-value < 0.05; ***: p-value < 0.001.
Independent cohort validation results confirmed the CRLncSig has stable prognostic power
In Supplementary Figures S1B, C, the upper parts showed the patients sorted by increasing risk score, the scatter plot in the middle showed the survival status of the LUADs, and the heatmap in the lower part showed the relative expression levels of nine hub lncRNAs. Each LUAD was assigned a risk score based on the risk score calculator we developed. Patients were categorized into high-risk and low-risk groups according to their median scores. Figure 5A shows the Kaplan-Meier estimator suggested that LUADs at high risk had poorer survival prospects than LUADs at low risk. Additionally, the results of the validation cohort (Figure 5B) indicate that LUAD patients who fall into high-risk groups have a less favorable prognosis. In addition, in Supplementary Figure S2A, we displayed each nine lncRNA’s prognosis ability in the form of Kaplan-Meier curves using the two cohorts’ data, showing that the LINC0183, LINC01711, and LINC00862 performed stable unfavorable impact on LUAD patients, while AC107464.3, LINC00324, AC009120.2, COLCA1, PRKAG2-AS1, and AC093010.2 helped the prognosis improvement of LUADs. Next, we identified whether the risk score can be adopted to independently predict the outcomes of LUAD patients. The risk score and clinical parameters, including age, gender, race, ethnicity, tumor stage, tumor origin, etc. were included in univariate and multivariable analysis (Figure 5C). In the univariate Cox regression, the risk score showed significant prognostic ability. In the multivariate Cox analysis, it was recognized as an independent prognostic factor with a hazard ratio of 3.89 (95% CI: 2.49–6.07, p = 2.25e-09) in the training cohort, 2.45 (95% CI: 1.59–3.77, p = 4.81e-05) in the validation cohort. In addition, the age of validation cohort also showed independent prognostic value, however, its significance did not show a consistence in the two cohorts. Moreover, each gene’s univariate Cox regression was shown in the chart exhibited in Supplementary Figure S2B. We then tested our lncRNA signature’s prognostic ability using ROC analysis (Figure 5D) and time-dependent AUC (Figure 5E). The CRLncSig AUC in the training cohort was 0.737 at 1-year, 0.661 at 3-year, and 0.651 at 5-year, while in the validation cohort was 0.613 at 1-year, 0.647 at 3-year, and 0.635 at 5-year. According to the time-dependent AUC performed in the training cohort (Figure 5E), our risk score was close compared to tumor stage, which is regarded as the prognosis gold standard. Interestingly, when we combined the risk score and tumor stage, their predictive ability AUC could generally reach more than 0.7. Consistently, when we took the risk score and tumor stage combined in the validation cohort (Figure 5E), the predictive ability was stable beyond all factors at all time points, hinting that our risk score is an excellent addition to the tumor stage. As a result of PCA results (Figure 5F), significant heterogeneity was observed between high-risk and low-risk patients in training and validation cohorts, which demonstrated that the risk score model was superior at discriminating between these groups. Additionally, a nomogram using seven factors, including age, grade, tumor, etc., was constructed to predict 1-, 3-, and 5-year overall survival in the TCGA-LUAD cohort (Figure 5G). According to the calibration curves, the nomogram accurately predicted the 1-, 3-, and 5-year overall survival of LUAD patients in the TCGA cohort (Figure 5H).
FIGURE 5. Validation of the CRLncSig in the involved studied cohorts. (A, B) Kaplan-Meier analysis was performed in the training and validation cohorts. Patients in each cohort and subtypes were divided into low-risk groups and high-risk groups based on their median risk score. The log-rank test with a p-value < 0.05 suggests the survival difference is significant. (C) Univariate and multivariable Cox proportional hazards models. *: the variables involved in the studied cohorts, explains as follows: Gender: male vs. female; Race: white vs. non-white; Ethnicity: Hispanic or Latino vs. non-Hispanic or Latino; Prior malignancy: yes vs. no; Tumor origin: upper lobe lung vs. non-upper lobe lung; Smoking history: ever vs. never. (D) ROC curves. The ROC curves value the accuracy for LUAD outcome prediction of our signature at 1-, 3-, and 5-year, respectively. (E) The tAUC analyses compare our signature’s prognostic ability with other available clinical characteristics. The larger the AUC, the stronger the model’s predictive ability. (F) Principal component analysis scatter plot. (G) The nomogram, a quantitative model for predicting clinical prognosis, to predict 1-year, 3-year, and 5-year OS in the LUAD patients of the TCGA-LUAD cohort using seven factors, including age, grade, tumor stage, risk score, smoking history, tissue origin, and prior malignancy. *: p-value < 0.05; ***: p-value < 0.001. (H) The calibration curves indicates that the nomogram accurately predicted the 1-, 3-, and 5-year OS of LUAD patients in the TCGA cohort. CRLncSig, cuproptosis-regulated lncRNA signature; HR, hazard ratio; L95%, 95% confidence interval lower; H95%, 95% confidence interval higher; HR, hazard ratio; ROC, receiver operating characteristic; AUC, area under the ROC curve; tAUC, time-dependent AUC; TCGA, The Cancer Genome Atlas; LUAD, lung adenocarcinoma; PCA, Principal components analysis; OS, overall survival; Exp, relative expression; p-value < 0.05 was considered statistically significant.
The CRLncSig’s relationships with apoptosis, necroptosis, pyroptosis, and ferroptosis
We got apoptosis, necroptosis, pyroptosis, and ferroptosis genes following our criteria, shown in Supplementary Table S6. The Pearson coefficient examined the relationships between our prognosis model and apoptosis, necroptosis, pyroptosis, and ferroptosis -related genes, respectively (Supplementary Table S6; Supplementary Figure S3). The analysis showed that ANLN, SELENBP1, TPX2, METTL7A, CDC20, CENPA, CAPN3, TUBA1C, CEP55, and CDC6 were the top ten correlated apoptosis-related genes, and overall, 2,469/3,681 (67.07%) genes significantly linked with the CRLncSig. CYLD, RIPK3, FADD, TLR3, TRPM7, ZBP1, IPMK, FAS, TNF, and FASLG were discovered as the top necroptosis-related that correlated with the lncRNA signature. As a whole, there were 13/20 (65.00%) necroptosis-related genes correlated with the signature pronouncedly. Moreover, the Pearson test found the top pyroptosis-related genes that correlated with our signature are NLRP1, CARD8, CYCS, CASP1, IRF2, NLRC4, NLRP3, GSDMB, NAIP, and GSDMD. In total, 35/50 (70.00%) pyroptosis-related genes correlated with our signature. The examination found RRM2, SLC2A1, AURKA, DUOX1, CDCA3, IL33, ACADSB, EPT1, SLC7A5, and HILPDA were top ferroptosis-related genes that correlated with our signature. To sum up, there were 238/380 (62.63%) ferroptosis-related genes correlating with our signature.
CRLncSig’s mechanisms were identified by GSEA
As a result of the GSEA analysis based on risk scores, significantly enriched gene sets were identified in the CRLncSig, which the top ten KEGG items were primarily related to cardiac muscle contraction, leukocyte transendothelial migration, arachidonic acid metabolism, ATP-binding cassette transporter, taurine and hypotaurine metabolism, T cell receptor signaling, aldosterone-regulated sodium reabsorption, VEGF signaling pathway, calcium signaling, Proximal tubule bicarbonate reclamation (Figure 6).
FIGURE 6. GSEA analysis with the KEGG gene set as the background identified relevant pathways of the CRLncSig. The significance threshold of this analysis was set as: p-value < 0.05, and FDR < 0.25. GSEA: Gene Set Enrichment Analysis; KEGG: Kyoto Encyclopedia of Genes and Genomes; FDR: false discovery rate.
CRLncSig’s potential links to the LUADs immunological status
There is widespread acceptance that cancer is a dynamic ecosystem in which subclonal populations, mainly cancer cells and non-malignant cells in the tumor microenvironment, cooperate to promote disease progression. A general examination of the tumor microenvironment is therefore necessary. Here, using data from the TCGA cohort, we quantify the scores, including immune score, stromal score, and ESTIMATE score, using the R package “ESTIMATE”. Then through the analysis of the Wilcoxon rank sum test and Pearson correlation coefficient test, it was found that all the scores of categories were lower in the high-risk group, and all scores were negatively correlated with the CRLncSig (Figures 7A–C). Eight mainstream immunoinformatic algorithms, including IPS, TIMER, CIBERSORT, CIBERSORT−ABS, QUANTISEQ, MCPCOUNTER, XCELL, and EPIC, were integrated into the process of screening immune cell types involved in tumor infiltration in high- and low-risk groups. These processes were carried out with Wilcoxon rank-sum test and Pearson correlation coefficient test and displayed in the shapes of the heatmap (Figure 7D) and lollipop (Figure 7E) plots, respectively. In the plots, only significant factors were outlined, and much detail was shown in Supplementary Table S7. The Venn diagram shown in Figure 7F demonstrated the intersection of the results from the heatmap and the lollipop plots, detailed showing that the CD4 T cell, B cell, CD8 T cell, and macrophage were the top cells that our risk correlated in the outline of the word cloud. As for the related immune functions, the scores for the APC_co_stimulation, CCR, Check−point, HLA, Inflammation−promoting, T_cell_co−inhibition, T_cell_co−stimulation, and Type_II_IFN_Reponse were lower in the high-risk group than in the low-risk group (Figure 7G). These findings revealed that the lncRNA signature might be associated with the immunological status of LUADs.
FIGURE 7. In-depth analytics on the relationship between the CRLncSig and the tumor microenvironment condition, immune cell infiltration, immune function. (A) A boxplot and a scatterplot display the relationship between the CRLncSig and immune score. (B) A boxplot and a scatterplot display the relationship between the CRLncSig and stromal score. (C) A boxplot and a scatterplot display the relationship between the CRLncSig and ESTIMATE score. (D) The heatmap demonstrated the infiltration distribution of eight types of immune cells generated from eight mainstream immunoinformatic algorithms in different CRLncSig scores. Only significantly distributed immune cells are exhibited. (E) The lollipop plot describes the infiltration correlations of eight types of immune cells with CRLncSig scores. Only significantly immune cells are shown. (F) The Venn diagram shows the intersection of the results from the heatmap and the lollipop plots, detailed showing that the CD4 T cell, B cell, CD8 T cell, and macrophage were the top cells that our risk correlated in the outline of the word cloud. (G) The violin plot shows the immune functions distributions in the high- and low-risk groups. CRLncSig, cuproptosis-regulated lncRNA signature; p-value < 0.05 was considered statistically significant; ns, not significant; *p-value < 0.05; **p-value < 0.01; ***p-value < 0.001; ****p-value < 0.0001.
The CRLncSig’s participation in immunotherapy and targeting potential immune checkpoint
We thoroughly explored the mutation characteristics of all tumor samples in TCGA-LUAD cohort. Only genes with significant mutation differences between low- and high-risk groups (p < 0.05) were selected and were arranged according to mutation frequency. As shown in the Figure 8A, we listed the top 20 mutated genes showing TP53 mutated most frequently approximately accounting for 53.3% in the cohort, followed by TTN (50.6%) and MUC16 (43.8%). Among the alterations, missense mutation was the most common variant classification. Wilcoxon tests revealed that groups with higher risk scores had higher TMB than groups with lower risk scores. Among risk scores and TMB, Pearson coefficients showed a positive correlation (Figure 8B). According to clinical trials and preclinical studies, the immune checkpoint blockade offers more durable clinical benefits, including treatment responses and long-term survival, to patients with higher TMB (Samstein et al., 2019; Stenzinger et al., 2019). Our results here demonstrated that our high risk LUADs might expect more responses from immunotherapy to a certain extent. Following that, we examined the potential clinical efficacy of immunotherapy based on the risk score subgroups using the TIDE. As a surrogate biomarker, TIDE scores can provide insight into whether a NSCLC patient will respond to immune checkpoint blockades, including anti-PD1 and anti-CTLA4, if therapy is initiated. Immune evasion is more likely to occur in patients with higher TIDE prediction scores, suggesting that immunotherapy has a lesser chance of being effective (Jiang et al., 2018; Fu et al., 2020; Chen et al., 2021a). Based on our results, high-risk patients had lower TIDE scores, meaning that immunotherapy was more beneficial for them (Figure 8C), which corresponded with our “TMB difference” findings above.
FIGURE 8. Determination of the relationship between the CRLncSig and immunotherapy. (A) The waterfall plot shows the top 20 genes mutated in the LUAD, and their difference in the high and low risk groups. (B) The TMB difference in the high and low-risk patients tested by the Wilcoxon rank-sum showing in form of boxplot. The correlation between TMB and the signature tested by the Pearson coefficient showing in form of scatterplot. (C) The TIDE difference in the high and low-risk patients tested by the Wilcoxon rank-sum showing in form of boxplot. The correlation between TIDE and the signature tested by the Pearson coefficient showing in form of scatterplot. (D) The lollipop plot describes the correlations between CRLncSig and immune checkpoints explained by the Pearson coefficient. (E) The violin plot shows the distribution of immune checkpoints in low and high- CRLncSig score tested by the Wilcoxon rank-sum. (F) The Kaplan–Meier estimator found that 8 checkpoints affected the prognosis of LUAD, namely SELP, CD40LG, IL10, IL2, KIR2DL3, CD28, BTLA, and TLR4. (G) The Cox Proportional-Hazards Model revealed 15 prognosis-related genes among 60 checkpoints. (H) The Venn diagram shows the intersection of the above results and found that SELP, CD40LG, IL10, IL2, KIR2DL3, CD28, and BTLA not only closely related to our CRLncSig, but also affected the LUAD prognosis. (I) The correlations between seven found checkpoint genes are shown as a heatmap and performed via the Pearson coefficient. All displayed correlations were statistically significant. (J) Heatmap exhibits these seven checkpoint genes’ roles in the immune system and immunotherapy from multiple immune datasets. In the immunotherapy cohorts (black module), the seven checkpoints were ranked descendingly based on the average score of a group with the immunotherapy cohorts, and they were ranked from high to low as KIR2DL3, IL10, IL2, CD40LG, SELP, BTLA, and CD28. CRLncSig: cuproptosis-regulated lncRNA signature; ns: not significant; *p-value < 0.05; **p-value < 0.01; ***p-value < 0.001; ****p-value < 0.0001; p-value < 0.05 was considered statistically significant; TMB, Tumor mutational burden; TIDE, Tumor Immune Dysfunction and Exclusion.
We curated 60 immune checkpoint genes based on past studies, and we found that 48 of the 60 genes were significantly associated with risk score using the Pearson correlation coefficient (Figure 8D). The top 5 were CD40LG (coefficient = −0.55), SELP (coefficient = −0.48), BTLA (coefficient = −0.48), TNFRSF14 (coefficient = −0.47), and EDNRB (coefficient = −0.44) (Figure 8D). As shown in Figure 8E, 46 checkpoint genes were differentially distributed between high-risk and low-risk groups based on the Wilcoxon signed-rank test. The Kaplan–Meier estimator found that 8 checkpoints affected the prognosis of LUAD, namely, SELP, CD40LG, IL10, IL2, KIR2DL3, CD28, BTLA, and TLR4 (Figure 8F). The Cox Proportional-Hazards Model revealed 15 prognosis-related genes among 60 checkpoints (Figure 8G). In order to simplify the results and balance the outputs of each analysis, we used Venn diagrams to intersect the above results and found that SELP, CD40LG, IL10, IL2, KIR2DL3, CD28, and BTLA not only closely related to our CRLncSig, but also affected the LUAD prognosis, which deserves more attention (Figure 8H). Interestingly, the relationships between each two of the seven checkpoints were calculated, showing all of them have positive correlations and are statistically significant (Figure 8I).
To further evaluate the potential impact of these seven checkpoint genes on the immune system and immunotherapy, we collected data from multiple immune datasets and presented them in a heatmap format (Figure 8J). It can be seen that these seven checkpoint genes exhibited higher T dysfunction values in the GSE12417_GPL570, METABRIC, TCGA Endometrial, and TCGA Melanoma cohorts than in the E-MTAB-179 cohort (grey module). In the immunotherapy cohorts (black module), the seven checkpoints were ranked descendingly based on the average score of a group with the immunotherapy cohorts, and they were ranked from high to low as KIR2DL3, IL10, IL2, CD40LG, SELP, BTLA, and CD28. These findings potentially hinted at the direction of future crosstalk research between our CRLncSig and immunotherapy.
Identification and validation of potential therapeutic agents for high risk score LUADs
CTRP and PRISM datasets contain gene expression profiles and drug sensitivity profiles of hundreds of CCLs, which can be used to build a drug response prediction model. There were 160 compounds shared between the two datasets. In total, 1770 compounds were found in the two datasets after removing duplication (Figure 9A; Supplementary Table S8). We identified candidate agents with greater drug sensitivity for patients with high-risk scores using two different approaches (Figure 9B). Based on CTRP and PRISM data, the analyses were conducted. The first step was to identify compounds with lower estimated AUC values in the high risk score group by performing a differential drug response analysis (log2FC > 0.9) between the highest and lowest risk score decile groups. Following that, having examined the correlation between AUC values and risk scores, compounds with a negative Spearman coefficient of 0.09 were selected. These analyses yielded six CTRP-derived compounds (including leptomycin B, paclitaxel, parbendazole, PHA−793887, triazolothiadiazine, and gemcitabine) (Figure 9C) and seven PRISM-derived compounds (including decitabine, docetaxel, NVP−AUY922, ganetespib, daunorubicin, nobiletin, and gemcitabine) (Figure 9D). AUC values for all of these compounds were lower in the risk score-high group, and there was a negative correlation between risk scores and AUC values.
FIGURE 9. Identification of candidate agents with higher drug sensitivity in high CRLncSig risk score patients. (A) A venn diagram for summarizing included compounds from CTRP and PRISM datasets. (B) Schematic outlining the strategy to identify agents with higher drug sensitivity in high risk score patients. (C) The results of Spearman’s correlation analysis and differential drug response analysis of six CTRP-derived compounds. The lower the value of the y-axis, the greater the drug sensitivity. (D) The results of Spearman’s correlation analysis and differential drug response analysis of seven PRISM-derived compounds. The lower the value of the y-axis, the greater the drug sensitivity. (E) Identification of most promising therapeutic agents for high risk score patients according to the evidence from multiple sources). Six CTRP-derived agents and seven PRISM-derived agents were shown on the left and right of the diagram, respectively. CRLncSig, cuproptosis-regulated lncRNA signature; Pos/Neg, Positively/Negatively; **p-value < 0.01; ***p-value < 0.001; p-value < 0.05 was considered statistically significant.
In spite of the 13 candidate compounds showing a higher drug sensitivity in high-risk patients, the above analyses alone cannot support the conclusion that these compounds are good enough. Further multi-angle analysis was carried out in order to assess their therapeutic potential in LUADs. First, we used CMap analysis to identify compounds whose gene expression patterns were opposed to the LUAD-specific expression patterns (i.e., gene expression increased in tumor tissues but decreased by treatment of certain compounds). Four compounds, including PHA-793887, gemcitabine, daunorubicin, and nobiletin, had CMap scores <−95, representing that these compounds might have potential therapeutic effects in LUADs (Figure 9E; Supplementary Table S9). Second, we calculated the fold-change value representing the difference between the candidate drug target in tumor and normal tissues. An increased fold-change value indicates that the candidate drug is more likely to be effective at treating LUAD (Figure 9E; Supplementary Table S9). The third step was to conduct a comprehensive literature review in PubMed (https://www.ncbi.nlm.nih. gov/pubmed/) to determine the experimental and clinical evidence that candidate compounds are beneficial in treating LUAD (Figure 9E; Supplementary Table S9) (Burkes and Shepherd, 1995; Fossella, 2002; Ramalingam and Belani, 2004; Shimamura et al., 2012; Garon et al., 2013; Momparler, 2013; Liu and Gao, 2017; Wang et al., 2019; Han et al., 2021; Sever et al., 2021). On a whole, gemcitabine, daunorubicin, and nobiletin are considered promising for treating high risk score LUADs based on our robust in silico evidence and confirmed in vitro findings. In contrast, PHA-793887, which showed an excellent CMap score but lacked vitro experiment validation, also deserved attention in further research.
Confirming the nine lncRNAs expression patterns in human tissues using real-time PCR and their potential in pan-cancer
To better understand the real-world expression pattern of each gene in the gene signature, we applied the real-time PCR to detect the above lncRNAs in human LUAD tissues (n = 9) and human normal lung tissues (n = 9) difference in expression. Table 3 shows the primer sequences for the nine lncRNAs, AC009120.2, AC093010.2, AC107464.3, COLCA1, LINC00324, LINC00862, LINC01711, LINC01833, and PRKAG2-AS1. Notably, all lncRNAs had different expression in LUAD and normal lung tissues (Figure 10A). It was found that LINC00324, LINC00862, LINC01711, and LINC01833 lncRNAs were upregulated in LUAD tissues, while the remaining lncRNAs were underexpressed. The upregulated genes in LUAD like LINC00862, LINC01711, and LINC01833 were also showed having unfavorable prognosis powerful in Supplementary Figure S2, and the downregulated genes in LUAD like AC009120.2, AC093010.2, AC107464.3, COLCA1, and PRKAG2-AS1 displayed owning protectable function for the LUAD outcomes, which further proved the validity of the gene signatures we found and provided clues for upcoming deeper research.
FIGURE 10. Confirming the CRLncSig’s nine lncRNAs expression patterns in human tissues using real-time PCR and their potential in pan-cancer. (A) The expression levels of the 9 signature lncRNAs in normal lung (n = 9) and LUAD (n = 9) tissues detected by real-time PCR. Data were means ± standard deviation. (B) The diagrams show the ability of AC107464.3, LINC01711, and COLCA1 in tumor and normal differential expression, prognosis, and tumor staging. CRLncSig: cuproptosis-regulated lncRNA signature; *p-value < 0.05; **p-value < 0.01; ***p-value < 0.001; ****p-value < 0.0001; p-value < 0.05 was considered statistically significant; BLCA, Bladder Urothelial Carcinoma; BRCA, Breast invasive carcinoma; COADREAD, Colon adenocarcinoma/Rectum adenocarcinoma Esophageal carcinoma; KICH, Kidney Chromophobe; KIPAN, Pan-kidney cohort; KIRC, Kidney renal clear cell carcinoma; KIRP, Kidney renal papillary cell carcinoma; LUAD, Lung adenocarcinoma; STAD, Stomach adenocarcinoma; STES, Stomach and Esophageal carcinoma.
Using pan-cancer expression patterns as a starting point, we assessed the potential of these nine lncRNAs. We analyzed the expression data of 34 cancer types, finding six of the nine lncRNAs, AC009120.2, LINC01833, PRKAG2-AS1, AC107464.3, COLCA1, and LINC01711, showed significant differences (tumor vs. normal) in 30 or more cancer types (Supplementary Figure S4; Table 3). Taking a step further, we then explored the prognostic ability of nine lncRNAs in pan-cancer. For this reason, we curated 44 cancer types according to the screening criteria and built the Cox models. The results displayed that the top five lncRNAs were AC107464.3, LINC01711, LINC00862, LINC00324, and AC009120.2; they affected 17, 16, 14, 14, and 12 cancer types, respectively (Supplementary Figure S5; Table 3). Finally, we tested the expression distribution of nine lncRNAs in different cancer stages, and according to the criteria we set, we found data from 30 cancer types for this analysis. The results showed that LINC01711, AC107464.3, LINC00324, COLCA1, and PRKAG2-AS1 were the top five lncRNAs, and they were differentially expressed in cancer stages of 12, 8, 8, 8, and 8 cancer types, respectively (Supplementary Figure S6; Table 3). According to our above results, we tried to find which lncRNA was significant in tumor and normal differential expression, prognostic ability, and tumor staging power and found that AC107464.3, LINC01711, and COLCA1 occupied the top 3 most counts of cancer types, which was 7, 6, and 5, respectively (Figure 10B; Table 3). Our brief exploration of the nine lncRNAs and pan-cancer prove the importance of our CRLncSig, which may shed light on future research in other cancers.
Discussion
In the past few years, copper has become increasingly recognized as an essential mineral nutrient for cell proliferation and death due to its intrinsic redox properties (Lelievre et al., 2020; Ge et al., 2022). This field has brought together researchers from various disciplines due to its rapid development. Including researchers could help translate basic copper chemistry and biology research into clinical applications, especially for cancer treatment (Ge et al., 2022). Interestingly, Tsvetkov found that copper can bind directly to the lipid acylation component of the tricarboxylic acid (TCA) cycle, causing it to aggregate lipid acylated proteins and iron-sulfur clusters (Tsvetkov et al., 2022). Unlike known forms of cell death such as pyroptosis, apoptosis, ferroptosis, and necroptosis, cuproptosis is completely new. Considering that the prognostic outcomes of LUAD vary widely, the development of a robust classifier is crucial to maximizing the benefits of personalized treatment and timely follow-up for patients of varying risk and prognosis. In the present research, we innovatively adopt the novel cuproptosis concept to establish CRLncSig predicting LUAD outcomes by digging TCGA and GEO databases (with a total human sample size exceeding 1,000 cases). Specifically, our novelty lay in adopting comprehensive bioinformatics analysis, including consensus clustering, LASSO regression, Kaplan-Meier curves, Cox models, ROC curves, and tAUC, and the validations in a large cohort. Our analysis revealed that the CRLncSig is potentially linked to the LUADs immunological status and participates in immunotherapy and targeting potential immune checkpoints like SELP, CD40LG, IL10, IL2, KIR2DL3, CD28, and BTLA. Above all, the in silico screening of 1770 compounds identified three drugs with potential therapeutic implications for high-risk LUADs, gemcitabine, daunorubicin, and nobiletin. Lastly, we describe signature lncRNAs’ real-world expression patterns using real-time PCR and their potential in pan-cancer by analyzing multi-cohort.
Many studies have revealed relevant risk parameters for its progression in malignant hepatocellular neoplasms (Moriyama et al., 2022; Park et al., 2022; Yang et al., 2022). However, the increased incidence of hepatocellular carcinoma in Wilson’s patients and in animal models has raised researchers’ concerns that abnormal copper accumulation may be related to an unidentified mechanism of tumor progression (Ge et al., 2022). Many studies have shown that tumors require a greater amount of copper than healthy tissue, leading to a link between copper and cancer (Ge et al., 2022). The balance between gastrointestinal absorption and biliary excretion is maintained for copper homeostasis once mammalian growth and development has been completed. Nevertheless, anabolic status requires significantly more copper in the diet to meet metabolic demands (Ge et al., 2022). Copper deficiency during embryonic and fetal development results in numerous structural and biochemical abnormalities (Ge et al., 2022). Since copper is necessary for mitochondrial cytochrome c oxidase, which is required for rapid cell division, cancer cells have an increased requirement for copper (Ge et al., 2022). Various cancers, including lung cancer, have been linked to elevated copper concentrations in animal models and patients (Ge et al., 2022). Copper imbalance affects mitochondrial respiration and interferes with glycolysis, insulin resistance, and lipid metabolism (Ge et al., 2022). Additionally, copper pathways, such as the ATOX-ATP7A-LOX pathway, contribute to metastatic spread (Ge et al., 2022). Aside from affecting tumor angiogenesis, copper regulates autophagy through ULK1 and ULK2, as well as promoting tumor formation, growth, and metastasis (Ge et al., 2022). As a result of this metal nutrient, copper, a number of pro-angiogenic factors are activated, including vascular endothelial growth factor, fibroblast growth factor 2, tumor necrosis factor, and interleukin 156. Many recent studies link copper signaling to cancer cell proliferation, tumor growth, and metastasis (Ge et al., 2022). Copper and cancer metabolism must be studied intensively in order to understand all aspects of cancer biology, such as tumor initiation, growth, and metastasis (Ge et al., 2022). It is urgently needed to find meaningful cuproptosis-related signatures for predicting tumor prognosis, which will be useful for understanding cuproptosis’ mechanism in cancer. According to Zhang and collogues, a prognostic lncRNA profile related to cuproptosis may offer new perspectives on hepatocellular carcinoma therapy (Zhang et al., 2022b). Bian’s team constructed a cuproptosis-related gene signature as a potential prognostic predictor in patients with clear cell renal cell carcinoma, offering novel insights into cancer treatment (Bian et al., 2022). For the moment, in LUAD, there is still no study showing that any cuproptosis-related signature can predict prognosis. Our present study attempts to fill this gap and construct a cuproptosis related lncRNA signature for predicting the outcomes of LUAD, providing more clues for this “virgin land” in follow-up studies.
Our signature contains nine lncRNAs (Table 2), which were AC009120.2, AC093010.2, AC107464.3, COLCA1, LINC00324, LINC00862, LINC01711, LINC01833, and PRKAG2-AS1. Furthermore, we performed real-time PCR validation and found that all lncRNAs were significantly expressed differently between LUAD and normal human lung tissues. A negative effect was demonstrated by LINC00862, LINC01711, and LINC01833 on LUAD outcomes, while the remaining lncRNAs detected a positive impact (Supplementary Figure S2). There are rare studies on LINC00862, which were only mentioned in one recent study. According to Yu and colleagues, LINC00862 levels were elevated in HCC tumor nodules and cirrhotic livers, suggesting it might be involved in the progression of HCC (Yu et al., 2021). Moreover, these researchers also examined the expression of LINC00862 in human HCC cell lines and found that it is upregulated in several HCC cell lines in comparison to THLE-2, a normal human liver cell line (Yu et al., 2021). Elevated expression of LINC01711 in bladder cancer correlates with reduced survival (Du et al., 2021). Through its interaction with miR-326 and FSCN1, LINC01711 has been shown to promote esophageal squamous cell carcinoma initiation and progression (Xu et al., 2021). Additionally, LINC01711 expression was positively correlated with TGF-β1, a critical factor in the TGF-β signaling pathway (Lee et al., 2005). Unfortunately, the research on LINC01711 and lung cancer is still lacking, and more efforts are needed. A study showed that LINC01833 sponging miR-519e-3p and regulating S100A4 expression promote LUAD invasion (Zhang et al., 2020). Overexpression of LINC01833 improves the proliferation and invasion ability of lung cancer cells, as well as promotes the transition from epithelial to mesenchymal state (Zhang et al., 2020).
In their study, Mitra et al. (2012) pointed out that excessive copper exposure can lead to apoptosis and cell cycle arrest, potentially causing immunotoxicity, and found that copper-induced immunosuppression appears to be regulated by the apoptotic pathway. In several human diseases, including cancer, necroptosis causes caspase-independent programmed cell death mediated by the MLKL signaling cascade (Zhang et al., 2022c). A high level of inflammation results in pyroptosis, a form of lysis-programmed cell death that occurs mainly during intracellular pathogens within the cell and may form part of an antimicrobial response (Tan et al., 2021). Ferroptosis is characterized by the inactivation of GPX4, followed by the accumulation of ROS and the binding of free Fe2+ to membrane lipids (Li et al., 2022c). Copper depletion limited GPX4 protein expression, the only enzyme capable of protecting against lipid peroxidation, according to Li et al. (2022c). In the presence of copper depletion, the dermal papilla cells (DPCs) were less sensitive to erastin (an inducer of ferroptosis), while the ferroptosis inhibitor ferrostatin-1 (Fer-1) partly prevented bathocuproinedisulfonic (BCS)-induced cell death (Li et al., 2022c). While the way of regulated cell death of cuproptosis, apoptosis, necroptosis, pyroptosis, and ferroptosis varies (Chen et al., 2021b; Tsvetkov et al., 2022), from our research, they seem to be somewhat related, such as our cuproptosis-related signature correlated with 2,321/3,681 (63.05%) apoptosis-related genes, 11/20 (55.00%) necroptosis-related genes, 34/50 (68.00%) pyroptosis-related genes, and 222/380 (58.42%) ferroptosis-related genes, providing potential explanations and inspirations for further research of cell death-related tumor mechanism.
Cancer immunotherapy prolongs the survival of deadly cancer patients and as more cancer patients become eligible for immune-based cancer treatments, reflecting that this approach is revolutionizing the field of oncology (Vanneman and Dranoff, 2012; Waldman et al., 2020). New therapeutic combinations and newly discovered drug targets are expanding the use of immunotherapy in cancer treatment (Vanneman and Dranoff, 2012; Waldman et al., 2020). Targeted strategies inhibit tumor progression by interfering with crucial molecular pathways, while immunotherapy produces durable and effective tumor destruction by stimulating the host’s own response (Vanneman and Dranoff, 2012; Waldman et al., 2020). The main challenge of immunotherapy is how to determine whether a certain biomarker is suitable for the host, and how to adjust its application strategy to maximize the benefits (Suresh et al., 2018). This study gives hints about which immunotherapy targets to use or under what circumstances to apply. We first found that our risk score was associated with TMB and TIDE, suggesting that our signature appeared to guide immunotherapy. Next, we followed the trail and found seven checkpoints, including SELP, CD40LG, IL10, IL2, KIR2DL3, CD28, and BTLA elated to our signature score. In the immunotherapy cohorts we chosen, the seven checkpoints were ranked from high to low as KIR2DL3, IL10, IL2, CD40LG, SELP, BTLA, and CD28. KIR2DL3, a transmembrane glycoprotein expressed by natural killer cells and T cell subsets, is responsible for sending inhibitory signals throughout the cell (Sim et al., 2016). KIR2DL3’s structure is important for cancer development (Sim et al., 2016). Neuroblastoma development is influenced by KIR2DL3, according to Sezgin et al. (2022). IL-10 is a cytokine with potent anti-inflammatory properties. The function of IL-10 is to prevent damage to the host and maintain normal tissue homeostasis by limiting the host’s immune response to pathogens (Oft, 2014). IL-10 suppresses inflammatory Th17 T cells and macrophages, which can trigger or promote tumor formation (Oft, 2014). The research by Vahl and colleagues showed that IL-10 competes with IFN-γ on the PD1/PDL1 pathway, potentially contributing to resistance to PD1/PDL1 immunotherapy in lung cancer patients (Vahl et al., 2017). A key function of IL-2 is activating the immune system, which could potentially eradicate cancer (Jiang et al., 2016). The FDA has approved the use of IL-2 as a monotherapy for metastatic renal cell carcinoma and metastatic melanoma (Jiang et al., 2016). A decrease in IL-2 levels and a marked increase in soluble IL-2 receptor concentrations have been observed in advanced NSCLC, and these are associated with poor outcomes (Jiang et al., 2016). The role of IL-2 activation in restoring lymphocyte immunocompetence against lung cancer has also been demonstrated (Jiang et al., 2016).
Individuals with LUAD exhibit high levels of heterogeneity, making it nearly impossible to find a treatment that is effective for everyone (Zito Marino et al., 2019). All current treatments for advanced LUAD lack corresponding biomarkers, so they do not achieve satisfactory therapeutic effects because they are population-based (Zito Marino et al., 2019). One of the main goals of this study is to find a tailored treatment strategy for a specific population, which is crucial to maximizing treatment effectiveness. As well as providing information on prognosis, the CRLncSig risk score can be used in precision oncology to guide targeted therapy. In particular, our study discovered three potential therapeutic agents for high-risk LUAD patients: gemcitabine, daunorubicin, and nobiletin. Non-small cell lung cancer is commonly treated with gemcitabine, a synthetic antimetabolite tumor drug (Sederholm et al., 2005). A synthetic version of gemcitabine was developed by Larry Hertel in the early 1980s for use as an antiviral medication. However, preclinical tests showed it could kill leukemia cells in vitro as well (Burkes and Shepherd, 1995). Gemcitabine was approved for the treatment of non-small cell lung cancer by the FDA in 1998 (Barton-Burke, 1999). In NSCLC, gemcitabine monotherapy has shown a response rate greater than 20%, a median survival of 7–9 months, and favorable side effects in more than 500 patients in six phase II studies (Sederholm et al., 2005). Although gemcitabine, which has been extensively studied, is effective for most lung cancers, some patients cannot get effective drug responses due to the heterogeneity of lung cancers (Zito Marino et al., 2019). The CRLncSig score of our research can potentially be an excellent indicator to guide clinical gemcitabine medication. Known as daunorubicin, daunomycin is an anthracycline antibiotic that binds to DNA and causes helical unwinding, ultimately inhibiting DNA synthesis and DNA-dependent RNA synthesis (Murphy and Yee, 2017). This drug slows or stops the growth of cancer cells and can be used to treat acute myeloid leukemia, acute lymphoblastic leukemia, chronic myelogenous leukemia, and Kaposi’s sarcoma (Murphy and Yee, 2017). Potapov et al. demonstrated that compared to doxorubicin, daunomycin lowered the proportion of the DNA-synthesizing cells in the drug-sensitive xenografts of lung cancer more significantly (Potapov Iu et al., 1986). The experimental data suggested that daunomycin would be highly efficient in the chemotherapy of patients with lung cancer (Potapov Iu et al., 1986). Using co-delivery liposomes targeting daunorubicin and dioscin, Wang and colleagues demonstrated significant antitumor effects in tumor-bearing mice, which may provide an effective strategy for the treatment of NSCLC (Wang et al., 2019). Nobiletin is a naturally occurring compound with potential therapeutic effects against a variety of cancer types (Ashrafizadeh et al., 2020). Compound 29d, a derivative of nobiletin, has been shown to increase paclitaxel accumulation in lung cancer cells by reducing P-gp activity, thereby enhancing its antitumor activity (Feng et al., 2020). In terms of enhancing the antitumor activity of adriamycin, nobiletin inhibits Akt and Wnt/β-catenin signaling pathway by increasing GSK-3β activity, and leads to decreased lung cancer cells viability (Moon et al., 2018).
There are some limitations to this study. We generated this CRLncSig from publicly accessible data. Although it has been confirmed to have stable prognosis ability through applied to another large independent cohort and have differential expression patterns in tumor and normal tissues via silico and real-time PCR approaches, its clinical applicability needs further confirmation with more parameters. Furthermore, there are still no wet laboratory facts to hold up the nine lncRNAs’ parts in cuproptosis-related mechanisms. Therefore, more research, which focuses in vivo and in vitro, is urgently needed to reveal more clues that support the signature’s potential future.
Conclusion
The present research constructed a novel and capable cuproptosis-related lncRNA signature, CRLncSig, for LUAD. Applying the signature to a large independent cohort and assessing its human tissue expression pattern using real-time PCR validated its stability and broad applicability. The signature owns the potential ability to undertake the role of precise immunotherapy. Our study identified potential immunotherapy targets and agents which might improve patients’ prognoses most effectively for those with high CRLncSig scores. In summary, this study has introduced new insights into personalized prognostication approaches and shed light on the integration of tailored prognosis prediction and precision therapy.
Data availability statement
Publicly available datasets were analyzed in this study. This data can be found here: This study used publicly available datasets. Their names and where to get them are as follows: TCGA and pan-cancer TCGA TARGET GTEx, https://xenabrowser.net; GSE29013, GSE30219, GSE31210, GSE37745, and GSE50081, https://www.ncbi.nlm.nih.gov/geo; CCLs database, https://depmap.org/portal; CTRP database, https://portals.broadinstitute.org/ctrp; PRISM database, https://depmap.org/portal/prism.
Ethics statement
Approval for this research was granted by the ethics committee of the First Affiliated Hospital of Zhengzhou University (2022-KY-0109-004), and all patients who participated in the research project gave their written informed consent. Every experiment was undertaken in compliance with the institutional ethical standards that Zhengzhou University had adopted.
Author contributions
CM contributed to the study design, data analysis, and draft manuscript. FL and ZH analyzed the data. ZH conducted a literature review. SZ, YY, and ZG revised the entire manuscript. ZG supervised the study. All authors approved the final submitted version.
Acknowledgments
CM thanks the technical support from the post-doctoral research station of the First Affiliated Hospital of Zhengzhou University.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar.2023.1113808/full#supplementary-material
References
Ashrafizadeh, M., Zarrabi, A., Saberifar, S., Hashemi, F., Hushmandi, K., Hashemi, F., et al. (2020). Nobiletin in cancer therapy: How this plant derived-natural compound targets various oncogene and onco-suppressor pathways. Biomedicines 8 (5), 110. doi:10.3390/biomedicines8050110
Barton-Burke, M. (1999). Gemcitabine: A pharmacologic and clinical overview. Cancer Nurs. 22 (2), 176–183.
Bian, Z., Fan, R., and Xie, L. (2022). A novel cuproptosis-related prognostic gene signature and validation of differential expression in clear cell renal cell carcinoma. Genes (Basel). 13 (5), 851. doi:10.3390/genes13050851
Burkes, R. L., and Shepherd, F. A. (1995). Gemcitabine in the treatment of non-small-cell lung cancer. Ann. Oncol. 6 (3), S57–S60. doi:10.1093/annonc/6.suppl_3.s57
Cao, R., and Lopez-de-Ullibarri, I. (2019). ROC curves for the statistical analysis of microarray data. Methods Mol. Biol. 1986, 245–253. doi:10.1007/978-1-4939-9442-7_11
Chalmers, Z. R., Connelly, C. F., Fabrizio, D., Gay, L., Ali, S. M., Ennis, R., et al. (2017). Analysis of 100,000 human cancer genomes reveals the landscape of tumor mutational burden. Genome Med. 9 (1), 34. doi:10.1186/s13073-017-0424-2
Chen, H., Carrot-Zhang, J., Zhao, Y., Hu, H., Freeman, S. S., Yu, S., et al. (2019). Genomic and immune profiling of pre-invasive lung adenocarcinoma. Nat. Commun. 10 (1), 5472. doi:10.1038/s41467-019-13460-3
Chen, X., Kang, R., Kroemer, G., and Tang, D. (2021). Broadening horizons: The role of ferroptosis in cancer. Nat. Rev. Clin. Oncol. 18 (5), 280–296. doi:10.1038/s41571-020-00462-0
Chen, Y., Li, Z. Y., Zhou, G. Q., and Sun, Y. (2021). An immune-related gene prognostic index for head and neck squamous cell carcinoma. Clin. Cancer Res. 27 (1), 330–341. doi:10.1158/1078-0432.ccr-20-2166
Clough, E., and Barrett, T. (2016). The gene expression Omnibus database. Methods Mol. Biol. 1418, 93–110. doi:10.1007/978-1-4939-3578-9_5
Du, Y., Wang, B., Jiang, X., Cao, J., Yu, J., Wang, Y., et al. (2021). Identification and validation of a stromal EMT related LncRNA signature as a potential marker to predict bladder cancer outcome. Front. Oncol. 11, 620674. doi:10.3389/fonc.2021.620674
Feng, S., Zhou, H., Wu, D., Zheng, D., Qu, B., Liu, R., et al. (2020). Nobiletin and its derivatives overcome multidrug resistance (MDR) in cancer: Total synthesis and discovery of potent MDR reversal agents. Acta Pharm. Sin. B 10 (2), 327–343. doi:10.1016/j.apsb.2019.07.007
Fossella, F. V. (2002). Docetaxel for previously treated non-small-cell lung cancer. Oncol. Willist. Park) 16 (6), 45–51.
Friedman, J., Hastie, T., and Tibshirani, R. (2010). Regularization paths for generalized linear models via coordinate descent. J. Stat. Softw. 33 (1), 1–22. doi:10.18637/jss.v033.i01
Fu, J., Li, K., Zhang, W., Wan, C., Zhang, J., Jiang, P., et al. (2020). Large-scale public data reuse to model immunotherapy response and resistance. Genome Med. 12 (1), 21. doi:10.1186/s13073-020-0721-z
Garon, E. B., Finn, R. S., Hamidi, H., Dering, J., Pitts, S., Kamranpour, N., et al. (2013). The HSP90 inhibitor NVP-AUY922 potently inhibits non-small cell lung cancer growth. Mol. Cancer Ther. 12 (6), 890–900. doi:10.1158/1535-7163.MCT-12-0998
Ge, E. J., Bush, A. I., Casini, A., Cobine, P. A., Cross, J. R., DeNicola, G. M., et al. (2022). Connecting copper and cancer: From transition metal signalling to metalloplasia. Nat. Rev. Cancer 22 (2), 102–113. doi:10.1038/s41568-021-00417-2
Ghandi, M., Huang, F. W., Jane-Valbuena, J., Kryukov, G. V., Lo, C. C., McDonald, E. R., et al. (2019). Next-generation characterization of the cancer cell line Encyclopedia. Nature 569 (7757), 503–508. doi:10.1038/s41586-019-1186-3
Goeman, J. J. (2010). L1 penalized estimation in the Cox proportional hazards model. Biom J. 52 (1), 70–84. doi:10.1002/bimj.200900028
Han, S. H., Han, J. H., Chun, W. J., Lee, S. S., Kim, H. S., and Lee, J. W. (2021). Nobiletin inhibits non-small-cell lung cancer by inactivating WNT/β-Catenin signaling through downregulating miR-15-5p. Evid. Based Complement. Altern. Med. 2021, 7782963. doi:10.1155/2021/7782963
Hu, H., Xu, Q., Mo, Z., Hu, X., He, Q., Zhang, Z., et al. (2022). New anti-cancer explorations based on metal ions. J. Nanobiotechnology 20 (1), 457. doi:10.1186/s12951-022-01661-w
Ji, Z. H., Ren, W. Z., Wang, H. Q., Gao, W., and Yuan, B. (2022). Molecular subtyping based on cuproptosis-related genes and characterization of tumor microenvironment infiltration in kidney renal clear cell carcinoma. Front. Oncol. 12, 919083. doi:10.3389/fonc.2022.919083
Jiang, P., Gu, S., Pan, D., Fu, J., Sahu, A., Hu, X., et al. (2018). Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat. Med. 24 (10), 1550–1558. doi:10.1038/s41591-018-0136-1
Jiang, T., Zhou, C., and Ren, S. (2016). Role of IL-2 in cancer immunotherapy. Oncoimmunology 5 (6), e1163462. doi:10.1080/2162402X.2016.1163462
Johnson, W. E., Li, C., and Rabinovic, A. (2007). Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics 8 (1), 118–127. doi:10.1093/biostatistics/kxj037
Lee, F. T., Mountain, A. J., Kelly, M. P., Hall, C., Rigopoulos, A., Johns, T. G., et al. (2005). Enhanced efficacy of radioimmunotherapy with 90Y-CHX-A''-DTPA-hu3S193 by inhibition of epidermal growth factor receptor (EGFR) signaling with EGFR tyrosine kinase inhibitor AG1478. Clin. Cancer Res. 11 (19), 7080s–7086s. doi:10.1158/1078-0432.CCR-1004-0019
Lelievre, P., Sancey, L., Coll, J. L., Deniaud, A., and Busser, B. (2020). 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). 12 (12), 3594. doi:10.3390/cancers12123594
Li, F., Niu, Y., Zhao, W., Yan, C., and Qi, Y. (2022). Construction and validation of a prognostic model for lung adenocarcinoma based on endoplasmic reticulum stress-related genes. Sci. Rep. 12 (1), 19857. doi:10.1038/s41598-022-23852-z
Li, F., Wu, X., Liu, H., Liu, M., Yue, Z., Wu, Z., et al. (2022). Copper depletion strongly enhances ferroptosis via mitochondrial perturbation and reduction in antioxidative mechanisms. Antioxidants (Basel). 11 (11), 2084. doi:10.3390/antiox11112084
Li, J., Chen, S., Liao, Y., Wang, H., Zhou, D., and Zhang, B. (2022). Arecoline is associated with inhibition of cuproptosis and proliferation of cancer-associated fibroblasts in oral squamous cell carcinoma: A potential mechanism for tumor metastasis. Front. Oncol. 12, 925743. doi:10.3389/fonc.2022.925743
Li, X., Ma, C., Luo, H., Zhang, J., Wang, J., and Guo, H. (2020). Identification of the differential expression of genes and upstream microRNAs in small cell lung cancer compared with normal lung based on bioinformatics analysis. Med. Baltim. 99 (11), e19086. doi:10.1097/MD.0000000000019086
Liberzon, A., Subramanian, A., Pinchback, R., Thorvaldsdottir, H., Tamayo, P., and Mesirov, J. P. (2011). Molecular signatures database (MSigDB) 3.0. Bioinformatics 27 (12), 1739–1740. doi:10.1093/bioinformatics/btr260
Lin, W., Chen, Y., Wu, B., Chen, Y., and Li, Z. (2021). Identification of the pyroptosisrelated prognostic gene signature and the associated regulation axis in lung adenocarcinoma. Cell Death Discov. 7 (1), 161. doi:10.1038/s41420-021-00557-2
Liu, J., Lichtenberg, T., Hoadley, K. A., Poisson, L. M., Lazar, A. J., Cherniack, A. D., et al. (2018). An integrated TCGA pan-cancer clinical data resource to drive high-quality survival outcome analytics. Cell 173 (2), 400–416.e11. doi:10.1016/j.cell.2018.02.052
Liu, Z., and Gao, W. (2017). Leptomycin B reduces primary and acquired resistance of gefitinib in lung cancer cells. Toxicol. Appl. Pharmacol. 335, 16–27. doi:10.1016/j.taap.2017.09.017
Lu, Y., Luo, X., Wang, Q., Chen, J., Zhang, X., Li, Y., et al. (2022). A novel necroptosis-related lncRNA signature predicts the prognosis of lung adenocarcinoma. Front. Genet. 13, 862741. doi:10.3389/fgene.2022.862741
Ma, C., Li, F., He, Z., and Zhao, S. (2022). A more novel and powerful prognostic gene signature of lung adenocarcinoma determined from the immune cell infiltration landscape. Front. Surg. 9, 1015263. doi:10.3389/fsurg.2022.1015263
Ma, C., Li, F., and Luo, H. (2021). Prognostic and immune implications of a novel ferroptosis-related ten-gene signature in lung adenocarcinoma. Ann. Transl. Med. 9 (13), 1058. doi:10.21037/atm-20-7936
Ma, C., Li, F., Wang, Z., and Luo, H. (2022). A novel immune-related gene signature predicts prognosis of lung adenocarcinoma. Biomed. Res. Int. 2022, 4995874. doi:10.1155/2022/4995874
Ma, C., Luo, H., Cao, J., Gao, C., Fa, X., and Wang, G. (2020). Independent prognostic implications of RRM2 in lung adenocarcinoma. J. Cancer 11 (23), 7009–7022. doi:10.7150/jca.47895
Ma, C., Luo, H., Cao, J., Zheng, X., Zhang, J., Zhang, Y., et al. (2020). Identification of a novel tumor microenvironment-associated eight-gene signature for prognosis prediction in lung adenocarcinoma. Front. Mol. Biosci. 7, 571641. doi:10.3389/fmolb.2020.571641
Ma, C., Luo, H., Liu, B., Li, F., Tschope, C., and Fa, X. (2018). Long noncoding RNAs: A new player in the prevention and treatment of diabetic cardiomyopathy? Diabetes Metab. Res. Rev. 34 (8), e3056. doi:10.1002/dmrr.3056
Mitra, S., Keswani, T., Dey, M., Bhattacharya, S., Sarkar, S., Goswami, S., et al. (2012). Copper-induced immunotoxicity involves cell cycle arrest and cell death in the spleen and thymus. Toxicology 293 (1-3), 78–88. doi:10.1016/j.tox.2011.12.013
Momparler, R. L. (2013). Epigenetic therapy of non-small cell lung cancer using decitabine (5-aza-2'-deoxycytidine). Front. Oncol. 3, 188. doi:10.3389/fonc.2013.00188
Moon, J. Y., Manh Hung, L. V., Unno, T., and Cho, S. K. (2018). Nobiletin enhances chemosensitivity to adriamycin through modulation of the akt/gsk3β/β⁻Catenin/MYCN/MRP1 signaling pathway in A549 human non-small-cell lung cancer cells. Nutrients 10 (12), 1829. doi:10.3390/nu10121829
Moriyama, M., Kanda, T., Midorikawa, Y., Matsumura, H., Masuzaki, R., Nakamura, H., et al. (2022). The proliferation of atypical hepatocytes and CDT1 expression in noncancerous tissue are associated with the postoperative recurrence of hepatocellular carcinoma. Sci. Rep. 12 (1), 20508. doi:10.1038/s41598-022-25201-6
Murphy, T., and Yee, K. W. L. (2017). Cytarabine and daunorubicin for the treatment of acute myeloid leukemia. Expert Opin. Pharmacother. 18 (16), 1765–1780. doi:10.1080/14656566.2017.1391216
Oft, M. (2014). IL-10: Master switch from tumor-promoting inflammation to antitumor immunity. Cancer Immunol. Res. 2 (3), 194–199. doi:10.1158/2326-6066.CIR-13-0214
Park, I., Kim, N., Lee, S., Park, K., Son, M. Y., Cho, H. S., et al. (2022). Characterization of signature trends across the spectrum of non-alcoholic fatty liver disease using deep learning method. Life Sci. 314, 121195. doi:10.1016/j.lfs.2022.121195
Potapov Iu, N., Krutova, T. V., Korman, D. B., and Pashkova, V. S. (1986). [Cytotoxic action of doxorubicin and daunomycin on human lung cancer cells]. Antibiot. Med. Biotekhnol 31 (8), 614–617.
Ramalingam, S., and Belani, C. P. (2004). Paclitaxel for non-small cell lung cancer. Expert Opin. Pharmacother. 5 (8), 1771–1780. doi:10.1517/14656566.5.8.1771
Samstein, R. M., Lee, C. H., Shoushtari, A. N., Hellmann, M. D., Shen, R., Janjigian, Y. Y., et al. (2019). Tumor mutational load predicts survival after immunotherapy across multiple cancer types. Nat. Genet. 51 (2), 202–206. doi:10.1038/s41588-018-0312-8
Sauerbrei, W., Royston, P., and Binder, H. (2007). Selection of important variables and determination of functional form for continuous predictors in multivariable model building. Stat. Med. 26 (30), 5512–5528. doi:10.1002/sim.3148
Sederholm, C., Hillerdal, G., Lamberg, K., Kolbeck, K., Dufmats, M., Westberg, R., et al. (2005). Phase III trial of gemcitabine plus carboplatin versus single-agent gemcitabine in the treatment of locally advanced or metastatic non-small-cell lung cancer: The Swedish lung cancer study group. J. Clin. Oncol. 23 (33), 8380–8388. doi:10.1200/JCO.2005.01.2781
Sever, B., Altintop, M. D., Ciftci, G. A., and Ozdemir, A. (2021). A new series of triazolothiadiazines as potential anticancer agents for targeted therapy of non-small cell lung and colorectal cancers: Design, synthesis, in silico and in vitro studies providing mechanistic insight into their anticancer potencies. Med. Chem. 17 (10), 1104–1128. doi:10.2174/1573406416666201021142832
Sezgin, G., Goruroglu Ozturk, O., Ozkan, A., Kupeli, S., and Bayram, I. (2022). Clinical impact of KIR2DS3 and KIR2DL3 genes in neuroblastoma patients. Med. Princ. Pract. 31, 532–539. doi:10.1159/000524656
Shimamura, T., Perera, S. A., Foley, K. P., Sang, J., Rodig, S. J., Inoue, T., et al. (2012). Ganetespib (STA-9090), a nongeldanamycin HSP90 inhibitor, has potent antitumor activity in in vitro and in vivo models of non-small cell lung cancer. Clin. Cancer Res. 18 (18), 4973–4985. doi:10.1158/1078-0432.CCR-11-2967
Siegel, R. L., Miller, K. D., Fuchs, H. E., and Jemal, A. (2021). Cancer statistics, 2021. CA Cancer J. Clin. 71 (1), 7–33. doi:10.3322/caac.21654
Sim, M. J., Stowell, J., Sergeant, R., Altmann, D. M., Long, E. O., and Boyton, R. J. (2016). KIR2DL3 and KIR2DL1 show similar impact on licensing of human NK cells. Eur. J. Immunol. 46 (1), 185–191. doi:10.1002/eji.201545757
Stenzinger, A., Allen, J. D., Maas, J., Stewart, M. D., Merino, D. M., Wempe, M. M., et al. (2019). Tumor mutational burden standardization initiatives: Recommendations for consistent tumor mutational burden assessment in clinical samples to guide immunotherapy treatment decisions. Genes Chromosom. Cancer 58 (8), 578–588. doi:10.1002/gcc.22733
Strasser, A., and Vaux, D. L. (2020). Cell death in the origin and treatment of cancer. Mol. Cell 78 (6), 1045–1054. doi:10.1016/j.molcel.2020.05.014
Sung, H., Ferlay, J., Siegel, R. L., Laversanne, M., Soerjomataram, I., Jemal, A., et al. (2021). Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J. Clin. 71 (3), 209–249. doi:10.3322/caac.21660
Suresh, K., Naidoo, J., Lin, C. T., and Danoff, S. (2018). Immune checkpoint immunotherapy for non-small cell lung cancer: Benefits and pulmonary toxicities. Chest 154 (6), 1416–1423. doi:10.1016/j.chest.2018.08.1048
Taminau, J., Meganck, S., Lazar, C., Steenhoff, D., Coletta, A., Molter, C., et al. (2012). Unlocking the potential of publicly available microarray data using inSilicoDb and inSilicoMerging R/Bioconductor packages. BMC Bioinforma. 13, 335. doi:10.1186/1471-2105-13-335
Tan, Y., Chen, Q., Li, X., Zeng, Z., Xiong, W., Li, G., et al. (2021). Pyroptosis: A new paradigm of cell death for fighting against cancer. J. Exp. Clin. Cancer Res. 40 (1), 153. doi:10.1186/s13046-021-01959-x
Thorsson, V., Gibbs, D. L., Brown, S. D., Wolf, D., Bortone, D. S., Ou Yang, T. H., et al. (2018). The immune landscape of cancer. Immunity 48 (4), 812–830.e14. doi:10.1016/j.immuni.2018.03.023
Tibshirani, R. (1997). The lasso method for variable selection in the Cox model. Stat. Med. 16 (4), 385–395. doi:10.1002/(sici)1097-0258(19970228)16:4<385::aid-sim380>3.0.co;2-3
Tsvetkov, P., Coy, S., Petrova, B., Dreishpoon, M., Verma, A., Abdusamad, M., et al. (2022). Copper induces cell death by targeting lipoylated TCA cycle proteins. Science 375 (6586), 1254–1261. doi:10.1126/science.abf0529
Vahl, J. M., Friedrich, J., Mittler, S., Trump, S., Heim, L., Kachler, K., et al. (2017). Interleukin-10-regulated tumour tolerance in non-small cell lung cancer. Br. J. Cancer 117 (11), 1644–1655. doi:10.1038/bjc.2017.336
Vanneman, M., and Dranoff, G. (2012). Combining immunotherapy and targeted therapies in cancer treatment. Nat. Rev. Cancer 12 (4), 237–251. doi:10.1038/nrc3237
Waldman, A. D., Fritz, J. M., and Lenardo, M. J. (2020). A guide to cancer immunotherapy: from T cell basic science to clinical practice. Nat. Rev. Immunol. 20 (11), 651–668. doi:10.1038/s41577-020-0306-5
Wang, Y., Fu, M., Liu, J., Yang, Y., Yu, Y., Li, J., et al. (2019). Inhibition of tumor metastasis by targeted daunorubicin and dioscin codelivery liposomes modified with PFV for the treatment of non-small-cell lung cancer. Int. J. Nanomedicine 14, 4071–4090. doi:10.2147/IJN.S194304
Wu, X., Kong, W., Qi, X., Wang, S., Chen, Y., Zhao, Z., et al. (2019). Icariin induces apoptosis of human lung adenocarcinoma cells by activating the mitochondrial apoptotic pathway. Life Sci. 239, 116879. doi:10.1016/j.lfs.2019.116879
Xu, M. L., Liu, T. C., Dong, F. X., Meng, L. X., Ling, A. X., and Liu, S. (2021). Exosomal lncRNA LINC01711 facilitates metastasis of esophageal squamous cell carcinoma via the miR-326/FSCN1 axis. Aging (Albany NY) 13 (15), 19776–19788. doi:10.18632/aging.203389
Yang, C., Guo, Y., Wu, Z., Huang, J., and Xiang, B. (2022). Comprehensive analysis of cuproptosis-related genes in prognosis and immune infiltration of hepatocellular carcinoma based on bulk and single-cell RNA sequencing data. Cancers (Basel). 14 (22), 5713. doi:10.3390/cancers14225713
Yang, C., Huang, X., Li, Y., Chen, J., Lv, Y., Dai, S., et al. (2021). Prognosis and personalized treatment prediction in TP53-mutant hepatocellular carcinoma: An in silico strategy towards precision oncology. Brief. Bioinform 22 (3), bbaa295. doi:10.1093/bib/bbaa295
Yoshihara, K., Shahmoradgoli, M., Martinez, E., Vegesna, R., Kim, H., Torres-Garcia, W., et al. (2013). Inferring tumour purity and stromal and immune cell admixture from expression data. Nat. Commun. 4, 2612. doi:10.1038/ncomms3612
Yu, A. T., Berasain, C., Bhatia, S., Rivera, K., Liu, B., Rigo, F., et al. (2021). PHAROH lncRNA regulates Myc translation in hepatocellular carcinoma via sequestering TIAR. Elife 10, e68263. doi:10.7554/eLife.68263
Zhang, A., Yang, J., Ma, C., Li, F., and Luo, H. (2021). Development and validation of a robust ferroptosis-related prognostic signature in lung adenocarcinoma. Front. Cell Dev. Biol. 9, 616271. doi:10.3389/fcell.2021.616271
Zhang, G., Sun, J., and Zhang, X. (2022). A novel Cuproptosis-related LncRNA signature to predict prognosis in hepatocellular carcinoma. Sci. Rep. 12 (1), 11325. doi:10.1038/s41598-022-15251-1
Zhang, T., Wang, Y., Inuzuka, H., and Wei, W. (2022). Necroptosis pathways in tumorigenesis. Semin. Cancer Biol. 86 (3), 32–40. doi:10.1016/j.semcancer.2022.07.007
Zhang, Y., Li, W., Lin, Z., Hu, J., Wang, J., Ren, Y., et al. (2020). The long noncoding RNA Linc01833 enhances lung adenocarcinoma progression via MiR-519e-3p/S100A4 Axis. Cancer Manag. Res. 12, 11157–11167. doi:10.2147/CMAR.S279623
Zhang, Z., Zeng, X., Wu, Y., Liu, Y., Zhang, X., and Song, Z. (2022). Cuproptosis-related risk score predicts prognosis and characterizes the tumor microenvironment in hepatocellular carcinoma. Front. Immunol. 13, 925618. doi:10.3389/fimmu.2022.925618
Zheng, X., Li, Y., Ma, C., Zhang, J., Zhang, Y., Fu, Z., et al. (2020). Independent prognostic potential of GNPNAT1 in lung adenocarcinoma. Biomed. Res. Int. 2020, 8851437. doi:10.1155/2020/8851437
Zhou, N., and Bao, J. (2020). FerrDb: A manually curated resource for regulators and markers of ferroptosis and ferroptosis-disease associations. Database 2020, baaa021. doi:10.1093/database/baaa021
Keywords: lung adenocarcinoma, lncRNA, signature, cuproptosis, targets, therapeutic agent, prognosis
Citation: Ma C, Li F, He Z, Zhao S, Yang Y and Gu Z (2023) Prognosis and personalized treatment prediction in lung adenocarcinoma: An in silico and in vitro strategy adopting cuproptosis related lncRNA towards precision oncology. Front. Pharmacol. 14:1113808. doi: 10.3389/fphar.2023.1113808
Received: 01 December 2022; Accepted: 06 February 2023;
Published: 15 February 2023.
Edited by:
Yuhua Yao, Hainan Normal University, ChinaReviewed by:
Run-Ze Li, Guangdong Provincial Hospital of Chinese Medicine, ChinaLi Shanbao, Shanghai Jiao Tong University School of Medicine, China
Copyright © 2023 Ma, Li, He, Zhao, Yang and Gu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Chao Ma, Y2hhby5tYUBhbHVtbmkuY2hhcml0ZS5kZQ==; Zhuoyu Gu, ZmNjZ3V6eUB6enUuZWR1LmNu
†ORCID: Chao Ma, http://orcid.org/0000-0003-1444-4668