Identification of disulfidptosis in esophageal squamous cell carcinoma based on single-cell and bulk RNA-seq data to predict prognosis and treatment response

We obtained the GSE53625↗ dataset (gene expression profiles and clinical data) from the Gene Expression Omnibus (GEO). After filtering, 179 tumor samples with complete expression and survival (overall survival [OS]) data were included as the training cohort. The R package TCGAbiolinks was used to download FPKM expression profiles (for log transformation), survival data, and clinical data from the TCGA-ESCA EC and kept 80 squamous cancer samples with both expression and survival data for validation. The expression profiles (UMI-count) of 60 squamous cancer samples from the GSE160269↗ single-cell dataset were selected for analysis in this project. A dataset of cutaneous melanomas treated with anti-CTLA-4 was downloaded and used to assess the predictive efficacy of signature in the immunotherapy cohort (14). Interaction relationships of disulfidptosis-related genes were obtained based on the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (https://www.string-db.org/↗), and Protein-Protein Interaction Networks (PPI) were constructed based on the interactions (10).

Based on the expression profile data of 10 disulfidptosis-related genes, unsupervised cluster analysis was applied to identify different expression patterns of disulfidptosis. The ConsensusClusterPlus package was used for the operation, the distance used for clustering was Pearson, the clustering method was Pam, and 1000 replications were performed to ensure the stability of the classification. The R packages survminer and survival were used to generate survival curves for prognostic analysis using the Kaplan-Meier method. The log-rank test was used to determine the significance of differences and resolve the correlation between samples with different expression patterns and OS. Additionally, the R package limma was used to identify differentially expressed genes between different subgroups. We identified differentially expressed genes with |log2FC| ≥ 0.585 (fold change ≥1.5) and FDR < 0.05.

The hazard ratio (HR) and prognostic significance of differential genes and screened genes with p < 0.05 as prognostic correlates were determined using univariate Cox regression analysis.

We applied LASSO regression (glmnet R package) to select key prognostic genes. The risk score was calculated by summing the product of each gene’s expression level (exp) and its LASSO coefficient (coef): Score = Σexp*coef. The samples were divided into high and low groups according to the score. We generated Kaplan-Meier survival curves to compare patient outcomes. The significance of the differences was determined using the log-rank test to further resolve the correlation between these two types of samples and OS. The predictions of the scoring system on scoring were evaluated using the receiver operating characteristic curve (ROC), and the area under the curve (AUC) was visualized using the R package timeROC. Moreover, univariate and multivariate Cox analysis were performed to explore the independent prognostic value of the score.

We performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses using the clusterProfiler R package (parameters pvalueCutoff = 0.05, pAdjustMethod = “BH”), and the R package GSVA was used to perform genes differentially expressed in the two groups of samples were screened using the difference multiplicity |log2FC| ≥ 0.585 (difference multiplicity great ne set variation analysis to annotate the potential functions of key genes and predict their potential molecular functions. Gene sets were downloaded from the Molecular Signatures Database (MSigDB) for the HALLMARK, KEGG, and Gene Ontology Biological Processes (GOBP) sublibraries and used to conduct GSVA analysis.

The immune score, stromal score, and tumor purity were calculated for each tumor sample using the ESTIMATE algorithm. The Spearman correlation was calculated between score and immune score, stromal score, and tumor purity. According to previous studies, the single-sample gene-set enrichment analysis (ssGSEA) algorithm was used to estimate the relative abundance of each cell infiltrate in TME. The distribution of immune cell infiltrates were compared between different grouped samples using the Wilcoxon test.

Based on the Genomics of Drug Sensitivity in Cancer (GDSC) (https://www.cancerrxgene.org/↗) and Chemical and Systems Biology Program (CTRP) (https://portals.broadinstitute.org/ctrp/↗) cancer genomics drug susceptibility databases, R package oncoPredict’s calcPhenotype algorithm was used to evaluate the drug IC₅₀ values for each sample in the training set. The spearman correlation was calculated between score and drug IC₅₀ to assess the correlation between drug sensitivity and signature. The difference between drugs IC₅₀ in high and low groups was compared.

After the original authors’ quality control, 60 ESCC samples were analyzed from the GSE160269↗ dataset using the R package Seurat (v4.1.0) and normalized using the NormalizeData function. The FindVariableFeatures function was used to identify highly variable genes based on their average expression values (greater than 0.1 and less than 8) and dispersion (greater than 1) to identify highly variable genes. Batch correction between samples was performed using the R package Harmony to avoid batch effects interfering with downstream analysis. Then, the data were scale-transformed, principal component analysis was used for dimensionality reduction, the top 50 principal components were selected for downstream analysis, and the RunUMAP function was used for visualization.

The cell type annotation information provided by the original authors of the dataset was extracted. Among them, the malignant epithelial cells were identified using the copykat method. The CellScore was calculated using the AddModuleScore function of the Seurat package based on the signature-containing genes. Malignant epithelial cells were classified into high and low groups based on the median CellScore of malignant epithelial cells. Cell stemness scores were calculated using the AddModuleScore function and stemness gene sets from previously published studies.

The R package monocle2 was used for trajectory analysis, time series analysis of malignant epithelial cells was proposed, and the R package CellChat was used for cell-to-cell communication analysis.

All analyses were performed using R version 4.1.2. For significance analysis between various values (expression, infiltration ratio, and various eigenvalues), the Wilcoxon rank sum test was used to compare differences between two groups of samples, whereas the Kruskal-Wallis was used to compare differences between multiple groups of samples. For plot presentation, ns indicates p > 0.05, * indicates p < 0.05, ** indicates p < 0.01, *** indicates p < 0.001, and **** indicates p < 0.0001. Survival curves were generated for the prognostic analysis by the Kaplan-Meier method, and determined the significance of differences using the log-rank test.

We compared their expression differences between normal and tumor samples to assess whether the expression of disulfidptosis-related genes affects tumor progression in ESCC. The results revealed that six genes were significantly highly expressed in tumor samples, whereas two genes were significantly low (Figure 1A). NDUFS1 showed the highest connectivity in the protein interaction network (Figure 1B), suggesting its potential as a hub gene regulating disulfidptosis.

We performed unsupervised cluster analysis using the R package ConsensusClusterPlus based on 10 disulfidptosis-related genes expression profiles of ESCC samples from the GSE53625↗ dataset. We identified two subgroups, named Cluster1 and Cluster2 (N = 85/94, Figures 1C–E). The Kaplan-Meier curves revealed that the two subgroups had significantly different prognoses, with Cluster1 having a worse OS (Figure 1F). We also compared the clinicopathological characteristics of different subgroups of ESCC using a heat plot. We discovered significant differences in the distribution of esophageal carcinoma sites (locations) between subgroups of patients (p < 0.05, Figure 1G).

We constructed a disulfidptosis signature for predicting the prognosis of ESCC patients based on differential genes between the disulfidptosis expression patterns. Initially, we screened 493 differentially expressed genes between the two disulfidptosis expression patterns using the limma package (Figures 2A, B). GO and KEGG enrichment analysis of these differentially expressed genes using the clusterProfiler package. The KEGG results indicated that these genes were primarily enriched in pathways related to cytochrome P450, DNA adducts, reactive oxygen species (ROS), suggesting potential associations with drug toxicity, drug resistance, or metabolism of environmental carcinogens, involving mechanisms of carcinogenesis driven by DNA damage or oxidative stress. For the GO analysis: Biological Process (BP) terms were predominantly linked to xenobiotic substance metabolism and stress response, lipid metabolism and inflammatory regulation. Molecular Function (MF) terms were mainly associated with oxidation-reduction and detoxification functions, monooxygenase and lipid metabolism, and heme binding and oxygen metabolism. Cellular Component (CC) terms were enriched in processes related to cornified envelope formation, neural synapses, and regulation of cell polarity. (Figure 2C).

Next, the univariate Cox regression analysis exposed that 47 of these differential genes were significantly associated with patient OS in the GSE53625↗ esophageal squamous cancer cohort (Supplementary Table 1), and the forest plot visualizes the Cox regression analysis results of the 20 genes with the smallest p-values (Figure 3A). Supplementary Figure 1 displays the KM curves of the six genes with the smallest to the greatest P-values.

Furthermore, LASSO-Cox regression analysis was performed based on these 47 genes. We performed 10-fold cross-validation under optimal conditions to determine the model’s penalty parameter (λ), eliminating 19 key prognostic factors affecting patient OS (Figures 3B–D). Based on the key prognostic factor expression levels and the linear combination of the corresponding weights, we construct the signature that can assess each patient’s prognosis. Table 1 illustrates the coefficients of each factor.

We calculate the risk score for each patient in the training set based on the constructed prognostic signature and divide them into high- and low-risk groups based on the median value. KM curve analysis and the log-rank test indicated that high-risk group patients had significantly shorter OS (p-value < 0.05, Figures 3E). The AUC for the predicted outcome of the sample was 0.712, 0.726, and 0.756 at one, two, and three years respectively (Figures 3F), indicating that the Score can provide a good characterization of sample OS. Then, we explored the independence of prognostic signature in the training set. Univariate and multivariate Cox regression models were constructed based on prognostic signature and clinical characteristics, and the results revealed that prognostic signature was an independent prognostic factor (HR = 2.31, p-value < 0.05, Figure 3I).

We used TCGA-ESCC as an independent validation set to assess the reliability of the prognostic signature. Patients were divided into high- and low-risk groups according to the prognostic signature risk score. Patient OS was also significantly lower in the high-risk group than in the low-risk group (Figure 3G), with predicted outcome AUCs of 0.649, 0.803, and 0.625 for the validation set samples at one, two, and three years, respectively (Figure 3H). Univariate and multivariate Cox regression models were constructed based on prognostic signature and clinical characteristics, and the results were consistent with the test set, again supporting prognostic signature as an independent prognostic factor (HR = 3.22, p-value < 0.05, Figure 3J).

The original authors of the dataset quality-controlled scRNA-seq data of 192,078 cells (containing 20,335 B cells, 10,346 Endothelial, 44,547 Epithelial, 27,881 Fibroblast, 1,136 Fibroblastic reticular cells (FRC), 18,514 Myeloid, 3,023 Pericytes, and 66,296 T cell). Moreover, 17,986 genes were detected. The PCA results presented a significant batch effect between samples (Figure 4A), and the batch effect between samples was removed using Harmony (Figure 4B). UMAP demonstrated the distribution of different cell types (Figure 4C), and the proportion of cell distribution in each sample was heterogeneous (Figure 4D).

We extracted epithelial cells identified by copykat as malignant and with at least one model gene detected (N = 19159, such as at least one model gene with expression greater than 0), calculated Stemness and CellScore for each cell, and classified CellGroup into high and low groups according to CellScore, using UMAP for presentation (Figures 4F, J, M). Trajectory analysis of malignant epithelial cells revealed three differentiation states (Figure 4E). State1 cells, characterized by low stemness (Figures 4G, H), represented the least aggressive phenotype. Tumor progression was associated with a gradual shift from State1 to State2/3, marked by increased stemness and loss of differentiation (Figure 4I). State1 → State2 and State1 → State3 trajectory routes have an increased CellScore (Figures 4K, L) and an increased proportion of high in the CellGroup (Figure 4N), indicating an increase in tumor malignancy.

Based on the bulk sequencing data, we calculated the pathway/biological process activities of KEGG and GOBP using the GSVA algorithm. We compared the differences between the activities of different groups using the rank sum test. The results demonstrated that the immune-related biological processes activities, such as T cell activation and T cell-mediated tumor cell immune response, were significantly higher in the high-score group than in the low-score group. High-score tumors exhibited stronger activation of pro-tumor pathways, including p53 signaling, chemokine production, and immune cell recruitment (Figure 5A).

Based on scRNA-seq data, we calculated differentially characterized genes between high and low groups of CellGroup and performed GO and KEGG enrichment analysis using clusterProfiler. The results presented that genes were significantly enriched in immune-related biological processes, such as regulation of T-cell activation and proliferation, and tumor progression-related pathways, such as cell cycle, P53 signaling pathway, and antigen processing and presentation (Figures 5B, C), consistent with the results of bulk sequencing data.

Next, we analyzed the enrichment differences between high and low groups using GSEA. The results disclosed that the high-risk group was significantly enriched in immune-related pathways, such as B-cell activation and lymphocytes, in both bulk and scRNA-seq data (Figures 5D, E).

We calculated the percentage of immune cell infiltration in each tumor sample using the bulk sequencing data. The results revealed that the score of prognostic signature was significantly positively correlated with a stromal score, immune score, and ESTIMATE score, and significantly negatively correlated with tumor purity (Figure 6A). The percentage of infiltrated cells, such as Natural killer cells, Type1 T helper cells, and regulatory T cells (Tregs), was significantly lower in the Low group samples than in the High group (Figure 6B).

Next, we performed intercellular communication analysis using the CellChat package, with extensive cellular communication among the cell populations (Figures 7A, B). When distinguishing between incoming and outgoing signaling, Malignant Epithelial Cells, Non-Malignant Epithelial Cells, Pericytes, Low neoplastic, Endothelial, and FRC are outgoing signaling dominant senders, while B cell, T cell, Fibroblast, High neoplastic, and Myeloid are incoming signaling dominant receivers (Figure 7C). Low and high neoplastic are in multiple tumor malignant progression-related pathways such as MIF, ITGB2, and FN1 with different cells for messaging (Figures 7D–F).

We explored the predictive efficacy of prognostic signature on sample prognosis in the immunotherapy cohort (14). We discovered that patients with the same low score in the immunotherapy cohort had a better prognosis (Figure 8A). However, there was no significant difference in risk scores among patients in the immunotherapy with/without response group (Figure 8B), nor was there a significant difference in the proportion of samples responding to immunotherapy with/without immunotherapy in the different risk groups (Figure 8C).

We predicted the IC₅₀ values of the drugs in the training set sample using the R package oncoPredict, GDSC, and CTRP database drug information combined with the expression profile of the training set. We compared the spearman correlation between the score and the IC₅₀ of each drug, ranked the drugs according to the absolute value of the correlation coefficient from the largest to the smallest, and selected the top six drugs with significant positive and negative correlations, respectively (Figures 8D, F, H, J, p < 0.05), with significant differences in drug IC₅₀ across Score groupings (Figures 8E, G, I, K).