The expression and prognostic value of disulfidptosis progress in lung adenocarcinoma

The expression profile of 23 DRGs

As shown in Figure 1A, among 23 DRGs, 3 DRGs (GYS1, NDUFA11, ACTN4) were located at Chromosome 20. In addition, 3 DRGs (LRPPRC, NCKAP1, NDUFS1) were situated at Chromosome 2. 2 DRGs (FLNB, RPN1) were found at Chromosome 3. 2 DRGs (NUBPL, INF2) were located at Chromosome 14. Others (CAPZB, SLC7A11, CD2AP, ACTB, TLN1, PDLIM1, SLC3A2, MYL6, IQGAP1, MYH10, DSTN, MYH9, FLNA) were located at Chromosome 4, 6, 7, 9, 10, 11, 12, 15, 17, 21, 22, and X, respectively. The Wilcoxon rank-sum test showed significantly different expressions in 22 DRGs, except for ACTN4, between normal and tumor tissues (Figure 1B). The correlation among 23 DRGs was presented in Figure 1C. Among 23 DRGs, 15 DRGs had gene mutation, although the mutation frequency for each DRGs was relatively low (Figure 1D).

Figure 1. The expression of 23 DRGs in LUAD. (A) Location of 23 DRGs in chromosomes. (B) The different expression of DRGs between normal and LUAD tissues. (C) The correlation network of 23 DRGs. (D) The mutation frequency of 23 DRGs. Abbreviations: DRGs: Disulfidptosis-related genes; LUAD: Lung adenocarcinoma; N: normal tissues. ^*P < 0.05, ^**P < 0.01, ^***P < 0.001.

Establishment and evaluation of a DRGs model

We utilized the least absolute shrinkage and selection operator (LASSO) regression (LASSO) analysis to analyze 23 DRGs and determined 8 DRGs highly related to the overall survival (OS) of LUAD patients in the TCGA set (Figure 2A). Then, those 8 DRGs were subjected to the multivariate Cox regression analysis. 5 candidate DRGs (ACTB, FLNB, NCKAP1, SLC3A2, SLC7A11) were screened out and used to build a prognostic model (Figure 2B) and a risk score formula to calculate the risk score of all patients. The risk score formula was presented as follows. The Kaplan–Meier plots demonstrated that all 5 DRGs were unfavorable factors for the clinical outcomes of LUAD patients (Figure 2C–2G). Figure 2H showed the correlation among 5 DRGs. According to the optimal cut-off risk score: 0.987, patients were divided into low- and high-risk groups. The principal component analysis (PCA) revealed that the samples in two risk groups were distributed in different areas (Figure 2I).

$$\begin{array}{l}\text{Risk score}=(0.430\times \exp \text{of ACTB})\\ +(0.248\times \exp \text{of FLNB})\\ +(0.408\times \exp \text{of}\text{NCKAP1})\\ +(0.336\times \exp \text{of}\text{SLC3A11})\\ +(0.153\times \exp \text{of}\text{SLC7A11})\end{array}$$

Figure 2. Construction of a DRGs model in TCGA set. (A) The LASSO regression analysis to filter out prognosis-related DRGs. (B) The multivariate Cox regression analysis to develop a DRGs model. (C–G) Survival analysis of 5 candidate DRGs (ACTB, FLNB, NCKAP1, SLC3A2, SLC7A11). (H) The correlation network of 5 candidate DRGs. (I) Patients in different risk groups gathering in two areas in PCA analysis. Abbreviations: LASSO: Least absolute shrinkage and selection operator; PCA: Principal component analysis.

Next, we calculated the area under the curve (AUC) value and the c-index in each set to estimate the predictive accuracy of the DRGs model. The results showed the AUC values and c-index based on risk score in all sets were higher than that based on clinical indexes (age, gender, smoking, TNM stage), hinting that the DRGs model performed well in predicting clinical outcomes of LUAD patients (Figure 3A, 3B). In addition, the calibration plot in each set revealed that when the line angle is 45°, it represents the best prediction result, indicating that the DRGs model had a good predictive ability (Figure 3C).

Figure 3. Evaluation of the DRGs model. (A) The time-dependent AUC value in each set. (B) The time-dependent C-index in each set. (C) The calibration plots of each set. Abbreviations: AUC: The area under the curve; C-index: concordance index.

Correlation between DRGs model and clinical characteristics

The distributions of samples with different TNM stage (P = 0.003), tumor size (P < 0.003), and lymph node metastasis (P = 0.025) were significantly diverse between different risk groups (Figure 4A). Moreover, there was a significant difference in the risk scores among patients with varying tumor sizes (P < 0.001, Figure 4E). Similar results were observed among patients with different TNM stages (Figure 4H). The risk score significantly increased in female patients (P = 0.042, Figure 4C) or with lymph node metastasis (P = 0.009, Figure 4F). However, we did not find the effect of age (Figure 4B), smoking (Figure 4D) and distant metastasis (Figure 4G) on the risk score.

