GSE98793 dataset, the expression profile of whole blood samples, was downloaded from the Gene Expression Omnibus (GEO) database. This dataset totally included 128 MDD cases (64 with anxiety symptoms and 64 without) and 64 healthy controls. 64 MDD cases without anxiety symptoms and 64 healthy controls were included in our analysis. All the participants of GSE98793 were from the GlaxoSmithKline–High-Throughput Disease-specific target Identification Program (GSK-HiTDiP) study. The MDD patients were evaluated by the semi-structured Schedule for Clinical Assessment in Neuropsychiatry (SCAN) (Wing et al., 1990), which was administered by trained staff. Furthermore, patients had a diagnosis of recurrent MDD (at least two episodes of depression satisfying DSM-IV or ICD10 criteria) were included as well. The exclusion criteria were as follows: if they had experienced mood incongruent psychotic symptoms, a lifetime history of intravenous drug use or diagnosis of drug dependency, depression secondary to alcohol or substance abuse or depression as clear consequence of medical illnesses or use of medications. Patients with co-morbid anxiety disorders, with the exception of obsessive compulsive and post traumatic stress disorders, were included. Patients with diagnosis of schizophrenia, schizoaffective disorders and other axis I disorders were excluded from the study (Leday et al., 2018). The detail information of participants included in the present study was showed in Table 1. GPL570 (Affymetrix Human Genome U133 Plus 2.0 Array) was detection platforms for GSE98793. Gene symbols were used to annotate the downloaded gene probes, eliminate probes without matching, and retain any gene probes with multiple matching.

“Limma” package (R Foundation for Statistical Computing) was used for gene differential expression analysis with the processed gene expression matrix (Diboun et al., 2006). Before the bioinformatic analysis, all the samples were tested in two batches, and batch information could be extracted from phenotypic data. Thus, removeBatchEffect from the limma package was used to remove the batch effect (Ritchie et al., 2015). We conducted the gene differential expression analysis with set threshold: |log₂(FC)|>0.1 and Benjamini–Hochberg-adjusted p < 0.05. Finally, we identified 2,216 DEGs in GSE98793. A total of 184 PS genes were obtained by inputting the keyword “pyroptosis” in the GeneCards database. The overlap of 2,216 DEGs in GSE98793 dataset and 184 PS genes were defined as pyroptosis-related DEGs, totally 23 pyroptosis-related DEGs acquired.

We applied random forest (RF) and least absolute shrinkage and selection operator (LASSO) to establish a diagnostic model for MDD. First, we used RF to screen candidate MDD-specific genes from pyroptosis-related DEGs. RF is a general technique for the training and prediction of samples based on the classification tree. The number of decision trees (ntree) and the value of MTRY in this study were 300 and 4, respectively. The RF was performed by R package “randomforest” (Kursa, 2014). Subsequently, to reduce the number of genes in the model and to solve the multicollinearity problem in regression analysis, we used LASSO logistic regression to screen feature genes and then construct a diagnostic model for MDD. The “glmnet” package was applied for LASSO algorithm (Friedman et al., 2010). Finally, a receiver operating characteristic (ROC) curve was created to investigate whether the built model could effectively predict MDD.

10-fold cross-validation as internal validation method was performed to confirm the predictive performance of our diagnostic model for MDD (Martinez et al., 2011). We chosed 10-fold cross-validation because 10-fold cross-validation techniques could test all data in the dataset and produce stable predictive accuracy. Therefore, the 10-fold cross-validation method with 2000 iterations of resampling was used for internal validatIon.

GSE76826 was used as an external validation dataset to examine the universality and reliability of the diagnostic model. GSE76826 dataset (expression profiling by array) included 20 MDD patients and 12 healthy controls (Miyata et al., 2016).

In addition to this, we examined the effectiveness of this diagnostic model in other mental illness as well. The GSE38484 dataset was used to analyze the diagnostic model. This dataset, based on the microarray platform of the Illumina HumanHT-12 V3.0 expression beadchip (GPL6947), included 96 healthy controls and 106 schizophrenia patients (Van Eijk et al., 2015).

The same model and the same coefficient were conducted for GSE76826 and GSE38484. The characteristic information of participants in GSE76826 and GSE38484 datasets were exhibited in Tables 2, 3.

