We searched and selected two microarray datasets in the GEO database. GSE189993 and GSE141022 were downloaded, and the gene expression data of the middle cerebral artery were analyzed. GSE189993 consists of 21 MMD and 11 non-MMD patients, and GSE141022 consists of four MMD and four non-MMD patients. After normalization, we used the SVA package in R script for background correction, and the ComBat function was used to remove the batch effects. After correction, the distribution of the data was visualized by principal component analysis (PCA), and PCA illustrated that different batches were at the same level. We used the mean as the unique expression value for duplicate gene symbols and the limma package to identify DEGs. p < 0.05 was considered significant. |log2 FC (fold change) | >1 was the criterion for identifying DEGs.

GSEA is a bioinformatics tool used to analyze enriched gene sets, KEGG pathways, and immunological signatures, as well as for analyzing genome-wide expression data. Genes in the gene set share the same function and other features. GSEA could screen the functional enrichment score for all expression profiles rather than different expression genes alone. A good score was presented in terms of rank-based and p-value-based pathway enrichment, and no false prediction was found (Nguyen et al., 2019). GSEA software was obtained from the official website (Subramanian et al., 2005). All gene expression profiles were divided into two groups, and the c2.cp.kegg.v7.4.symbols.gmt subset obtained from MsigDB (Liberzon et al., 2011) was used to evaluate molecular mechanisms. The minimum and maximum number of gene sets ranged from 5 to 5000, and 1000 was set as the threshold of permutations. p < 0.05 and FDR<0.25 were considered significant. To further explore the immune pathways, the c7.immunesigdb.v7.4.symbols.gmt subset was also obtained from MsigDB to evaluate immune-related pathways. The minimum and maximum number of gene sets ranged from 5 to 5000, and 1000 was set as the threshold of permutations. p < 0.01 and FDR<0.1 were considered significant.

We used the STRING online database to construct the PPI network (von Mering et al., 2003). The network was drawn by Cytoscape based on the results of the STRING analysis, and the most important module was identified by MCODE plugin (Bader and Hogue, 2003), and hub genes were identified by the hub gene plugin. The top 50 genes were selected as hub genes from five algorithms, namely, Closeness, Degree, MCC, MNC, and EPC. The real core genes were obtained from the intersection of the result. The default parameters for MCODE were selected as the criterion.

WGCNA is usually used to analyze the gene expression patterns. According to the WGCNA algorithm, we performed Pearson’s correlation and average linkage first. β was known as the soft-thresholding parameter, and it emphasized strong correlations between genes and penalized weak correlations. Next, we transformed the adjacency matrix into a topological overlap matrix (TOM) and calculated the corresponding dissimilarity (1-TOM) after confirming the soft-thresholding power. A gene’s network connectivity could be measured by TOM. To classify the genes sharing similarity into gene modules, we then carried out average linkage hierarchical clustering in accordance with TOM-based dissimilarity, and gene dendrograms should have had a minimum size of 30. To further explore these modules, the dissimilarity of the module eigengenes was computed, and we merged some modules in the module dendrogram after choosing a cut line. Last, the relationship between specific clinical features and the modules was analyzed.

Gene Ontology (GO) and KEGG are the most important bioinformatics tools for the analysis of gene annotation and biological processes. First, the identified PPI hub genes and WGCNA gene set were merged together to obtain the combined core genes. Based on these core genes, GO terms and KEGG pathway enrichment were performed, respectively. Next, we obtained immune gene symbol lists from the website. In order to identify immune DEGs, ImmPort (The Immunology Database and Analysis Portal database) was selected to identify immune-related genes. The ImmPort gene set contains 2,489 immune-related genes. The overlapping genes were obtained from the intersection of three datasets, hub genes of DEGs, and hub genes of the most important molecular and immune-related gene set. GO terms and KEGG pathways were analyzed as well based on these immune-related core genes. p < 0.05 was considered significant for GO terms and KEGG pathway analysis.