Figure 4. Correlation between risk score and clinical characteristics. (A) The distribution of clinical characteristics between high- and low-risk groups. (B) The difference of risk score between Age ≤65 and Age >65, (C) different gender, (D) smoking and non-smoking, (E) tumor size, (F) with and without lymph node metastasis, (G) with and without distance metastasis, (H) different TNM stage. ^**P < 0.01.

The DRGs model is an independent prognostic factor for the prognosis of LUAD

The Kaplan–Meier plot in the TCGA set illustrated that the high-risk score predicted poorer clinical outcomes (P < 0.001, Figure 5A). It was confirmed in GEO set 1 (P < 0.001, Figure 5B), GEO set 2 (P < 0.001, Figure 5C), and GEO set 3 (P < 0.001, Figure 5D). Furthermore, stratification analyses demonstrated the low-risk patients had a better prognosis in each subgroup (Supplementary Figure 1).

Figure 5. The DRGs model is an independent prognostic factor for the prognosis of LUAD. Survival difference between high- and low-risk groups (A) in TCGA set, (B) GEO set 1, (C) GEO set 2, (D) GEO set 3. (E) The univariate Cox regression analyses of risk score and clinical characteristics in TCGA set. (F) The multivariate Cox regression analyses of risk score and clinical characteristics in TCGA set.

To compare the DRGs model with clinical parameters (age, gender, stage, tumor size, lymph node metastasis, and distance metastasis), we performed the univariate and multivariate Cox regression model. The univariate Cox regression model showed that the DRGs model was an essential prognosis-related influence factor (HR = 3.423, 95% CI (2.431, 4.820), P < 0.001, Figure 5E). The multivariate Cox regression model demonstrated that the DRGs model was an independent prognostic factor for the OS of LUAD patients (HR = 3.316, 95% CI (2.300, 4.782), P < 0.001, Figure 5F). Similar results were observed in GEO set 1, 2, and 3 (Supplementary Figure 2).

Gene set enrichment analysis (GSEA)

Then, we performed the GSEA to explore the altered biological roles and signaling pathways between the low- and high-risk group. A total of 27 pathways were enriched, including Pathways in cancer, P53 signaling pathway, and so on, that were highly associated with tumorigenesis and development (Figure 6A). Figure 6B, 6C presented the top 5 pathways enriched in high- and low-risk groups.

Figure 6. Gene set enrichment analysis (GSEA) between high- and low-risk groups. (A) 27 pathways enriched in GSEA. (B) Top 5 pathways enriched in high-risk group. (C) Top 5 pathways enriched in low-risk group. Abbreviations: FDR: False discovery rate; NES: Normalized enrichment score.

Tumor mutational burden (TMB)

The mutational landscape showed that the frequent mutation events in the low-risk group (96.53%) were more than that in the high-risk group (90.99%) (Figure 7A, 7B). In the high-risk group, the most frequently mutated gene was TTN (50%), followed by TP53 (50%) and MUC56 (44%). In the low-risk group, the top 3 frequently mutated genes were TP53 (45%), TTN (43%), and CSMD3 (38%). However, there was no difference in TMB between the two risk groups (P = 0.180, Figure 7C). Survival analysis revealed that high TMB had a better prognosis in LUAD patients (P = 0.008, Figure 7D). Furthermore, it could better predict prognosis when combining the TMB and the risk score (Figure 7E).

Figure 7. The tumor mutation burden characteristics in low- and high-risk group. (A) Mutational landscape in the low-risk group. (B) Mutational landscape in the high-risk group. (C) The difference in TMB between two groups. (D) Survival difference between high- and low-TMB. (E) Survival analysis of TMB along with risk score. Abbreviation: TMB: Tumor mutational burden.

Investigation of the relation of the DRGs model to the immune statue

We applied the single-sample gene set enrichment analysis (ssGSEA) to analyze the level of 15 tumor-infiltrating immune cells (TIICs) and the activity of 13 immune-related functions. The results demonstrated that the enrichment scores of 9 TIICs (aDCs, B cells, CD8 T cells, iDCs, Mast cells, Neutrophils, Tfh, Th1 cells, TIL) in the low-risk group were significantly enhanced (Figure 8A). Similar results were observed in 4 immune-related functions (Check-point, HLA, T cell co−stimulation, Type II IFN response, Figure 8B). The immune score in the high-risk group was significantly higher than in the low-risk group (Figure 8C). The TIDE score and Exclusion score in the high-risk group were increased (Figure 8D), whereas the IPS values of PD-1/PD-L1 and CTLA4 were decreased (Figure 8E). Drug sensitivity analysis manifested a significant difference in the sensitivity of 47 anti-tumor drugs between high- and low-risk groups (Supplementary Figure 3).

