Combining machine learning with external validation to explore necroptosis and immune response in moyamoya disease

Abstract

Moyamoya disease (MMD) is a rare chronic vascular disease leads to cognitive impairment and stroke with its etiology unknown. The relationship between necroptosis or necroinflammation and MMD pathogenesis was poorly understood. Differentially expressed necroinflammation and necroptosis related genes (DE-NiNRGs) were selected based on the public gene expression data from Gene Expression Omnibus (GEO) and validated by our self-test data of MMD patients and control group. Functional enrichment analysis, PPI network and multi-factors regulation network construction of DE-NiNRGs were employed to discover the connections between these genes. DE-NiNRGs and immune cells correlation analysis provided evidence for the relationship between DE-NiNRGs and necroinflammation in MMD patients. We then established an MMD prediction model using support vector machine (SVM) and selected DE-NiNRGs as features. The DE-NiNRGs based MMD prediction model had excellent performance on test set with the area under the curve (AUC) higher than 0.9. Four genes, PTGER3, ANXA1, ID1, and IL1R1, that contributed significantly to the SVM model and passed the test of validation set are key genes in DE-NiNRGs. The upregulation of PTGER3 expression indicated that necroptosis and angiogenesis were promoted in MMD patients, whereas the downregulation of ANXA1 expression indicated that the migration and differentiation of immune cells are closely related to MMD pathogenesis. These findings provided new inspiration for our study of the immune-related pathogenesis and therapeutic targets of MMD.

Introduction

Moyamoya disease (MMD) is a rare cerebrovascular disease characterized by progressive stenosis or occlusion of the intracranial part of internal carotid artery and the smoke-like associated network of abnormally dilated collateral vessels on angiography [1]. MMD has been observed throughout the world and is now the most common pediatric cerebrovascular disease in East Asia [2, 3]. The typical complications include hemiparesis, dysarthria, aphasia, cognitive impairment and intracranial hemorrhage [4, 5], resulting in a high incidence of disability and even death. There are currently proteomics and genomics studies targeting MMD, but the pathogenesis of it remains unclear [6, 7]. Currently, surgical revascularization has been used as primary treatment by adjusting the blood flow to the affected cerebral hemisphere [8], but no treatment can now reverse the disease progress due to the lack of understanding of the detailed mechanism of MMD occurrence and development.

Necroptosis is a gene regulated inflammatory cell death pathways which has been a focus in recent years due to its significant contribution to the pathogenesis of many diseases [9, 10]. It utilizes a signal transduction pathway involves receptor interacting protein kinase (RIPK1/RIPK3), mixed lineage kinase domain-like protein (MLKL) and upon stimulation of death receptors (DR) and can potentially activate necroinflammation via the release of damage associated molecules (DAMPs) [11]. Necroptosis and necroinflammation are strongly associated with various pathologies of ischemic vascular diseases including stroke and myocardial infarction (MI) [12, 13]. MI involves irreversible cell death of cardiomyocytes due to coronary artery blockage. Increased level of RIPK3 expression was found in ischemic portions of MI mouse model, which indicates necroptosis. It was also found that necroptosis inhibitor, Nec-1 had an apparent cardioprotective effect for mice undergoing cardiac ischemia. Necroptosis played a similar role in pathology of stroke and blocking the expression of necroptosis related genes not only prevents cell death but also guides the inflammatory response into a neuroprotective type [14]. MMD is also a typical ischemic vascular disease and its pathological characteristics includes stroke. This suggests that necroptosis and necroinflammation is closely related to the mechanism of MMD.

Recently, the role played by inflammatory response and immune mediators in MMD has garnered increasing attention. Many articles had taken close observation on the relationship between immune cell infiltration and MMD. With the help of histopathological and immunohistochemical findings methods, aberrant expression of IgG and S100A4 protein was found in intracranial vascular wall of MMD patients, indicating that the thickened intima of MMD patients was related to immune cell infiltration [15]. Involvement of intrinsic immune reaction in MMD pathogenesis had also been suggested by the evidence that CD163 + M2-polarized macrophages mediated autoimmune led to tissue remodeling and angiogenesis [16].