Finally, another dataset (GSE157628) was utilized to validate immune-related core genes.

Expression values from the original combined expression matrix in samples of GSE189993 and GSE141022 were extracted. Expression values of the same genes from microarray data in samples of GSE157628 were extracted as well. The expression pattern was compared between the two groups. To evaluate diagnostic accuracy, the receiver operating characteristic (ROC) curve was applied, and we calculated the area under the curve (AUC). Ideally, the AUC should be 1, and values over 0.5 are considered predictive.

The cell-type identification by estimating relative subsets of the RNA transcript algorithm, also called CIBERSORT (Newman et al., 2015), is widely used in evaluating cell fractions. Based on our gene expression profile, we selected the computational method to calculate the score of immune infiltrating cells. The fractions of 22 types of infiltrating immune cells in each sample were visualized by a bar graph, and the correlation between immune cell subtypes and gene expression was visualized by a correlation heat map.

To further confirm immune-related genes as potential diagnostic immune biomarkers, we divided the 60 samples into the MMD and non-MMD groups; the MMD group had 35 patients and the non-MMD group had 25 patients. The immune scores of 22 infiltrating immune cell subtypes were compared. After analyzing immune infiltration, we performed Spearman’s correlation to determine the relationship between the expression levels of identified immune biomarkers and 22 subtypes of infiltrating immune cells.

Data analysis and statistical analysis were conducted by R Studio (version 4.1.2). R package “ggplot2” and “pheatmap” and other packages were used to generate the pictures of volcano plots, heat maps, and principal component analysis. Box plots were drawn by GraphPad Prism 7.0 software using the Mann–Whitney U test. The ROC curve was performed with SPSS 17.0 software. The results were imported into Sangerbox (Shen et al., 2022), and the data were visualized. In the statistical analysis, we selected False Discovery Rate (FDR) for multiple corrections, and p < 0.05 was considered statistically significant.

The raw data were processed, and the results showed that the samples from each dataset were comparable (Figures 2A, B). After normalization, we found 348 DEGs between MMD and control groups, which consisted of 89 upregulated genes and 259 downregulated genes. The analysis results of the expression matrix are depicted in the volcano plot (Figure 2C). The top 50 and last 50 DEG genes are depicted in the heat map (Figure 2D) (Supplementary Table S1).

As expected, the exported results of GSEA-based enrichment analysis illustrated a strong correlation between immune cells and immune activities. According to the KEGG enrichment analysis, there were 15 enriched terms, and they were all upregulated (p < 0.05 and FDR<0.25), mainly involving primary immune-related pathways, including natural killer cell-mediated cytotoxicity, antigen processing and presentation, primary immunodeficiency, and the Fc epsilon RI signaling pathway; pathways related to immune-associated diseases include autoimmune thyroid disease, type-I diabetes mellitus, graft versus host disease, asthma, and allograft rejection; pathways related to cytokine activation include cytokine–cytokine receptor interaction and neuroactive ligand receptor interaction. The immune-related pathways results are visualized in Figure 3 (Supplementary Table S2).

To further characterize the signature of immune activities related to the pathology of MMD, GSEA was used again to identify the immunological pathways enriched in immune signature gene sets. The criteria p < 0.01 and FDR <0.15 showed that 20 prominent immune pathways were enriched totally. There were 6 activated pathways (Figure 4A) and 14 suppressed pathways (Figure 4B) among them. The gene sets in CD4 T cell, CD8 T cell, Th1 cells, monocytes, mast cells, dendritic cells (DCs), and peripheral blood mononuclear cells were mainly enriched (Supplementary Table 3).

