A genetic signature linked to disulfidptosis for predicting lung adenocarcinoma outcomes and immune environment

The transcriptomic data, somatic mutation data, and clinical information of LUAD were obtained from the Cancer Genome Atlas (TCGA) database (http://portal.gdc.cancer.gov/↗) and then under pretreatment for subsequent analyzes. The copy number variation (CNV) information was downloaded from Xena database (http://xena.ucsc.edu/↗)¹⁴. Ultimately, a cohort consisted of 507 patients was included in our investigation and was randomly divided into the training set (355 patients) and testing set (152 patients) at a ratio of 7:3 using the R package “caret”. The comparative analysis of clinical information between these two sets was performed utilizing the chi-square test methodology.

According to present research⁸, 16 genes which implicated in disulfidptosis pathways were selected and denoted as disulfidptosis-related genes (DRGs) (Table S1). Utilizing a lncRNA annotation data from the GENCODE database (http://www.gencodegenes.org/↗), a comprehensive cohort of 16,876 lncRNAs linked to disulfidptosis processes was extracted by the use of the Perl programming language. Subsequently, a comprehensive set of 104 lncRNAs associated with disulfidoptosis was identified and validated through Pearson’s correlation analysis (|Pearson R| > 0.4, P < 0.001).

The training set served as the foundational basis for constructing the DRLS, while the testing and entire sets were employed to assess and validate the prediction accuracy of the DRLS. Prognostic lncRNAs were screened and examined through a series of analytical methodologies, encompassing univariate and multivariate Cox regression analysis¹⁵, LASSO regression analysis¹⁶, by utilizing the “glmnet” package. Additionally, the DRLS score for each LUAD patient was computed via the formula as follow:

Coef and Expr represent the coefficient and the DRlncRNAs expression respectively.

Patients in the training, testing as well as entire sets were respectively categorized to the low and high risk groups based on the DRLS score by using the “survminer” package. Subsequently, the Kaplan–Meier (K–M) analysis and log-rank test were employed for the assessment of disparities in overall survival (OS) among these two groups, along with using the “survival” package for statistical analyses. Receiver operating characteristic (ROC) was applied for the accuracy evaluation of the developed signature in predicting OS of LUAD patients using the package “timeROC”. Principal Component Analysis (PCA)¹⁷ and t-distributed Stochastic Neighbor Embedding (t-SNE) were applied for further facilitating dimensionality reduction and visualize the distinctions between these two risk groups. Furthermore, univariate and multivariate Cox regression analyses were employed to evaluate the independent prognostic significance of the DRLS score.

Nomograms were constructed through the combination of the risk scores and clinical variables including gender, age, and stage, with the aim of predicting the 1-, 3-, and 5-year OS of LUAD patients. The nomograms were developed using the “RMS” package. Furthermore, the Hosmer–Lemeshow test was applied to generate calibration plots for depicting the assistance between actual outcomes and predictive results. Additionally, we generated Decision Curve Analysis (DCA) curves and ROC curves for further validation of the nomogram.

The “LIMMA” R package was employed to discern the differentially expressed genes (DEGs) in the high and low DRLS groups (|Log2FC| > 2.0, P < 0.05). Subsequently, functional analyses encompassing Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG)^18–20 were carried out utilizing the ‘ClusterProfiler’, ‘org.Hs.eg.db’, ‘enrichplot’, ‘ggplot2’, ‘RColorBrewer’, ‘dplyr’, ‘ggpubr’, and ‘ComplexHeatmap’ packages. The “GSVA” R package was used to conduct Gene Set Variation Analysis (GSVA) for exploring potential variations in biological functions. The gene set files utilized for these analyses were acquired from the Molecular Signatures Database or previously published research. Additionally, the “c2.cp.kegg.Hs.symbols.gmt” gene set file was employed as a reference dataset for performing Gene Set Enrichment Analysis (GSEA).

The “maftools” package was employed to quantify tumor mutation burden (TMB). Following this, all LUAD patients were stratified into two groups based on distinct TMB levels determined by the median TMB score. Additionally, the association between DRLS and TMB was investigated using Spearman correlation analysis.

Firstly, the “ESTIMATE” R package was applied for calculation of the immune scores, stromal scores, and ESTIMATE scores of each patient to evaluate the population of immune cells and stromal cells between distinct groups. And then, the relationship between the three scores and DRLS were further evaluated using Spearman correlation analysis. Subsequently, 8 distinct algorithms, including TIMER, xCELL, CIBERSORT, CIBERSORT-ABS, quanTIseq, MCP-counter, EPIC, and ssGSEA were employed to access the infiltration of immune cells^21–28. Finally, the comparative analysis focused on the immune checkpoints (ICs) expression was conducted to explore the potential utility of DRLS in predicting responses to immunotherapy.

The tumor immune dysfunction and exclusion (TIDE) algorithm was accessed and utilized online (http://tide.dfci.harvard.edu/↗)²⁹. The lower TIDE scores tend to associated with more favorable responses toward Immune Checkpoint Inhibitors (ICIs) therapy. Immunophenoscores (IPS) for LUAD were acquired from the Cancer Immunome Atlas (TCIA) database (http://tcia.at/home↗)³⁰. The Cancer Drug Sensitivity Genomics of Cancer Cell Lines (GDSC) database (https://www.ancerrxgene.org/↗) were applied for assessment of potential clinical application of DRLS, along with the calculation of half-maximal inhibitory concentration (IC50) for commonly used anti-tumor drugs using the R package “pRRophetic”. The predictive results were estimated applying ridge regression analysis, as for the prediction precision was measured via tenfold cross-validation.

All patients diagnosed with LUAD were stratified into two clusters based on the expression levels of DRlncRNAs, employing the “ConsensusClusterPlus” R package. This classification aimed to investigate potential molecular subtypes. Following this, distinctions between the two clusters regarding survival rates, time-related parameters, immune infiltration, and response to immune therapy were assessed using the methodologies previously outlined.

The human bronchial epithelial cell line BEAS-2B and two human lung adenocarcinoma cell lines (A549 and H1975) were procured from the Institute of Biochemistry and Cell Biology of the Chinese Academy of Sciences. All the cells were cultured in RPMI-1640 medium (Gibco; Thermo Fisher Scientific, Inc.) supplemented with 10% fetal bovine serum (FBS) (Gibco; Thermo Fisher Scientific, Inc.) and 1% penicillin–streptomycin. The culture conditions were maintained at a temperature of 37 °C, a CO₂ concentration of 5%, and a humidity level of 95%.

Short interfering RNAs (siRNAs) were synthesized by KeyGEN (Nanjing, China). Transfection procedures were performed employing Lipofectamine 2000 (Invitrogen, USA), following the manufacturer’s instructions. The assessment of transfection efficiency was conducted 48 h after transfection, utilizing quantitative real-time polymerase chain reaction (qRT-PCR). The target sequences are as follows:

Total RNA was extracted from the cells using TRIzol® reagent (Invitrogen; Thermo Fisher Scientific, Inc.) in accordance with the manufacturer’s instructions. Subsequently, reverse transcription was carried out utilizing a PrimeScript Reverse Transcription Kit (Takara Bio, Inc.). For RT-qPCR, an StepOnePlus™ Real-Time PCR System (Applied Biosystems; Thermo Fisher Scientific, Inc.) was employed. The conditions for the RT-qPCR reaction involved an initial step at 95 °C for 5 min, following 40 cycles comprising denaturation at 95 °C for 15 s, annealing at 60 °C for 30 s, and extension at 72 °C for 45 s. The 2^−ΔΔCt computing method was employed to calculate the relative mRNA expression, which was normalized using β-actin as a reference. All reactions were conducted in triplicate. The gene primers applied within the study are as follows:

The CCK-8 assay was conducted employing CCK-8 solution (Dojindo Laboratories, Inc.), in strict accordance with the manufacturer’s instructions. Cells were seeded into 96-well plates at a density of 2 × 10³ cells/100 µl per well. Following incubation periods of 0, 24, 48, 72, or 96 h, 10 µl CCK-8 per well solution was added to the 96-well plates to measure the relative number of viable cells, with the absorbance at 450 nm was evaluated by a microplate reader.

In the colony formation assay, LUAD cells were seeded in 6-well plates at a density of 1 × 10³ cells/well. Subsequently, these plates were incubated at 37 °C with a 5% CO₂ for a duration of 14 days. Following the incubation period, colonies were fixed using 4% paraformaldehyde and stained with 0.1% crystal violet obtained from Beyotime Institute of Biotechnology. The population of colonies in each well was counted, with each colony comprising more than 50 cells. The experiments were performed in triplicate.

LUAD cells were seeded in 6-well plates when the concentration reached 90%. A linear scratch wound was generated within the monolayers of cells using a 20 µl pipette tip. Subsequently, cells were cultured within an FBS-free medium and followed by carefully washing with PBS to remove debris and suspended cells. The images of each wound were recorded at 0 and 36 h by the use of an inverted microscope.

Transwell chambers (8-µm pore size; Corning, Inc.) were applied for the migratory capabilities evaluation of LUAD cells. In each upper Transwell chamber, 4 × 10⁴ cells were seeded with 300 µl of serum-free medium, while the lower chamber was filled with 700 µl of RPMI-1640 medium supplemented with 10% FBS. After a 2-day incubation period, cells were fixed using 4% paraformaldehyde and subsequently stained with 0.1% crystal violet obtained from Beyotime Institute of Biotechnology. Non-migratory cells residing on the upper insert chamber membrane were carefully eliminated using a cotton swab. Following this, images of the stained cells were acquired using an inverted microscope. Cell counts were then performed in five randomly selected fields for the evaluation of the LUAD cells migratory potential. All experiments were replicated in triplicate.

Statistical analyses were conducted using the R, Perl, and GraphPad Prism 8 platforms. Group differences were assessed using Student’s t-test and the Wilcoxon rank sum test. The association between two variables was evaluated using Pearson’s and Spearman’s correlation methodologies. A significance threshold of P < 0.05 was considered indicative of statistical significance in this study.

According to the research focused on disulfidptosis, there are 16 key genes which take an important part in the process of disulfidptosis, therefore we selected these 16 disulfidptosis-related genes (DRGs) for further study, Fig. 1 depicted the flowchart of this research.

Firstly, we studied the expression levels of 16 DRGs in 58 normal samples and 524 LUAD tumor samples from TCGA database (Fig. 2A). It was observed that, except for FLNB, the remaining 15 DRGs displayed differential expression between tumor and normal tissues. Subsequently, protein–protein interaction (PPI) networks for the 16 DRGs were constructed applying the String website and Cytoscape software. And the findings revealed tight interactions among these 16 DRGs (Fig. 2B). In addition, we also investigated the frequency of copy number variation (CNV) mutations, and found that all of the 16 DRGs contained a large number of CNV mutations. The decreased extensive CNVs were observed in OXSM, FLNB, SLC7A11, NDUFA11, GYS1, while the overall enhancement of CNVs were presented in ACTB, TLN1, NUBPL, FLNA (Fig. 2C). We illustrated the CNV alteration locations on the chromosomes for these 16 DRGs (Fig. 2D). Next, we extracted 16,876 lncRNAs from the GENCODE database using Perl-based methods. Subsequently, 104 lncRNAs associated with disulfidptosis were identified by the use of Pearson’s correlation analysis (|Pearson R| > 0.4, P < 0.001, Table S2), and their co-expression network with DRGs was visualized in a Sankey’s plot (Fig. 2E). Furthermore, LUAD patients were involved in our study, which categorized into the training set and testing set. The training set was applied for identifying prognostically relevant LncRNAs and constructing DRLS, while the testing set served as the validation of the predictive efficacy of DRLS. There was no significantly statistic difference observed in the comparative analysis in age, stage, gender, and TNM between the two groups (P > 0.05, Table S3). Subsequently, 17 DRlncRNAs highly correlated with prognosis through univariate COX regression analysis were identified (Fig. 2F). Finally, we generated a prognostic correlation network diagram for these 17 DRlncRNAs (Fig. 2G), identifying all DRlncRNAs as risk factors, except for PTPRN2-AS1, AC009065.4↗, AC012615.1↗, and AC090559.1↗, which were considered protective factors.

According to the univariate COX regression analysis, 17 DRlncRNAs were identified and then applied in LASSO regression analysis. Subsequently, 12 DRlncRNAs were identified as significantly linked to the LUAD patients prognosis using the optimal lambda value (Fig. 3A,B). Finally, 5 DRlncRNAs were determined via multivariate COX regression analysis (Table S4). Among these, AC012615.1↗ was a prognostic protective factor, while AL365181.2↗, GSEC, AC093673.1↗, AL606834.1↗ were prognostic risk factors (Fig. 3C). In addition, we examined the correlation of the five DRlncRNAs with clinical data and DRGs (Fig. 3D,E). Our findings revealed a negative correlation between AC012615.1↗ and the majority of DRGs, whereas the remaining four DRlncRNAs exhibited a positive correlation with most DRGs (Fig. 3E). Following that, the disulfidptosis-related lncRNA-based signature (DRLS) score of each patient was computed applying the formula below: DRLS Score (Risk Score) = 0.152323368591861 × expr AL365181.2↗ + 0.228911037835446 × expr GSEC + 0.203628306081756 × expr AC093673.1↗ + (− 0.44972537859839) × expr AC012615.1↗ + 0.292888245903172 × expr AL606834.1↗.

The patient cohort was stratified into high-risk and low-risk groups according to the optimal cut-off value of DRLS score. Subsequently, K–M analysis was performed among the training set, and the results revealed a significant improvement in OS of LUAD patients with low DRLS score compared to those with high DRLS score (Fig. 3F, P < 0.001). Further investigation involved assessing the distribution of survival times, survival statuses, and risk scores between the two groups. Simultaneously, the expression profiles of five DRlncRNAs were compared within these two sets (Fig. 3G). Our results highlighted remarkably increased expression levels of AL365181.2↗, GSEC, AC093673.1↗, and AL606834.1↗ in the high-risk group, consistent with their respective positive correlation coefficients in the DRLS score calculation formula, indicating their role as risk factors. Additionally, in the time-based ROC curves, the area under the curve (AUC) values of 1, 3, and 5 years were 0.701, 0.668, and 0.618, respectively. As for the combined ROC curves, the AUC values considering clinical features were 0.668 (Risk), 0.658 (Stage), 0.523 (Gender), and 0.519 (Age) (Fig. 3H). Ultimately, the t-SNE and PCA were utilized to reduce dimensionality and visualize the distinct characteristics between the two groups. The results conclusively showed significant and consistent differences in the distributions of the high-DRLS and low-DRLS groups (Fig. 3I).

We further examined the reliability and predictive capacity of established DRLS by calculating DRLS scores (risk scores) for LUAD patients among the testing and entire sets utilizing the same computational formula. Similar to the training set, patients in the two sets were stratified into high and low risk groups. Subsequent K–M analyses were conducted among the two sets, and the findings consistently revealed the low-risk groups displayed markedly superior OS compared to the high-risk groups (Fig. 4A,B, P < 0.001). Additionally, the progression-free survival (PFS) was demonstrated notably enhanced in the groups with low DRLS score in comparison to the high-DRLS groups (P < 0.001) among the entire set. An evaluation for the distribution of survival time and survival status, based on DRLS scores, was also conducted in the two groups among both the testing and entire sets. The expression profiles of these 5 DRlncRNAs in the two groups displayed consistent outcomes with those observed in the training sets (Fig. 4D,E). Subsequently, ROC curves of the two sets reflected the same trends observed in the training sets. As for the AUC values at 1, 3, and 5 years in the testing set were 0.669, 0.681, and 0.671, respectively. While the corresponding AUC values were 0.690, 0.671, and 0.631 within the entire set (Fig. 4F). Furthermore, PCA and t-SNE were also applied in the two sets to visualize the sample distributions of high-DRLS and low-DRLS groups (Fig. 4G,H). These findings implied the DRLS based on 5 DRlncRNAs could consistently and accurately predict the outcomes for LUAD patients with a high level of dependability.

Considering the widespread clinical application of TNM staging and the reliable predictive ability of DRLS, we constructed a nomogram through combining TNM staging with DRLS in the training set. Subsequently, this nomogram was used to compute scores of each patient, enabling better prediction of their 1, 3, and 5-year prognoses (Fig. 5A). We also generated a forest plot based on the nomogram, which highlighted T stage, N stage, and risk score as the primary prognostic factors (P < 0.05, Fig. 5B). Next, we plotted calibration curves to further assess the predictive accuracy of this nomogram among the training and testing sets. These curves illustrate a concordance between the observed and predicted values (Fig. 5C,D). DCA was employed to assess the utility and effectiveness of this prediction model. The findings indicated the potential of this nomogram for accurate prediction of the LUAD patients’ survival probabilities at various time intervals. (Fig. 5E). Furthermore, time-dependent ROC curves were constructed using this nomogram among the training and testing sets. In the training set, the AUC values for 1, 3, and 5 years were 0.715, 0.7, and 0.705, respectively (Fig. 5F). While in the testing set, the corresponding AUC values were 0.831, 0.74, and 0.677 (Fig. 5G). In summary, this nomogram demonstrated the ability for accurate prognosis prediction and the potential in refining clinical treatment decisions for LUAD patients. Finally, patients were stratified into early-stage (I–II) and advanced-stage (III–IV) groups, and K–M analysis was applied for validating the disparities in survival of the high and low-risk groups among patients in each stage. We found that regardless of the stage, LUAD Patients categorized within the low-risk group exhibited notably superior survival rates in comparison to those in the high-risk group (P < 0.05, Fig. 5H,I).

An investigation employing 513 DEGs that differentiated between the high-risk and low-risk groups (|Log2FC| > 1.0, P-value < 0.05) were conducted to unravel the underlying molecular mechanisms responsible for the marked distinctions observed in our multifaceted analyses between the two groups (Table). The GO analysis unveiled significant enrichment of these DEGs in cellular components (e.g., cytoplasmic region and plasma membrane bounded cell projection cytoplasm) and molecular functions (e.g., oxidoreductase activity and serine-type endopeptidase inhibitor activity) (Supplementary Fig.A). Furthermore, the KEGG analysis revealed a predominant association of these DEGs with pathways pertinent to biological process of tumors, such as the PI3K-Akt signaling pathway, and pathways linked to metabolic processes, including arachidonic acid metabolism (Supplementary Fig.B). Subsequently, we explored the differences in pathways enriched with DEGs in the low- and high-risk groups through GSVA analysis. Notably, we observed significant enrichment of DEGs in the high-risk group across a spectrum of metabolic pathways, encompassing glycolysis gluconeogenesis, pentose phosphate pathway, and galactose metabolism. Furthermore, these DEGs exhibited enrichment in pathways associated with p53 signaling, cell cycle, and DNA replication, which could potentially promote tumor progression in the high-risk cohort (Supplementary Fig.C). S5 S1 S1 S1

Subsequently, correlation analysis was performed to examine the association between five DRlncRNAs employed in constructing the DRLS and the risk score with regard to tumor-related pathways. Our analysis revealed positive correlations between AC093673.1↗, AC606834.1, and the risk score with tumor-related pathways. Conversely, AC012615.1↗ exhibited negative associations with all tumor-related pathways, thereby reinforcing our earlier conclusion that AC012615.1↗ functions as a protective factor (Supplementary Fig. S1D). Furthermore, the results of GSEA revealed notable pathway enrichments in the high-DRLS group, including cell cycle, focal adhesion, metabolism of xenobiotics by cytochrome P450, the P53 signaling pathway, and steroid hormone biosynthesis (Supplementary Fig. S1E). On the contrary, the low-risk group demonstrated significant pathway enrichments in areas such as allograft rejection, asthma, and the immune network for IgA production (Supplementary Fig. S1F). In summary, the DEGs identified between the two groups with distinct risk score may exert a significant influence on the pathogenesis of LUAD.

Tumor mutational burden (TMB) was the measure used to assess mutation frequency within tumor cells, and we conducted an analysis of somatic mutations of the high and low risk groups. This analysis unveiled 15 genes with the highest mutation frequencies in both groups, including TP53, TTN, MUC16, CSMD3, RYR2, LRP1B, ZFHX4, USH2A, and KRAS (Fig. 6A,B). The specifics of these mutations are visually presented in Fig. 6C. The most prevalent mutation type was missense mutations, characterized by a greater frequency of single nucleotide polymorphism (SNP) compared to insertions (INS) and deletions (DEL). Unlike C > T and C > G transitions, C > A mutations emerged as more frequent. The top 10 mutated genes of LUAD were displayed in the horizontal bar graph. Additionally, a quantitative assessment of TMB indicated the patients with high-DRLS scores displayed notably elevated TMB scores in comparison to those with low-risk scores (Fig. 6D, P < 0.05). To investigated the impact of TMB on prognosis, whole LUAD patients were categorized into high and low TMB groups, with the median TMB score serving as the cutoff point. The results of K–M analysis revealed the patients with high TMB scores exhibited improved OS in comparison to those with low TMB scores (Fig. 6E, P < 0.05). Interestingly, our subsequent analysis, which sought to predict the outcomes of LUAD patients by combining TMB level and risk scores, revealed that patients with low-risk score and high TMB score exhibited markedly improved prognoses in comparison to patients in all other groups. This observation suggested the risk score surpassed TMB in predicting individual outcomes (Fig. 6F, P < 0.001). Finally, we examined the relation between the risk score, TMB, and immune cell infiltration, visually representing the results through a correlation chord diagram (Fig. 6G).

Over the past few decades, a large number of studies has highlighted the predominance of TME in the development and progression of various types of tumors. Immune cells, as essential components of the TME, engage in complex interactions and regulatory processes with tumor cells, subsequently influencing tumor progression and treatment outcomes. Therefore, the immune score, stromal score, and ESTIMATE score were computed for each patient employing the ESTIMATE algorithm, with the aim of exploring potential distinctions in immune cell infiltration within the high- and low-risk groups. The results indicated an increased level of immune cell infiltration within the low-risk group, characterized by significantly improved immune and ESTIMATE scores when compared to those of the high-risk group (Fig. 7A). In correlation analysis, it was observed that the risk score exhibited a negative correlation with the immune score, stromal score, and ESTIMATE score, while positively related to tumor purity (Fig. 7B). Following this, the association of risk score and tumor-infiltrating immune cells was investigated using 7 distinct algorithms, in order to ensure a more dependable assessment (Fig. 7C). Our findings indicated that infiltration levels of diverse immune cell types, including CD8+ T cells, B cells, NK cells, macrophages, and mast cells, exhibiting an inverse correlation with the risk score. This observation suggested the possibility of increased immune cell infiltration within the TME of patients categorized by low-DRLS score. Next, the ssGSEA algorithm was utilized for evaluating the immune status of TME in the two groups with distinct risk scores. Our results implied the low-risk group manifested higher ssGSEA scores in diverse immune cells, such as activated CD8 T cells, mast cells, NK cells, and effector memory CD4 T cells, in comparison to the high-risk group (Fig. 7D, P < 0.05). Additionally, the enhanced immune function scores were observed among the low-risk group, such as activated Dendritic Cells, B cells, Human Leukocyte Antigen, inactive Dendritic Cells. These findings suggested that patients categorized by low DRLS exhibited a more proactive immune response against tumor progression, which in turn may lead to a more favorable outcome (Fig. 7E, P< 0.05). Furthermore, the evaluation of 50 immune checkpoints (ICs) expressions was conducted within the two groups with different level of risk score (Fig. 7F). We identified differential expression in 18 ICs within the two groups (P < 0.05). Notably, immune checkpoints such as CD27, CD28, CD48, CD160, CD244, CD40LG, ADORA2A, BTLA, BTNL2, IDO2, and CD200R1 exhibited significantly higher expression levels in patients of the low-risk group. This suggests that patients categorized by low DRLS scores may potentially benefit from immunotherapies targeting these specific ICs. Finally, analyses were conducted for exploring the correlation of 5 DRlncRNAs and risk scores with 38 immune checkpoints. The findings revealed a tightly positive correlation between AC093673.1↗ and AC012615.1↗ with most ICs, while AL365181.2↗, GSEC, AL606834.1↗, and risk scores exhibited a negative correlation with the most ICs (Fig. 7G). These findings indicated the elevated levels of immune infiltration within low-risk group, which could potentially contribute to an enhanced prognosis. Additionally, patients with elevated expression levels of most ICs demonstrated the potential to benefit from specific immunotherapies.

The TIDE website was utilized to assess tumor immune evasion and the impact of immunotherapy on both the low- and high-risk groups. The considerably higher scores were observed in the high-risk group, suggesting that patients characterized by high DRLS score could be more susceptible to immune evasion, potentially resulting in a less favorable response to immunotherapy (Fig. 8A). Following this, we obtained Immunophenoscore (IPS) for LUAD patients from the TCIA database with the aim of evaluating the predictive value of DRLS in the context of immunotherapy. This assessment aimed to estimate patients’ immunogenicity, predicting their potential responsiveness to immune checkpoint blockades (ICBs). The findings showed elevated IPS scores in the low-risk group, suggesting the potential of patients characterized by low DRLS scores might more favorably respond to ICBs (Fig. 8B). Given that chemotherapy and targeted therapies are standard treatments for intermediate and advanced stage LUAD, we conducted an additional investigation to assess how patients categorized by DRLS respond to commonly used antitumor agents. Interestingly, when comparing the IC50 values of these agents between the two groups, a significant improvement was observed in response to targeted therapies (e.g., Dasatinib, Erlotinib, Oretinib, Savolitinib, Trametinib, and Ulixertinib) and chemotherapeutic treatments (e.g., 5-Fluorouracil, Cisplatin, and Cytarabine) in the low-risk group. These findings indicated that patients characterized by elevated risk scores might respond more favorably to these drugs (Fig. 8C).

We employed an unsupervised consensus clustering method to investigate distinct molecular subtypes according to the DRlncRNAs expression. Two clusters were generated for patients with LUAD, and we confirmed an optimal cluster stability at k = 2 (Fig. 9A). We employed PCA to evaluate the sample distribution within cluster A and cluster B (Fig. 9B). Subsequently, we explored the associations between the clusters and two groups. The findings revealed that cluster A primarily included patients categorized to the low-risk group, whereas cluster B was predominantly composed of patients characterized with high-risk score (Fig. 9C). Subsequent survival analysis indicated a significantly better prognosis of cluster A in comparison to cluster B (Fig. 9D). We generated a clinically relevant heat map by integrating clinical data, the expression levels of 5 DRlncRNAs, and the clustering results (Fig. 9E). Similarly, we assessed immune cell infiltration within the two clusters using 7 algorithms utilized in previous investigation. Our findings revealed higher levels of immune cell infiltration, such as B cells and CD4+ T cells within the TME of cluster A (Fig. 9F). This implied that the TME of cluster A recruited a greater number of immune cells and initiated a more adaptive immune response, thereby fostering an inflammatory tumor microenvironment. Analyzing the TME scores for the various clusters revealed the enhanced immune scores, stromal scores, and ESTIMATE scores among the cluster A when compared to cluster B. Nonetheless, cluster A exhibited lower tumor purity than cluster B (Fig. 9G).

In this study, we assessed the expression of five DRlncRNAs, which were utilized in our model construction, in LUAD cell lines including A549, H1975, and BEAS-2B using RT-qPCR (Fig. 10A). Among these, four DRlncRNAs exhibited upregulation in both LUAD cell lines (AL365181.2↗, GSEC, AC093673.1↗, AL606834.1↗), while AC012615.1↗ showed significant downregulation in both LUAD cell lines. According to our research, GSEC is one of the four highly expressed DRlncRNAs in tumor samples and related to a relatively higher hazard ratio (HR 1.257), indicating an unfavorable prognosis for LUAD patients. Furthermore, we observed an association between GSEC and ferroptosis, cuproptosis and disulfidptosis. Consequently, we conducted in vitro experiments for investigating the role of GSEC in LUAD progression. Firstly, GSEC was effectively silenced in A549 cells (Fig. 10B). Subsequently, we evaluated the impact of GSEC downregulation on LUAD cell proliferation through CCK-8 experiments. We found that the proliferation capacity of A549 cells decreased upon GSEC knockdown in comparison to the control group (Fig. 10C), implying that GSEC may promote LUAD cell proliferation. Similarly, colony formation experiments demonstrated a marked reduction in the number of colonies within GSEC downregulation group in comparison to the control group (Fig. 10D). Transwell assays also implied that GSEC knockdown significantly inhibited the migration of A549 cells (Fig. 10E). Wound healing experiments revealed that GSEC downregulation significantly delayed wound healing (Fig. 10F). These results implied the involvement of GSEC in the proliferation and migration of LUAD cells, highlighting its potential in promoting LUAD progression.

RNA-sequencing data and clinical data of LUAD were downloaded from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/↗). Further inquiries can be directed to the corresponding author.