Construction and Validation of a Novel Disulfidptosis-Related lncRNA Signature to Predict Immune Landscape and Precise Treatment in Stomach Adenocarcinoma

1. Introduction

Gastric cancer is one of the most commonly diagnosed cancers in the world with nearly 1 million new cases and more than 720,000 deaths worldwide each year[1]. Worldwide gastric cancer statistics from 2020 showed 770,000 deaths and 1.1 million newly diagnosed cases. These numbers are projected to climb significantly, with forecasts indicating 1.8 million new cases and 1.3 million deaths by 2040[2]. Early gastric cancer typically presents with subtle and nonspecific symptoms, resulting in more than 60% of patients being diagnosed after metastasis[3]. This late diagnosis contributes to poor prognosis, with only 5% of patients with metastatic GC surviving beyond five years[4,5].

Luckily, the development of immunotherapy such as immunocheckpoint inhibitors (ICIs) shed light on the treatment of advanced and recurrent gastric cancer[6]. However, a substantial proportion (40-60%) of patients treated with immune checkpoint inhibitors demonstrated resistance to therapy[7]. Thus, it is of great significance to find reliable biomarkers to predict the survival of stomach adenocarcinoma (STAD) patients, identify high-response individuals and drugs, and ultimately improving the survival of GC patients.

Disulfidptosis is a kind of original regulated cell death that is different from ferroptosis and cuproptosis. Liu et al. first defined it as a unique form of cell death induced by abnormal accumulation of disulfide bonds in high-SLC7A11 cells. Here we focus on drlncRNAs related to 24 DRGs (disulfidptosis-related genes) mentioned in the literature including GYS1, NDUFS1, OXSM, LRPPRC, NDUFA11, NUBPL, NCKAP1, RPN1, SLC3A2, SLC7A11,INF2, CD2AP, ACTN4, PDLIM1, IQGAP1, DSTN, CAPZB, ACTB, MYL6, MYH9, MYH10, TLN1, FLNA, FLNB. Among these 24 DRGs, the last 14 are genes that encode actin-related proteins and upregulate following the glucose starvation[8]. Given that upregulated SLC7A11 is common in human cancers, the discovery of disulfidptosis expands the framework of programmed cell death and may lead to new treatment and therapeutic targets for multiple cancers [9,10,11].

Accumulating research have reported that lncRNAs serving as an "oncogene" could promote the proliferation, invasion and metastasis of gastric cancer cells, and can be used to predict survival, guide precise clinical therapy of STAD patients as potential biomarkers[12,13,14,15]. However, disulfidptosis-related lncRNA in STAD have barely been reported so far.

The present study aims to develop a novel disulfidptosis-related lncRNA signature to predict the prognosis and immune landscape of STAD. We believe that this research will provide valuable reference in identifying effective drugs and individuals with high response to immunotherapy.

2. Materials and methods

2.1. Date source

