The Gene Expression Omnibus (GEO) is a public repository that archives and distributes high-throughput gene expression and other functional genomics data sets, with web-based tools for data visualization and analysis [22]. We downloaded the gene expression datasets, GSE98793 [23] and GSE32280 [24], which are associated with major depressive disorder (MDD) patients from the GEO database using the R package GEOquery [25] The species source for both the GSE98793 and GSE32280 datasets was Homo sapiens and the data platform was GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array. We also downloaded the GSE24982 [26] and GSE30691 [27] gene expression profiles of Neuropathic Pain (NP) patients, in which the species source for both was Rattus norvegicus. The data platform for the GSE24982 dataset was GPL1355 [Rat230_2] Affymetrix Rat Genome 230 2.0 Array, whereas the data platform for the GSE30691 dataset was the GPL85 [RG_U34A] Affymetrix Rat Genome U34 Array.

The chip GPL platform file was used for all dataset probe name annotations. We selected data from 128 patients with major depressive disorder (MDD) (MDD group, group: MDD) and 64 healthy controls (control group, group: Control) from the GSE98793 datasets of whole blood samples for inclusion in the analysis. A total of 16 subjects from the GSE32280 datasets were analyzed, including 8 examples of peripheral blood lymphocytes from matched healthy controls (control group, group: Control) and 8 peripheral blood lymphocyte samples from MDD patients. In addition, in the Neuropathic Pain (NP) dataset, GSE24982, we used a total of 40 mRNA data samples based on the L4 and L5 Dorsal Root Ganglion (DRG), including 20 control (group: Control) mRNA samples and 20 mRNA samples from the spinal nerve ligation model of neuropathic pain in rats (NP group, group: NP). The GSE30691 dataset contains mRNA samples from 11 rat neuropathic pain spinal nerve ligation (Ch) models (NP group, group: NP) and 9 Sham (SH) control (group: Control) mRNA data samples. The specific grouping of the information from these datasets is listed in Table 1. The datasets we selected are all expression profiling by array, and the sample size is sufficient to meet our analysis requirements. In addition, to ensure that our control samples match the experimental conditions, the healthy control samples and diseased samples in these datasets are as consistent as possible in terms of sample collection methods, time, and other experimental conditions to minimize the impact of confounding factors.

The MDD (GSE98793, GSE32280) and NP (GSE24982, GSE30691) datasets were combined for analysis. To minimize the variance in sample combinations across batches, we first standardized the datasets internally using the ControlizeBetweenArrays function of the R limma package [28]. We corrected the combined data for batch effects using the remove batch effect function, which enabled us to obtain the combined MDD and NP datasets. The MDD datasets contained 136 cases (disease group, group: Case/MDD) and 72 controls (control group, group: Control). The NP datasets contained 31 cases (disease group, group: Case/NP) and 29 controls (control group, group: Control). The expression values for the samples in the MDD and NP datasets were then analyzed using Principal Component Analysis (PCA) before and after correction [29].

Differentially expressed genes (DEGs) are a subset of genes that express differently among experimental conditions, used to determine biological functions or predict therapeutic outcomes. To identify the potential mechanism of action of DEGs in MDD and NP and the associated biological features and pathways, we performed a differential expression analysis on the case and control groups using the R limma package. Genes that met the criteria |logFC| >0 and p-value <0.05 were considered DEGs for subsequent studies. |logFC| represents the absolute log2 value of the fold change in gene expression. To obtain the common differentially expressed genes (Co-DEGs) associated with MDD and NP, we selected the DEGs from the differential analysis of the MDD and NP datasets DEGs and constructed a Venn diagram. A volcano plot was generated using ggplot2 of the R package and a heat map was used to display the results.