Figure 8. The relation of the DRGs model to immune statue. Comparison of the enrichment score of (A) TIICs, (B) immune-related functions, (C) TME, (D) TIDE score, (E) The immunophenoscore value between low- and high-risk groups. ^*P < 0.05, ^**P < 0.01, ^***P < 0.001.

scRNA-seq analysis

Before filtering, there were 45,632 features for 23,747 cells in the 8 LUAD samples. Then, we performed data standardization and quality control and finally selected 18,977 cells and the top 2000 highly expressed and variable genes for further analysis. The PCA for reducing data dimensionality reduction revealed no noticeable separation trend of cells. Nonlinear dimension reduction was performed with the t-Stochastic Neighbor Embedding (t-SNE) algorithm, which successfully clustered the cells into 14 clusters (Figure 9A). Then, we annotated all clusters, and 11 cell types were identified (Figure 9B). In addition, the expression level of SLC7A11 and SLC3A2 were most abundant in Epithelial cells (Figure 9D). NCKAP1 and FLNB had the highest expression in Endothelial cells. ACTB was expressed in most cell types. However, the expression level in Naive CD4 T was the richest. In this study, cluster 2, 8, and 13 were annotated into Epithelial cells, which mainly contained cancer cells and cancer stem cells. Thence, we included the Epithelial cells in the pseudo-time cell differentiation trajectory analysis. The results were presented in Figure 9F and demonstrated the evolutionary pattern of Epithelial cells. Based on the ordering of pseudotime, Epithelial cells (cluster 2) were divided into Epithelial cells (cluster 8) and Epithelial cells (cluster 13). Figure 9G showed the 5 candidate DRGs expression in different developmental states. Finally, we further analyzed the Epithelial cells with the t-SNE algorithm and clustered Epithelial cells into 6 sub-clusters (Figure 9C). The expression of SLC7A11, SLC3A2, NCKAP1, and ACTB in sub-cluster 4 was highest, whereas FLNB in sub-cluster 2 was highest (Figure 9E).

Figure 9. scRNA-seq data analysis. (A) The t-SNE algorithm divided the cells into 14 clusters by 20 principal components. (B) The tSNE plot revealing 14 clusters was annotated into 11 different cell types. (C) The t-SNE algorithm divided Epithelial cells into 6 sub-clusters. (D) The expression of 5 candidate DRGs in 11 cell types, (E) in 6 Epithelial cells sub-clusters. (F) The trajectory analysis of Epithelial cells with top 10 marker genes. (G) 5 candidate DRGs expression in different developmental states. Abbreviations: ScRNA-seq: Single-cell RNA sequencing; t-SNE: t-Stochastic Neighbor Embedding.

Construction of a prognostic nomogram

Using clinical characteristics (age, gender, smoking, TNM stage) and risk score, we develop a prognostic nomogram in TCGA set for predicting 1-, 3-, and 5-year OS of LUAD patients (Figure 10A). The AUCs and c-index in each set were more than 0.7 (Figure 10B, 10C). In addition, the calibration plot showed an optimally fit with the ideal model (Figure 10D).

Figure 10. Construction of a prognostic nomogram. (A) A nomogram for predicting the 1-, 3-, and 5-year overall survival rate in LUAD patients. (B) The time-dependent AUC value. (C) The time-dependent C-index. (D) The calibration plot.

Patients set

Four independent sets, including 1306 samples, were collected from TCGA and GEO databases in the study. TCGA set (n = 465) was from TCGA-LUAD and level 3 RNA sequencing. GEO set 1 (n = 432), GEO set 2 (n = 83), and GEO set 3 (n = 226) were from GSE68465 [14], GSE30219 [15], and GSE31210 [16], respectively. The gene expression profile in GSE68465, GSE30219, GSE31210 was raw data and normalized with the robust multi-array average algorithm via the R package Affy (3.17). The corresponding clinical data of all patients were also collected. Patients lacking pathologic diagnosis, survival time, or survival status would be removed. All patients’ detailed demographic and baseline information was presented in Table 1. ScRNA-seq data from 8 human LUAD samples were gathered from the GSE171145 [17] dataset in the GEO database. Disulfidptosis-related genes were obtained from previous literature [8, 9, 18, 19].

Establishment and evaluation of a disulfidptosis-related genes model (DRGs) model