The flowchart of the present study was showed in Figure 1. We downloaded the transcriptome data of 448 STAD patients (412 tumor samples and 36 normal samples) and clinical data of 443 STAD patients from the Cancer Genome Atlas (TCGA) database (https://tcga-data.nci.nih.gov/tcga/). To separate mRNAs and lncRNAs, Strawberry Perl software (version 5.30.0-64bit) was downloaded from https://strawberryperl.com/. 24 DRGs were collected from the literature by Liu et al. mentioned before and Pearson correlation analysis was then performed to filter out drlncRNAs. In order to minimize errors, patients with missing survival information were excluded and 407 STAD patients remained eventually. The simple nucleotide variation (SNV) data of 434 STAD patients were also obtained from TCGA database and then used to calculate tumor mutational burden (TMB). The data of STAD microsatellite status were obtained from TCIA website (https://www.tcia.at/ ). Finally, we acquired the tumor immune dysfunction and exclusion (TIDE) score from http://tide.dfci.harvard.edu.

2.2. The development and verification of the arlncRNAs signature

The entire set containing 407 STAD patients were randomly divided into train set and test set at 1:1 ratio. Then, statistical analysis was performed to confirm the reliability of the grouping. We first developed the drlncRNA prognostic model in the train set. By conducting univariate cox regression analysis, we obtained 29 drlncRNAs that were significantly associated with overall survival (OS) (p<0.05). To avoid overfitting, we next executed LASSO (least absolute shrinkage and selection operator)-penalized Cox analysis with 10-fold cross-validation in the condition of p<0.05 and screened out 7 drlncRNAs. Finally, a prognostic model composed of 3 drlncRNAs was constructed after performing multivariate Cox regression analysis. The test set and the whole set were then used to validate the performance of the model. The risk score was calculated according to following formula:

Risk score= lncRNA1exp × coef1 + lncRNA2exp × coef2 +… + lncRNAnexp × coefn.

lncRNAnexp: expression of lncRNAn; coefn: risk coefficient of lncRNAn based on the model

All STAD patients were then segmented into high and low risk groups (high and low disulfidptosis score groups) based on the median calculated risk score as the cutoff value. To assess the accuracy of our prognostic model, we performed independent prognostic analysis to validate whether our drlncRNA signature can predict the prognosis of GC patients independent of age, gender, grade and stage. Subsequently, we used R packages "survival", "survminer" and "timeROC" to plot ROC curves of the train set, test set, entire set and different clinicopathological characteristics respectively, and corresponding AUC (area under curve) value were also calculated. In addition, to further verify the performance of our prognostic model, Kaplan–Meier curves of different clinicopathological variables such as age, grade, stage, size of the primary tumor (T stage), metastasis of regional lymph node (N stage) and distant metastasis of tumor (M stage) were obtained

2.3. Independent prognostic analysis and ROC curve plotting

Independent prognostic analysis was utilized to verify whether our model can predict the outcomes of patients independent of age, gender, stage and grade. The accuracy of the prognostic signature was assessed through ROC curves and AUC (area under curve) values. (R packages "limma", "scatterplot3d", "survival", "survminer" and "timeROC" were used in this progress)

2.4. Nomagram and Calibration

Considering clinicopathological factors and risk score, a predictive nomogram which was designed to predict 1-, 3-, and 5-year OS was created utilizing R packages "survival", "regplot" and "rms". In the meantime, 1-, 3-, and 5-year calibration curves were also plotted to evaluate the predictive precision of the nomogram.

2.5. Gene set enrichment analysis

The Gene Ontology (GO) database classified the functions of genes into three categories, including biological process (BP), cellular component (CC) and molecular function (MF). Using R package "clusterProfiler", we conducted GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis to identify potential enriched functions and signaling pathways, and visualized with R packages "GOplot" and "ggplot2".

2.6. Investigation of TME, immune infiltration and immune checkpoints

Based on the result of GO and KEGG, we conducted immune-related analyses. ESTIMATE algorithm is usually used to calculate immune, stromal, and comprehensive scores of TME[16]. Here we executed ESTIMATE algorithm with R packages "limma" and "ggpubr" to explore the potential correlation between risk score and TME score. The infiltration difference of distinct immune cell subpopulations between high and low disulfidptosis score groups was investigated via 7 platforms including XCELL, TIMER, QUANTISEQ, MCPCOUNTER, EPIC, CIBERSORT and CIBERSORT-ABS, and visualized in a bubble chart. R packages "limma", "scales", "ggplot2", "ggtext", "reshape2", "tidyverse" and "ggpubr" were utilized in this progress. Subsequently, we calculated the scores of 29 immune cells and functions in 448 STAD patients via single-sample gene set enrichment analysis (ssGSEA) using R packages "GSVA", "limma" and "GSEABase", and visualized in two boxplots with R packages "limma", "ggpubr" and "reshape2". Finally, we conducted immune checkpoint analysis to explore differentially expressed immune checkpoints between high and low disulfidptosis score group.

2.7. TMB, MSI and TIDE

We first downloaded the SNV data of 434 STAD patients from TCGA database and calculated the TMB of all patients utilizing Strawberry Perl software. Then we plotted the waterfall diagrams of 15 genes with the highest mutation rates in high and low risk groups through "maftools" package. Kaplan–Meier curves of “H-TMB + high risk”, “H-TMB + low risk”, “L-TMB + high risk” and “L-TMB + low risk” were obtained by R packages "survival" and "survminer". In addition, according to the extracted STAD microsatellite status data, we drew a barplot and a boxplot to illustrate the relationship between MSI and risk score. Finally, we compared the difference of TIDE score between high and low disulfidptosis score groups.

2.8. Drug sensitivity analysis

The half maximal inhibitory concentration (IC50) values of multiple drugs were calculated and compared between high and low disulfidptosis score groups utilizing R package "pRRpphetic", and the results were presented as boxplots.

3. Results

3.1. Data source and processing

Using Strawberry Perl software, we processed the TCGA transcriptome data and gained the expression data of 19938 mRNAs and 16876 lncRNAs. By combining the mRNAs expression profile with 24 disulfidptosis-related genes, we obtained the expression data of 24 DRGs. Then, based on the criteria of |Pearson R|>0.4 and p < 0.001, Pearson correlation analysis between 16876 lncRNAs and 24 DRGs was conducted and 367 drlncRNAs were finally identified (Figure 2A).

3.2. Construction and validation of the drlncRNA prognostic model

By combining 367 drlncRNAs expression data with 407 GC patients’ clinical information, we acquired the drlncRNAs expression data profile and survival data for 407 samples, which was used for prognostic model building. We first performed univariate cox regression analysis and preliminary screening out 29 drlncRNAs (Figure 2B, p<0.05). Subsequently, LASSO-penalized Cox analysis was executed to avoid overfitting and 7 drlncRNAs significantly related to overall survival (OS) were obtained (Figure 2C,D, p<0.05). Eventually, the drlncRNA risk signature comprising 3 lncRNAs (AC107021.2, AC016394.2 and AC129507.1) and corresponding risk coefficient was acquired (Figure 2E, p<0.05). In the meantime, the relationship between 24 DRGs and 3 drlncRNAs is presented in Figure 2F. According to the previous calculation formula, risk score = AC107021.2exp × 0.66000986856416 + AC016394.2exp × -0.435804339506525 + AC129507.1exp × 0.636605050383588.

A total of 407 STAD samples were then allocated into two groups randomly at a 1:1 ratio with 204 samples in the train set and 203 samples in the test set. Based on the median risk score as cutoff value, patients were divided into high and low risk group. According to the risk score and survival status of train set patients manifested in Figure 3A,B, we can conclude that the mortality of patients was positively correlated with risk score. Obviously, the risk heatmap of 3 drlncRNAs constructed the model (Figure 3C) showed that the upregulated AC107021.2 and AC129507.1 were significantly correlated with higher disulfidptosis score, while AC016394.2 was highly expressed in low disulfidptosis score group. The Kaplan–Meier curve of train set (Figure 3D, p<0.001) indicated that the survival time of patients in the high risk group were notably lower than that in the low risk group, which was consistent with our expectations. To demonstrate the predictability of the arlncRNA model, we conducted the same analysis in the test set and whole set. The results of the test group (Figure 3E-H) and whole group (Figure 3I-L) were consistent with the TCGA training group.

Independent prognostic analysis showed that the 3 drlncRNAs model was a reliable prognosis predictor, independent of clinical features such as age, stage (Figure 3M, N). Subsequently, the area under the ROC curves (AUC values) was calculated in train group (1-year AUC = 0.741, 3-year AUC = 0.654, and 5-year AUC = 0.781; Figure 3O), test group (1-year AUC = 0.638, 3-year AUC =0.640, and 5-year AUC = 0.707; Figure 3P) and entire group (1-year AUC = 0.692, 3-year AUC = 0.646, and 5-year AUC = 0.723; Figure 3Q). The AUC values of ROC curves for different clinical features were also calculated, and the ROC curve based on our risk model showed the greatest AUC value (Figure 3R), which illustrated that our risk model as a prognosis predictor was considerably accurate.

In order to further verify the accuracy of the model, we plotted K-M curves of patients in high and low risk groups based on different clinical features including age , grade, stage, T stage, N stage and M stage.We can recognize that our prognostic model can accurately predict the prognosis of patients independent of clinical variables (Figure 4). For patients in M1 stage, there was no significant difference in survival between our high and low risk groups, which may be due to the lack of sufficient samples.

3.3. Nomogram and calibration curves

For the purpose of predicting the survival time of patients utilizing our drlncRNA signature, we plotted a nomogram that can predict 1-, 3-, and 5-year survival rates. (Figure 5A). In the meantime, calibration curves were plotted to confirm the accuracy of our nomogram (Figure 5B).

3.4. Functional analysis of the model

The circle chart and corresponding notes (Figure 6A,B), column chart (Figure 6C) and bubble chart (Figure 6D) of GO and KEGG enrichment analysis showed that these three drlncRNAs were closely related to sulfur compound binding, canonical WNT signaling pathway, cell-substrate adherens junction and cAMP signaling pathway.

3.5. Immune infiltration status

In order to explore the difference of tumor microenvironment between high and low disulfidptosis score groups, we utilized ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumour tissues using Expression data) algorithm to score the purity of tumors. According to the boxplots, we recognized that the score of the high risk group was higher than those of the low risk group in StromalScore(p=2.5e-10), ImmuneScore(p=0.039) and ESTIMATEScore(p=6.1e-06) (Figure 7A). Subsequently, we analyzed the differences in immune cell subpopulations in TME between high and low risk groups. The bubble chart showed that the high risk group corresponded to a higher degree of immune cells infiltration, such as CAF(Cancer associated fibroblast) in XCELL, MCPCOUNTER and EPIC, Macrophage 2 in QUANTISEQ, CIBERSORT and CIBERSORT-ABS (Figure 7B). In addition, the boxplot of 29 immune cells and functions between high and low risk groups showed that high risk group had higher levels of MAST cells infiltration and lower levels of MHC class 1 expression (Figure 7C,D). ICP analysis revealed that the expression levels of ICPs such as PDCD-1, CD274 (PD-L1) and CTLA-4 were upregulated in the low risk group (Figure 7E).

