A genetic marker model related to disulfidptosis to predict survival and immunotherapy response in lung adenocarcinoma

The RNA-sequencing (RNA-seq)-based transcriptome profiling data, clinical information and somatic mutation data of over 500 LUAD patients were downloaded from The Cancer Genome Atlas (TCGA) database. Normal control samples were excluded for further analysis. LUAD cases with insufficient information about survival time, age, and tumor stage were also removed.

We obtained 25 disulfidptosis-related genes based on previous studies (Liu et al., 2023; Yang et al., 2023; Zhao et al., 2023). Pearson correlation analysis was performed to identify lncRNAs that exhibit co-expression patterns with the disulfidptosis-related genes, with the absolute value of correlation coefficient >0.4 and p < 0.001 as the screening threshold. These lncRNAs were defined as disulfidptosis-related lncRNAs.

A total of 507 LUAD samples with survival information were randomly divided into two groups, one for model construction (the training group, n = 254) and one for model validation (the test group, n = 253). In the training group, disulfidptosis-related lncRNAs that were associated with patients’ overall survival were obtained by performing univariate Cox regression analysis. After LASSO regression analysis to determine lncRNAs with minimum deviation, a prognostic model based on eight disulfidptosis-related lncRNAs was established through multivariate Cox regression analysis. The risk score was the sum of products of the expression value of each of the eight lncRNAs and its regression coefficient, risk score=∑i=0nβiExpi (Zhu et al., 2020). Based on the median risk score, patients were grouped into high- and low-risk subgroups, and survival analysis was carried out to evaluate the significance of the prognostic model. Samples in the test group were used to validate the reliability of this prognostic model. Multivariate Cox regression analysis was conducted to evaluate whether the risk score derived from the model is an independent prognostic factor of LUAD patients.

Differentially expressed genes (DEGs) were screened between the high and low risk groups, according to the screening criteria: |log2 fold change| > 1 and false discovery rate (FDR) < 0.05. After that, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were carried out to gain insights into possible molecular events that distinguish between high- and low-risk groups. GO terms or KEGG pathways were considered significantly enriched when FDR was less than 0.05. With a focus on Gene Ontology gene sets, gene set enrichment analysis (GSEA) was also performed based on gene expression profiles between the two groups. A gene set was considered enriched when p-value and FDR were less than 0.05 and 0.25, respectively.

The ESTIMATE R package was used to analyze the abundances of infiltrating stromal and immune cells in LUAD tissues using gene expression data (Yoshihara et al., 2013). The ESTIMATE algorithm generates three scores based on single sample GSEA, including stromal, immune and estimate scores. Their differences between high and low risk LUAD groups were compared.

The CIBERSORT tool (Chen et al., 2018) was further employed to estimate the abundances of 22 immune cell types in each of the LUAD samples. In addition, single sample GSEA was performed using immune-related gene sets to evaluate multiple immune functions of each sample, and the activities of these immune functions were compared between two risk groups.

According to the total number of somatic base substitutions, the tumor mutational burden (TMB) and mutation frequencies in each sample were calculated. Differences in TMB between the high- and low-risk groups of patients were analyzed. According to the median TMB score, LUAD patients were divided into two groups and survival analysis was performed to explore the influence of TMB on patients’ overall survival. The combined effect of TMB and risk score on patient prognosis was also investigated.

We exploited the Tumor Immune Dysfunction and Exclusion (TIDE) platform (Jiang et al., 2018) to predict LUAD patient response to anti-PD-1 and anti-CTLA4 immunotherapy. TIDE prediction scores are negatively associated with immunotherapy response. Differences in response to immunotherapy between the high- and low-risk groups of patients were analyzed by comparing the TIDE scores.