Despite inflammation and immune response in MMD had been studied at cellular level and some important gene expression in immune cells had been characterized in the pathogenesis of MMD, there is not yet a report that specifically analyzed the role of necroinflammation and necroptosis related genes (NiNRGs) played in the disease.

In this study, we focused on the NiNRGs differentially expressed in MMD and non-MMD patients. Both public MMD datasets from Gene Expression Omnibus (GEO) and our self-test MMD dataset measured from patient samples were used in the analysis to make the result more robust and reliable, and enhance the depth of biological interpretation. Machine learning methods of random forest (RF), support vector machine (SVM), generalized linear model (GLM), and extreme gradient boosting (xgboost) were utilized to select the key NiNRGs related to MMD and create a prediction model for MMD based on the expression of key NiNRGs. Function enrichment analysis was implied to provide information on the underlying mechanisms of necroptosis and necroinflammation in MMD. And by identifying the key genes and the correlation between differential gene expression and immune infiltration, we mapped the roads to further validation and functional studies to understand the precise role of necroptosis and necroinflammation genes in MMD etiology and provided ideas for future drug development targeting necroptosis and necroinflammation.

Methods

Gene expression data and preprocessing

The discovery cohort (training set) used in this study comes from the raw data of GSE189993 in the GEO database with 21 disease samples (MMD) and 11 control samples (Control). Most of MMD samples included in the training set came from patients aged from 30 to 60 years old, and only 2 of them belonged to patients under 18 years old. Also, atherosclerotic disease was excluded from the dataset. The Combat function in SVA package was used to perform batch correction on this dataset, which utilizes empirical Bayesian methods to estimate the parameters of batch effects. By estimating parameters such as means and variances for batches based on batch information and other covariates in the data, the Combat function adjusts the data based on these parameters. For each gene, an adjustment factor based on the variance of the gene and the size of the batch effect was calculated in different batches, and then the expression value of the gene in each batch was multiplied by the corresponding adjustment factor to obtain the adjusted expression value. After the removal of batch effect, the dataset were used as our training set for subsequent analysis.

The data in our validation set came from self-test data. It included gene expression data from 13 patients, including 3 non-MMD patients and 10 MMD patients. STA samples were obtained from patients during surgery. The detailed information including clinical features of included participants were provided (Table 1). All participants in either MMD group or non-MMD group were adults aged, and atherosclerotic disease was totally excluded. Also, all MMD patients suffered from bilateral MMD. Total RNA of each STA sample was used as input material for the RNA sample preparations of RNA sequencing. Then mRNA was purified from total RNA by using poly-T oligo-attached magnetic beads. Fragmentation, cDNA synthesis, adaptor ligation, purification and PCR was carried out to construct a library for sequencing.

Table 1 Detailed clinical information for the participants of discover cohort

Necroptosis related genes (NRGs) and necroinflammation related genes (NiRGs) were selected by searching for “necroptosis” based on specific criteria (category protein coding and score > 0.8) in the GeneCrads database and in previous reports.

Differential expression analysis

Limma analysis was used to detect differentially expressed genes (DEGs) between individuals with MMD and healthy controls in our training set. The threshold for identifying significantly DEGs was set to P Value < 0.05 and|logFC|>1. Based on the results of differential expression analysis, volcano maps and heatmaps were generated to visualize the identified DEGs.

Correlation analysis of MMD with necroptosis and necroinflammation

Expression matrix of NRGs and NiRGs in the MMD samples of training set was extracted to perform Spearman correlation analysis. And a threshold of correlation coefficient| R|>0.5 and P value ≤ 0.001 were set to screen for necroinflammation and necroptosis related genes (NiNRGs). The intersection of differentially DEGs and NiNRGs was used to generate a Venn plot.

Functional enrichment analysis and construction of DE-NiNRGs PPI network

According to the differential expression analysis results between MMD and Control, the genes were sorted in descending order of logFC value. Then, clusterProfiler was used to perform KEGG pathway enrichment analysis on the genes, and adjusted P < 0.05 enrichment results were selected. GO and KEGG enrichment were performed on DE-NiNRGs, and pathways with P value < 0.05 were viewed as significantly enriched.