3.6. Immunotherapy response analysis

In this part, TMB, MSI, TIDE analyses were conducted. First of all, we compared the TMB of the high and low risk groups and found that the patients in the low risk group had higher TMB (Figure 8A). In order to further explore the difference of mutant genes between high and low risk groups, we plotted the corresponding waterfall diagram. Among the 15 genes with the highest mutation rate, only TP53 gene had a higher mutation rate in the high risk group than in the low risk group, and we noticed that the mutation rate of ARID1A gene in the high risk group was significantly lower than that in the low risk group (Figure 8B). Next, we explored the relationship between TMB, patient survival, and risk score (Figure 8C). We found that the survival of GC patients with “H-TMB” was better than that of patients with “L-TMB”(p=0.009), and the level of TMB could predict the prognosis of patients independent of risk score. In addition, when we explored the relationship between MSI and risk score, we found that the proportion of patients with “MSI-H” in the low risk group (32%) was much higher than that in the high risk group (9%)(Figure 8D). The risk score of patients with “MSI-H” was much lower than patients with “MSS”(p=4.7e-10) and “MSI-L” (p=1.7e-05) (Figure 8E). Finally, TIDE (Tumor Immune Dysfunction and Exclusion) analysis revealed that patients in the high risk group had higher TIDE scores (Figure 8F).