In the TCGA set, we investigated the difference in the expression of 23 DRGs between normal and LUAD tissues with the Wilcoxon rank-sum test. Then, we preliminary identified prognosis-related genes among 23 DRGs using the least absolute shrinkage and selection operator (LASSO) regression analysis with R package “glmnet” (4.0.2). Under the minimum lambda, genes with coefficient ≠ 0 were chosen for screening genes. The iteration was set as 1,000 and 10-fold cross-validations to avoid overfitting. Then, the screened genes were subjected to the multivariate Cox regression analysis to filter candidate genes and develop a DRGs model for predicting the prognosis of LUAD. The risk score was calculated as Risk score = h0(t) × exp (β1X1 + β2X2 + … + βnXn), where β refers to the regression coefficient; X represented the gene expression level; h0(t) is the benchmark risk function. Also, a risk score formula was formed, and patients’ risk scores were calculated. Finally, we used the “surv cutpoint” functions of the package “Survminer” (0.4.9) in R to determine the optimal cut-off value of risk score, which divided patients into low- and high-risk groups. Based on the cut-off value, patients in GEO set 1, 2, and 3 were also divided into two groups.

The AUC value and the c-index in four sets were calculated to assess the predictive capacity of the DRGs model. The calibration plot with a boot-strapping set of 1,000 resamples was also pictured to investigate the predictive accuracy of the DRGs model.

GSEA

The GSEA was performed to explore the altered biological roles and signaling pathways between the low- and high-risk group with the R package “clusterProfiler” (4.2.2). The “c2.cp.v7.2.symbols.gmt (Curated)” was selected as annotated gene set. The false discovery rate (FDR) <0.05 was set as the threshold.

Assessment of the relation of the DRGs model to the immune statute and drug sensitivity

The R package: GSVA (3.1.7) was applied to carry out the ssGSEA reflecting the level of TIICs and the activity of immune-related functions. The R package: “ESTIMATE” (1.0.13) was used to calculate the immune and stromal score in the tumor microenvironment. The immunophenoscore (IPS) value of samples was gathered from The Cancer Immunome Atlas (TCIA, https://tcia.at) to present the response to immune checkpoint inhibitors treatment [20]. The TIDE score was collected from Tumor Immune Dysfunction and Exclusion (TIDE, http://tide.dfci.harvard.edu). The higher TIDE score represents the worse response to immune therapy.

We took advantage of the R package “pRRophetic” (version 6) [21] to calculate the semi-inhibitory concentration values of commonly used anti-tumor drugs for LUAD to look for variations in the efficacy of medicines between the two groups. P < 0.01 was set as the cut-off value.

TMB

The DNA somatic mutation data of corresponding LUAD patients were also downloaded from TCGA and further analyzed with the “maftools” R package (3.17). The ChromPlot function in the “maftools” R package was used to visualize the output results. The results were presented with waterfall diagram to show the variation distribution of genes with high somatic mutation frequency in LUAD samples and two different groups. Further, it was used to evaluate the difference in survival of the patients classified in the high-TMB and low-TMB groups.

scRNA-seq data processing

The R package: “Seurat” (3.4) [22] was used to analyze the transcript count matrix for quality control and preliminary data exploration. The filtering threshold was set as follows:

• Excluding genes detected in less than 3 cells

• Excluding cells with <50 genes detected

• Excluding cells with >10% mitochondrial gene expression

Then, the expression profiles were normalized with the Log Normalization algorithm and were subsequently normalized using a linear regression model. The top 2000 highly expressed and variable genes were selected for PCA to determine significant and influential dimensions. The t-SNE algorithm was used to reduce the dimension of the top 20 principal components and gather major cell clusters. The marker genes between difference clusters were identified with |log2 (fold change) |>1 and adjusted P value < 0.05 as the threshold. Cell annotation was carried out with the “SingleR” package (3.17) [23] and reports from the literature [24–26]. Finally, single-cell trajectory analysis was performed with the “Monocle 2 algorithm” [27].

Statistical analysis

The categorical data were presented as Numbers and compared with the chi-square test. The measurement data were presented as Mean ± standard deviation (SD) and compared with the Wilcoxon rank-sum test. Correlation analysis was performed with the Spearman correlation test. Survival analysis was performed with the Kaplan–Meier plot and compared by the log-rank method. Finally, the univariate and multivariate Cox regression analysis determined the independent prognostic predictors. P < 0.05 was set as the cut-off value. All statistical analyses were conducted with R 4.1.1 (https://www.r-project.org).

Availability of data and materials

Bulk RNA-seq data were analyzed in this study, this data can be found at: The Cancer Genome Atlas (TCGA-LUAD, https://portal.gdc.cancer.gov) and Gene Expression Omnibus (GSE68465, GSE31210, and GSE 30219; https://www.ncbi.nlm.nih.gov). ScRNA-seq data were collected from GEO-GSE171145.