A protein-protein interaction (PPI) network is made up of different proteins that engage with one another to function in a variety of biological processes, including signaling, regulation of gene expression, substance metabolism, energy production, and cell cycle regulation. A database for identifying known proteins and predicting protein interactions is the STRING database [30]. We constructed a PPI network of Co-DEGs linked to both NP and MDD disease [minimum required interaction score: middle confidence (0.400)], which was visualized with Cytoscape [31] (version 3.9.1). Using the MCODE plugin, we mined the hub nodes with connections to other PPI network nodes (K score: 2, Cutoff degree: 2, Cutoff node score: 0.2) [32]. These nodes were highly interconnected with one another and may play a role in regulating various biological processes associated with NP and MDD.

A molecular complex with a particular biological activity may be represented by closely connected local areas in the PPI network. We also used four algorithms to mine the scores of Co-DEGs in PPI networks that are connected to other PPI network nodes, which included MCC (Matthews Correlation Coefficient metric) [33], MNC (the maximal neighborhood coefficient), EPC (edge percolated component), and Closeness. According to our rankings of the Co-DEGs. The top 20 Co-DEGs across the four algorithms were considered as hub genes (hub genes, mRNA).

Large-scale functional enrichment studies, including biological process (BP), molecular function (MF), and cellular component (CC), are frequently carried out using the Gene Ontology (GO) database [34]. The Kyoto Encyclopedia of Genes and Genomes (KEGG) is a widely used database that genomic data for biological pathways, diseases, and drugs [35]. We conducted GO and KEGG enrichment analyses for MDD and NP disease-associated Co-DEGs, PPI network hub nodes, and hub genes in the PPI network using the clusterProfiler [36]. p values < 0.05 and FDR values (q value) < 0.05, which were considered statistically significant as a criteria for selection and Benjamini-Hochberg was used to correct p values.

To more precisely measure the percentage of various immune cells in samples associated with MDD and NP. Single-sample gene-set enrichment analysis (ssGSEA) was used to identify highly sensitive and specific differentiation of the various human immune cells in the tumor microenvironment (TME). The algorithm generated a set of 28 genes to mark different tumor-infiltrating resistant cell types from a study of published tumor immune infiltration articles [37, 38]. The degree of immune cell infiltration in each sample was represented by an enrichment score computed by ssGSEA in the GSVA tool of R. The varying abundance of infiltrating immune cells between the Case (MDD/NP) group and the Control group are shown using a heat map and a boxplot for the MDD and NP datasets. We then calculated the correlation between immune cells and Co-hub genes in different disease samples by combining the gene expression matrix of both disease samples. The correlation heat map was displayed using the R package heatmap.

Finally, the expression of Co-hub genes in the MDD and NP datasets was examined using the Spearman method. The results of the correlation analysis were presented by plotting the correlation heat map. Next, we selected the results for the Co-hub genes with the same trend and displayed them by plotting the correlation scatter graph with the R package ggplot2.

We then established a group comparison graph for the Co-hub genes in various groups (Case/Control) of the MDD and NP datasets. The Receiver operating characteristic (ROC) Curve [39] is a composite indicator that represents continuous variables of sensitivity and specificity. We used the R package pROC program to plot the ROC curves of the Co-hub genes in the MDD and NP datasets and calculated the Area Under the Curve (AUC) of the ROC curve to determine the diagnostic significance of the Co-hub genes. The Receiver operating characteristic (ROC) curve AUC typically falls between 0.5 and 1. The diagnostic impact increases as the AUC gets near 1. The AUC exhibits low accuracy in the range of 0.5–0.7, some accuracy in the range of 0.7–0.9, and high accuracy in the range of 0.9 or above.

R software was used for data processing and analysis in this paper (Version 4.2.2). We calculated the normally distributed variables using an independent Student t-test to compare two sets of continuous variables. We used the Mann-Whitney U test to examine differences among factors with non-normal distributions (i.e., Wilcoxon rank sum test). If not explicitly indicated, a Spearman correlation was used to calculate the correlation coefficients between different molecules. The p values were two-sided with p < 0.05 being the threshold for statistical significance.

