Gene expression information for SCLC and normal samples were gathered from the gene expression omnibus (GEO) database. GSE149507 was selected as the training cohort, with 18 SCLC samples and 18 normal samples. GSE108055 was selected as the test cohort, with 12 SCLC samples and 10 normal samples. FRGs (n = 355) utilized for this research were identified from the ferroptosis database (FerrDb). Supplementary Table 1 shows the gene list. Potential candidate-targeted drugs for marker genes were predicted through the drug-gene interaction database (DGIdb, https://www.dgidb.org/).

Expression data on 204 FRGs from SCLC and normal samples were obtained from GSE149507. Differentially Expressed Ferroptosis-Related Genes (DE-FRGs) between SCLC and normal samples were assessed by Student’s t-test in R (Supplementary Table 2). Genes with p 0.05 were deemed to be statistically significant, and relationships between the DE-FRGs were examined using Pearson’s correlation analysis.

The R language software package V4.2.0 (https://www.r-project.org/) and https://www.reactome.org were used to perform gene ontology (GO) and Reactome enrichment, respectively, for analysis of the functions and pathways of DE-FRGs.

The best diagnostic gene markers for SCLC were identified using support vector machine recursive feature elimination (SVM-RFE) together with the least absolute shrinkage and selection operator (LASSO). LASSO regression algorithm improves the prediction accuracy by using regularization, which was implemented using the “glmnet” R package. The SVM-RFE model was then set up by the SVM package and compared with 10-fold cross-validations by using mean misjudgment rates. The overlapping biomarkers revealed by these two algorithms were determined as the best gene biomarkers for SCLC. The diagnostic efficacy of the best marker genes was assessed from the plot of the receiver operating characteristic (ROC) curve and the computation of the area under the curve (AUC). Meanwhile, a logistic regression model with all marker genes was established using “glm” R package to distinguish SCLC samples from normal samples in GSE149507. The ROC curve was used to assess the diagnostic capability of this logistic regression model.

To study functional enrichment and pathway enrichment, we analyzed the relationship between the 6 candidate marker genes identified previously and the remaining genes in GSE149507 using R package gene sets enrichment analysis (GSEA) (V.4.1.0). All genes from high to low were then sorted according to their correlation with marker genes. The sorted genes were used as gene sets in additional studies. Supplementary Tables 3,4 summarize the GO and KEGG enrichment data for each marker gene, respectively.

Gene sets variation analysis (GSVA) analysis was performed using the R package GSVA (V.1.38.0). The difference in GSVA scores between high-expression and low-expression marker genes was examined using the “Limma” R package. The filter condition was set to $|$t$|$$>$ 2, p 0.05. In the group with high expression, the pathway was deemed active when t $>$ 0, while in the group with low expression, it was deemed active when t 0.

Immune cells in the tumor microenvironment were evaluated in the GSE149507 dataset using CIBERSORT. The sum of the proportions of 22 infiltrating immune cells determined by CIBERSORT was 1. The association between immune cells and marker genes was subsequently examined using Pearson’s correlation analysis. Specific information on the percentage of 22 types of infiltrating immune cells is shown in Supplementary Table 5.

The microRNAs (miRNAs) of targeted marker genes were predicted by the miRanda, miRDB, and Targetscan databases. miRNAs selected by all three databases were selected for the subsequent study. Long non-coding RNA (lncRNA) binding to selected miRNAs were screened using SpongeScan Version1 (http://spongescan.rc.ufl.edu/). Cytoscape V3.9.1 (https://cytoscape.org/) was subsequently implemented to map the competing endogenous RNA (ceRNA) regulatory network.

The gene expression level of all samples was compared by Student’s t-test in R V4.2.0 (https://www.r-project.org/). All statistical analyses and the production of Venn diagrams were performed using R packages. Statistical significance was considered to be p 0.05 (*p 0.05; ** p 0.01; *** p 0.001).

A total of 204 DE-FRGs were found between SCLC and normal samples in the GSE149507 dataset, comprising 85 up-regulated genes and 119 down-regulated genes (Supplementary Table 1). A clustering heatmap shows how the expression of DE-FRGs varied across all samples (Fig. 1), while a correlation heatmap shows how these DE-FRGs were related (Supplementary Fig. 1).

Fig. 1.

Clustering heatmap for the expression of 204 DE-FRGs.

GO function enrichment and Reactome enrichment analyses were performed to investigate DE-FRG-related functions and pathways. GO analysis results showed that BP (Biological Process) was enriched for the functions of: ‘response to oxidative stress’, ‘cellular response to chemical stress’, and ‘response to nutrient levels’. Enrichment of CC (Cell Component) included ‘transcription regulator complex’. Enrichment of MF (Molecular Function) contained ‘DNA-binding transcription factor binding’ and ‘RNA polymerase II-specific DNA-binding transcription factor binding’ (Fig. 2A,B). Reactome results showed that DE-FRGs were significantly enriched for ‘signaling by Interleukins’, ‘Interleukin-4 and Interleukin-13 signaling’, and ‘fatty acid metabolism’ (Fig. 2C). The ‘autophagy’ pathway and ‘cellular response to chemical stress’ pathway were also enriched. These results suggest that the effect of DE-FRGs on the progression of SCLC may be related to immune regulation, oxidative stress, transcription factors, and autophagy.

Fig. 2.

Function and pathway analysis of DE-FRGs. (A) Circle diagram of GO analysis. (B) GO enrichment analysis revealed that DE-FRGs were significantly related to the functions of “oxidative stress response”, “chemical stress cell response”, “transcriptional regulatory complex”, and “DNA-binding transcription factor binding”. (C) Reactome pathway analyses showed that DE-FRGs were significantly enriched in “interleukin signal transduction” and “fatty acid metabolism”. GO, gene ontology.

The aim of this work was to determine whether DE-FRGs could be used for the diagnosis of SCLC. LASSO regression and the SVM-RFE algorithm were performed to find the optimal DE-FRGs for distinguishing SCLC from normal tissue in the GSE149507 dataset. For LASSO regression analysis, the penalty parameters were tuned through ten-fold cross-validation, with eight SCLC characteristic genes identified (Fig. 3A,B). DE-FRGs were then filtered using the SVM-RFE algorithm, and 64 were considered to be characteristic genes (maximal accuracy = 0.942, minimal Root Mean Square Error (RMSE) = 0.0583) (Fig. 3C,D). In the ensuing analysis, six marker genes were identified as common genes by LASSO regression in conjunction with the SVM-RFE algorithm. These were the ATG3, MUC1, RRM2, IDH2, PARP1, and EZH2 genes (Fig. 3E).

Fig. 3.

Six DE-FRGs were identified as characteristic genes for SCLC diagnosis. (A,B) LASSO logistic regression identified 8 characteristic genes related to SCLC. (C,D) The SVM-RFE algorithm selected 64 characteristic genes. Six marker genes were ultimately determined to be the best genes for SCLC (maximal accuracy = 0.942, minimal RMSE = 0.0583). (E) Overlapping genes were determined using LASSO regression and the SVM-RFE algorithm. (F) ROC curves of 6 marker genes. (G) ROC curve for the logistic regression model established by the 6 marker genes. SCLC, small cell lung cancer; ROC, receiver operating characteristic; SVM, support vector machine; LASSO, least absolute shrinkage and selection operator.

Using the “glm” R package, the six marker genes were subsequently used to develop a logistic regression model. Results from ROC curve analysis showed that ATG3, MUC1, RRM2, IDH2, PARP1, and EZH2 have good diagnostic efficacy in the GSE149507 cohort, with AUC values of 0.796, 0.951, 0.997, 0.994, 0.997 and 0.994, respectively (Fig. 3F). When the six genes were combined into one regression model, the AUC value was 1.000 (Fig. 3G), indicating this model could more accurately distinguish SCLC from normal tissue than single marker genes.

Next, GSE108055 was used as the test cohort to confirm the expression level of the six candidate marker genes. The gene expression patterns for EZH2, IDH2, MUC1, PARP1, and RRM2 observed in the test cohort were consistent with those observed in the training cohort. EZH2 (p = 3.1 $\times$ 10${}^{-6}$), IDH2 (p = 2.2 $\times$ 10${}^{-5}$), PARP1 (p = 3.1 $\times$ 10${}^{-6}$) and RRM2 (p = 3.1 $\times$ 10${}^{-6}$) had higher expression levels in SCLC samples than in normal samples, whereas MUC1 (p = 0.025) had lower expression in SCLC samples than in normal samples (Fig. 4A).

Fig. 4.

Validation of marker genes in the test cohort. (A) Expression of marker genes in the GSE108055 cohort. (B,C) ROC curves for marker genes in the GSE108055 cohort.

We then used the test cohort GSE108055 to confirm the accuracy of our logistic regression model. The ROC curve is shown in Fig. 4B,C. According to this result, the AUCs of EZH2, IDH2, MUC1, PARP1, and RRM2 were all $>$0.7, while the AUC for the sum of these genes was 1. This confirms the logistic regression model developed earlier can specifically distinguish between SCLC and normal samples.

We performed single-gene GSEA-KEGG pathway analysis and GSEA-GO function analysis to further evaluate the potential for the six marker genes to differentiate SCLC samples from normal samples. The top six results with the most significant enrichment are shown in the figures. From the GSEA-KEGG analysis, the six marker genes had strong associations with the cell cycle, cytokine-cytokine receptor interactions, lysosome, hematopoietic cell lineage, complement and coagulation cascades, as well as some disease pathways (asthma, autoimmune thyroid disease, viral myocarditis, leishmaniasis infection). Moreover, the results showed that EZH2 and IDH2 were associated with the JAK-STAT signaling pathway, while ATG3 and PARP1 were associated with the intestinal immune network for IgA production (Supplementary Fig. 2A–F).

GSEA-GO function analysis revealed that with the exception of MUC1, the other five marker genes were all related to immune response (adaptive immune response, lymphocyte-mediated immunity, T-cell-mediated immunity, immune receptor activity, humoral immune response, and immune response regulation signal pathway) (Supplementary Fig. 3A–F).

Depending on the expression level of each marker gene, the pathways for different activation levels between SCLC and normal samples were analyzed by GSVA. GSVA-KEGG pathway analysis revealed that EZH2 overexpression was associated with ‘ARACHIDONIC ACID METABOLISM’, ‘JAK STAT SIGNALING PATHWAY’, ‘PPAR SIGNALING PATHWAY’, and ‘TYROSINE METABOLISM’. IDH2 up-regulation was associated with ‘PPAR SIGNALING PATHWAY’ and ‘TYROSINE METABOLISM’, while PARP up-regulation was linked to the ‘JAK STAT SIGNALING PATHWAY’ and ‘PPAR SIGNALING PATHWAY’. Furthermore, down-regulation of RRM2 was associated with the ‘PPAR SIGNALING PATHWAY’, ‘JAK STAT SIGNALING PATHWAY’, and ‘ARACHIDONIC ACID METABOLISM’. These pathways were found to have a strong connection to the occurrence and development of lung cancer. Moreover, IDH2 and PARP overexpression could activate the ‘NOD-LIKE RECEPTOR SIGNALING PATHWAY’ and the ‘TOLL-LIKE RECEPTOR SIGNALING PATHWAY’, both of which are related to the immune reaction. Down-regulation of RRM2 could also activate these two pathways. All of the above results are shown in Supplementary Fig. 4A–E.

GSVA-GO functional analysis showed that overexpression of EZH2, IDH2, PARP1, and RRM2, and down-regulation of MUC1 were all partially related to the immune reaction, including the ‘IGA_IMMUNOGLOBULIN_COMPLEX’, ‘INTERLEUKIN-21_PRODUCTION’ and the ‘IPAF_INFLAMMASOME_COMPLEX’ (Supplementary Fig. 5A–E).

GSEA and GSVA analysis showed a close connection between the candidate marker genes and immune response. The CIBERSORT algorithm was therefore used to investigate differences in infiltrating immune cells between SCLC and normal samples. As shown in Fig. 5A, SCLC samples had lower proportions of memory-resting CD4 T cells, naive B cells, monocytes, activated dendritic cells, resting mast cells, and neutrophils compared with normal samples. However, the proportion of gamma delta T cells, T follicular helper cells, and M1 macrophages were higher in SCLC than in normal samples.

Fig. 5.

Immune correlation analysis. (A) Differences in the immune microenvironment between SCLC and normal samples were analyzed with the CIBERSORT algorithm. (B) Pearson’s correlation analysis showed that MUC1 was negatively correlated with activated mast cells and positively correlated with resting mast cells. EZH2 was negatively correlated with neutrophils (*p 0.05, **p 0.01, ***p 0.001). SCLC, Small cell lung cancer.

Associations between immune cells and the six marker genes were subsequently investigated using Pearson’s correlation analysis. MUC1 showed a negative correlation with activated mast cells and a positive correlation with resting mast cells. EZH2 was negatively correlated with neutrophils (Fig. 5B). These findings imply that MUC1 and EZH2 may be involved in determining the immune microenvironment of SCLC.

Targeted drugs for the marker genes were predicted by the DGIdb database, and two-parameter interaction relationships were set to defaults (Supplementary Table 6). Cytoscape software was used to visualize the results (Fig. 6). This analysis found 40 drugs that target the marker genes. Amongst these, four were for MUC1, twelve for RRM2, six for IDH2, twelve for PARP1, and seven for EZH2. Unfortunately, none were found that targeted ATG3. Inhibitors of EZH2 included CPI-1205 and TAZEMETOSTAT, while for IDH2, they included ENASIDENIB. Inhibitors for PARP1 included NIRAPARIB, OLAPARIB, VELIPARIB, INIPARIB, TALAZOPARIB, RUCAPARIB CAMSYLATE, TALAZOPARIB TOSYLATE, with the antagonist including RUCAPARIB and the binder including NIACINAMIDE. Inhibitors of RRM2 included HYDROXYUREA, FLUDARABINE PHOSPHATE, GALLIUM NITRATE, CLOFARABINE, GEMCITABINE HYDROCHLORIDE, GEMCITABINE, TEZACITABINE, and CLADRIBINE.

Fig. 6.

Targeted drugs for the marker genes as predicted by the DGIdb database.

In order to construct a ceRNA network, we first used miRanda, miRDB, and Targetscan databases to predict miRNAs that bind to the six marker genes. We then predicted lncRNAs that bind to the miRNAs using the spongeScan database. Finally, a lncRNA-miRNA-messenger RNA (mRNA) ceRNA regulatory network containing 217 nodes and 245 edges was constructed (Supplementary Fig. 6).

A total of nineteen lncRNAs were found that competitively bind with five miRNAs to regulate ATG3. Amongst these, 12 lncRNAs bound to hsa-miR-338-3p. For the regulation of EZH2, eighteen lncRNAs were found to bind seven miRNAs, with hsa-miR-1207-3p shared by nine lncRNAs, and hsa-miR-20a-3p shared by four lncRNAs. A total of sixteen lncRNAs regulated the expression of IDH2 by binding to four miRNAs, including five lncRNAs that bound to hsa-miR-767-5p and nine that bound to hsa-miR-423-5p. A total of twenty-eight lncRNAs could regulate PARP1 expression by binding to seven miRNAs, of which hsa-miR-539-5p was bound by seven lncRNAs, hsa-miR-7-5p by nine lncRNAs, and hsa-miR-221-5p and hsa-miR-20a-3p by four lncRNAs each. Finally, a total of 39 lncRNAs regulated the expression of RRM2 by binding 12 miRNAs. Amongst these, hsa-let-7a-3p, hsa-miR-150-5p, and hsa-miR-197-3p were bound by four miRNAs each, hsa-miR-342-3p by six miRNAs, and hsa-miR-539-5p by seven miRNAs. Of note, AC091153.4 could simultaneously regulate EZH2 and IDH2 by combining with hsa-miR-101-3p. Moreover, four lncRNAs could simultaneously regulate EZH2 and PARP1 expression by binding to hsa-miR-20a-3p. In addition, RP11-94C24.13 could regulate IDH2 and RRM2 simultaneously by combining with hsa-let-7a-5p. Detailed information regarding this ceRNA network is shown in Supplementary Table 7.