In order to identify the hub genes, 384 DEGs were analyzed with the PPI network (Figure 5A). The three most significant modules were clustered via MCODE (Figure 5B), module 1 was made up of six DEGs, and module 2 and module 3 consisted of three DEGs. According to five commonly used classification algorithms in cytoHubba, the top 50 genes were identified by each algorithm. Finally, we took the intersection of Closeness (top 50), Degree (top 50), MCC (top 50), MNC (top 50), and EPC (top 50) to acquire 30 hub genes (Figure 5C) (Supplementary Table S4).

The co-expression modules were generated after WGCNA. Gene sets contained in co-expression modules frequently shared high topological overlap similarity. To guarantee a scale-free network, we chose the soft-threshold β = 4 (scale-free R2 = 0.87) in the present study. The adjacency matrix was converted into a TOM matrix (Figure 6A), which showed the similarity between different nodes because the weighted correlation was considered in the matrix. Ultimately, we obtained 34 co-expression modules by WGCNA (Figure 6B). In order to select the most important clinical module, we performed correlation analysis between the identified modules and clinical traits (Figure 6C). Among these modules, the thistle^l module was selected, and it showed a strong correlation with the MMD group (cor = 0.59 and p = 4.2 e–4). Automatically, all the 31 genes contained in this module were identified as the hub genes (Figure 6D) (Supplementary Table S5).

First, we obtained 61 core genes with the combination of the aforementioned two hub gene sets (Figure 7A) (Supplementary Table S6). In order to study the biological process of the core genes, GO and KEGG functional analyses were conducted. The biological process analysis of the hub genes was mainly related to the immune response, immune system process, regulation of the immune response and immune system process, vesicle-mediated transport, exocytosis, secretion, and so on (Figure 8A). Cell component analysis of the hub genes was primarily enriched in the extracellular matrix, cytoplasmic vesicle, intracellular vesicle, secretory vesicle, secretory granule, plasma membrane part, and so on (Figure 8B). Molecular function analysis was primarily enriched in the protein tyrosine kinase activity, phosphatidylinositol phosphate binding, phosphatidylinositol binding, phosphorylated amino acid binding, ion channel binding, titin binding, and so on (Figure 8C). KEGG analysis revealed that the proteoglycans in cancer, human T-cell leukemia virus I infection, phospholipase D signaling pathway, Jak-STAT signaling pathway, Fc epsilon RI signaling pathway, and ErbB signaling pathway were mainly enriched (Figure 8D). The results of the functional analysis revealed that immune activities may be mainly involved in the progression of MMD. The detailed results are illustrated in Figure 8.

Subsequently, we identified 13 immune-related genes overlapped in the PPI hub genes, WGCNA hub genes, and the immune genes set, as shown in the Venn plot (Figure 7B), including CTSS, NPPA, PTPN11, FGR, SYK, TLR8, IL13, CSF2RB, CCR1, MAP2K1, CXCR1, and TGFB2. To screen biological differences of the immune-related genes between MMD and the control group, GO and KEGG functional analyses for the identified immune genes were performed as well and showed results similar to those of the core genes. GO analysis showed that these immune-related hub genes mainly involved the immune response and immune system process and their regulation, vesicle-mediated transport, exocytosis, secretion in BP (Figure 9A), extracellular matrix, cytoplasmic vesicle, intracellular vesicle, lysosome, mast cell granule and immunoglobulin complex in CC (Figure 9B), protein tyrosine kinase activity, signaling receptor binding, cytokine receptor activity/binding, type-III transforming growth factor beta receptor binding, and G protein-coupled peptide receptor activity in MF (Figure 9C). KEGG analysis revealed that these immune hub genes primarily involved apoptosis, cytokine–cytokine receptor interaction, phospholipase D signaling pathway, Jak-STAT signaling pathway, NF-kappa B signaling pathway, HIF-I signaling pathway, Fc epsilon RI signaling pathway, chemokine signaling pathway, B-cell receptor signaling pathway, NK cell-mediated cytotoxicity, platelet activation, and vascular smooth muscle contraction (Figure 9D). The consistent results revealed a close relationship between MMD and the immune mechanism. The detailed results are illustrated in Figure 9.