3.7. Drug sensitive analysis

By performing drug sensitivity analysis, we identified a series of drugs (e.g., gemcitabine, ABT.888) that exhibited lower IC50 values in GC patients with low disulfidptosis scores (Figure 9).

4. Discussion

A literature search revealed 24 genes involved in the disulfidptosis pathway. Co-expression analysis enabled us to identify lncRNAs associated with disulfidptosis. Subsequently, using Cox regression analysis combined with LASSO penalty analysis, we constructed a prognostic model comprising three disulfidptosis-related lncRNAs. AC107021.2 has emerged as a diagnostic marker for gastric and lung adenocarcinomas[17,18], while AC016394.2 and AC129507.1 show diagnostic potential in both gastric adenocarcinoma and prostate carcinoma[19,20]. However, the relationship between disulfidptosis-associated lncRNAs and immune landscape in STAD, along with their prognostic significance, remains to be elucidated. Our study investigated drlncRNAs in gastric cancer and established a three-drlncRNA prognostic signature to predict the prognosis and provide precise and individual clinical treatment guidance for patients with GC.

Through independent prognostic analysis, we demonstrated that our risk model's ability to predict patient outcomes remains significant when accounting for other clinical parameters. By drawing the ROC curve, we calculated AUC (Area Under the Curve) values, which validate our risk model's accuracy in predicting patient survival outcomes across 1-, 3-, and 5-year timeframes. Kaplan-Meier curves of different clinical features, including age, grade, stage, T stage, N stage, and M stage, illustrated that our risk model can predict patient outcomes regardless of clinical variables. Through nomogram and calibration curves, our model enables individualized prediction of patient survival at 1-, 3-, and 5-year timepoints.