By combining the STRING database to predict and analyze whether there are interactions between DE-NiNRGs encoded proteins, only those interactions that have been empirically validated and have a total score higher than 0.4 are used to construct the PPI network.

DE-NiNRGs and immune cells correlation analysis

This study used the single-sample Gene Set Enrichment Analysis (ssGSEA) algorithm and the expression data of the training set to calculate the immune infiltration score of each sample. Box plots were drawn based on the immune cell scores of samples in different groups, with a significance threshold set at P value < 0.05.

According to the ssGSEA results, the correlation between 28 immune cell immune infiltration scores in disease samples was calculated using the Spearman method and displayed by a heatmap. The correlation and P value between DE-NiNRGs genes and various immune cells were calculated by using the Spearman method, with a significance threshold set at P < 0.05, and a heatmap was plotted for display.

Identification and characterization of necroinflammation and necroptosis related subtypes in MMD patients

In order to determine the different subtypes of necrotic apoptosis associated with necrotic inflammation, we performed unsupervised clustering on samples of smoke disease patients in the training set based on DE-NiNRG gene expression data. We set the number of clustering clusters k from 2 to 5 and used the “pam” clustering algorithm. The final number of subtype groups, was considered based on the clustering results heatmap and cumulative distribution map.

The immune infiltration characteristics of different MMD subtypes were evaluated by the abundance of immune cells in the sample calculated using ssGSEA and xcell. A boxplot was generated for display. The significantly enriched functions and pathways were selected based on FDR < 0.05 criteria using Gene Set Variation Analysis (GSVA).

The construction of multi-factors regulation network of DE-NiNRGs

