Disulfidptosis-related long non-coding RNAs predict prognosis and indicate therapeutic response in non-small cell lung carcinoma

Data source and processing

A total of 1,096 patients diagnosed with NSCLC were enrolled in this study. The transcriptome profile and somatic mutation datasets were obtained from the Genomic Data Commons data portal (https://gdc.cancer.gov/access-data/gdc-data-portal), along with the clinical information, including age, sex, survival status, follow-up time, smoking history, pathological stage, and TNM stage [16]. The annotations for lncRNAs were retrieved from the National Center for Biotechnology Information (NCBI, https://www.ncbi.nlm.nih.gov/) [17]. Genes known to be markers for disulfidptosis (CAPZB, DSTN, ACTB, MYL6, TLN1, FLNB, FLNA, IQGAP1, MYH10, MYH9, ACTN4, PDLIM1, CD2AP, INF2, SLC7A11, NOX4, NOXA1, and ENOX2) were acquired from previous literature [7]. The overview for this study is depicted in Figure 1.

Figure 1:

Overview of this study.

Screening disulfidptosis-related lncRNAs

The co-expression analysis between and disulfidptosis marker genes and lncRNAs was conducted using Pearson correlation to identify DRLs. The threshold was established where the correlation coefficient’s absolute value >0.4 and p-value <0.01. The Sankey plot based on these findings was visualized using the “ggalluvial” R package (version 0.12.5) [18].

Construction of DRLs prognostic signature (DRLPS)

The R package “caret” (version 6.0-94) was employed to randomly separate NSCLC samples into the training and validation cohorts. The prognostic significance of DRLs was assessed through a univariate Cox proportional hazards regression analysis, with a significance level set at p <0.05. Afterward, the LASSO algorithm was executed using the “glmnet” package (version 4.1-8) in the R software. A penalized function was developed, incorporating 10-fold cross-validation, to identify the optimal prognostic DRLs. Subsequently, the DRLs prognostic signature was developed using a multiple stepwise Cox regression analysis [19]. Meanwhile, the scoring system for disulfidptosis was established based on this prognostic signature, and the formula was as follows:

where βi represents the hazard ratio coefficient of Cox regression analysis and EXP (lncRNAi) means the expression of each lncRNA of DRLPS. Consequently, the patient was labelled with a DS.

Evaluation of prognostic ability

NSCLC patients were categorized into high- and low-DS groups based on the median DS. The survival analysis was conducted to assess the prognostic ability of DRLPS in NSCLC. ROC analysis using R package “survivalROC” (version 1.0.3.1) was employed to evaluate the predictive accuracy of overall survival (OS) at 1-, 3-, and 5-years. Validating the reliability of the prognostic signature was done using PCA. The “scatterplot3D” R package (version 0.3-44) was applied for visualization [20].

The independent prognostic performance of DRLPS was assessed using univariate and multivariate Cox regression methods among clinical features. A nomogram was constructed, incorporating the DS, age, sex, smoking history, as well as the TNM and pathological stages. The calibration curve was used to detect the accuracy of the nomogram.

Functional enrichment assessment

The gene set enrichment analysis (GSEA) was used to assess the differences in biological function between the low- and high-DS cohorts. Gene ontology (GO) and REACTOME pathways were involved in the analysis using R package “clusterProfiler” (version 4.6.2) [21]. The analysis of GO terms reveals the enriched functions of gene clusters in relation to the biological processes. The REACTOME analysis uncovers the potential signaling pathways that the genes interacted through. A p-value below 0.05 was deemed to represent statistical significance.

Association between DRLPS and immune microenvironment

ESTIMATE (https://bioinformatics.mdanderson.org/estimate/) is a commonly employed tool for evaluating the levels of stromal and immune cells in tumor tissues. The potential impact of DRLs on the tumor microenvironment was investigated by determining the stromal and immune score using the transcriptome profile of NSCLC [22]. The R package “immunodeconv” (version 2.1.0) was applied to calculate immune cell abundance between low- and high-DS groups [23]. The immune-related functions were examined using the “ssGSEA” algorithm (version 1.46.0) [24].

Prediction of immunotherapy and drug sensitivity

A practical web application known as Tumor Immune Dysfunction and Exclusion (TIDE, http://tide.dfci.harvard.edu/) has been developed to assess the efficacy of immunotherapy for individuals [25]. The TIDE score reflects the effectiveness of treatment with anti-PD1 and anti-CTLA4 drugs. The higher the TIDE scores, the greater the likelihood that cancer patients would benefit from immunotherapy. To evaluate the potential utility of DRLPS for NSCLC treatment, the response of the patients in the low- and high-DS groups to the commonly used medicines in NSCLC therapy was assessed. The half-maximal inhibitory concentration (IC₅₀) for these drugs was calculated using the R package “oncoPredict” (version 0.2) [26].

Genomic alteration evaluation

The R package ‘maftools’ (version 2.14.0) was employed for the visualization of the landscape of genomic variations and calculation of tumor mutation burden (TMB) for NSCLC patients in high- and low-DS groups [27]. The joint analysis of DS and TMB was performed using Kaplan–Meier survival analysis to further evaluate the prognostic impact on the NSCLC patients.

RT-PCR assay

The expression levels of six prognostic DRLs in lung cancer cells were determined using Reverse Transcription-Polymerase Chain Reaction (RT-PCR). The primer sequences are provided in Table S1. The TRizol reagent (TIANGEN, Beijing, China) was utilized to extract total RNA from A549, NCI-H1975, and BEAS-2B cells. Subsequently, the obtained RNA was reverse transcribed into cDNA according to the manufacturer’s guidelines (Thermo Scientific, Waltham, USA). For PCR amplification, a reaction mixture (20 μL) was prepared with 10 μL of 2× PerfectStart Green qPCR SuperMix (TransGen, AQ601-04, Beijing, China), 2 μL of cDNA diluted in nuclease-free water, and 1 μL of a 10 μM primer with nuclease-free water. ABI Q1 sequencers were used for amplification and detection procedures. GAPDH was utilized as the reference gene. Quantifying the transcript was conducted using the 2^−ΔΔCT method.

Statistical analysis

The R software (version 4.2.3) was used for performing the statistical analyses. The distribution of samples between the training and validation cohorts was assessed by the Chi-square test. Pearson correlation coefficient with a cutoff of |Cor| >0.4 and p-value <0.01 was used to identify DRLs. The quantitative comparison between the low- and high-DS groups was conducted using the Wilcoxon rank-sum test. The Kaplan–Meier survival analysis was implemented to estimate the predictive ability for NSCLC prognosis, followed by the log-rank p-value. The accuracy of prognostic prediction was assessed through ROC analysis. Univariate and multivariate regression analyses were performed using Cox proportional hazard models. Statistical significance was considered when the p-value fell below 0.05.

Identifying disulfidptosis-related lncRNAs

The boxplot showed a significant differential expression of all disulfidptosis marker genes between the NSCLC tumor and adjacent normal tissue (Figure 2A). The cumulative rate at which these genes were altered reached 23.28 % (Figure 2B), and this highlights the critical role that disulfidptosis plays in tumorigenesis and cancer development. The Pearson correlation analysis revealed a strong association between 232 lncRNAs and 16 markers that are related to disulfidptosis. The relationship between these candidate DRLs and disulfidptosis marker genes was illustrated in the Sankey diagram (Figure 2C).

Figure 2:

Identification of disulfidptosis-related lnRNAs in NSCLC. (A) The expression of disulfidoptosis marker genes in the tumor and adjacent normal tissue of NSCLC. (B) Mutational status of disulfidptosis marker genes of NSCLC. (C) Sankey diagram between 16 disulfidptosis marker genes and 232 disulfidptosis-associated lncRNAs. (D) The correlation between disulfidptosis marker genes and six lncRNAs of DRLPS.

Establishing the DRLs prognostic signature

A joint approach was employed in this analysis. The univariate Cox analysis screened 21 lncRNAs based on the extent to which they were associated with disulfidptosis, as well as their correlation with prognosis in NSCLC (Figure 3A). Subsequently, the LASSO regression analysis identified 14 optimal features among these lncRNAs (Figure 3B and C). Finally, a prognostic signature that consisted of six DRLs (AL606489.1, LINC00857, AP003555.1, AP000695.1, AC113346.1, and LINC01615) was established using the multiple stepwise Cox regression analysis. The survival analysis revealed a significant association between each lncRNA of DRLPS and the prognosis of NSCLC (Figure S1A–F). The RT-PCR assay demonstrated an obvious up-regulation of AL606489.1, LINC00857, AP000695.1, and AC113346.1 in NSCLC cells, while the downregulation of LINC01615 was observed in NCI-H1975 cells (Figure S2). Based on the heatmap of Pearson correlation coefficient, correlation was noted between disulfidptosis marker genes and these prognostic lncRNAs (Figure 2D). DS was calculated for each NSCLC patient based on the regression coefficient and expression of lncRNAs. DS=0.6578*EXP_AL606489.1+0.8820 *EXP_LINC00857+0.5217*EXP_AP003555.1+0.8724*EXP_AP000695.1+0.1019 *EXP_AC113346.1−0.4471*EXP_LINC01615.

Figure 3:

Construction of disulfidptosis-related lncRNAs prognostic signature in NSCLC. (A) Identification of prognostic DRLs using univariate Cox regression analysis. (B, C) The LASSO algorithm was utilized to determine the optimal features among 21 prognostic DRLs. (D) Kaplan–Meier survival curve of low- and high-DS cohorts in the training set. (E) The survival statistics distribution with disulfidptosis scores in the training set. (F) Heatmap of prognostic six lncRNAs of DRLPS in the training set. (G) Kaplan–Meier survival curve of low- and high-DS groups in the validation set. (H) The survival statistics distribution with disulfidptosis scores in the validation set. (I) heatmap of DRLs of six lncRNAs of DRLPS in the validation set.

There were no significant differences (p>0.05, Table 1) observed in the clinical characteristics between the training and validation cohorts. This ensured the confidence of the classified groups in the subsequent analysis. Patients in the training cohort were categorized into low- and high-DS groups based on the median DS. The prognostic performance of DRLPS in NSCLC was evaluated through survival analysis. The survival curves indicated that the patients with high-DS were remarkably associated with poor overall survival (Figure 3D). The survival statistics distribution with disulfidptosis scores is illustrated in Figure 3E, which indicates a rise in deaths with increasing DS. Moreover, the expression level of DRLPS lncRNAs was obviously elevated in high-DS group (Figure 3F). These findings were consistently observed across the validation and entire cohorts (Figures 3G–I and 4A–C), thus confirming the accuracy of the established prognostic signature. Compared to disulfidptosis-related genes and lncRNAs, the prognostic DRLs exhibited superior efficacy in distinguishing between high- and low-DS populations (Figure 4D–G). This further validated the reliability of DRLPS.

Figure 4:

Validation of the prognostic signature using the entire cohort and PCA. (A) Kaplan–Meier survival curve of the low- and high-DS groups in the entire cohort. (B) The survival statistics distribution with disulfidptosis scores in the entire set. (C) Heatmap of the DRLPS lncRNAs in the entire set. Scatter plot of PCA between low- and high-DS subgroups based on (D) all genes, (E) disulfidptosis-associated genes, (F) DRLs, and (G) prognostic DRLs.

Findings from the analysis of the association between DS and clinical factors revealed that DS significantly affects age, pathological stage, histopathology, as well as stages T and M (Table 2).

Disulfidptosis score serves as an independent risk factor

The ROC curves demonstrated the accuracy of DRLPS in predicting the overall survival at 1-, 3-, and 5-year periods, with the AUC values of 0.679, 0.711, and 0.740, respectively (Figure 5A). Meanwhile, the DS demonstrated a superior prognostic value than other clinical characteristics (Figure 5B). To further evaluate the prognostic significance of DRLPS in comparison with clinical characteristics, both univariate and multivariate Cox regression analyses were conducted. The univariate Cox analysis demonstrated a significant correlation between the DS (HR=1.36, Cox p=0.0018), pathological stage (HR=1.91, Cox p<0.001), N stage (HR=1.67, Cox p<0.001), T stage (HR=1.60, Cox p<0.001), age (HR=1.24, Cox p=0.029) and the prognosis of NSCLC (Figure 5C). The findings from the multivariate Cox regression approach were as follows: DS (HR=1.32, Cox p=0.0049), pathological stage (HR=1.59, Cox p=0.0025), T stage (HR=1.60, Cox p=0.041) and age (HR=1.32, Cox p=0.0055). These results also showed the predictive ability of NSCLC prognosis (Figure 5D). These findings revealed the significance of high DS as an independent risk indicator for unfavorable outcomes. Additionally, the nomogram (Figure 5E) was employed to quantify the prognostic performance of DRLPS based on the DS and clinical factors. The reliability of this nomogram was confirmed by the calibration plot (Figure 5F).

Figure 5:

Prognostic performance of DRLPS in NSCLC. (A) Time-dependent ROC curves depicting overall survival at 1, 3, and 5 years in NSCLC. (B) ROC curves of overall survival for disulfidptosis score and clinical factors in NSCLC. (C) Univariate Cox regression analysis of disulfidptosis score and clinical characteristics in NSCLC. (D) Multivariate Cox regression analysis of disulfidptosis score and clinical characteristics in NSCLC. (E) establishment of a nomogram utilizing disulfidptosis score and clinical features. (F) The calibration curves of the nomogram.

DRLPS is related to the immune mediations and metabolism in NSCLC

The GO analysis revealed that the gene cluster of the low-DS group was primarily involved in the biological processes that are relevant for activating immune responses and T-cells, leukocyte-mediated immunity, as well as the regulation of lymphocytes. On the contrary, genes with high DS were found to be involved in the functions such as oxidative phosphorylation, nuclear chromosome segregation, mitochondrial respiratory chain complex assembly, meiotic cell cycle, and cytokinetic process (Figure 6A). With respect to the REACTOME pathways, T cell receptor (TCR) signaling, programmed death-1 receptor (PD-1) signaling, CTLA4 inhibitory signaling, interleukin-10 signaling, and signaling by interleukin were the prominent pathway terms in the low-DS cohort. It’s important to note that the enrichment of several pathways linked to metabolism was noted in the high-DS group (Figure 6B). These include selenoamino acid metabolism, O-like glycosylation of mucins, mitochondrial translation, respiratory electron transport, the citric acid TCA cycle, and termination of O-glycan biosynthesis. These results unveiled the correlation between DRLs and immune regulation as well as metabolism in the microenvironment of NSCLC. Therefore, it can be extrapolated that these DRLs may serve as the therapeutic targets for NSCLC treatment.

Figure 6:

Functional enrichment analysis of DRLPS. (A) Enriched GO terms of biological process in the high- and low-DS groups. (B) Enriched REACTOME pathways in the high- and low-DS groups.

A significant association between DRLs and immune infiltration

Given the potential role of DRLs function in NSCLC immunity, the ESTIMATE and ssGSEA were performed to evaluate the intra-tumoral immune cell infiltration between the high- and low-DS groups. The ESTIMATE analysis revealed higher scores of immune and stromal cells in the low-DS group compared to the high-DS group (Figure 7A). ssGSEA results uncovered the immune cell populations of NSCLC patients (Figure 7B). A remarkable difference was observed in the two NSCLC subgroups with regard to naive B-cells, activated natural killer (NK) cells, activated/resting mast cells, macrophages M0/M1, activated dendritic cells, and various functional T-cells (Figure 7C). These findings demonstrated the significant correlation of DRLs with immune cell infiltration. Additionally, patients in the low-DS group exhibited heightened immune function activity, as evidenced by increased MHC class I level, type II interferon (IFN) response, cytolytic activity, inflammation promotion, CCR signaling, and co-stimulation/suppression of antigen-presenting cells (APCs) and T-cells (Figure 7D). This may contribute to anti-tumor potency in the NSCLC microenvironment.

Figure 7:

Exploration of the association between DRLPS and immune microenvironment of NSCLC. (A) Comparison of stromal and immune infiltration levels between the high- and low-DS groups. (B) Immune cell population of NSCLC patients in the high- and low-DS groups. (C) Comparison of tumor-infiltrated immune cells between the high- and low-DS groups. (D) Heatmap of immune-related activity in the high- and low-DS groups of NSCLC patients. (E) Expression difference of the immune checkpoint markers between the high- and low-DS groups. (F) Comparison of TIDE score in the high- and low-DS cohorts. *p<0.05, **p<0.01, and ***p<0.001. ns, not significant.

Promising role of DRLPS in predicting immunotherapy and drug sensitivity

Comparisons on the expression levels of immune checkpoint genes in groups with low- and high-DS were done (Figure 7E). Within the cohort with a low-DS, the expression levels of 24 genes were significantly higher compared to those in the group with a high DS. Notably, PDCD1, CD274, and CTLA4, which are commonly used in clinical treatments, were among these 24 genes. To investigate the potential predictive significance of DRLPS in immunotherapy, the TIDE score was used to evaluate the impact of immunosuppressive factors on patients with NSCLC. The higher TIDE score observed in the group with a lower DS suggests that these patients exhibited heightened responsiveness to immunotherapy that involved anti-PD1 and anti-CTLA4 drugs (Figure 7F), suggesting the potential role of DRLPS in predicting the immunotherapeutic efficacy. Afterward, the relationship between the DS and drug sensitivity for the NSCLC treatment was also explored. The IC₅₀s of 16 medicines, including bleomycin, doxorubicin, etoposide, gemcitabine, methotrexate, temozolomide, alectinib, dabrafenib, trametinib, cabozantinib, afatinib, panatinib, AZD6244, MS-275, CGP43231, and ruxolitinib exhibited significant negative correlation with the DS (Figure S3). Patients with a high DS exhibited enhanced drug susceptibility, as evidenced by lower IC₅₀ values, in contrast to patients with low DS (Figure 8). The findings demonstrated the possibility of using DRLPS in guiding the clinical decisions of these medicines for NSCLC patients.

Figure 8:

The significant association of drug sensitivity (IC₅₀) with disulfidptosis score in NSCLC treatment. Notes: ** p value<0.01; *** p value <0.001.

Tumor mutation burden and genomic variation landscape

Increasing studies pointed out that genomic alterations are closely related to patients’ sensitivity to tumor treatment. Patients with a high tumor mutation burden may contain large amounts of neoantigens, which activate the immune system to target cancerous cells, thereby promoting apoptosis. In this study, the TMB levels and genomic variation landscape were involved in the analyses to investigate the relationship between DRLs and genomic alterations in NSCLC. A higher TMB level was noted among patients exhibiting low DS (p=0.029, Figure 9A). Apparently, higher TMB levels significantly correlated with an increased survival probability for NSCLC patients (p=0.02, Figure 9B). Additionally, the joint survival analysis combined with the DS and TMB level was also performed and the result revealed that these two indicators had a strong synergistic effect in predicting NSCLC prognosis (p=0.00071, Figure 9C). The waterfall plot of genetic variation displayed the mutation landscape for the top 20 altered-frequency genes in the low- and high-DS cohorts (Figure 9D and E). A higher alteration frequency was found in the group with a low DS, and TP53 exhibited the highest proportion of mutations in both subgroups of NSCLC.

Figure 9:

The correlation between the disulfidptosis score and tumor mutation burden in NSCLC. (A) Comparison of TMB between the low- and high-DS groups in NSCLC. (B) Survival curve of the low- and high-TMB groups. (C) Survival analysis combined with DS and TMB level for NSCLC patients. (D, E) Genomic alteration landscape of the top 20 mutated genes in the low- and high-DS cohorts.