As anticipated, enrichment analysis revealed a strong association between disulfidptosis-related lncRNAs and sulfur compound binding. The identification of disulfidptosis revealed a promising therapeutic strategy for cancer treatment, wherein sulfur compound binding orchestrates this cell death process through the activation of specific signaling cascades[10]. GO and KEGG analysis also revealed other functions and pathways that are closely related to drlncRNA, such as GPCR(G protein coupled receptor) binding, cAMP signaling pathway. A recent study revealed that activation of the Gαs-PKA signaling pathway drives CD8+ T cell dysfunction and confers resistance to immunotherapy[21]. cAMP exhibits dual regulatory effects on tumor cell survival and proliferation through its interactions with multiple immune components of the TME, particularly T cells and tumor-associated macrophages (TAMs)[22,23,24,25,26,27].

Given that GO and KEGG enrichment analyses revealed associations between drlncRNAs and immune function, we proceeded to investigate the immune landscape by comparing immune cell infiltration patterns between high- and low-risk groups. The elevated TME scores observed in high-risk patients indicate enhanced immune cell infiltration and enriched stromal cell content in their tumor microenvironment. Macrophages exhibit two distinct polarization states: M1 and M2[28]. In our analysis of immune cell infiltration patterns within the TME, we observed contrasting infiltration levels between M1 and M2 across the high- and low-risk groups: high levels of M1 infiltration were observed in the low-risk group, while high levels of M2 infiltration were observed in the high-risk group. While M1 play a crucial role in anti-tumor immunity, M2 facilitate tumor progression by promoting immune escape, angiogenesis, and extracellular matrix remodeling in tumor cells[29,30]. This may explain the poor survival of patients in the high disulfidptosis score group. Besides, the distinctions in 16 immune cells and 13 immune-related pathways across two subgroups were investigated using single-sample Gene Set Enrichment Analysis (ssGSEA). In the high-risk group, we observed elevated levels of mast cells accompanied by downregulation of MHC class I (Major Histocompatibility Complex I) molecules. A recent study shows that mast cells contribute to resistance against anti-PD-1 immunotherapy, and targeted depletion of mast cells enhances the therapeutic efficacy of immune checkpoint blockade[31]. Additional research shows that cancer cells can escape immune detection by reducing the expression of MHC-I. This decrease in MHC-I expression represents a key pathway for both inherent and adaptive resistance to immunotherapeutic interventions in cancer patients[32]. Unlike traditional therapy, this approach doesn't directly target tumor cells. Instead, it works by alleviating immunosuppression and stimulating the body's natural anti-tumor immune response, demonstrating remarkable efficacy across various treatment-resistant tumors. When combined with chemotherapy, targeted therapy, radiotherapy, and other treatment modalities, it enhances overall therapeutic outcomes[33,34,35,36]. Elevated expression of immune checkpoints, including PDCD1, PD-L1, and CTLA-4, in the low-risk group suggests a potentially enhanced response to immunotherapy.