The transcription factors that regulate DE-NiNRGs genes were analyzed by using the transcription factor database TRRUST database (https://www.grnpedia.org/trrust/). At the same time, through the miRWalk database (http://129.206.7.150/), we predicted the miRNA that regulate DE-NiNRGs, setting the screening criteria as follows: score value > 0.8, the binding region is located in 3UTR, and the miRNA appears in the TargetScan and miRDB databases. Then we constructed an mRNA miRNA network using Cytoscape software.

Prediction and analysis of small molecular drugs based on DE-NiNRGs

DGIdb (drug gene interaction database, https://digdb.org/) was used to predict the targeting relationship between DE-NiNRGs and drugs. Genes were selected with interaction scores > 1 and drug action pairs, and Cytoscape(version 3.9.1, https://cytoscape.org/) was used to build a gene drug interaction network.

Building MMD predictive models using machine learning

The expression values of DE-NiNRGs were used as explanatory variables to construct RF model, SVM model, GLM model, and xgboost model separately. Further analysis of these four models, and comprehensive consideration on their residual reverse cumulative distribution map, residual box plot, and ROC curve, led us to select the best model, which later appeared to be SVM.The SVM function of R package e1071 is mainly used to construct Support Vector Machine (SVM) models, optimize the model parameters using gene expression data, visualize the SVM model using the explanatory properties of DALEX package, create residual reverse cumulative distribution map, residual boxplot, and Receiver Operating Characteristic (ROC) curve to evaluate the performance of SVM models. After selecting the best model, we extracted the top 20 explanatory variables that contributed to the classification, and generated a nomogram by selecting genes that have been validated through the validation set.

Results

The 37 DE-NiNRGs of moyamoya disease

We used data from GSE189993 as our training dataset. We removed the batch effects to eliminate systematic bias between MMD data and HC data (Fig. 1B). We filtered differentially expressed genes(DEGs) based on a P value < 0.05 and|logFC|>1 threshold and got 5320 DEGs (Fig. 1C). Intersection of 5320 DEGs and 238 necroptosis and necroinflammation related genes (NiNRGs) was performed using Venn analysis, resulting in 37 DE-NiNRGs(Fig. 1D). Limma difference test showed there were 7 DE-NiNRGs that showed consistent trends in both training and validation sets and P < 0.05: PCSK9, TM6SF2, ID1, RB1, ANXA1, PTGER3, IL1R1.

Fig. 1

The extraction of DE-NiNRGs and functional enrichment analysis. A: STA sample acquisition method. MMD: moyamoya disease. STA: superficial temporal artery. MCA: middle cerebral artery. HC: healthy control, non-MMD patients. B: PCA plots of training set (GSE189993) data after batch correction. Blue: Control samples. Yellow: MMD samples. C: Volcano plot of differentially expressed genes (DEGs) in the MMD training dataset. Red: upregulated genes, Blue: downregulated genes. D: Venn map of DEGs and necroinflammation and necroptosis related genes (NiNRGs). Green: DEGs of GSE189993. Blue: NiNRGs. Intersection: DE-NiNRGs. E: Bubble chart of GO-BP enrichment result with top 10 significance level. BP: Biological Process. F: Bubble chart of GO-CC enrichment result with top 10 significance level. CC: Cellular Component

Functional enrichment analysis

GO-BP function analysis result showed that DE-NiNRGs mainly participated in the regulation of cytokine production pathway, inflammatory pathway and immune effector process (Figs. 1E and 2D). 2 DE-NiNRGs (PTGER3 /IL1R1) that past the limma difference test in both training and validation datasets were involved in the regulation of inflammatory pathway.

Fig. 2

Functional enrichment analysis. A: Bubble chart of GO-MF enrichment result with top 10 significance level. MF: Molecular Function. B: Bubble chart of KEGG pathway enrichment result with top 10 significance level. C: Histogram of KEGG pathway enrichment result with top 10 significance level. D: Connection map of GO-BP enrichment pathways and DEGs involved in these pathways. E: The GSEA enrichment results. It included the top 10 enrichment pathways

The gene expression products of DE-NiNRGs were enriched in cytoplasmic side of plasma membrane(Fig. 1F). According to the GO-MF analysis result, function of DE-NiNRGs expression products mainly linked to ubiquitin-protein ligase binding or ubiquitin-like protein ligase binding, which played an important role in modulating inflammatory signaling and immunity(Fig. 2A). KEGG analysis result also showed that DE-NiNRGs are enriched in necroptosis(Fig. 2B).

GSEA enrichment considered all differentially expressed genes of MMD but not just 37 DE-NiNRGs. We selected the top 10 pathways with high enrichment scores of GSEA enrichment result. The top 10 pathways include some immune-related pathways like autophagy and ubiquitin mediated proteolysis(Fig. 2E).

Immune infiltration related to DE-NiNRGs

Through examination of the difference between the immune infiltration level of MMD group and control group, we selected 6 kinds of immune cells that have distinct difference between groups: CD56 bright natural killer cell, immature B cell, macrophage, mast cell, type 17 T helper cell and type 2 T helper cell. Except type 2 T helper cell, all other 5 kinds of immune cells’ infiltration level are higher in MMD group (Fig. 3A).

Fig. 3

Immune infiltration analysis. A: Box plot of immune infiltration score. Blue: control group. Red: MMD group. Horizontal axis: 28 kinds of immune cells. Vertical axis: Immune infiltration score. ‘*’: P value between 0.05 and 0.01. ‘**’: P value < 0.01. B: Heatmap of correlation between DE-NiNRGs and immune cells in disease samples of training dataset. ‘ ’: P value > 0.05. ‘*’: P value between 0.05 and 0.01. ‘**’: P value < 0.01

The correlation between DE-NiNRGs genes and immune cells in disease samples was evaluated by ‘Spearman’ correlation (Fig. 3B). The interaction between DE-NiNRGs and immune cells were analyzed: ANXA1 significantly positively correlated with CD56 bright natural killer cell, mast cell and type 2 T helper cell, PTGER3 strongly positively correlated with type 17 T helper cell, and L1R1 strongly correlated with type 2 T helper cell. Besides, ID1 positively correlated with immature B cell and mast cell. In addition, CD47, FPR1, TLR2, S100A9, RIPK3 and CXCL5 showed strong positive correlation with immune cells.

The multi-factors regulation network

We selected miRNAs that regulate DE NiNRG based on the following criteria: score value > 0.8, the binding region is located in 3UTR, and the miRNA appears in the TargetScan and miRDB databases. The association between the DE-NiNRGs and miRNAs was visualized by the network (Fig. 4A). The network is mostly scattered than fully connected, suggesting that there wasn’t a large and intensely regulation network for DE-NiNRGs. In the contrary, the regulatory pathways of these genes are relatively independent.

Fig. 4

Multifactor regulatory network of DE-NiNRGs and identification of necroinflammation and necroptosis related subtypes in patients with MMD. A: mRNA-miRNA regulatory network of DE-NiNRGs. Pink diamonds: DE-NiNRGs. Green polygons: miRNAs. B: Unsupervised clustering result of MMD samples in training set, k = 2. Darker color in the consensus matrix represented higher possibility that two samples belonged to the same cluster. C: CDF (Cumulative Distribution Function) curves of different k value in clustering. D: PCA plot of 2 MMD subtypes. Red and green dots represented samples belonged to the different clusters. E: Curve of the change of area under CDF curves. F: Tracking plot of unsupervised clustering of MMD samples with different k values. The horizontal axis represents different samples, and different colors represent the different clusters to which each sample belonged under different k value on the vertical axis. G: Box plot of immune infiltration score. Blue: MMD cluster (1) Red: MMD cluster (2) Horizontal axis: 28 kinds of immune cells. Vertical axis: Immune infiltration score. ‘*’: P value between 0.05 and 0.01. ‘**’: P value < 0.01

The prediction of small molecular drugs

DGIdb database was used to predict the targeted relationship between DE-NiNRGs genes and drugs, and a gene drug interaction network diagram was constructed by using Cytoscape(Fig. 5C). None of these potential small molecular drugs was showed to have ability to interact with more than 2 target genes, and most of these small molecular drugs targeted only one gene. Among the 7 hub DE-NiNRGs, PCSK9, PTGER3, IL1R1 and ANXA1 were found to have multiple potential small molecular drugs that targeted them.

Fig. 5

The ROC curves and prediction and analysis of small molecule drugs. A: ROC curve for SVM on validation dataset (test dataset). B: ROC curve for SVM on training dataset. C: Gene-drug interaction network. Pink diamonds: genes among DE-NiNRGs. Blue triangles: small molecule drugs

The identification and characterization of 2 MMD subtypes

Based on the DE-NiNRGs expression dataset, unsupervised clustering analysis was performed on all MMD samples in training dataset. PAC was used to search the optimal number of clusters, which was 2 in this study (Fig. 4B, C, E and F). According to the clustering results, it can be seen that when K = 2, the heatmap is most clearly divided without too many cluttered lines(Fig. 4B). At the same time, the cumulative contribution of CDF is closest to parallel to the X-axis(Fig. 4C). Therefore, the optimal number of clusters is set to 2, and MMD patient samples are divided into two subtypes: cluster1 contains 28 samples and cluster2 contains 4 samples(Fig. 4B). According to the PCA result, cluster1 and cluster2 had distinct distribution(Fig. 4D).

2 MMD subtypes had different characteristics of immune infiltration(Fig. 4G). The levels of immune cell infiltration of basophils, immature dendritic cells(iDC), mesenchymal stem cells(MSC), plasmacytoid dendritic cells(pDC) in cluster2 were significantly higher than in cluster1.

The construction of SVM model for MMD prediction

To select the optimal MMD prediction model, four prediction models were constructed using the expression data of DE NiNRGs, including Random Forest Model (RF), Support Vector Machine Model (SVM), Generalized Linear Model (GLM), and Extreme Gradient Boosting Model (xgboost). According to the residual reverse cumulative distribution map, residual box plot, and receiver operating characteristic (ROC) curve, the SVM model has the smallest residual and an AUC value greater than 0.8 (Fig. 6A, B and C). Considering the SVM model as the best model, the top 20 explanatory variables with contribution values in the SVM model are selected for subsequent analysis (Fig. 6D).

Fig. 6

MMD prediction model using machine learning methods. A: Box plot distribution of the absolute residual of different machine learning models. Horizontal axis: the absolute value of residual. Vertical axis: machine learning models. B: Reverse cumulative distribution of the absolute residual of machine learning methods. C: Receiver Operating Characteristic Curve (ROC) for different machine learning models. D: Feature importance of SVM model. 52 lines correspond to 52 genes in DE-NiNRGs. E: Nomogram for evaluate the degree of influence of 4 crucial genes on MMD prediction. F: Calibration curve of column chart model. G: The DCA curve

Four genes, ID1, ANXA1, PTGER3, and IL1R1, among the top 20 explanatory variables in the SVM model, passed the validation set validation. The feature importance of these four genes, from high to low, is as follows: PTGER3, ANXA1, ID1, IL1R1 (Fig. 6E). Calibration curves, decision curves, and ROC curves were also plotted to evaluate the model’s performance (Fig. 6F and G). According to the results, it can be concluded that the AUC values in both the training and validation sets are greater than 0.8, which indirectly indicates that the predictive performance of the column chart is good.

Discussion

MMD is a progressive cerebrovascular disease with its etiology remained unclear, and yet no effective drugs has now been developed to prevent or reverse the progression of MMD. Currently, very few of studies have focused on the role necroptosis and necroinflammation played in MMD. In this study, we utilized various machine learning methods including SVM, RF and GLM to identify the most important necroptosis-related and necroinflammation-associated genes in MMD, and analyzed the relationship between hub genes and various immune cells through immune infiltration analysis. The abnormal expression of PTGER3 and ANXA1 in MMD patients has been found to be closely related to the promotion of necroptosis, migration and differentiation of immune cells, and angiogenesis. These genes may serve as a breakthrough for further research on the pathogenesis of MMD.

Previous studies on the pathophysiological features of MMD patients suggested caspase-dependent apoptosis as a contributory mechanism in the associated degradation of the arterial wall [17]. Meanwhile, infiltration of macrophages and T cells were found in the occlusive blood vessels in MMD patients [18]. Research on cell death and inflammation mediators in MMD have been done to explicit the mechanism of MMD progression, but none of them can fully explain the immune-related mechanisms underlying MMD. Necroptosis and necroinflammation were discovered to be closely related to stroke and cardiovascular disease and may be an important mechanism of MMD to be discovered.

In this study, we selected 52 DE-NiNRGs and 4 most important DE-NiNRGs in machine learning MMD prediction model. Prostaglandin E2 receptor 3 (PTGER3 or EP3) belongs to the receptor of prostaglandin E2 (PGE2) family which is a subtype of G-protein-coupled receptors [19]. PGE2 stimulation on EP3 can directly decrease intracellular cAMP concentration and rise Ca²⁺. The activation of EP3 receptors by PGE2 can influence necroptosis. By promoting necroptosis through EP3 activation, PGE2 may serve as a mediator that leads cell death towards necroptosis when apoptosis is inhibited. The upregulation of EP3 expression we found in MMD patients indicates that necroptosis was involved and may play a significant role in MMD pathogenesis.

The upregulation of PTGER3 expression is also involved in various vascular diseases. The EP3 receptor can regulate the contraction and relaxation of blood vessels by affecting the function of smooth muscle cells. Its activation usually leads to vasoconstriction, affecting local blood flow and blood pressure, thus playing a certain role in cardiovascular disease and hypertension [20]. Also, EP3 can facilitate the cerebral small artery remodeling, which leads blockade of the EP3 receptor a promising treatment option of cerebral small vessel disease (CSVD) [21]. In this study, PTGER3 was upregulated in the disease group compared to the control group. This upregulation suggests that PTGER3 promote the pathogenesis connected with internal carotid artery contractile and vascular network generation in MMD.

Annexin A1 (ANXA1) also has great importance in the immune filtration of MMD. ANXA1 is part of the annexin family of proteins, which are characterized by their ability to bind phospholipids in a calcium-dependent manner. ANXA1 can modulate the migration and activation of immune cells, promoting the resolution of inflammation in brain microenvironment [22, 23]. Tumor-derived ANXA1 is related to recruitment and polarization of bone marrow-derived macrophages in order to suppress the immunoenvironment [24]. Also, ANXA1-FPR1 pathway was closely related to microvascular invasion (MVI) in tumor development, specifically related to MVI-associated intercellular communication [25]. In this study, ANXA1 was down-regulated in MMD samples and it significantly positively correlated with CD56 bright natural killer cell, mast cell and type 2 T helper cell. This suggests that in MMD patients, the suppression of ANXA1 expression leads to an enhanced immune response in the cerebrovascular system and further cell damage and apoptosis which potentially exacerbating the disease.

Inhibitor of DNA binding 1(ID1) encodes a helix-loop-helix (HLH) protein that acts as a transcription factor that participated in multiple biological processes. It plays an essential role in T helper type 9 cell (Th9) and T regulatory cell (Treg) differentiation [26, 27]. In this study, ID1 was upregulated in MMD group and is commonly positively correlated with immune cells. This may be a hint on ID1 participated inflammation and T cells differentiation process in MMD development. Interleukin-1 receptor type 1 (IL-1R1) is a specific receptor for interleukin-1, and activates various intracellular pathways such as the MAPK and NF-κB pathways like ID1. IL-1R1 can directly interact with RIP1/RIP3, leading to the formation of necrosome complex and then initiate hemin-induced neuronal necroptosis [28]. In this study, IL1R1 gene was down regulated in MMD samples which could either be a signal of activation or repression of necroptosis and that need further validation.

Most of DE-NiNRGs are enriched in pathways for regulation of inflammatory or necroptosis. The enrichment result was not specific enough but confirmed that DE-NiNRGs we selected did significantly related to necropotosis and necroinflammation.

Immune-mediated inflammation of the vasculature is a common cause of vascular narrowing and occlusion in vascular diseases without atherosclerosis including MMD [29, 30]. As we discussed above, the four major DE-NiNRGs all had a tendency to increase the immune infiltration level in general and even lead to necroptosis. Necroptosis was commonly related to cerebral ischemia disease [12] and it has been proven that ischemic brain injury can activate necroptosis pathways [31]. MMD is also occlusive vascular disease and involves ischemic brain injury but lack of research on its necroptosis and necroinflammation condition. The correlation analysis we did between DE-NiNRGs and immune infiltration revealed the most important genes and 6 kinds of crucial immune cells which participated. However, immune infiltration based on transcriptomic data has its inherent defect, including the disability to capture this spatial heterogeneity due to the bulk nature of the samples, potentially leading to misestimations of immune cell populations. Because of that, further research like spatial transcriptomics analysis was necessary to validate the significance of these immune cells in the development of MMD.

We also dedicated to identify and characterize 2 subtypes of MMD according to the differentially differential gene expression and immune infiltration. Previous study has used systemic immune-inflammatory markers to predict MMD and its subtypes and selected NLR and SII level as criteria [32]. This has significant clinical implications, but the genetic differences behind different subtypes of MMD have not been elucidated. In this study, 2 MMD subtypes were identified based on DE-NiNRGs expression and they showed great immune infiltration difference of basophils, iDC, pDC and MSC. This may indicate distinct inflammation activation pathways in subtypes related to the differential expression DE-NiNRGs.

Limitation

Though we used self-test data for validation, the limited size of samples still made it difficult to further interpret the experimental results. The sample size included in this study was limited due to the difficulties in collecting samples from MMD patients. 32 samples data in GSE189993 and 20 samples data in GSE157628 from GEO database were used in this study as discovery cohort. All these samples in discovery cohort came from the same organization and single ethnic group, which may limit the representativeness of the study. More samples from multiple ethnic groups were needed to further validate our discovery.

Besides, the lack of animal models restricted us to validate the significance of DE-NiNRGs and the related pathways in MMD development in vitro and weakens the level of evidence to draw conclusion. Future studies to validate our conclusion in vivo or on MMD organoid models were needed.

Conclusion

In this study, we obtained MMD data from Gene Expression Omnibus (GEO) dataset as training set and our self-test data as validation set. Differential analysis and correlation analysis identified 8104 MMD differentially expressed genes (DEGs) and 52 DE-NiNRGs. Immune correlation analysis identified 7 immune cells whose immune infiltration level differ significantly between MMD and control group, and the relationship between they and DE-NiNRGs was analyzed. Besides, we identified 2 MMD subtypes based on DE-NiNRGs expression and characterized them by their immune infiltration level. Furthermore, we built a MMD prediction model and selected 4 significant features which are also hub genes among DE-NiNRGs. This provides us with a new perspective on understanding MMD from the standpoint of necroptosis and necroinflammation related genes.