The technical route of this study is shown in Figure 1. Using the ControlizeBetweenArrays function of the limma package, we normalized the two major MDD datasets (GSE98793 and GSE32280) and the two NP datasets (GSE24982 and GSE30691), respectively. The MDD datasets (Figures 2A, B) and the NP datasets (Figures 2E, F) were obtained through a batch effect correction of the combined data using the “remove batch effect” function. The MDD datasets include 136 MDD samples (group: MDD) and 72 control samples (group: Control). The NP datasets consist of 31 NP samples (group: NP) and 29 control samples (group: Control). In addition, to convert mouse to human genes for subsequent analysis, the R package homologene was used to conduct an ID transformation of the NP datasets.

To verify the effect of removing the batch effect (Figures 2C, D, G, H), we grouped the MDD, and NP datasets according to the source of the samples. For the dataset expression matrix, before, and after the batch effect was eliminated, we performed a Principal Component Analysis (PCA). The results indicated that after the batch removal process, the batch effect was essentially eliminated from the MDD and NP datasets.

Following Principal Component Analysis (PCA), we obtained DEGs between different groups of the MDD and NP datasets. The MDD datasets yielded 21,655 DEGs, of which 1708 met the criteria. For the Case/MDD group, there were 824 upregulated DEGs and 884 downregulated DEGs. The NP datasets yielded 4076 DEGs, of which 2,136 met the criteria. For the Case/NP group, 978 exhibited high expression (low expression in the Control group with positive logFC), and 1,158 genes exhibited low expression (high expression in the Control group with negative logFC). The results for the MDD and NP datasets are depicted in volcano plots (Figures 3A, B).

To obtain Co-DEGs for MDD and NP, we focused on the intersection of 48 upregulated and 45 downregulated Co-DEGs that were presented in a Venn diagram (Figures 3C, D). The data and annotation for these Co-DEGs is listed in Table 2. We examined differential expression of 93 Co-DEGs from the MDD and NP datasets in various groups. Heat maps were used to display the results of the differential analysis using the R package heatmap (Figures 3E, F). As shown in Figure 3, 93 Co-DEGs showed significant differences in expression between the groups based on the MDD and NP datasets.

For the 93 Co-DEGs, we performed a functional enrichment analysis using the Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway databases (Table 3). The results indicated that 93 Co-DEGs were primarily enriched in biological processes (BP), such as fibrinolysis, regulation of inflammatory response, and developmental maturation. They were also enriched in cellular component (CC) tertiary granules and complement and coagulation cascade from the KEGG pathway database. The results are shown in the form of bubble plots and a network diagram (Figures 4A, B). We then combined the logFC values with the enrichment analysis, which generates a z score by providing the logFC values for the 93 Co-DEGs. The effects of GO and KEGG enrichment analysis by joint logFC are shown as circle plots (Figure 4C) and bubble plots (Figure 4D). The results indicated that the 93 Co-DEGs from the MDD dataset were primarily located in the BP pathway (Figures 4C, D).

After excluding the Co-DEGs that did not have a connection with other nodes, we constructed a PPI network (Figure 5A) consisting of 54 Co-DEGs using the STRING database. We analyzed the nodes that have connections with other nodes in the PPI network using the MCODE plugin. We then used the genes in cluster1 and cluster2 of the results as hub nodes for the Co-DEG PPI network, in which we obtained an MCODE cluster network consisting of 6 Co-DEGs (PLAU, TIMP2, HGF, ANGPT, MMP9, EPO) (Score = 4.8) (Figure 5B) and an MCODE cluster network (Score = 3) (Figure 5C) consisting of 3 Co-DEGs (RPL30, TTC37, WDR61). The 9 (hub node) genes included PLAU, TIMP2, HGF, ANGPT, MMP9, EPO, RPL30, TTC37, and WDR61. These nodes warrant additional study as they may be important in regulating the entire BP.