The oncoPredict R package was used for predicting drug sensitivity in LUAD patients based on gene expression data. The required training sets were derived from the Genomics of Drug Sensitivity in Cancer database (GDSC) and downloaded from oncoPredict’s Open Science Framework (https://osf.io/c6tfx/↗) (Maeser et al., 2021). We used the calcPhenotype function to obtain drug sensitivity scores of each patient. Differences in response to multiple drugs between the high- and low-risk groups were compared based on the drug sensitivity scores.

Data analysis was performed in R (4.2.2). Two-tailed Student’s t-test was used to compare statistical differences between two groups. The Kaplan-Meier estimate and log-rank test were used for survival analysis. Unless otherwise indicated, differences were considered statistically significant when p < 0.05.

RNA-seq data of LUAD patients were downloaded from TCGA. According to the annotation of gene type, protein coding mRNAs and lncRNAs were distinguished. To identify lncRNAs implicated in disulfidptosis, Pearson correlation analysis was conducted based on expression levels of lncRNAs and 25 disulfidptosis-related genes. Following stringent screening criteria (|Pearson R| > 0.4 and p < 0.001), 127 lncRNAs were screened out and their expression were correlated with 20 of the 25 disulfidptosis-related genes (Figure 1A; Supplementary Table S1).

We randomly divided 507 LUAD patients with survival information into two groups, the training group was for model construction and the test group for model validation. Among the 127 lncRNAs, only those associated with LUAD patient survival were considered to be enrolled in model construction. We performed univariate Cox regression analysis to exclude lncRNAs that showed no significant effect on survival, and 14 disulfidptosis-related lncRNAs were left. Of them, three were prognostically favorable lncRNAs and eleven were prognostically unfavorable ones (Figure 1B). Based on these 14 prognosis-associated disulfidptosis-related lncRNAs and using LASSO Cox regression analysis, a prognostic model comprised of 8 disulfidptosis-related lncRNAs was further established (Figures 1C, D). Then we assigned each patient a risk score per the formula of the prognostic model: risk score = (0.433 * KTN1-AS1 expression value (EV)) + (0.099 * AL365181.3↗ EV) + (0.274 * MANCR EV)—(0.604 * LINC01352 EV)—(0.404 * AC090559.1↗ EV) + (0.425 * AC093673.1↗ EV) + (0.374 * AP001094.3↗ EV)—(0.340 * MHENCR EV). The expression correlations between the 8 lncRNAs and 25 disulfidptosis-related genes were shown in Figure 1E, AC093673.1↗ and AL365181.3↗ showed positive associations while AP001094.3↗ and MHENCR showed negative correlations with most disulfidptosis-related genes. According to the median risk score, the training group was divided into high-risk and low-risk groups. As expected, patients in the high-risk group had shorter overall survival time, demonstrating the prognostic significance of the 8 lncRNAs-based model (Figure 2A). Similar survival analysis results were observed in the test group and after combination of the two groups (Figures 2B, C), which indicates that our prognostic model is reliable. Moreover, the 8 disulfidptosis-related lncRNAs exhibited consistent expression patterns per the risk scores between the training and test groups (Figure 2D).

We merged the training and test groups into one and divided it into two groups. Each patient’s risk score and survival state were shown in Figures 2E, F a high risk score was positively correlated with an increased probability of death. In addition to predicting dismal overall survival, a high risk score also indicated poor progression-free survival (Figure 2G).

We next asked whether our prognostic model was interfered by other clinical factors. We enrolled four clinical features of LUAD patients, including age, gender, tumor stage and risk score, for Cox regression analysis. According to univariate analysis, tumor stage (HR: 1.639, 95% CI:1.426–1.884, p < 0.001) and the 8 disulfidptosis-related lncRNAs-based risk score (HR: 1.225, 95% CI: 1.155–1.299, p < 0.001) are two hazardous factors that affect prognosis (Figure 3A). Moreover, the multivariate Cox regression analysis suggested that risk score (HR: 1.245, 95% CI: 1.167–1.328, p < 0.001), together with tumor stage (HR: 1.647, 95% CI: 1.428–1.900, p < 0.001), are independent prognostic factors (Figure 3B). To further evaluate the predictive accuracy of the lncRNA-based prognostic model, ROC curve analysis was performed. The AUC values for 1-, 3-, and 5-year survival are 0.703, 0.673, and 0.654, respectively, indicating high accuracy of our prognostic model (Figure 3C). Of note, our prognostic model is almost as accurate as tumor stage in prognosis prediction, as reflected by similar AUC values (0.673 versus 0.687) (Figure 3D). These results suggest that our disulfidptosis-related lncRNAs-based model can serve as an independent and accurate prognostic indictor.

The practicability of the lncRNA-based model in prognosis prediction of LUAD patients with the same disease stage.

Since both our prognostic model and tumor stage are good survival predictors, we wondered what are the advantages of our model as compared with tumor stage in terms of prognosis prediction. To that end, LUAD patients were divided into early stage group (stage I or II) and advanced stage group (stage III or IV) according to the disease stage. It was found that in both groups, LUAD patients with high risk scores had a poorer overall survival rate than patients with low risk scores (Figures 3E, F). These results suggest that our prognostic model can distinguish between patients at high and low risk, even at the same disease stage.

To further gain insights into the biological differences between the high- and low-risk groups, we performed differentially expressed gene analysis and identified 643 DEGs between the two groups. GO enrichment analysis results indicated that these DEGs are associated with microtubule-based movement, humoral immune response, cilium movement, and other biological processes (Figure 4A). KEGG pathway analysis showed that the DEGs are involved in systemic lupus erythematosus and neutrophil extracellular trap formation (Figure 4B). Furthermore, GSEA that incorporates transcriptome data was carried out. Our analysis showed that nucleosome assembly, DNA packaging complex, nucleosome, protein DNA complex, and structural constituent of chromatin are the top five significantly enriched terms in the high-risk group, while in the low-risk group, B cell receptor signaling pathway, complement activation, immunoglobulin complex, T cell receptor complex, and immunoglobulin receptor binding are the top five significantly enriched cellular processes (Figures 4C, D).

Tumor immune microenvironment plays key roles in determining tumor progression. Considering the above GSEA result that showed enrichment of immune related processes in the low-risk group, we speculated that the tumor immune microenvironment is different between the high- and low-risk LUAD groups. According to the ESTIMATE algorithm, the immune scores are significantly lower in the high-risk group than in the low-risk group (Figure 5A), indicating less infiltration of immune cells in high-risk LUAD. We next used CIBERSORT approach to investigate the abundances of diverse immune cell types in the LUAD tissues. As shown in Figure 5B, high-risk LUADs have less infiltration of monocytes, resting Dendritic cells and resting mast cells, but increased infiltration of M0 macrophages. In addition, we analyzed multiple immune functions between the two LUAD groups. Strikingly, among the 29 kinds of immune functions, 25 showed lower function scores in high-risk LUADs than in low-risk LUADs, such as B cells, CD8⁺ T cells, and cytolytic activity (Figure 5C). Together, these results suggest that high-risk LUADs, as classified by our disulfidptosis-related lncRNA model, may have compromised immune responses in their tumor microenvironment, resulting in tumor progression and worse overall survival.

TMB, the number of somatic mutations per megabase of genomic sequence, is a potential predictive biomarker in many solid tumors (Sha et al., 2020). We analyzed and compared gene mutation frequency and TMB between the two LUAD groups. High-risk LUADs exhibited significant higher TMB than low-risk LUADs (Figure 6A). The TMB score of each LUAD sample and the top 20 most frequently mutated genes as well as their mutation types were shown in Figures 6B, 6C. Mutation frequencies of almost all of these genes were higher in high-risk LUADs than in low-risk LUADs. TP53 and TTN were mutated in over half of the high-risk LUADs. High TMB can bring benefits for LUAD patients in their survival, and the overall survival rate is much higher in high TMB patients than in low TMB patients (Figure 6D). We next investigated the effects of risk score and TMB on LUAD patients’ overall survival, patients were divided into four subgroups based on these two factors. We found that patients with high TMB and low risk scores exhibit the best prognosis, their 10-year survival rate is around 60%. In contrast, patients with low TMB and high risk scores have the poorest prognosis, the 5-year survival rate is merely about 25%. There are no significant differences in overall survival between high-risk patients with high TMB and low-risk patients with low TMB, and the survival rate of these patients is between the other two subgroups (Figure 6E).

Drug resistance is a major cause of cancer relapse and cancer-related death. Immune checkpoint inhibitors (ICBs) have exhibited impressive therapeutic effects in certain cases of NSCLC. To explore the role of our disulfidptosis-related lncRNA model in predicting response to immunotherapy, we analyzed the correlations between LUAD risk score derived from the model and TIDE score. High-risk LUAD patients have significant lower TIDE scores (Figure 7A), suggesting that immune checkpoint inhibitors are more effective in these patients. Since our analysis showed that high-risk LUADs have high TMB (Figure 6A), our results are in line with previous finding that higher TMB was associated with clinical efficacy of anti-PD-1 therapy (Rizvi et al., 2015). In addition, the associations between LUAD risk score and sensitivity to other antitumor agents were also investigated. Drug sensitivity scores were generated by the calcPhenotype function in the oncoPredict package, based on gene expression data of LUAD patients and preprovided training datasets. As compared with patients in the low-risk group, patients in the high-risk group are less sensitive to diverse types of antineoplastic drugs, including EGFR tyrosine kinase inhibitors (gefitinib, erlotinib, and AZD3759) (Figure 7B), MEK and ERK inhibitors (Trametinib, PD0325901, Ulixertinib and ERK_6604) (Figure 7C), inhibitors of cell cycle-related kinases (AZD7762, BI-2536 and MK-1775) (Figure 7D), MET inhibitors (Savolitinib, Foretinib and Crizotinib) (Figure 7E), and drugs that disturb genome integrity (Talazoparib, AZD6738, VE821 and GDC0810) (Figure 7F). These results suggest that our disulfidptosis-related lncRNA model is a potential tool to predict response of LUAD patients to ICBs and other common antineoplastic drugs.

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/. Supplementary Material