TMB refers to the total number of nonsynonymous somatic mutations found in the cancer genome. As a biomarker for ICI therapy response, MSI-H demonstrates predictive value across diverse tumor types, potentially due to its role as a mechanistic mediator of immunotherapy outcomes[37]. Accumulating evidence suggests that TMB levels serve as a potential predictor of patient responsiveness to immune checkpoint inhibitor therapy, particularly when combined with PD-1 expression and MSI status for enhanced prognostic accuracy[38,39,40,41]. Our analysis demonstrated that patients with low disulfidptosis scores exhibited significantly higher tumor mutational burden (TMB) (p=1.2e-09) and a greater proportion of MSI-H cases (32%) compared to the high-score group. These findings suggest that patients with low disulfidptosis scores may achieve enhanced therapeutic responses to immunotherapy. Moreover, we found patients with elevated TMB showed improved survival outcomes regardless of disulfidptosis scores. To comprehensively analyze the mutational landscape differences between high- and low-risk groups, we generated waterfall plots depicting the 15 most frequently mutated genes. Notably, among the 15 most frequently mutated genes, only the tumor suppressor gene TP53 showed a higher mutation frequency in the high-risk group compared to the low-risk group. In contrast, we observed that ARID1A mutations were significantly less prevalent in the high-risk group, establishing a distinct mutational pattern between the two risk groups. TP53 represents the most frequently mutated gene across human cancers[42]. Its mutations not only compromise its tumor-suppressive functions but also confer oncogenic properties to the mutant p53 protein[43]. Studies have demonstrated that TP53 mutations correlate with poor clinical outcomes in cancer patients[44]. This observation, on the other hand, provides a mechanistic explanation for the poor survival outcomes observed in patients with high disulfidptosis scores. Multiple clinical trials have demonstrated that ARID1A-mutated solid tumors exhibit enhanced responsiveness to immune checkpoint blockade (ICB) therapy across various cancer types, independent of MSI status or TMB[45,46,47]. These findings suggest that patients with high disulfidptosis scores may exhibit resistance to immune checkpoint inhibitor therapy.

TIDE (Tumor Immune Dysfunction and Exclusion) analysis further revealed that patients in the high-risk group exhibited a higher likelihood of developing resistance to immunotherapy, suggesting limited therapeutic benefit from immune checkpoint blockade in this subgroup. Moreover, comprehensive drug sensitivity analysis identified multiple potential therapeutic agents that showed enhanced efficacy in the low-risk group, while also revealing drugs prone to resistance in the high-risk group, thereby providing valuable insights for personalized treatment strategies. Recent clinical trials have demonstrated that the combination of gemcitabine plus cisplatin with PD-L1 inhibitors exhibits promising therapeutic efficacy in patients with advanced biliary tract cancer[48]. In another phase I dose-escalation study, the combination of ABT.888(veliparib) with PD-1 inhibitor plus platinum-based doublet chemotherapy demonstrated promising therapeutic efficacy in patients with metastatic or advanced non-small cell lung cancer[49]. Based on these findings, the combination of immune checkpoint inhibitors with conventional chemotherapeutic agents represents a promising therapeutic strategy for gastric cancer.

Additionally, we realized that this study still had limitations. Firstly, the 3 drlncRNAs prognostic model was developed and verified on the basis of the results of these retrospective analyses of TCGA. It is necessary to conduct prospective cohort studies to further validate this risk model. Secondly, more independent immunotherapy cohorts were needed to validate the predictive value of the risk model for immunotherapy response. Lastly, we should not only verify the accuracy and practicality of the above lncRNA signature by collecting a large number of clinical data, but also further explore the mechanism of these lncRNAs in GC.

5. Conclusions

In conclusion, our established disulfidptosis-related lncRNA prognostic model enables accurate prediction of patient survival and facilitates the identification of potentially effective drugs and treatment-sensitive individuals, ultimately contributing to improved survival outcomes in STAD patients.