On the 9 hub nodes (PLAU, TIMP2, HGF, ANGPT, MMP9, EPO, RPL30, TTC37, WDR61), we performed GO and KEGG enrichment analyses (Supplementary Table S1). The results indicated that the BP, negative regulation of the apoptotic signaling pathway, was largely abundant in the 9 hub nodes. Enrichment in CC, such as tertiary granules, and transcriptionally active chromatin as well as MF, such as serine-type endopeptidase activity, was observed. The KEGG pathways, including proteoglycans in cancer, PI3K-Akt signaling pathway, and RNA degradation, were also enriched. The results are presented using bar graphs in Figure 5D. In addition, the GO analysis outcomes are shown for the BP pathway (Figure 5E), CC pathway (Figure 5F), and MF pathway (Figure 5G) as a network diagram, whereas the results of KEGG pathway enrichment analysis are displayed as a circular network diagram (Figure 5H).

We calculated the PPI network using the cytoHubba plugin for Cytoscape using four algorithms: MCC (Matthews Correlation Coefficient metric), MNC, EPC (edge percolated component), and Closeness. The top 20 Co-DEGs with the best scores were selected to further identify the hub genes in the Co-DEGs PPI network (Figures 6A–D). The color of the dotted blocks in the graph, from yellow to red, represents a gradual increase in rating. Next, we focused on the intersection of the top 20 Co-DEGs obtained by each of the four algorithms (MCC, MNC, EPC, Closeness) to obtain the hub genes and drew Venn diagrams to display the results (Figure 6E). We obtained a total of 8 hub genes, which included ANGPT2, BMP2, CEBPA, EPO, HGF, MMP9, PLAU, and TIMP2.

We performed GO and KEGG enrichment analyses for the hub genes (Supplementary Table S2). The results indicated that the 8 hub genes were primarily enriched in BP, such as response to hypoxia. They were also enriched in CC, such as tertiary granules, and MF, such as serine-type endopeptidase activity. The KEGG pathways identified included transcriptional misregulation in cancer. We displayed the findings of the GO and KEGG enrichment analyses using bubble plots (Figure 6F). We also displayed the GO gene functional enrichment analysis findings (Figure 6G) as a network diagram and the KEGG pathway enrichment analysis results (Figure 6H) as a circular network diagram.

We first obtained Co-hub genes by taking the intersection of the 9 hub nodes (PLAU, TIMP2, HGF, ANGPT, MMP9, EPO, RPL30, TTC37, WDR61) and the 8 hub genes (ANGPT2, BMP2, CEBPA, EPO, HGF, MMP9, PLAU, TIMP2) (Figure 7A). As shown in Figure 7A, we obtained 6 Co-hub genes, which included ANGPT2, EPO, HGF, MMP9, PLAU, and TIMP2. We analyzed their interaction using the STRING database and visualized them with Cytoscape (Figure 7B). Next, the 6 Co-hub genes (ANGPT2, EPO, HGF, MMP9, PLAU, TIMP2) were analyzed for functional similarity. We calculated the semantic similarity among GO terms, sets of GO terms, gene products, and gene clusters using the R package GOSemSim. The results of this functional similarity analysis between the 6 Co-hub genes are presented as a boxplot in Figure 7C. HGF had the highest value of functional similarity with the other Co-hub genes among the 6 Co-hub genes.

Using the GeneMANIA database, we constructed a Co-hub gene PPI network after retaining the nodes with links to the 6 Co-hub genes (Figure 7D). As shown in Figure 7D, there are five types of interactions between the nodes of our constructed Co-hub genes PPI network and the 6 Co-hub genes, including Co-expression, Pathway, Physical Interactions, Shared protein domains, and Predicted.

The relative infiltration of 28 immune cell types in the disease group (Case/MDD/NP) and control group samples in the MDD and NP dataset expression matrix were determined using the ssGSEA algorithm. We presented the findings for both datasets using a complex heat map (Figures 8A, B). The results indicated that among the disease and control group samples of MDD (Figure 8A) and NP (Figure 8B), the infiltration by the 28 cell types varied significantly. Next, we used the Mann-Whitney U test to analyze the differential degree of infiltration of these cells among the different groups (Case/Control) in the two datasets and the results are grouped in comparison plots (Figures 9A, B). The control group exhibited a higher abundance of activated B cells, effector memory CD8 T cells, memory B cells, and type 1 T helper cells compared with the disease group. In contrast, the disease group contained a more activated dendritic cells, macrophages, and regulatory T cells compared with the control group (Figure 9A).

