The LUAD data were obtained from The Cancer Genome Atlas (TCGA)-LUAD (https://portal.gdc.cancer.gov/), GSE31210, GSE68465, and GSE72094 in Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/), which included RNA sequencing data, somatic mutation data, copy number variation (CNV) data, and corresponding clinical information. The TCGA-LUAD cohort, which consisted of 539 LUAD tissues and 59 cancer-adjacent tissue samples, served as a training cohort. Meanwhile, GSE31210, GSE68465, and GSE72094 were utilized for the validation cohort, which consisted of 1,150 LUAD patients. Additionally, the “sva” package was used to correct the batch effect between the different datasets by adopting the “combat” algorithm (17). Moreover, patients who lacked OS time were filtered out. Finally, 1569 eligible patients were encompassed in the study. The detailed clinical characteristics of all LUAD patients is listed in Supplementary Table 1 . Additionally, 18 DRs were collected from the previous study ( Supplementary Table 2 ) (12).

A consensus clustering algorithm was employed to discern optimal subtypes founded on the expression of 18 DRs using the R package “ConsensusClusterPlus”. The number of clusters (K) and their stability (with 1,000 repeats for mast k = 9) were determined by the consensus clustering algorithm (18). The clustering was based on dividing centromeres with “Euclidean” distances (the most common and familiar distance measurement methods and correlation of K-Means clustering). Additionally, a T-distributed Stochastic Neighbor Embedding (tSNE) analysis was conducted to decrease the dimensions and differentiate the subtypes of information (tSNE can preserve local similarities between data points and is one of the most used unsupervised clustering visualization methods).

DEGs among DRs subtypes were determined using the “limma” R package with the filtering criteria of log (fold change) >1 and False-discovery rate (FDR) <0.05 (19). To further investigate the potential functions and enriched pathways of DEGs, functional enrichment analyses were conducted on DEGs employing the “clusterprofiler” R package (20).

The DS was calculated to quantify the disulfidptosis patterns of the LUAD. Accordingly, a univariate cox regression analysis was conducted to identify the DEGs are related to prognosis. Subsequently, the patients were segregated into distinct gene cluster groups according to the expression of prognostic DEGs using the unsupervised clustering method. Based on prognostic DEGs, the 10 machine learning algorithms, including ‘Least Absolute Shrinkage and Selection Operator (LASSO, “glmnet” package)’, ‘Ridge (“glmnet” package)’, ‘Elastic network (“glmnet” package)’, ‘StepCox (“survival” package)’, ‘Survival support vector machine (survival-SVM, “survivalsvm” package)’, ‘CoxBoost (“CoxBoost” package)’ (21), ‘Supervised principal components, “superpc” package’ (22), ‘partial least squares regression for COX, “plsRcox” package’, ‘random survival forest (RSF, “randomForestSRC” package)’, ‘generalized boosted regression modeling, “gbm” package’ were used to constructed the models. Briefly, 101 combinations of 10 machine learning algorithms were used to build the models based on a leave-one-out cross-validation (LOOCV) framework. Models with <3 genes were excluded. Simultaneously, each patient’s linear score and concordance index (C-index) are calculated based on various models, and the optimal model is selected based on the average highest C-index in the training and testing cohorts. Utilizing the optimal model, the DS for each patient was determined. See the supplementary methods table for details. The patients were then separated into high- and low-DS groups using the median DS value. Subsequently, Kaplan–Meier survival analysis was utilized to compare the OS rates of patients in various DS groups.

The “maftools” package was utilized to evaluate and visualize the mutation type and frequency of the genes (23). Correspondingly, the tumor mutation burden (TMB) of each LUAD sample was calculated based on the total count of somatic mutations per megabase (MB) in the exon coding region of the human genome. Different mutation types were classified as either synonymous or nonsynonymous mutations. The nonsynonymous variants included Frame_Shift_Del, Frame_Shift_Ins, In_Frame_Del, In_Frame_Ins, Missense, Nonsense, Nonstop, Splice_Site, and Translation_Start_Site. The maftools analysis focused on identifying significantly mutated genes (P <0.05) between the two groups and assessing the interaction effect of gene mutations. Only genes with mutations occurring 30 times or more in at least one group were considered for both analyses. GISTIC 2.0 was used to identify significant regions within CNV data (24). To quantify and compare CNVs, we calculated the fraction of altered genome (FGA), fraction of genome gained (FGG), and fraction of genome lost (FGL) for each sample. FGA reveals the proportion of genomic segments that have been altered. FGG/FGL considers specifically the genomic segments that have undergone gain or loss, respectively.

TME at LUAD were evaluated under four aspects. First, the immune score, stromal score, ESTIMATE score, and tumor purity were calculated using the ESTIMATE algorithm (25). Secondly, three algorithms, Single Sample Gene Set Enrichment Analysis (ssGSEA), Tumor Immune Estimation Resource (TIMER, https://cistrome.shinyapps.io/timer/), and “MCPcounter” were used to quantify the relative infiltration of immune cells in the entire cohort (26–28). Thirdly, the seven steps of cancer immunity cycle were analyzed using the Tracking Tumor Immunophenotype (TIP) website (http://biocc.hrbmu.edu.cn/TIP/) (29). In the fourth step, 35 inhibitory immune checkpoints were extracted from a prior study. Subsequently, Gene Set Enrichment Analysis (GSEA) was utilized to identify underlying mechanisms in hallmark gene sets with the recommended criteria (FDR <0.25 and NES >1), in order to determine the underlying hallmark pathways associated with DS (30). Meanwhile, “Gene Set Variation Analysis” (GSVA) package was applied to the two DS groups with an adjusted p value <0.01 (31). The “h.all.v7.4.symbols.gmt” hallmark gene sets from the MSigDB database were employed for GSVA implementation.

To explore the predictive value of DS in LUAD patients after immunotherapy, we compared the immune dysfunction and exclusion (TIDE, http://tide.dfci.harvard.edu/) score; additionally, the submap algorithm was applied to compare the efficacy of immunotherapy among various DS subtypes (32, 33). In addition, four immunotherapy-treated cohorts, IMvigor210, GSE35640, GSE79671, and GSE173839, were collected to investigate the immunotherapy response ability of DS. The sensitivity of tumor cell lines to potential drugs was obtained from the Cancer Therapeutics Response Portal (CTRP, https://portals.broadinstitute.org/ctrp) and Profiling Relative Inhibition Simultaneously in Mixtures (PRISM, https://depmap.org/portal/prism/). The more sensitive a cell line is to a potential drug, the lower its area under the curve (AUC).

ScRNA-seq data was downloaded from GSE131907. Subsequently, the “Seurat” R package was utilized to measure the gene expression levels by processing the raw data from each sample (34). Cells with fewer than 200 detected genes were removed. Accordingly, the top 2000 highly variable genes were selected for subsequent clustering analysis. Following this, single cells were classified into distinct subgroups via the application of the FindNeighbors and FindClusters functions (dim = 15 and resolution =0.2). In addition, the tSNE was constructed utilizing the top 15 primary components. Subsequently, immune cells and tumor cells were identified using the “single R” and “copyKAT” packages (35, 36). Pseudotime trajectory analysis is a method used to study the temporal order of cells during development, differentiation, or other biological processes. Its goal is to reduce the dimensionality of single-cell data from high-dimensional to one-dimensional, thus representing the temporal changes of cells in a pseudotime manner. This aids in understanding the timeline of processes like differentiation, development, and transformation in tumor cells. Hence, the cell trajectory cancer cell populations were ordered in pseudotime using the “Monocle” package (37).

Two human LUAD cells (PC-9 and H838) as well as a normal bronchial epithelial cell (BEAS-2B) were obtained from the American Type Culture Collection (Manassas, VA, USA). The PC-9 and H838 cells were grown in RPMI 1640 supplemented with 10% FBS and 1% penicillin-streptomycin, while BEAS-2B cells were grown in DMEM supplemented with 10% FBS and 1% penicillin-streptomycin in humidified air at 37°C and 5% CO₂ (This condition simulates the physiological environment in the human body, aiding in maintaining normal cellular growth, metabolism, and function). To construct KYNU knockdown and negative control (NC) cell lines, H838 and PC-9 cells were seeded in 6-well plates at a density of 5 × 10⁴ cells/well and transfected with 50 nM siRNAs-KYNU and siRNA negative control (siRNA-NC) using lipofectamine 3000 (This reagent possesses high efficiency, broad spectrum, and low toxicity) following the manufacturer’s guidelines (Hanheng, Shanghai, China). After 48 h of transfection, subsequent experiments were conducted (Typically, after 48 hours of transfecting siRNA, the target gene’s expression is effectively disrupted in the vast majority of cases, while cells continue to maintain a relatively healthy and appropriate growth state). The siRNA sequences are provided in Supplementary Table 3 .

For this study, a total of 15 TM tissue samples (HLugA030PG04-1) were utilized, including 15 LUAD and 15 adjacent non-tumor tissues. All the tissues were procured from Shanghai Outdo Biotech Co., Ltd. (Shanghai, China). TM was stained using IHC with KYNU. In brief, the antibodies were diluted to the suitable concentration and incubated overnight with the sections at 4°C. The avidin–biotin and streptavidin complex were then incubated with the biotinylated goat anti-rabbit IgG secondary antibody for 0.5 h. The cell nuclei were counterstained blue with Hematoxylin. Each specimen was graded based on the intensity (0: absent, 1: mild, 2: moderate, and 3: pronounced) and the proportion of positively stained cells (0: 0%, 1: 1%–25%, 2: 26%–50%, 3: 51%–75%, and 4: 76%–100%). The final IHC scores were then calculated as the product of the intensity score and the percentage score (38).

RNA extraction was performed employing RNA extraction kit (AG, Changsha, China) in accordance with the manufacturer’s instructions. The complementary DNA (cDNA) was synthesized from total RNA for each sample by using a reverse transcription kit (Promega, Madison, Wisconsin). Subsequently, the mRNA expression was quantitatively analyzed via RT-PCR (Roche Light Cycler^® 480II System). The relative gene expression was normalized to β-actin, acting as the control gene. The primer sequences for the genes are outlined in Supplementary Table 4 . The fold change of the target gene was computed utilizing the 2^−ΔΔCt method.

Cells in logarithmic phase were seeded into 96-well plates at density of 2,000 cells per 100 µl. They were then cultured for 0, 24, 48, 72, and 96 h. At the end of the incubation cycle, CCK-8 reagent (Saiku, Shanghai, China) was added to each well for 2 h at 37°C. Subsequently, the absorbance values of each well were measured at 450 nm. For colony formation assay, 500 cells were seeded into each well of a 12-well plate and incubated. After 10 days, the cells were stained with crystal violet solution, and fixed with 4% paraformaldehyde for 0.5 h. Finally, colonies consisting of more than 50 cells were counted.

The transwell experiment was executed using a transwell chamber (Transwell, Corning Costar, USA). For migration assays, a total of 4 × 10⁴ cells were resuspended in 200 μl serum-free medium and placed to the upper chamber. The lower chamber was filled with 700 μl of medium containing 10% FBS. After 18 h, cells that migrated through the membrane were fixed and stained with hematoxylin. Following a 24 h-incubation period, cells that traversed the membrane were fixed with 4% paraformaldehyde for 30 minutes and stained with a crystal violet solution for 15 min. Subsequently, they were observed with a microscope, and five random fields of view were selected to count the cells. The protocol for the invasion assay was similar to the migration assay, with the exception that the upper chambers were coated with 70 μl of diluted Matrigel. In addition, the wound healing assay was utilized to examine cell migration. Serum-starved LUAD cells (2 × 10⁵) were seeded in 6 well plates and then transfected with siRNA-NC or siRNA-KYNU. Upon reaching >95% cell confluence, a sterile 20 μl pipette tip was used to create a scratch in the monolayer of cells. After a PBS wash, the cells were incubated for 0 and 24 h in the corresponding basic culture medium. The monolayer cells were examined microscopically, and the gap distance was measured quantitatively to ascertain LUAD cell migration.

All statistical analyses and representations were conducted using R (version 4.2.1) and GraphPad Prism (version 9.00). The Chi-squared test was used to compare the proportion of individuals within two groups. Additionally, continuous variables in two or more groups were compared using Wilcoxon rank-sum test or Kruskal–Walli’s test. For correlation analysis, the Pearson correlation test was employed. A p-value less than 0.05 was deemed statistically significant.

Gene expression analysis using bulk RNA-seq demonstrated that the majority of DRs displayed relative higher expression levels in LUAD tissues compared to para-carcinoma tissues ( Figure 1A ). As depicted in Figure 1B , 104 of 616 (16.88%) LUAD samples possessed genetic mutations. A total of 13 of 18 DRs were found to be mutated, with CNOT1 exhibiting the highest rate of mutation. Among them, missense mutations were found to be the most frequent ( Figure 1B ). In order to unmask the genetic modifications in DRs, we presented an overview of the frequency of somatic and copy number mutations with malignancies. Moreover, analysis of these 18 DRs revealed that CNV alterations were common. NDUFS2, NUBPL, PPM1F, EPAS1, and LRPPRC displayed widespread CNV amplification, whereas CCNC, NDUFA11, OXSM, and GYS1 showed widespread CNV deletions ( Figure 1C ). Figure 1D depicts the locations of CNV alterations in LUAD DRs. The correlation network composed of 18 DRs is illustrated in Figure 1E (correlation coefficient >0.4; the positive correlation is represented by the red line, negative correlation is represented by the blue line).

In order to delve further into understanding the DRs’ expression pattern in the oncogenesis of LUAD, 1569 patients from 4 independent LUAD cohorts were included. Accordingly, we performed unsupervised clustering and classification on the combined LUAD cohort based on DRs. Our results showed that k = 2 was the optimal choice ( Figures S1 ). In addition, the tSNE results revealed significant differences between the two clusters in terms of DRs expression ( Figure 2A ). Moreover, the Kaplan–Meier survival analysis demonstrated that the DRcluster A had a greater survival advantage than DRcluster B ( Figure 2B ). In addition, the clinicopathological characteristics of the various DR clusters also revealed significant differences ( Figure 2C ). To further explore the biological behavioral difference between these two clusters, we also conducted a GSVA enrichment analysis ( Figure 2D ). The results demonstrated that DRcluster A was primarily enriched for carcinogenic pathways like focal adhesion, EMC receptor interaction, and others. Moreover, Figure 2E data also revealed significant differences in the relative expression of immune infiltration cells across two DRclusters.

A total of 51 DEGs were identified from two DRclusters using the “limma” package. These DRs subtypes-related genes were significantly enriched in cellular metabolism ( Figures 3A, B ). Subsequently, we conducted a univariate Cox regression analysis for DEGs and identified 39 DEGs with prognostic importance for LUAD. To further investigate the heterogeneity of DRs subtypes, we performed unsupervised clustering on 39 DEGs. The 1569 patients with LUAD were also separated into two gene clusters, designated geneCluster A and B ( Figure S2 ). Similar to DRcluster, these genes accurately differentiate LUAD patients, with distinct clusters of genes exhibiting variations in DRs, survival rates, and immune cell infiltration ( Figures 3C–F ). Following this, the LOOCV framework was used to fit 101 prediction models to both training and testing sets. The C-index for each model was calculated, and based on the model with the highest average C-index (0.713), “Lasso+RSF” was deemed the optimal model ( Figures 4A–C , Supplementary Table 5 ). A DS was calculated for each patient based on the expression of 7 genes (KRT6A, NEIL3, KYNU, ABCC2, SFTPC, CPS1, and INSL4) weighted by their regression coefficients in the model, and patients were divided to high- or low‐DS groups based on the median cutoff point of DS ( Supplementary Table 6 ). As evident from the K–M survival analysis, OS rates were significantly diminished in the high-DS group compared to the low-DS group ( Figures 4D–G ). Moreover, the relationship between different types of patients and their prognoses was analyzed ( Figure 5 ), with results suggesting that a low DS was related to a better prognosis in all patient categories.

A time-dependent receiver operating characteristic (ROC) curve was employed to assess the validity of DS, as the AUC value for TCGA (0.93–0.96), GSE31210 (0.68–0.79), GSE68465 (0.63–0.66), GSE72094 (0.63–0.72), and meta-cohort (0.71–0.78) ( Figures 6A–E ). In addition, the C-index of clinical factors in patients with LUAD was determined ( Figures 6F–J ). Notably, the DS had a higher predictive efficacy than the vast majority of clinical indicators. Subsequently, both DS and clinical indicators were subjected to univariate and multivariate Cox analyses. In all cohorts, the DS was determined to be an independent indicator of OS prognosis ( Tables 1 – 4 ). In order to determine the prognostic efficacy of DS, we combined 56 previously published LUAD prognostic models and conducted a comparative analysis of each model’s C-index. These models were developed using a variety of biologically relevant features, including autophagy, EMT, ferroptosis, hypoxia, necroptosis, glycolysis, and m6A methylation. Accordingly, DS was found to exhibit superior performance relative to the vast majority of models across all cohorts ( Figure 7 , Supplementary Table 7 ). Cumulatively, these results, therefore, demonstrate that the DS would be a valuable LUAD prognostic model.

Using the “maftool” package, a comparison was made between the distribution differences of somatic mutations observed in high- and low-DS groups. ( Figures 8A, B ). Comparing the frequency of mutants between the high- and low-DS groups, more somatic mutations, both synonymous and nonsynonymous, were observed in the high-DS group ( Figures 8C–E ). In addition, maftools analysis results showed that 17 genes mutated more frequently in LUAD patients in the high-DS group, including KEAP1, STK11, SMARCA4, TTN, NCKAP5, COL5A2, ANKRD30A, TEX15, PCDH15, GRIN2B, AHNAK, FAT4, FMN2, FAT1, ZNF804B, DOCK2 and COL22A1 ( Figure 8F ), and there was extensive co-mutation between these genes ( Figure 8G ). This is consistent with the mutation analysis described above; accordingly, TMB was found to be higher in the high-DS group compared to the low-DS group ( Figure 8H ). Subsequently, LUAD patients were classified into two mutation groups based on their TMB score. When combining DS and TMB, we discovered that patients with low TMB from the high-DS group had the worst prognosis ( Figure 8I ). Subsequently, we used the GISTIC 2.0 software to decipher the amplification and deletion of CNA on chromosome. Compared to the high-DS group, the low-DS group had a greater burden of amplification and deletion at both the arm and focal levels ( Figure 8J ). However, no significant disparities were observed between the high- and low-DS groups in terms of FGA, FGG, and FGL ( Figure 8K ). These results indicate a certain degree of correlation between DS and mutations.

In order to evaluate the discriminative potential of the DS subgroup for the TME and its applicability in immunotherapy, we simultaneously evaluated the abundance of immune cell infiltration across multiple samples using four distinct algorithms. Unsurprisingly, as the DS increased, the number of immune cells declined ( Figure 9A ). Apparently, the activation of key steps in the cancer immunity cycle, such as step 3 (priming_and_activation) and step 4 (CD4 T cell recruiting, Dendritic cell recruiting, Macrophage recruiting, Monocyte recruiting, and T cell recruiting) appeared to be significantly higher in the low-DS group than in the high-DS group ( Figure 9B ). Subsequently, the expression profile of immune checkpoints in the two DS groups were further evaluated. Accordingly, the analysis revealed that the low-DS group demonstrated elevated expression levels of immune checkpoints, including HHLA2 and CD48 ( Figure 9C ). Given the observed upregulation of immune-related characteristics in the group with low DS, its underlying biological mechanisms were investigated further. As evident from the findings, DS exhibited correlations with multiple metabolic pathways ( Figure 9D ). In addition, DS demonstrated a strong correlation with numerous immunotherapeutic strategies ( Figure 9D ). To investigate the cancer signaling pathways associated with DS, GSVA analysis was performed on high- and low-DS groups. Using a predetermined threshold, we discovered that 16 signature pathways were significantly upregulated in the high-DS group compared to the low-DS group ( Figure 9E ). GSEA validated that 12 of these pathways were upregulated in high-DS patients, the majority of which are known to be carcinogenic ( Figures 9F–I ).

The potential of an immunotherapy response was subsequently predicted for each immune cluster utilizing the TIDE algorithm and submap analysis. Accordingly, lower TIDE scores were observed in the low-DS group, implying a higher sensitivity to immunotherapy in these patients ( Figure 10A ). Moreover, the submap results indicated that the group with a low DS level was more sensitive to CTLA4 inhibitors ( Figure 10B ). Although we evaluated an individual’s immunotherapy efficacy using two algorithms, it remains critical to directly compare the curative efficacy of immunotherapy cohorts across various DS groups. As a result, four immunotherapy cohorts were included for further analysis. In the IMvigor210 cohort, patients with DR/PR had significantly longer OS compared to patients with SD/PD, whereas the influence of DS on patient prognosis was minimal ( Figure 10C ). However, patients who responded better to immunotherapy had lower DS levels across all cohorts ( Figures 10D–G ). Although there were no significant differences in patient survival and DS between the two groups for some cohorts, the propensity for these results was consistent for the other cohorts. In addition, these results demonstrate that the DS was able to predict the efficacy of ICBs and can provide direction for the deployment of immunotherapy. Thus, to pinpoint candidate drugs that may exhibit heightened sensitivity in LUAD patients, we conducted drug response predictions using CTRP- and PRISM-derived data. Finally, the cross-correlation of the two pharmacogenomics databases allowed us to predict four drugs or compounds (including SB−743921, ispinesib, cabazitaxel, and gemcitabine) with therapeutic potential in patients ( Figures 10H, I ).

To analyze the expression of DS in TME, we used the LUAD single-cell dataset GSE131907 from the GEO database. All the cells were partitioned into 16 clusters using the k-nearest neighbor (KNN) clustering algorithm ( Figure 11A ). Subsequently, using the “single R” and “copycat” packages to annotate all cells, we were able to identify 9 distinct cellular subtypes, including B cells, endothelial cells, epithelial cells, cancer cells, macrophage cells, monocyte cells, smooth muscle cells, NK cells, and T cells ( Figure 11B ). Most of these cells are important components of the TME mentioned in the above results. Subsequently, we investigated the single-cell transcriptome localization of 7 genes in DS ( Figure 11C ). Concurrently, the DS for each cell was calculated, which showed that cells with a high DS predominately resided in the region of cancer cells ( Figure 11D ). Additionally, the temporal sequence of cancer cellular differentiation was revealed by the analysis of the pseudotime trajectory. Accordingly, low DS cancer cells appear at an earlier pseudotime than high DS cancer cells, which are primarily found in the earliest stages of differentiation ( Figures 11E–H ).

The qRT-PCR experiments were performed on LUAD cell lines to confirm the expression levels of DS genes in LUAD. Our finding discovered that ABCC2, NEIL3, KYNU, and CPS1 exhibited elevated expression in LUAD cell lines, while KRT6A and SFTPC were found to be underexpressed ( Figure 12A ). In addition, KYNU exhibited the most significant correlation between high expression and unfavorable poor patient prognosis. No reports have documented the role of KYNU in LUAD to date. As a result of these considerations, KYNU was selected as the focus of further experiments. Correspondingly, IHC staining analysis showed that the protein level of KYNU expression was elevated in LUAD tumor tissues relative to paracancerous tissues ( Figure 12B ). Subsequently, we investigated the function of KYNU in PC-9 and H838 LUAD cell lines through a series of cell-based experiments. Initially, the effect of the siRNA was confirmed via RT-qPCR ( Figure 13A ). As demonstrated by CCK-8 and clone formation assays, KYNU knockdown inhibits LUAD cell growth and their clone formation capacity ( Figures 13B, C ). In addition, wound healing and transwell assays confirmed that KYNU knockdown inhibited the cell migration and invasion capabilities of LUAD cells ( Figures 13D, E ). The precision of these findings corroborated that the expression of KYNU mirrored the variations anticipated through bioinformatic prediction.