GSE53987 was based on the platform of the Affymetrix Human Genome U133 Plus 2.0 Array (Lanz et al., 2019). There were 17 subjects with MDD and 18 healthy controls in this dataset, and three brain regions (hippocampus, prefrontal cortex, and striatum) were included. We used the GSE53987 dataset to verify the diagnostic model. The characteristic information of participants in GSE53987 dataset were showed in Table 4.

Consensus clustering is an algorithm used to identify subgroup members and verify subgroups based on resampling. We performed consensus clustering with PAM mehod (Wilkerson and Hayes, 2010) to identify distinct pyroptosis regulation clusters according to the expression profiles of 23 PS genes. PCA was then used to further validate different regulation clusters.

Single-sample gene set enrichment analysis (ssGSEA) was used to quantify the relative abundance of 28 immune cell types related to immune response. In ssGSEA, the relative abundance of each immune cell was expressed as an enrichment score that was normalized to a uniform distribution of 0–1. A deconvolution approach CIBERSORT (http://cibersort.stanford.edu/) was used to evaluate the relative abundances of 22 distinct leukocyte subsets with gene expression profiles in the blood samples.

We utilized GSVA analysis by “GSVA” package (Hänzelmann et al., 2013) to explore the differentiation in biological processes between different pyroptosis regulation clusters. The well-defined biological signatures were derived from the Hallmark geneset (MSigDB database v7.1) (Denny et al., 2018). The GO annotation for different clusters was performed using the R package “clusterProfiler” (Kursa, 2014) with a false discovery rate (FDR) cutoff of <0.01.

The consensus clustering algorithm classified MDD patients into two distinct pyroptosis regulation subtypes. We next identified DEGs between two different clusters using the “limma” package. Specifically, gene expression data were normalized using “voom” function and then inputted to the “lmFit” and “eBayes” functions to calculate the differential expressed statistics. The selection criteria were an adjusted P of <0.01 and |FC| of >1.0.

To quantitatively analyze the pyroptosis subtypes, PCA was used to quantify the pyroptosis level of individual patients. First, PCA was used to distinguish pyroptosis subtypes. Then, the formula was performed to measure the pyroptosis scores (PS) as following: PS score = PC1_i, where PC1 represents principal component 1 and i represents DEG expression.

The framework and workflow are described as following. First, the characteristic information and gene expression profiles of MDD cases and healthy controls in GSE98793, GSE76826, and GSE53987 were obtained from the GEO database (https://www.st-va.ncbi.nlm.nih.gov/gds/?term=). Using GSE98793 dataset, we developed a novel diagnostic model for MDD based on PS genes by machine learning methods (RF and LASSO). GSE76826 and GSE53987 datasets were used as external validation and postmortem brain tissue samples valiadtion. Furthermore, GSE38484 dataset (96 healthy controls and 106 schizophrenia patients) was used to examine whether this diagnostic model was unique for MDD. Then, immune cell infiltration profiles and immune-related pathway enrichment were compared between healthy controls and MDD cases as well.

To discover the connections between PS genes and MDD subtypes, MDD cases were divided into two subtypes (A and B clusters) by consensus clustering analysis according to the PS genes expression profiles. The particular immune and biological characteristics of these two clusters were observed.

We developed a pyroptosis-related signature score, the “PS-score,” to quantify the pyroptosis phenotype subtype.

The GSE98793 dataset and Genecard database included 2,216 DEGs and 184 PS genes, respectively, among which 23 PS genes with significant expression differences were distributed on chromosomes as shown in Figure 1A. The upregulated genes included C-type lectin member 5 A (CLEC5A); cathelicidin antimicrobial peptide (CAMP); Toll-like receptor 2 (TLR2); adenosine A3 receptor (ADORA3); NOD-, LRR-, and pyrin domain-containing protein 3 (NLRP3); cluster of differentiation 14 (CD14); CD274; NLR family CARD domain containing 4 (NLRC4), NLR family apoptosis inhibitory protein (NAIP); calpain 1 (CAPN1), G-protein coupled estrogen receptor 1 (GPER1); Vitamin D receptor (VDR); fibroblast growth factor 21 (FGF21); interleukin-1 receptor antagonist (IL1RN); forkhead box O3 (FOXO3); serpin family B member 1 (SERPINB1); and activating transcription factor 6 (ATF6). The downregulated genes included high mobility group box 1 (HMGB1), processing of precursor 1 (POP1), immunity-related GTPase family (IRGM), baculoviral IAP repeat-containing protein3 (BIRC3), granzyme A (GZMA), and urotensin-II (UTS2) (Figures 1B,C).

Interactions among the 23 pyroptosis-related DEGs were observed, and each gene was mapped to the STRING database to show their interaction relationships (interaction minimum>0.9, highest confidence) and visualize in Cytoscape.

There were totally 54 interaction relationships among 23 pyroptosis-related DEGs. Figure 1D indicated the protein-protein interaction (PPI) network diagram.

The correlation analysis showed that CD14 was significantly positively correlated with NLRP3, while CD14 was significantly negatively correlated with HMGB1 in all healthy controls and MDD samples (Figure 2). These results indicated that expression imbalances of PS genes played important roles in the occurrence and development of MDD.

Firstly, the RF was used to screen MDD-specific genes that optimally differentiated MDD cases from healthy controls. When the number of decision trees reached 300, the error changes of the three kinds gradually decrease (Figure 3A). The 23 PS genes identified by MeanDecreaseAccuracy and MeanDecreaseGini were showed in Figure 3B. The top 10 of these—UTS2, NLRC4, GZMA, GPER1, IL1RN, CAPN1, NLRP3, HMGB1, ATF6, and ADORA3—were candidate genes.

Secondly, to establish an MDD-specific diagnostic model, LASSO regression was conducted based on the above 10 genes and basic phenotype information (age and sex), we contained eight genes (GPER1, GZMA, HMGB1, IL1RN, NLRC4, NLRP3, UTS2, CAPN1) according to lambda.1se (Figures 3C,D).

Furthermore, to improve the diagnostic efficiency of biomarkers, a novel diagnostic risk score was constructed by multiplying the gene expression. The total risk score was imputed as follows: (1.163 × GPER1 expression level) + (−0.474 × GZMA expression level) + (−0.588 × HMGB1 expression level) + (0.003 × IL1RN expression level) + (1.111 × NLRC4 expression level) + (0.479 × NLRP3 expression level) + (−0.294 × UTS2 expression level) + (0.644 × CAPN1 expression level). Receiver operating characteristic (ROC) analysis was used to evaluate the diagnostic ability of eight genes, which showed a favorable diagnostic value, with an AUC of 0.795 (95% CI 0.721–0.868) (Figure 3E). The coefficient of eight genes were showed in Table 5.

The internal validation of our diagnostic model was performed by 10-fold cross-validation (n = 2000). The model demonstrated good discrimination (bias-corrected AUC = 0.774, 95% CI = 0.679–0.864) (Figure 4A). The external analysis using independent dataset GSE76826 revealed a good performance of our diagnostic model for stratifying MDD patients (AUC = 0.891, 95% CI = 0.657–0.941) (Figure 4B).

Furthermore, the accuracy of our diagnostic model was validated using brain tissue GSE53987 dataset. The AUC was 0.86 (95% CI = 0.669–0.933), which also indicated good performance (Figure 4C).

In order to observe the effectiveness of this diagnostic model in other mental illness, the analysis of schizophrenia GSE38484 dataset was performed. The AUC of GSE38484 dataset was 0.593 (95% CI = 0.507–0.686) (Figure 4D), suggesting that this diagnostic model was not effective and accurate for schizophrenia. We speculated that our diagnostic model constituting of eight PS genes was more appropriate for MDD.

Many studies have indicated that MDD is accompanied by immune dysregulation, pyroptosis is closely related to immune response as well, thus we want to further explore the relationship between pyroptosis and immune response in MDD. The differences in immune cell infiltration profile and immune-related pathway enrichment between healthy controls and MDD cases were quantified.

Firstly, ssGSEA was used to calculate the relative abundance of immune cells in each sample. Eight of the immune cell types discovered significant changes in depression samples, with activated dendritic cells, immature dendritic cells, monocytes, and neutrophils showing upregulation and activated B cells, activated CD4 T-cells, effector memory CD8 T-cells, and immature B cells showing downregulation (Figure 5A). These results suggested significantly altered immune cell profiles in MDD cases comparing with healthy controls (Figure 5B). BIRC3 exhibited the most obvious positive correlation with activated CD4 cells, with low scores in MDD (Figures 5C–E). Meanwhile, HMGB1 showed a significant negative correlation with monocyte. The score for HMGB1 was lower in MDD, while monocytes showed the opposite phenomenon (Figures 5F–H).

Similarly, the enrichment fractions of immune-related pathways were also calculated by ssGSEA (Figure 6A). We observed differences in the antimicrobial pathway, B cell receptor (BCR) signaling pathway, chemokines, cytokines, TCR signaling pathway, and transforming growth factor (TGF)-β family member receptors between healthy controls and MDD cases (Figure 6B). CAPN1 was positively correlated with antimicrobials (Figures 6C–E), while FOXO3 was negatively correlated with the TCR signaling pathway (Figures 6F–H). We speculated that these changes immune cells and immune-related pathways played important roles in occurrence and development of MDD.

To explore the connections between the expression profiles of 23 PS genes (CLEC5A, CAMP, TLR2, ADORA3, NLRP3, CD14, CD274, NLRC4, NAIP, CAPN1, GPER1, VDR, FGF21, IL1RN, FOX O 3, SERPINB1, ATF6, HMGB1, POP1, IRGM, BIRC3, GZMA, and UTS2) and MDD subtypes, we performed unsupervised clustering using data from 64 MDD cases. By increasing the clustering variable (k) from 2 to 10, we observed the highest intra-group correlation and lowest inter-group correlation for k = 2, suggesting that the 64 MDD cases could be divided into two clusters according to these genes. Cluster A included 30 samples and cluster B contained 34 samples (Figures 7A–C).

PCA analysis indicated these two clusters differentiated significantly (Figure 7D). The heatmap and boxplot showed the expression differences of 23 PS genes between two clusters, in which the expressions of ATF6, BIRC3, CAMP, CAPN1, CD14, CLEC5A, GZMA, NLRP3, POP1, SEPRINB1, and TLR2 were significantly increased in cluster B, while the expressions of FGF21, FOX O 3, HMGB1, IL1RN, and VDR were significantly decreased. ADORA3, CD274, GPER1, IRGM, NAIP, NARC4, and UTS2 did not differ significantly (Figures 7E,F).

To explore the immune characteristics between these two clusters, we measured the immune cells enrichment fractions, immune pathway activity, and human leukocyte antigen (HLA) gene expression profiles. The two clusters revealed completely different immune characteristics: for example, memory B cells and cytokines were mainly enriched in cluster A, while HLA genes mostly concentrated in cluster B (Figures 8A–C).

In addition to immune characteristics, we further explored the biological functions. We analyzed the enrichment scores of the Hallmark pathways and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways in the two clusters by GSVA, which revealed enrichment of pathways in cluster A, including fatty acid metabolism and neuroactive receptor interaction (Figures 8D,E).

To further explore the pathological mechanisms of MDD related to pyroptosis, we developed a pyroptosis-related signature score, the “PS-score,” including the phenotype-related genes to quantify the pyroptosis regulation pattern of each MDD sample. The “PS-scores” for two distinct subtypes were calculated. The “PS-scores” of cluster B was significantly higher than that of cluster A (Figure 9A). Figure 9B exhibited that the relationships between gender, pyroptosis clusters, and “PS-scores”.

IFN-α is currently suggested to be an important link in the pathogenesis of MDD. Clinical studies have shown significantly higher serum IFN levels in depression patients than those of normal people. After 12 weeks of antidepressant treatment, IFN levels were significantly lower than those before treatment (Shelton et al., 2011). Thus, we explored the relationship between two subtypes and expression profiles of IFN-α isoforms. Figure 9C indicated that the expression profiles of IFN-α1, IFN-α2, IFN-α4, IFN-α5, IFN-α7, IFN-α8, IFN-α10, IFN-α14, IFN-α16, IFN-α17, and IFN-α21 were all lower in cluster B, although the difference was not significant for IFN-α17.