CD56dim natural killer cells, natural killer T cells, regulatory T cells, and type 17 T helper cells had higher infiltration in the NP dataset control group compared with that in the disease group. In contrast, the disease group had higher infiltration of activated CD8 T cells, activated dendritic cells, effector memory CD8 T cells, eosinophils, gamma delta T cells, mast cells, MDSCs, memory B cells, monocyte neutrophils, plasmacytoid dendritic cells, T follicular helper cells, and type 2 T helper cells (Figure 9B).

Figures 9A, B shows that there were statistically significant differences (p < 0.05) between the relative immune infiltration of the MDD and NP dataset samples compared with control group (Case/Control) samples for activated dendritic cells, effector memory CD8 T cells, memory B cells, and regulatory T cells.

We also calculated the correlation for these 4 immune cell types with the expression of the 6 Co-hub genes (ANGPT2, EPO, HGF, MMP9, PLAU, TIMP2) in the MDD and NP dataset disease samples (Figures 9C, D). The results indicated that the expression of these 6 genes and the relative abundance of the four immune cells tended to be significantly positive and less negatively correlated (p < 0.05) in the MDD dataset samples (Figure 9C). The results of the correlation analysis between the infiltration of the 4 immune cells and the expression of 6 the Co-hub genes in the NP dataset disease samples revealed that there was a significant positive correlation (p < 0.05) between these infiltrating cells and the 6 genes. Activated dendritic cells showed the highest correlation with the expression of these genes (Figure 9D).

Finally, the expression of activated dendritic immune cells, memory B cells and the 6 Co-hub genes had a significant positive correlation in the MDD and NP dataset disease samples (Figures 9C, D). In contrast, the expression of the 6 Co-hub genes was significantly negatively correlated with the infiltration abundance of effector memory CD8 T immune cells in these datasets (Figures 9C, D).

The correlation in expression between the 6 Co-hub genes (ANGPT2, EPO, HGF, MMP9, PLAU, TIMP2) was analyzed in the MDD and NP datasets by Spearman correlation analysis and presented as a correlation heat map (Figures 10A, B). There was a positive correlation among most Co-hub genes in both datasets; however, quite a few correlations of Co-hub genes in the NP datasets did not show a statistical difference (p > 0.05) (Figures 10A, B).

We selected interaction pairs that were statistically different (p < 0.05) with the same trend in the MDD and NP datasets for further analysis. Based on the Spearman algorithm, we generated scatter plots to show the correlation analysis results for the Co-hub genes, ANGPT2 and HGF, in the MDD and NP datasets as well as the correlation analysis result of PLAU and TIMP2 in the two datasets (Figures 10C–F). Among them, there were negative correlations between the expression of ANGPT2 and HGF (R = −0.21, Figure 10C), (R = 0.44, Figure 10D) in both datasets. Conversely, the expression of PLAU and TIMP2 (R = 0.23, Figure 10E), (R = 0.44, Figure 10F) had a positive correlation in both datasets.

For both datasets, we also examined the expression differences for the 6 Co-hub genes (ANGPT2, EPO, HGF, MMP9, PLAU, and TIMP2) between the disease group (group: Case/MDD/NP) and the control group (group: Control). The results are presented through grouped comparison plots (Figures 11A, B). Of the 6 common hub genes, only ANGPT2, MMP9, PLAU, and TIMP2 were statistically significantly different (p < 0.05) in the different groups of the MDD datasets. Although MMP9, PLAU, and TIMP2 expression were all significantly higher in the disease group of the MDD datasets compared with the control group, ANGPT2 expression was significantly lower in the disease group compared with the control group (Figure 11A). The 6 common hub genes in the NP datasets were all statistically significantly (p < 0.05) different in the different groups. Among them, ANGPT2 gene expression was upregulated in the normal group, whereas EPO, HGF, MMP9, PLAU, and TIMP2 were upregulated in the disease group (Figure 11B).