Among the 13 immune-related hub genes, the expression of six genes showed significant differences between the MMD and the control groups. The results indicated that the expressions of BTK, FGR, SYK, CSF2RB, and CXCR1 were upregulated in MMD, whereas the expression levels of PTPN11 were downregulated (Figure 10A). To verify the potential of immune-related hub genes as diagnostic biomarkers, we further downloaded the dataset GSE157628 as the external validation dataset. The expression pattern of BTK, FGR, SYK, and PTPN11 showed consistent results (Figure 10B). The results indicated that BTK, FGR, SYK, and PTPN11 genes showed strong stability and that they have the potential to be new diagnostic biomarkers. The correlation/anti-correlation among the 13 differentially expressed genes identified in the discovery cohorts was explored (Figure 10C). A positive correlation was demonstrated between CSF2RB and CCR1(r = 0.72) and CXCR1 and CCR1(r = 0.74), whereas TGFB2 has a negative correlation with FGR (r = −0.48). Moreover, GSEA was used to further explore the pre-defined immune-associated pathways of the four immune biomarkers. Single-gene GSEA results showed that several kinds of immune cell pathways were mainly enriched in MMD, including NK cells, monocyte cells, DC cells, and neutrophils (Figure 11).

In addition, we constructed the diagnostic model by using the logistic regression algorithm. ROC analysis was performed, and the four immune-related genes were visualized by ROC curves. ROC analysis confirmed that BTK, FGR, SYK, and PTPN11 could distinguish MMD from control patients, and the accuracy of AUC ranged from 0.699 to 0.859 (Figure 12A). For single genes, PTPN11 had the highest accuracy, with an AUC of 0.80. The results also revealed that the diagnostic model constructed by combined genes had better diagnostic performance, and the highest AUC was 0.859. Furthermore, in the validation set, the diagnostic power was much higher both in single and combined gene diagnostic models, with an AUC between 0.808 and 0.970 (Figure 12B). The results of the ROC analysis showed strong stability of the diagnostic performance between the different datasets, and PTPN11 exhibited an excellent diagnostic capability. These identified four immune-related genes might function as novel diagnostic biomarkers and potential immunotherapeutic targets.

According to the CIBERSORT algorithm, the state of immune cell infiltration was reconstructed in MMD and the control group. Then, the immune infiltration score of 22 immune cell subtypes was obtained, and the proportions of immune cells in each sample are illustrated in Figure 13. The immune score results indicated that the fractions for eosinophils in the MMD group were remarkably higher than those of the control (p = 0.03). Compared with the control group, monocytes and neutrophils were upregulated in the MMD group, whereas CD8⁺ T cells, resting memory CD4⁺ T cells, and M0 macrophages were downregulated; however, they were not significant (Figure 14A). In addition, the correlation among the differentially infiltrated types of immune cells was analyzed, and several pairs of immune cell subtypes were found to have different degrees of positive or negative correlation (Figure 14B). This result suggests that M1 macrophages and resting DCs, M2 macrophages and activated mast cells, and M0 macrophages and resting NK cells showed the most synergistic effect. Together, resting NK cells and neutrophils, resting NK cells and native NK cells, and M0 macrophages and neutrophils represented the most competitive effect. Furthermore, the results of the in-depth relationship between the core genes and immune cell subtypes reported that BTK had a positive relationship with neutrophil infiltration and a negative relationship with resting NK cell infiltration, PTPN11 had a positive relationship with neutrophil infiltration and a negative relationship with resting CD8 T cell infiltration, FGR showed a weak relationship with all the immune cell subtypes, and SYK had a positive relationship with neutrophil infiltration and a negative relationship with T regulatory cell (Tregs) infiltration (Figure 14C). These results fully indicated that the expression of immune-related genes had a strong correlation with immune infiltration.