Unlocking the Potential of Disulfidptosis‐Related LncRNAs in Lung Adenocarcinoma: A Promising Prognostic LncRNA Model for Survival and Immunotherapy Prediction

The datasets for this study were sourced from The Cancer Genome Atlas (TCGA) website (https://www.cancer.gov/ccg/research/genome‐sequencing/tcga↗), which provided the LUAD transcriptome RNA‐seq (STAR‐Counts), clinical, and simple nucleotide variation datasets. We gathered a total of 497 samples that were used for the analysis. To annotate the downloaded data, we utilized the annotation file from the GENCODE website (https://www.gencodegenes.org↗), which allowed for the identification of 16876 lncRNAs. And to validate the results of this study, GSE29013↗; GSE30219↗, GSE31210↗, AND GSE37745↗ were recruited. The use of these standardized and publicly available resources ensured the reliability and reproducibility of our results.

According to Liu et al. [1], the guide RNAs GYS1, LRPPRC, NCKAP1, NDUFA11, NDUFS1, NUBPL, OXSM, RPN1, SLC3A2, and SLC7A1 were relatively fold changed between the glucose‐replete and ‐starved groups, and these 10 RNAs were adopted in our study as disulfidptosis‐related genes for further analyses. The DRLs were identified using the Pearson correlation with a cor‐filter ≥ 0.4 and p < 0.001. Then, 91 lncRNAs were identified to be co‐expressed with the DRGs and retrieved for the next analyses. The limma, dplyr, ggalluvial, and ggplot2 R packages were recruited for analyses and visualization in this section.

Firstly, the univariable regression algorithm was applied to evaluate the candidate prognostic DRLs. Then, the LASSO and SVM‐RFE algorithms were applied to select prognostic lncRNAs from the candidate DRLs [9, 10]. By employing the multivariate Cox regression, a comprehensive risk model comprising the prognostic lncRNAs and the relevant correlation factors was derived. Then, we analyzed the correlation between the DRGs and prognostic DRLs using the relevant package. The process of classification of the high‐ and low‐risk groups was executed through the use of the median risk score from the 5‐DRLs signature model. To validate the model's ability to predict survival, a rigorous analysis was conducted, whereby all patient samples were randomly divided into training and test groups. Employing the Kaplan–Meier regression method and the receiver operating characteristic (ROC) curves coupled with the concordance index (C‐index) provided precise and detailed validation of the robustness of the risk model.

We constructed a nomogram of survival prediction based on the risk score and the clinical parameters in the whole TCGA‐LUAD cohort. The nomogram could be used for evaluating the time‐dependent predictive efficacy of the overall survival rates at 1, 3, and 5 year. Following the development of the predictive nomogram, we developed a calibration curve to demonstrate the nomogram model's prediction ability.

Usually, the dimensional features could be reduced by the principal component analysis (PCA), which is a normally used statistical technique. Thus, we applied the scatterplot3d R package to conduct PCA and judge the classification capability of varied gene clusters.

The relationship in the risk score, tumor mutant burden, and the tumor immune microenvironment was evaluated by the maftools R package, ESTIMATE R package, and CIBERSORT R script v1.03, respectively. The tumor immune dysfunction and exclusion (TIDE) analysis was carried out by uploading the expression matrix online (http://tide.dfci.harvard.edu/↗), and then the results for response prediction were downloaded for the following analyses. Finally, the sensitivities of chemotherapy and potential targeted therapy drugs were evaluated using data from the Genomics of Drug Sensitivity in Cancer (GDSC, https://www.cancerrxgene.org/↗) database and the oncoPredict R package.

Comparative analysis of gene expression between low‐ and high‐risk groups assessed different expression genes (DEGs), and the Gene Ontology (GO) analysis provided valuable insights into the molecular functions, cellular components, and biological roles that were enriched in the DEGs. Meanwhile, the Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis enabled us to identify the specific signaling pathways that were significantly affected by these DEGs. Going a step further, instead of concentrating on DEGs in the above analyses, Gene set enrichment analysis (GSEA) sorts all the genes according to their expression and then gives each gene an ES enrichment score. Thus, GSEA was also applied for further understanding of the involved functions and signaling pathways.

Human lung epithelial BEAS‐2B cells, human lung adenocarcinoma cell line NCI‐H1975, NCI‐1299, and human large cell lung cancer cell line NCI‐H460 were purchased from the Procell Life Science & Technology (Wuhan, China) and were cultured in DMEM and RPMI‐1640 medium (KeyGEN BioTECH, China) containing 10% fetal bovine serum (Gibco, USA) and 1% Penicillin–Streptomycin (KeyGEN BioTECH, China) at 37°C in 5% CO₂.

The total RNA of BEAS‐2B, NCI‐H1975, NCI‐H1299, and NCI‐H460 cells was extracted using the RNA extraction kit (Beyotime, China) and reverse‐transcribed into cDNA using the HiScript II Q SuperMix for qPCR (Vazyme, China) following the users' manual. Then, the qRT‐PCR was carried out using the Taq Pro Universal SYBR qPCR Master Mix (Vazyme, China). The primer sequences were as follows: GAPDH‐forward: TGT ACG CCA ACA CAG TGC TG, GAPDH‐reverse: TCA GGAGGA GCA ATG ATC TTG, AL563181.2↗ ‐forward: CCC AAG AAC CTT AAC CAC CTC A, AL5631812‐reverse: CGT CAC TGG TCT GGC TTT CA, EMSLR‐forward: GTT TCC ACC TAG GAC TAC AGG CT, EMSLR‐reverse: CCC CGC CGA TCC AAT TTC TC, AC092718.4‐forward: TGG AGG GAG GAA GCC ATT CT, AC092718.4‐reverse: GGT TGG TGC TGT TGA GGA GT, ENTPD3‐AS1‐forward: AA GAC ACT TGA GGG GGA GGT, ENTPD3‐AS1‐reverse: GGG TCT CGC TAA CAC CAG AG, AL606489.1↗ ‐forward: AGG AAA GAC ATC AGC AGA GAG C, and AL606489.1↗ ‐reverse: GTA GCT GGG TGG GTG CAT TC. The ΔΔCq method was applied to calculate the relative lncRNA expression.

The vector and pcDNA3.1‐AC092718.4 were purchased from OBiO Technology (Shanghai, China). NCI‐H1975 and NCI‐H1299 cells were seeded in 6‐well plates at a count of 1 × 10⁵ and then cultured overnight. The next day, the mix of 2.5 μg of plasmid DNA (vector or pcDNA3.1‐AC092718.4), Opti‐MEM Medium (Gibco, USA), and Lipo8000 (Beyotime, China) was added into the 6‐well plates, which was prepared according to the users' manual. Then the medium was replaced by the fresh complete medium 24 h after transfection. The RNA was extracted 48 h after transfection and the relative expression of lncRNA AC092718.4 was analyzed by qRT‐PCR.

The effect of lnc‐AC092718.4 on the proliferation of LUAD cells was analyzed by the colony‐forming assay. NCI‐H1975 and NCI‐H1299 cells transfected with vector or pcDNA3.1‐AC092718.4 were seeded in triplicate at a density of 4 × 10³ cells/well in 6‐well plates. After 7 days of culturing in a 37°C incubator with 5% CO2, the medium was removed. Then, the colonies were fixed with 4% polyformaldehyde and stained with 0.5% crystal violet. The plates were photographed using a smartphone and quantified using ImageJ.

The transwell chambers coated with or without Matrigel (Sigma, USA) were applied to explore the invasive ability or migration ability in the H1975 and H1299 cells, respectively. Briefly, cells were cultured in the serum‐free RPMI‐1640 for 24 h before the assay, then the Matrigel was diluted proportionally with precooled serum‐free 1640 medium. The transwell chambers coated with diluted Matrigel were used for detecting the invasive ability, while untreated chambers were used to assay migratory ability. Then the starvation preconditioning cells were counted and seeded in the upper chamber at a count of 1 × 10⁵. RPMI‐1640 medium containing 10% serum was added into the lower chamber and the plates were incubated for 24 h. Finally, Matrigel on the upper side of the insert was removed, and the cells on the lower side were fixed in methanol and stained using the crystal violet solution. The stained cells were counted as invasive/migrative cells.

The proliferation ability was also detected by CCK8 assay. Briefly, cells were cultured in 1 × 10³ quantities in 96‐well plates, and the cell viability was assessed per 24 h over 96 h according to the users' manual of the CCK8 kit (Abbkine, China).

We utilized the R software (version 4.2.2) as our primary analytical tool to perform all necessary statistical analyses. The relevant R packages (limma, dplyr, survival, ggalluvial, ggplot2, care, glmnet, survminer, timeROC, tidyverse, ggExtra, pheatmap, rms, pec, regplot, survcomp, scatterplot3d, colorspace, stringi, circlize, RColorBrewer, ggpubr, org.Hs.eg↗.db, enrichplot, estimate, clusterProfiler, reshape2, e1071, maftools, oncoPredict, and parallel) were recruited during the analyses. Statistical significance was set at p < 0.05.

According to Liu et al. [1], the guide RNAs GYS1, LRPPRC, NCKAP1, NDUFA11, NDUFS1, NUBPL, OXSM, RPN1, SLC3A2, and SLC7A1 were relatively fold changed between the glucose‐replete and ‐starved groups in the clear cell renal adenocarcinoma (ccRCC) cells, and these 10 DRGs were adopted for further analyses in this study. A total of 497 LUAD samples recruiting in the TCGA, database (https://www.cancer.gov/ccg/research/genome‐sequencing/tcga↗), of which 16876 lncRNAs were included. The clinical characteristics are shown in Table 1. The co‐expression analysis was screened between disulfidptosis‐related genes (DRGs) [1] and lncRNAs to identify the DRLs by co‐filter ≥ 0.4 and p < 0.001. Then we got 91 DRLs and most lncRNAs were related to the disulfidptosis‐related gene NDUFA11 (Figure 1A). After merging with survival data, univariable cox‐regression analysis accessed 11 candidate prognostic DRLs showed in the forest diagram (Figure 1B). Next, utilizing the methodology of LASSO regression analysis, a risk model was formulated based on candidate DRLs with the lowest error rates (Figure 1C,D). Subsequently, these candidate RNAs underwent thorough examination via multivariate Cox regression, ultimately culminating in the identification of five prognostic DRL signatures. The risk score was calculated as follows: (0.3223) × AC092718.4 + (0.1775) × AL365181.2↗ + (0.1567) × AL606489.1↗ + (0.3231) × EMSLR + (−0.4627) × ENTPD3‐AS1. The correlation between the 10 DRGs and these five prognostic DRL signatures is shown in Figure 1E. We noted that some correlation tendency was not unique. In detail, AC092718.4, AL365181.2↗, and EMSLR were only positively correlated with the disulfidptosis genes, while AL606489.1↗ and ENTPD3‐AS1 were both positively and negatively correlated with the disulfidptosis genes. The aforementioned results implied that a complex network might be involved in regulating the lncRNAs and genes.

The classification of the high‐ and low‐risk groups was executed through the use of the median risk score from the 5‐DRLs signature model. For validating the prognostic model, the whole cohort was randomly divided into a training group and a test group in a ratio of 1:1. Survival analysis results indicated that LUAD patients with high risk showed worse survival in both training (p < 0.001) and test (p = 0.004) groups (Figure 1F,G). The survival probabilities of the risk scores at 1, 3, and 5 year were denoted by an ROC value of 0.625, 0.666, and 0.666, respectively (Figure 1J). In the whole cohort, the Kaplan–Meier survival analysis curves effectively showcased that the low‐risk group had a better overall survival (p < 0.001, Figure 1H) and progression‐free survival (p = 0.019, Figure 1I) rate in contrast to their high‐risk counterparts. The Supplementary figure vividly presented the distribution of the risk score and survival status, both of which differ strikingly between the low‐ and high‐risk groups. It was revealed that patients with higher risk scores appeared shorter survival times and worse survival status (Figure S1A,C,E). Additionally, the expression of AL365181.2↗, EMSLR, and AC092718.4 was increased, whereas the level of ENTPD3‐AS1 was decreased in the high‐risk group, which might be a safeguarding element for LUAD patients (Figure 1K, Figure S1B,D).

Finally, datasets from GEO were recruited to explore the prognostic predicting value of the five signature DRLs in our prognostic model. However, because of the re‐annotation of lncRNAs in different microarrays and platforms, we did not observe data about lncRNA AL606489.1↗. Thus, we further analyzed the remaining four lncRNAs and performed the survival analysis of single lncRNA in GSE29013↗, GSE30219↗, GSE31210↗, and GSE37745↗ for AC092718.4, AL365181.2↗, ENTPD3‐AS1, and EMSLR (Figure 1L–O, Figure S1F–I). Interestingly, high expression of AC092718.4, AL365181.2↗, or EMSLR was an independent prognostic factor and related to worse OS, PFS, and RFS (recurrence‐free survival) in lung cancer patients. And lncRNA AC092718.4 was the only shared prognostic factor in all the above four GSE datasets. Meanwhile, consistently with the prognostic model, ENTPD3‐AS1 was a safeguarding element and high expression of ENTPD3‐AS1 was related to better OS in the GSE30219↗ dataset (Figure 1M).

In the whole TCGA‐LUAD cohort, the forest diagrams of univariable (Figure 2A) and multivariable (Figure 2B) Cox regression analyses indicated that the risk score could predict the survival independently as if the clinical stage (p < 0.001). ROC curves of clinical information suggested that the risk score just like the clinical stage had relatively good accuracy in predicting patient outcomes in terms of prognosis and survival (Figure 2C). C‐index curves also demonstrated that the risk‐scoring model had a good prediction power (Figure 2D). Together, the above outcomes revealed that the five DRL risk model possessed good stability in forecasting patient prognosis.

Using risk score, gender, age, and clinical stage, we constructed a nomogram excepting to predict the OS of LUAD patients at 1, 3, and 5 year. A randomly picked patient sample got 128 points in total, whose OS incidence was 90.1%, 66.1%, and 42.4% at 1, 3, and 5 years after the diagnosis of LUAD, respectively (Figure 2E). Calibration plots of the nomogram showed the excellent agreement between observed and OS predicted by the nomogram at 1, 3, and 5 years (Figure 2C,F‐index = 0.693).

To judge the classification capability of the five DRLs risk model, the PCA was applied in all genes, DRGs, DRLs, and prognostic DRLs (Figure 2G–J). The results demonstrated that the prognostic lncRNAs owned a more effective distributing ability for low‐ and high‐groups (Figure 2J), indicating that the 5‐DRL signatures were an essential prognostic element for LUAD patients. Finally, the risk model also exhibited great prognostic value in both early (p = 0.013) and advanced‐stage (p = 0.002) LUAD patients (Figure 2K,L). Consistent with the previous results in all patients, better OS appeared in patients with lower risk scores either in the earlier stage or relatively advanced stage, which demonstrated once again that the risk model based on the five prognostic DRLs has a nice value in the prognostic prediction of LUAD patients.

To examine the tumor microenvironment (TME) score disparities between the low‐ and high‐risk groups, we employed the assistive services of ESTIMATE. The consequential findings disseminated that samples belonging to the low‐risk group garnered a considerably augmented stromal score, immune score, and ESTIMATE score (Figure 3A), which often represented better therapeutic outcomes and survival. Then, we adopted the CIBERSORT algorithm for a comprehensive insight into the distribution model of 22 tumor‐infiltrating immune cells between the low‐ and high‐risk groups (Figure 3B). The subjects with high‐risk score presented with lower infiltration abundance of plasma cells, resting memory CD4⁺ T cells, monocytes, resting DCs and higher abundance of activated memory CD4⁺ T cells, regulatory T cells (Tregs), resting NK cells, M0 and M1 macrophages (Figure 3C). Moreover, A positive correlation of risk score was found with the score of APC co‐inhibition, MHC class I, para‐inflammation, and type I IFN response, while a negative correlation existed with the score of B cells, HLA, mast cells, neutrophils, T helper cells, and type II IFN response (Figure 3D).

As there exists a direct correlation between tumor mutation burden (TMB) and both the clinical efficacy of immunotherapy and the survival prognosis of patients afflicted by cancer [11], it was analyzed in high‐ and low‐risk groups. About 94.61% of samples mutated in types of a missense mutation, nonsense mutation, frame‐shift deletion, frame‐shift insertion, or in‐frame deletion in the high‐risk group (Figure 4A), and the frequency of mutations in the low‐risk group was 85.88% (Figure 4B). A significant difference was noted in the TMB between the high‐risk and low‐risk groups (Figure 4C, p = 0.008). Although a nonsignificant difference showed in TIDE (Tumor immune dysfunction and exclusion) score between those two groups (Figure 4D), better overall survival was shown in patients in the high mutant burden group (p = 0.024, Figure 4E). Inspiringly, survival results based on the integration of TMB and risk score indicated that significant differences between all the four groups (p < 0.001), and patients with high‐TMB and the low‐risk score had the best survival probability (p < 0.001, Figure 4F).

Finally, the possibility of assessing the pretreatment drug sensitivity was conducted based on the 5‐lncRNA prognostic model in LUAD. Patients with higher risk scores seemed more resistant to chemotherapy, such as gemcitabine, paclitaxel, docetaxel, cisplatin, and irinotecan (Figure 4G–K), which constitute the most commonly used chemotherapy protocols in LUAD patients with various stages. Compared to the small cell or squamous cell lung cancer patients, more LUAD patients could benefit from the target therapies including but not limited to EGFR TKIs, VEGFR inhibitors, or c‐MET inhibitors. Intriguingly, consistent with results from chemo‐drug analyses, LUAD patients in the high‐risk group were also resistant to c‐MET inhibitor Savolitinib, MEK, ERK, CHK, and ATR inhibitors (Figure 4L–P).

To gain a greater understanding of the possible biological processes and signaling pathways associated with the prognostic risk score signatures, we employed gene ontology (GO) enrichment as well as KEGG pathway analyses that were centered around the risk group based on the five prognostic DRLs. The outcomes of these analyses revealed that there was a notable enrichment of DEGs in gland development, collagen‐containing extracellular matrix, and receptor‐ligand activity (Figure 5A,B). The top three enriched signal pathways were cell cycle, complement and coagulation, and neuroactive ligand–receptor interaction (Figure 5C,D).

Additionally, we implemented GSEA in the high‐ and low‐risk groups, respectively. In GSEA‐KEGG analyses, the DEGs of the high‐risk group demonstrated significant enrichment in pathways about cell cycle, DNA replication, oocyte meiosis, spliceosome, and steroid hormone biosynthesis pathways (Figure 5E), whereas the DEGs in low‐risk group were mainly enriched pathways about alpha‐linolenic acid metabolism, asthma, the intestinal immune network for IgA production and systemic lupus erythematosus (Figure 5F). In GSEA‐GO analyses, the DEGs that characterized the high‐risk group exhibited notable concentration in processes relating the mitotic nuclear division, regulation of mitotic nuclear division, regulation of nuclear division, sister chromatid segregation and intermediate filament functions (Figure 5G), whereas the DEGs in low‐risk group were mainly concentrated in B‐cell receptor signaling pathway, cilium movement, ciliary plasm, immunoglobulin complex, and immunoglobulin complex circulating functions (Figure 5H).

The relative expressions of the five signature DRLs were examined by qRT‐PCR. Compared to benign BEAS‐2B cells, most signature DRLs were overexpressed in lung cancer cell lines (Figure 6A, black asterisks). The expressions of AC092718.4 and AL365181.2↗ were decreased in H1975 and H1299 cells compared to H460 cells known to be highly aggressive (Figure 6A, red asterisks) [12]. Combined with the aforementioned results that high expression of AC092718.4 was a poor prognosis factor in all four GSE datasets (Figure 1L–O), lncRNA AC092718.4 was chosen to further understand the involved functions and underlying mechanisms in LUAD cell lines. After transfection with pcDNA3.1‐AC092718.4, the expressions of AC092718.4 were significantly increased in H1975 (Figure 6B) and H1299 (Figure 6C). And results of transwell assays indicated that overexpression of AC092718.4 enhanced the invasion and migration abilities of H1975 (Figure 6D) and H1299 cells (Figure 6E). Moreover, we applied colony‐forming and CCK8 assays to explore the role of AC092718.4 in tumor proliferation. As the outcomes show, the upregulation of AC092718.4 contributed to the colony formation of LUAD cells (Figure 6F,G). Consistently, we observed that the absorbance of OD 450 nm was enhanced by the overexpression of AC092718.4, indicating an improved tumor proliferative capacity in pcDNA3.1‐AC092718.4‐transfected cells (Figure 6H). Experiments including but not limited to AC092718.4 downregulation, bulk, and single‐cell RNA‐sequencing will be completed in our future study.