Next, we plotted the ROC curves for the four Co-hub genes (ANGPT2, MMP9, PLAU, TIMP2) in both datasets (Figures 11C–J). As shown in Figures 11C–F, the expression of ANGPT2 (AUC = 0.594, Figure 11C), MMP9 (AUC = 0.595, Figure 11D), PLAU (AUC = 0.586, Figure 11E), and TIMP2 (AUC = 0.659, Figure 11F) all had low accuracy for the diagnosis of MDD disease, whereas the expression of the Co-hub genes, ANGPT2 (AUC = 0.734, Figure 11G) and MMP9 (AUC = 0.761, Figure 11H), exhibited a certain accuracy for the diagnosis of NP disease. In contrast, the expression of PLAU (AUC = 0.663, Figure 11I) and TIMP2 (AUC = 0.623, Figure 11J) all had low accuracy for the diagnosis of NP disease.

To verify the expression differences of 4 co-hub genes (ANGPT2, MMP9, PLAU, TIMP2) in MDD and NP datasets, we also used the Wilcoxon rank sum test to analyze the expression levels of these 4 co-hub genes in GSE98793 and GSE24982 datasets. The differences between the disease group (group: Case/MDD/NP) and the control group (group: Control) were shown through comparison charts (Figures 12A, B). As shown in Figures 12A, 11B, MMP9, PLAU, and TIMP2 all had statistically significant differences (p < 0.05) in different groups of the MDD datasets. The expression levels of MMP9, PLAU, and TIMP2 in the disease group of GSE98793 dataset were significantly higher than those in the control group (Figure 12A). In the NP dataset, ANGPT2, MMP9, and PLAU all had statistically significant differences (p < 0.05) between different groups. The expression level of ANGPT2 was upregulated in the control group while the expression levels of MMP9 and PLAU were upregulated in the disease group (Figure 12B).

We then plotted the ROC curves of the 4 co-hub genes (ANGPT2, MMP9, PLAU, TIMP2) in GSE98793, GSE32280, GSE24982 and GSE30691 datasets and presented the results (Figures 12C–R). As shown in Figures 12C–F, the expression of co-hub genes ANGPT2 (AUC = 0.582, Figure 12C), MMP9 (AUC = 0.600, Figure 12D), PLAU (AUC = 0.602, Figure 12E), and TIMP2 (AUC = 0.685, Figure 12F) all had low accuracy in diagnosing MDD in the GSE98793 dataset. Similarly, as shown in Figures 12G–J, the expression of co-hub genes ANGPT2 (AUC = 0.609, Figure 12G), MMP9 (AUC = 0.578, Figure 12H), PLAU (AUC = 0.641, Figure 12I), and TIMP2 (AUC = 0.547, Figure 12J) all had low accuracy in diagnosing MDD in the GSE32280 dataset.

As shown in Figures 12K–N, the expression of co-hub gene MMP9 (AUC = 0.948, Figure 12L) had high accuracy in diagnosing NP in the GSE24982 dataset. The expression of ANGPT2 (AUC = 0.735, Figure 12K) and PLAU (AUC = 0.857, Figure 12M) had some accuracy in diagnosing NP in the GSE24982 dataset while the expression of TIMP2 (AUC = 0.605, Figure 12N) had low accuracy in diagnosing NP in the GSE24982 dataset.

As shown in Figures 12O–R, the expression of common hub genes MMP9 (AUC = 0.717, Figure 12P), PLAU (AUC = 0.889, Figure 12Q), and TIMP2 (AUC = 0.758, Figure 12R) had some accuracy in diagnosing NP in the GSE30691 dataset while the expression of ANGPT2 (AUC = 0.535, Figure 12O) had low accuracy in diagnosing NP in the GSE30691 dataset. Finally, we summarized the clinical information of different groups in the GSE98793 dataset and presented the results in a clinical data table (Supplementary Table S3).