Developing a prognostic model using machine learning for disulfidptosis related lncRNA in lung adenocarcinoma

After excluding patients with incomplete clinical records, we analyzed data from 916 LUAD patients derived from The Cancer Genome Atlas (TCGA) datasets [https://portal.gdc.cancer.gov/↗] and Gene Expression Omnibus (GEO) [https://www.ncbi.nlm.nih.gov/geo/↗] datasets. The transcripts per million (TPM)-normalized expression dataset for TCGA-LUAD was retrieved using the “easyTCGA” package. To enhance the interpretability of expression values, we applied a logarithmic transformation with a base of 2. Data from the GEO were obtained from the Affymetrix® GPL570↗ platform, employing the Human Genome U133 Plus 2.0 Array for gene expression analysis. The raw data from Affymetrix® were processed with the robust multiarray averaging method using the "Affy" R package. We corrected for batch effects using the ComBat algorithm and merged three GEO datasets to form a comprehensive Meta-GEO cohort. Expression values for each gene were standardized as z-scores across all patient cohorts. A flowchart depicting the study methodology is presented in Fig. 1, and detailed data from the three GEO datasets are provided in Supplementary Table 1.

In our comprehensive literature review, we identified 24 genes associated with disulfide-triggered cell death⁶. Using the "limma" package, we examined the relationship between 24 disulfidptosis-related genes (DRGs) and lncRNAs, identifying significant lncRNAs associated with disulfidptosis (p-value < 0.010, |Correlation|> 0.3). These lncRNAs, combined with data from the GEO database, formed the basis of our model. To enhance model accuracy and robustness, we implemented a novel strategy starting with a univariate regression analysis of lncRNAs from TCGA-LUAD to identify potential prognostic lncRNAs. Firstly, we initiated a univariate regression analysis on the sequentially identified lncRNAs from TCGA-LUAD, pinpointing potential prognostic-linked lncRNAs. We then designated TCGA-LUAD as our training set and the three GEO datasets as validation. We engaged ten machine learning methods—Ridge Regression, Least Absolute Shrinkage and Selection Operator (Lasso) Regression, Stepwise Cox Regression, CoxBoost, Random Survival Forests (RSF), Elastic Net (Enet), plsRcox, Supervised Principal Components (SuperPC), survival-Support Vector Machines (SVM), and Gradient Boosting Machine (GBM). Each model was chosen for its unique strengths and ability to complement the analysis in various ways:

(1) Ridge Regression: Known for its ability to handle multicollinearity, Ridge regression applies L2 regularization to prevent overfitting, making it suitable for our dataset with a high dimensionality of genomic features.

(2) Lasso Regression: Lasso employs L1 regularization which not only helps in reducing overfitting but also in feature selection by shrinking some coefficients to zero. This is particularly useful in identifying the most relevant DRLs influencing LUAD progression.

(3) Stepwise Cox Regression: This method applies both forward and backward selection techniques to identify significant variables in survival analysis, ideal for modeling the impact of DRLs on patient survival.

(4) CoxBoost: CoxBoost adjusts for high-dimensional data by incrementally fitting the Cox model, enhancing its predictive accuracy in our setting where numerous potential predictors are involved.

(5) RSF: An extension of the random forest that handles censored data, RSF is adept at dealing with complex interactions and non-linear relationships, crucial for understanding intricate DRL interactions.

(6) Enet: Combining the properties of both Ridge and Lasso, Elastic Net is effective when dealing with datasets where numerous features are correlated, providing a balanced approach to regularization and variable selection.

(7) plsRcox: This model is useful for reducing dimensionality while maintaining the relationship with the survival outcome, helpful in elucidating the impact of DRLs on LUAD.

(8) SuperPC: SuperPC is employed to enhance predictive accuracy by focusing on principal components most associated with survival outcomes, refining our analysis of DRLs' prognostic power.

(9) survival-SVM: This model adapts the traditional SVM for survival analysis, which is beneficial for handling non-linear patterns in the data, providing a robust classification of risk groups.

(10) GBM: GBM builds an ensemble of decision trees sequentially, improving the model iteratively, which is key in accurately classifying patients into risk categories based on their DRL profiles.

These methods were amalgamated in 101 unique combinations for variable selection and model formulation. Finally, we evaluated the efficacy of our model on the training and validation datasets employing the concordance index (C-index), subsequently calculating an average score. The optimal selection was the Enet model with α = 0.2. Leveraging the survminer package, we determined the risk score threshold from the surv_cutpoint for the LUAD-TCGA dataset, subsequently segmenting both TCGA-LUAD and the three GEO datasets into low-risk and high-risk groups.

Utilizing the “survminer” package in R, we conducted Kaplan–Meier (KM) curve analyses on TCGA-LUAD and three GEO datasets, probing the disparities in LUAD patients' overall survival (OS) and disease-free survival (DFS) between the high-risk and low-risk cohorts. To evaluate the potential of the risk score as a standalone prognostic marker, we utilized univariate and multivariate cox analyses. With the “timeROC” package, we generated receiver operating characteristic (ROC) curves for one, two, and three-year durations for the risk score taken from TCGA-LUAD and the three GEO datasets, subsequently computing the area under the ROC curve (AUC). We delved into the relationship between risk score and patient attributes, including TNM staging, gender, smoking history, and mutations in the EGFR/KRAS/ALK genes. Additionally, we devised a nomogram that amalgamates age, gender, staging, and risk score to maximize predictive accuracy. The clinical relevance of this nomogram was then gauged using calibration and decision curve analysis (DCA).

The TMB was sourced from the TCGA-LUAD database. TMB quantification was carried out utilizing the “tmb” function from the “maftools” package. For visualization and in-depth analysis, a waterfall plot was employed to delineate the top 10 genes exhibiting the highest TMB in patients diagnosed with TCGA-LUAD. Moreover, a survival curve analysis was conducted to discern the disparities in OS among patient groups categorized based on distinct TMB levels.

We conducted a differential expression analysis with criteria set at |logFC|> 1 and p-value < 0.050 to identify significant differences between low-risk and high-risk cohorts. The differentially expressed genes (DEGs) identified through this analysis were then subjected to functional enrichment analysis. Specifically, we used Gene Ontology (GO)^13–15 enrichment to identify relevant biological processes, molecular functions, and cellular components. We utilized the GO database for Gene Set Enrichment Analysis (GSEA). Additionally, we performed pathway enrichment analysis using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database [https://www.genome.jp/kegg/↗] to uncover key signaling pathways associated with the DEGs.

We utilized tools like CIBERSORT, MCP-counter, EPIC, TIMER, xCell, and quanTIseq to probe differences in the TME between high-risk and low-risk cohorts. The analysis was streamlined using the R package “IOBR”, and the results were visualized as heatmaps via the “pheatmap” package. Further, we displayed CIBERSORT data using box plots, highlighting variations in immune cell infiltration between high-risk and low-risk cohorts. Drawing from our prior research¹⁶, we used a select set of 28 genes, recognized for representing immune cells, and carried out single sample genome enrichment analysis (ssGSEA). In our final steps, we utilized the “ESTIMATE” algorithm to quantify the components of immune and stromal cells in the TME of each LUAD patient.

We tapped into the Tumor Immune Dysfunction and Exclusion (TIDE) scores, noted for its predictive accuracy in immune checkpoint blockade (ICB) treatment outcomes¹⁷, alongside the Immune Prognostic Scores (IPS), a marker of immunogenicity. Data was sourced from the TIDE [http://tide.dfci.harvard.edu/↗] and The Cancer Immunome Atlas (TCIA) [https://tcia.at/home↗]. We juxtaposed the IPS and TIDE scores across our patient cohorts. Building on a study from Cancer Medicine¹⁸, we also identified genes linked to tumor-infiltrating lymphocytes (TLS), pivotal for immunotherapy success. To probe the relationship between risk scores and ICB treatment outcomes, we leveraged a compilation of 79 immune checkpoint genes (ICGs) drawn from prior research¹⁹. Individual samples received unique scores based on both TLS and ICG gene sets, and scores from the high-risk and low-risk cohorts were then compared. Exploring the signature's potential for targeted drug therapy in LUAD, we gauged the half maximal inhibitory concentration (IC50) of pertinent chemotherapeutic agents using the “pRRophetic” R package. We then employed Wilcoxon's test to discern IC50 variations between the risk groups.

We carried out a differential gene analysis using GTEx and TCGA-LUAD on the Xena Functional Genomics Explorer (https://xenabrowser.net/↗). We subsequently performed a single-cell analysis utilizing the TISCH database (http://tisch.comp-genomics.org/home/↗) and conducted a comprehensive ceRNA analysis. We began by using DIANA-LncBase v3 (www.microrna.gr↗) to identify miRNAs targeted by LINC00857. We then used these miRNAs to target the corresponding mRNAs in miRTarBase (http://miRTarBase.cuhk.edu.cn/↗).

We employed the R 4.3.1 software for statistical analysis. To assess the relationship between lncRNAs and DRGs, we applied the pearson correlation method. To construct prognostic models, we first conducted univariate cox regression to assess the prognostic significance of each lncRNA individually. LncRNAs with a p-value less than 0.05 in the univariate analysis were subsequently included in the machine learning methods. The model's predictive accuracy was gauged using DCA curves, C-index, ROC analysis, and the calibration curves, while KM analysis shed light on the survival dynamics of distinct groups. Our data underwent scrutiny using a mix of statistical tests such as the Chi-square, Wilcoxon, and Log-rank tests. The choice of test was based on data characteristics, distribution, and the number of groups in comparison. For all tests, statistical significance was denoted by a p-value below 0.05.

In our investigation of the LUAD landscape, we analyzed 16,882 lncRNAs derived from the TCGA-LUAD database. This comprehensive evaluation led to the identification of 708 DRLs, which demonstrate significant interactions with DRGs, as depicted in a sankey diagram (Fig. 2A). Through further analysis incorporating data from three GEO databases, we narrowed these DRLs down to 199 lncRNAs consistently present across datasets, suggesting a pivotal role in LUAD pathogenesis (Fig. 2B). Our prognostic assessment using univariate cox regression analysis revealed 37 lncRNAs with significant implications for LUAD patient outcomes (Fig. 2C). Leveraging these lncRNAs, we constructed a predictive model employing an ensemble of machine learning techniques, with the ensemble model (Supplementary Table 2) achieving a notably high C-index of 0.677[95% confidence interval (CI) 0.63 to 0.73], suggesting robust predictive performance (Fig. 2D). This model's effectiveness was further validated through a risk stratification system, categorizing patients into high and low-risk groups based on their lncRNA expression profiles. This stratification was substantiated by principal component analysis (PCA), which confirmed the distinct separation between the risk groups, underscoring the potential of our model in clinical risk assessment (Fig. 2E).

Our survival analysis using the TCGA-LUAD dataset revealed a significant distinction in OS between the high- and low-risk groups identified through our model (p < 0.001, log-rank test) (Fig. 3A). This finding was consistently replicated across three independent GEO datasets, demonstrating significant differences in both OS (GSE31210↗, p = 0.001; GSE30219↗, p = 0.019; GSE50081↗, p = 0.025) (Fig. 3B–D) and DFS (GSE31210↗, p < 0.001; GSE30219↗, p = 0.009; GSE50081↗, p = 0.023) (Supplementary Fig. S1A–C). The predictive power of the risk score was superior to that of traditional prognostic factors such as age, gender, and staging, as evidenced by the C-index comparison (Supplementary Fig. S1D). The risk score also emerged as an independent prognostic indicator in our univariate and multivariate cox analyses (p < 0.001) (Supplementary Table 3). Multicollinearity within the model was assessed using the variance inflation factor, which was below 10 for all variables (Supplementary Table 4). The AUC analysis further validated the robustness of our model, with one-year, two-year, and three-year AUCs of 0.76, 0.72, and 0.74, respectively, in the TCGA-LUAD dataset (Fig. 3F). The external validation using GEO datasets underscored the model's accuracy, particularly notable in GSE30219↗, GSE50081↗ and GSE31210↗ for the evaluated intervals (Fig. 3G,I).

Further analysis showed gender-specific differences in risk scores across various pathological stages. In early stages (I and II), men exhibited significantly higher risk scores compared to women (Stage I: p = 0.015; Stage II: p = 0.006; Wilcoxon test) (Supplementary Fig. S2A,B). However, these differences were not observed in later stages (III/IV) (p = 0.900, Wilcoxon test) (Supplementary Fig. S2C), suggesting stage-specific risk dynamics. In addition, our study uncovered notable disparities in risk scores among patients with mutations in EGFR, ALK, and KRAS genes in the GSE31210↗ dataset (p < 0.001, Kruskal–Wallis test) (Supplementary Fig. S2D). Patients harboring these mutations also exhibited better OS compared to those without (p = 0.018, log-rank test) (Fig. 3E), highlighting the potential prognostic relevance of genetic profiles in LUAD. The impact of smoking, a known risk factor for LUAD, was evident as significant differences in risk scores between smokers and non-smokers were observed in analyses of the GSE30210↗ and GSE50081↗ datasets (GSE31210↗, p = 0.003; GSE50081↗, p = 0.027; Wilcoxon test) (Supplementary Fig. S2E,F).

To enhance our model's utility in clinical decision-making, we developed a nomogram that incorporates the identified risk scores alongside essential clinical parameters—age, gender, and TNM staging. This integration aims to provide a more comprehensive tool for predicting the prognosis of LUAD patients (Fig. 4A). We rigorously validated the nomogram's predictive accuracy using calibration curves, which compare the predicted survival probabilities against the observed outcomes. The results demonstrated a high degree of concordance, indicating that our nomogram accurately reflects patient survival rates (Fig. 4B). Further assessment through DCA (Fig. 4C-E) confirmed that the nomogram provides substantial clinical benefit. Notably, the analysis showed that the nomogram significantly outperforms the predictive capabilities of the risk score alone, particularly in terms of net benefit across a wide range of threshold probabilities.

A marked difference in TMB was discerned between the high- and low-risk cohorts (p < 0.001 by wilcoxon test) (Fig. 4F). The waterfall plot delineates the mutational landscape of the ten most prevalent genes across both risk strata. In the low-risk cohort, approximately 84.53% of specimens exhibited gene mutations (Fig. 4G), whereas in the high-risk stratum, mutations were observed in roughly 95.33% of specimens (Fig. 4H). Predominant mutations within the high-risk category included TP53, TTN, and CSMD3.

The differential expression analysis revealed a total of 1474 DEGs between the low-risk and high-risk cohorts. Among these, 568 genes were upregulated and 906 genes were downregulated. The volcano plot (Supplementary Fig. S2G) illustrates the distribution of these DEGs. These results indicate that specific genes are significantly associated with risk stratification in our study cohort. In the GO analysis (Fig. 5A,D), DEGs showed predominant enrichment in terms of molecular functions such as organic anion transport, carboxylic acid transport. Regarding cellular components, the main enrichment was observed in the apical plasma membrane (Fig. 5C). Figure 5E demonstrates the GSEA results, highlighting significant enrichment of specific gene sets related to metabolic processes, DNA binding, and hyperkeratosis. The KEGG result highlighted a significant enrichment of DEGs in neuroactive ligand-receptor interaction and the cAMP signaling pathway (Fig. 5B).

To validate the precision of our results, we employed seven techniques: CIBERSORT, EPIC, MCP-counter, xCell, TIMER, quanTIseq, and ssGSEA, to assess immune cell penetration in both high-risk and low-risk categories (Fig. 6A). With the ssGSEA data, we explored the connection between TME and several characteristics of lung adenocarcinoma patients, such as age, gender, and disease stage (Fig. 6B). We then visualized this data with box plots for both CIBERSORT and ssGSEA (Fig. 6C,D). These plots showed that the infiltration levels of B cells memory, T cells CD4 memory resting, and Monocyte was notably lower in the high-risk group compared to the low-risk group. With the help of the “ESTIMATE” algorithm, we evaluated the stromal (Fig. 6F), immune (Fig. 6E), and ESTIMATE scores (Supplementary Fig. S3A) across the different risk groups. This allowed us to gauge tumor purity. Our study suggests that the high-risk group has reduced stromal, ESTIMATE, and immune scores. Conversely, the score of tumor purity in the low-risk group is less than that in the high-risk group (Supplementary Fig. S3B).

We calculated the TIDE score and forecasted the immunotherapy response in both groups of the high risk and low risk (Fig. 7A). Based on results from both datasets, patients in low-risk group seem more inclined to show a positive reaction to immunotherapy. Additionally, IPS for the combination of anti-CTLA4 and anti-PDL1 treatment, as well as for anti-CTLA4 alone, was consistently higher in the low-risk group (Fig. 7B,C). However, the analysis of anti-PDL1 treatment alone (P = 0.170) did not reach statistical significance (Fig. 7D). This suggests that low-risk patients may respond better to anti-CTLA4 and/or anti-PDL1 immunotherapy. Recently, research has found a link between tumor TLS and outcomes in several tumor types. In line with these discoveries, our review of TCGA-LUAD dataset showed that LUAD patients with high TLS scores had more favorable outcomes than those with low scores (Fig. 7F). We also noticed that the TLS score was higher in the low-risk group compared to the high-risk group (Fig. 7E).

In our assessment of the relationship between risk scores and sensitivity to chemotherapy, we measured the IC50 for some widely used chemotherapeutic medicine. Our findings showed that the high-risk group was more sensitive to drugs like Cisplatin, Vinblastine, Cytarabine, Vinorelbine, Bexarotene, Cetuximab, Docetaxel, and Doxorubicin than the low-risk group (Fig. 8A–P).

Through differential gene analysis of tumor tissues and normal tissues, 13,995 DEGs (|logFC|> 1.5, p-value < 0.050) (Fig. 8Q, Supplementary Fig. S3C) were identificated. By cross-referencing with the 27 lncRNAs that form our prognostic model, we pinpointed LINC01003. Supplementary Fig. S4A presents a heatmap demonstrating the expression levels of LINC01003 across different NSCLC datasets and cell types. The results indicate that LINC01003 is differentially expressed, with notable high expression in monocytes/macrophages and endothelial cells across several datasets, suggesting its potential involvement in these cell types within the NSCLC tumor microenvironment. Supplementary Figure S4B further illustrates the expression profile of LINC01003 in different cell populations from the GSE143423↗ dataset. The violin plot shows significant expression of LINC01003 in malignant cells, compared to other cell types, indicating its potential role in tumor progression.

To decipher the LINC00857 related regulatory mechanisms, we constructed a lncRNA-miRNA-mRNA network (Supplementary Fig.C). This network illustrates the intricate interactions between LINC00857 and various miRNAs and mRNAs. In this network, LINC00857 acts as a central regulatory hub, potentially influencing gene expression by sequestering multiple miRNAs, such as hsa-miR-4709-5p, hsa-miR-760, and hsa-miR-340-5p. These miRNAs, in turn, are connected to a wide array of target genes, including YWHAZ, BCL2L2, PTEN, and MYC, which are critical in cellular processes such as cell cycle regulation, apoptosis, and signal transduction. S4

The public datasets were obtained from TCGA (https://portal.gdc.cancer.gov/↗) and GEO (https://www.ncbi.nlm.nih.gov/geo/↗). GEO Accession Numbers: GSE31210↗, GSE30219↗, and GSE50081↗.