Identification of oxidative phosphorylation-related genes in moyamoya disease by combining bulk RNA-sequencing analysis and machine learning

Introduction: Moyamoya disease (MMD) is a chronic cerebrovascular disease that can lead to ischemia and hemorrhagic stroke. The relationship between oxidative phosphorylation (OXPHOS) and MMD pathogenesis remains unknown. Methods: The gene expression data of 60 participants were acquired from three Gene Expression Omnibus (GEO) datasets, including 36 and 24 in the MMD and control groups. Differentially expressed genes (DEGs) between MMD patients MMD and control groups were identified. Machine learning was used to select the key OXPHOS-related genes associated with MMD from the intersection of DEGs and OXPHOS-related gene sets. Gene ontology (GO), Kyoto encyclopedia of genes and genomes (KEGG), gene set enrichment analysis (GSEA), Immune infiltration and microenvironments analysis were used to analyze the function of key genes. Machine learning selected four key OXPHOS-related genes associated with MMD: CSK, NARS2, PTPN6 and SMAD2 (PTPN6 was upregulated and the other three were downregulated). Results: Functional enrichment analysis showed that these genes were mainly enriched in the Notch signaling pathway, GAP junction, and RNA degradation, which are related to several biological processes, including angiogenesis, proliferation of vascular smooth muscle and endothelial cells, and cytoskeleton regulation. Immune analysis revealed immune infiltration and microenvironment in these MMD samples and their relationships with four key OXPHOS-related genes. APC co-inhibition (p = 0.032), HLA (p = 0.001), MHC I (p = 0.013), T cellco- inhibition (p = 0.032) and Type I IFN responses (p < 0.001) were significantly higher in the MMD groups than those in the control groups. The CSK positively correlated with APC co-inhibition and T cell-co-inhibition. The NARS2 negatively correlated with Type I IFN response. The SMAD2 negatively correlated with APC co-inhibition and Type I IFN response. The PTPN6 positively correlated with HLA, MHC I and Type I IFN responses. Discussion: This study provides a comprehensive understanding of the role of OXPHOS in MMD and will contribute to the development of new treatment methods and exploration of MMD pathogenesis.


Introduction
Moyamoya Disease (MMD) is a cerebrovascular disease characterized by chronic progressive occlusion of the distal end of the internal carotid artery or the main branches of the Willis circle.This occlusion causes the compensatory proliferation of small blood vessels at the brain base, presented as a 'smoke-like' vascular network on digital subtraction angiography (DSA) (Gonzalez et al., 2023).MMD may have catastrophic consequences such as hemorrhagic or ischemic stroke.In addition, it can cause serious symptoms such as motor dysfunction, sensory dysfunction, cognitive disorders, and epilepsy (Liu et al., 2023;Fang et al., 2024).However, its pathogenesis remains unclear.Currently, intracranial revascularization surgery is the primary treatment for MMD (Ihara et al., 2022).Although several targets related to MMD pathogenesis, such as RNF213 and DIAPH1, have been discovered recently, corresponding drugs have not been developed based on these targets (He et al., 2023a;Fang et al., 2024).Few effective drugs for treating MMD are available because of the unknown molecular mechanisms underlying its occurrence (Ihara et al., 2022).Therefore, the mechanism of MMD should be urgently clarified to improve treatments.
Oxidative phosphorylation (OXPHOS), a crucial biochemical process in eukaryotic cells, is the final stage of aerobic respiration, followed by glycolysis and citric acid cycle (Wilson, 2017).It reduces oxygen to generate high-energy phosphate bonds in the form of adenosine triphosphate (ATP) and is the most efficient stage in the electron transport chain to generate ATP (Nolfi-Donegan et al., 2020).The ATP generated through OXPHOS is the primary energy source for cells and plays a crucial role in various physiological and pathological cellular processes (Lei et al., 2017).Recently, numerous studies have investigated the role of OXPHOS in vascular diseases.A previous study has reported that CHCHD4 orchestrates mitochondrial OXPHOS and antagonizes the aberrant growth and migration of pulmonary artery smooth muscle cells in pulmonary arterial hypertension (PAH) (Wang et al., 2023).The glutamine antagonist 6-diazo-5-oxo-L-norleucine (DON) inhibits vascular smooth muscle cell (VSMC) proliferation and migration by attenuating OXPHOS in atherosclerosis (Park et al., 2021).OXPHOS is essential for purinergic receptor-mediated angiogenic responses in vascular endothelial cells (ECs) (Lapel et al., 2017).Pathological mutations in the OXPHOS system subunits can increase vascular endothelial growth factor (VEGF) production and pathological angiogenesis (Bayona-Bafaluy et al., 2019).Thus, OXPHOS plays an important role in vascular disease pathogenesis.
Previous studies have demonstrated that OXPHOS is related to abnormal EC, VSMC, and angiogenesis, which might play important roles in vascular disease pathogenesis.As a cerebral vascular disease, the pathological characteristics of MMD also includes abnormal EC proliferation, abnormal VSMC migration, pathological angiogenesis, and vascular remodeling (Takagi et al., 2016;He et al., 2023b;He et al., 2023c).This generates evidence that OXPHOS disorder or regulation is closely related to MMD mechanism.However, studies on the relationship between MMD and OXPHOS remain scarce.Although a previous study has mentioned their relationship (Kanamori et al., 2021), further exploration is needed to study the role of OXPHOS in MMD, such as signaling pathways, related genes, and the immune response.Therefore, investigating the association between OXPHOS and MMD is urgently needed.
In this study, we collected MMD-related gene expression datasets and OXPHOS-related gene sets from Gene Expression Omnibus (GEO) and Genecards database.Machine learning algorithms of Least absolute shrinkage and selection operator (LASSO) and support vector machine-recursive feature elimination (SVM-RFE) were utilized to selected the key OXPHOS-genes related to MMD.We finally identified 4 key OXPHOS-related genes, CSK, NARS2, PTPN6, and SMAD2, and analyzed their function by KEGG and GO enrichment analysis, which showed that these genes were most relevant to MMD pathogenesis.GSEA analysis were also performed to provide the profile of immune infiltration and microenvironment about these genes.Ultimately, we identified the associations between these 4 key OXPHOS-related genes and the pathogenesis of MMD.

Data acquisition
Gene expression data were retrieved from RNA sequencing data from the GEO database (https://www.ncbi.nlm.nih.gov/geo/info/datasets.html).The series matrix data files GSE189993, GSE157628, and GSE141024 enrolled 32, 20, and 8 patients, respectively, including 11, 9, and 4 in the control group, respectively, and 21, 9, and 4 in the MMD group, respectively.All the series matrix files were annotated using the GPL16699 file.OXPHOS-related gene sets were downloaded in the GeneCards database (https://www.genecards.org/).Those OXPHOS-related genes include protein-coding genes, functional element genes and RNA genes.The relevance score of GeneCards is a comprehensive score that measures the correlation between genes and specific keywords or topics based on multiple data sources and algorithms.Choosing a higher reliability score threshold can reduce the interference of noise and unrelated genes, enhance the reliability of analysis results, and ensure that the number of genes is not too large.Referring to several previous studies about screening genes with specific functions in GeneCards dataset (Li et al., 2022;Liu et al., 2022), we set the threshold of relevance score >15 to ensure that the selected genes have sufficient correlation and potential biological significance of oxidative phosphorylation.

Gene differential expression analysis
After merging the three datasets, the correction of the microarray data was performed using a surrogate variable analysis (SVA) package to reduce batch effects among them.Subsequently, we used principal component analysis (PCA) to draw a sample distribution map before and after correction, visually demonstrating the changes in batch effects before and after correction.The R package 'limma' was utilized to identify the differentially expressed genes (DEGs) between MMD and control samples by analyzing the gene expression data from the GEO cohort.Due to the small sample size in this study, p-value <0.05 was set as the criteria for filtering the DEGs to reduce the omissions of some potential biologically significant genes, rather using padj or FDR value.Meanwhile, the |log2 Fold Change (FC)| >1 was also set as the criteria for DEGs identification combined with p < 0.05 to ensure that the screened genes have significant differences in expression levels.Therefore, the final criteria of filtering DEGs was set as p < 0.05 and |log2 FC| >1, which was also be utilized in previous studies (Liu et al., 2020;Wang et al., 2022).

GO and KEGG functional enrichment analysis
GO and KEGG enrichment analyses were performed to evaluating the biological function, biological process, molecular function, and cellular component of the intersection genes between DEGs and OXPHOS-related gene set.When performing the enrichment analysis, we utilized the R package 'clusterProfiler' to perform the GO and KEGG enrichment analysis.The default threshold for 'clusterProfiler' package to do enrichment analysis is a p-value <0.05 and a Q value <0.2.In this study, to make the enrichment results more accurate, we further limited the Q value as less than 0.05, which is a relative strict standard.Finally, we set the cutoffs of GO and KEGG enrichment analysis as a p-value <0.05 and a Q value <0.05, referring to several previous studies (Liu et al., 2023;Li et al., 2024).

Key genes and possible disease diagnostic indicator selection
LASSO logistic regression and SVM-RFE algorithm were used to select the potential disease diagnostic markers.LASSO logistic regression was conducted using 'glmnet' package.When performing LASSO logistic regression, we used the 'glmnet' package to fit gene expression levels based on the glmnet function, with the family set as gaussian distribution, and automatically selects the optimal lambda value.By using the 'cv.glmnet' function for 10 fold cross validation, the optimal regularization parameter lambda was selected, and the performance curves of the model were plotted for each lambda value during the cross validation process.Finally, we selected the optimal lambda value by observing the lowest point of the curve.SVM-RFE machine learning model was constructed using 'e1071' package to further identify the significance of indicators for disease diagnosis.When performing SVM-RFE, we first used the SVM-RFE method to perform feature selection on DEGs to obtain the ranking of each gene.SVM model was established using the "e1071"R package, and the SVM model is constructed based on the top 10 genes to evaluate the lowest error rate of the model.The gene sets with the error rate reaching the lowest point and the genes obtained by the LASSO were intersected to obtain the key genes.

Gene set enrichment pathway analysis
Gene set enrichment analysis (GSEA) was used to investigate the DEGs between the high-and low-expression groups based on the expression profiles of patients with MvMD.GSEA was performed using a predefined gene set 'c2. cp.Kegg' to rank them according to their differential expression levels in the two types of samples.The preset gene sets were examined to determine whether they were enriched at the top or bottom of the ranking table.This study used GSEA software (http://www.broadinstitute.org/gsea) to compare the differences in KEGG signaling pathways between high expression and low expression groups, and explored the molecular mechanisms of key genes in the two groups of patients.The times of substitutions was set as 1,000 and the type of substitution was set as 'phenotype.' Pathways with a p-value <0.05 were considered significantly enriched.

Immune cell infiltration analysis
Single-sample GSEA (ssGSEA) was performed to evaluate the immune cells types in the microenvironment using 'GSVA' R package.We performed ssGSEA to quantify the immune cells in the expression profile and calculated the relative proportions of 19 types of infiltrating immune cells, including activated Dendritic Cells (aDCs), Antigen-Presenting Cells co-inhibition (APC_co_ inhibition), Antigen-Presenting Cells co-stimulation (APC_co_ stimulation), B lymphocytes (B_cells), C-C Chemokine Receptor (CCR), Checkpoint, Human Leukocyte Antigen (HLA), Inflammation-promoting, Macrophages, Major Histocompatibility Complex class I (MHC_class_I), Neutrophils, Parainflammation, plasmacytoid Dendritic Cells (pDCs), T cell co-inhibition, T follicular helper cells (Tfh), Tumor-Infiltrating Lymphocytes (TIL), Regulatory T cells (Treg), Type I Interferon Response (Type_I_IFN_Reponse) and Type II Interferon Response (Type_ II_IFN_Reponse).Spearman's correlation analysis was used to analyze the gene expression and immune cell content.p < 0.05 was considered statistically significant.

Statistical analysis
All statistical analyses were performed using the R software (version 4.2.2) and the figures were plotted using the 'ggplot' R package (version 3.3.6).Statistical significance was set at p < 0.05.

DEG identification and functional enrichment analysis in the MMD cohort
The MMD-related datasets GSE189993, GSE157628, and GSE141024 were downloaded from the GEO database and used in this study, including 60 samples from patients with MMD (n = 36) and control participants (n = 24).The detailed clinical information for the participants was listed in Supplementary Table S1 in Supplementary Material.Correction was performed using the SVA algorithm to reduce the batch effect between microarrays.The PCA plot shows that this process reduced the batch effect (Figure 1A).
A total of 1,548 genes were differentially expressed between patients with MMD and control participants (850 were upregulated and 698 were downregulated in the MMD group), according to p < 0.05 and |logFC|>1.The volcano plot displays an overview of the DEGs between patients with MMD and control participants (Figure 1B).The detailed information of these 1,548 DEGs were listed in Supplementary Table S2 in Supplementary Material.A total of 389 OXPHOS-related genes with >15 relevance score was screened from OXPHOS-associated gene sets in the GeneCards database (https://www.genecards.org/).The detailed information of these 389 OXPHOS-related genes were listed in Supplementary Table S3 in Supplementary Material.Finally, 27 intersecting genes were selected from the screened OXPHOSassociated genes and genes that were differentially expressed between the patients with MMD and control participants (Figure 1C).The 27 intersecting genes were also listed in Supplementary Table S3 in Supplementary Material.
GO biological process enrichment analysis revealed that these 27 genes were mainly enriched in peptidyl tyrosine phosphorylation, peptidyl tyrosine modification, and gland development.GO molecular function enrichment analysis revealed that these genes were mainly enriched in protein tyrosine kinase activity, phosphotyrosine residue binding and phosphatase binding (Figure 1D).The GO cellular component enrichment was not significant.KEGG enrichment analysis showed that these genes were mainly enriched in the apelin signaling pathway, PI3K-Akt signaling pathway, and endocrine resistance (Figure 1E).

Identification of key phosphorylationassociated DEGs by machine learning
To further identify the key genes affecting MMD, LASSO logistic regression and SVM algorithms were performed based on the previously obtained intersecting genes.The results identified 13 genes as MMD characteristic genes using LASSO regression (Figures 2A,B).We also used the SVM algorithm to identify characteristic genes for MMD and screened the top seven genes with the lowest error rates (Figure 2C).The genes screened by LASSO regression and the SVM algorithm were then intersected, and four genes were selected, including three downregulated genes, CSK, NARS2, and PTPN6, and one upregulated gene, SMAD2 (Figure 2D).These four genes were considered key genes in MMD and further investigated.

Functional enrichment analysis for the key genes in MMD
Next, we investigated the signaling pathways enriched in these four genes and explored the potential molecular mechanisms affecting MMD progression.GSEA results showed that CSK was enriched in the adipocytokine signaling pathway, neurotrophin signaling pathway, proteasome, neuroactive ligand-receptor interaction, olfactory transduction, and RNA degradation (Figure 2E).NARS2 was enriched in aminoacyl-tRNA biosynthesis, inositol phosphate metabolism, neuroactive ligand-receptor interaction, olfactory transduction, PPAR signaling pathway, and RNA degradation (Figure 2F).PTPN6 was enriched in arachidonic acid metabolism, the B cell receptor signaling pathway, glycosylphosphatidylinositol (GPI) anchor biosynthesis, neuroactive ligand-receptor interaction, primary immunodeficiency, and proteasomes (Figure 2G).SMAD2 was enriched in the cytosolic NDA sensing pathway, GAP junction, glycerolipid metabolism, Notch signaling pathway, pentose phosphate pathway, and toll-like receptor signaling pathway (Figure 2H).

Analysis of immune infiltration and microenvironment
Next, the correlation between the four key genes and immune cell infiltration was analyzed.The heatmap shows the proportion of 19 immune cell types in each sample (Figure 3A).A thermodynamic chart displays the relationship with each immune cell type (Figure 3B).The results showed that APC co-inhibition (p = 0.032), CCR (p = 0.019), checkpoint (p = 0.001), HLA (p = 0.001), MHC I (p = 0.013), T cell-co-inhibition (p = 0.032), TFH cells (p < 0.001), TIL cells (p = 0.001), Treg cells (p = 0.012), and Type I IFN responses (p < 0.001) were significantly higher in the MMD groups than those in the control groups (Figure 3D).
Next, we analyzed the relationships between the four key genes and immune cells.The results showed that CSK significantly positively correlated with T cell co-inhibition and APC coinhibition, and significantly negatively correlated with B cells and inflammation promotion (Figure 3C).NARS2 significantly negatively correlated with pDCs and the Type I IFN response (Figure 3C).SMAD2 positively correlated with macrophages, but significantly negatively correlated with TFH cells, APC coinhibition, pDCs, and the Type I IFN response (Supplementary Figure S1F).PTPN6 positively correlated with TIL, CCR, macrophages, HLA, neutrophil MHC I, and inflammation (Supplementary Figure S1F).Bubble plots showed the relationship between the four key genes and the immune microenvironment, including chemokine, immunoinhibitor, receptor, immunostimulatory, and MHC (Supplementary Figures S1A-E).

In depth analysis of the correlation between the key genes and MMD
To further explore the potential molecular mechanisms of the four key genes in MMD, we analyzed other related genes.We first analyzed differences in the expression of MMD-related genes obtained from the GeneCards database between the MMD and control groups.The results showed that APOE, DNM2, IL10, TGFβ1, and VWF expression levels in MMD group were significantly higher than that in control group (Figure 4A).Next, we analyzed the correlation between MMD-related genes and these four key genes.Significant correlations were observed among these variables.Among these, PTPN6 was significantly and positively correlated with DNM2, whereas CSK was significantly and negatively correlated with APOE (Figure 4B).

Discussion
A considerable proportion of patients with MMD may suffer from serious cerebral hemorrhage or ischemic stroke; thus, timely intervention is required.Currently, intracranial revascularization surgery is the main treatment for treating MMD (Ihara et al., 2022).Due to possible postoperative complications, such as stroke or infection, surgery cannot achieve a satisfactory prognosis in some patients with MMD (Kuroda and Houkin, 2008;Liu et al., 2023).Therefore, novel treatment methods should be urgently developed.However, owing to the unclear mechanism of MMD, there is a lack of effective drugs to prevent its progression, as well as effective animal models (Rallo et al., 2021;Ihara et al., 2022).Although some biological targets have recently been identified by genomics, proteomics, or other methods, animal modeling and drug research based on these targets have achieved unsatisfactory results (Rallo et al., 2021;Ihara et al., 2022).This suggests that other mechanisms may be involved in MMD development.Several studies have reported the role of OXPHOS in vascular diseases (Lapel et al., 2017;Bayona-Bafaluy et al., 2019;Park et al., 2021;FIGURE 4 Analysis of the correlation between the key genes and MMD.(A) Box plot for the expression levels of DEGs for MMD related to the key phosphorylation-associated genes.(B) Bubble plot for the relationship between DEGs for MMD and the key phosphorylation-associated genes.MMD = Moyamoya disease; DEGs = Differentially expressed genes.
Frontiers in Genetics frontiersin.org08 Wang et al., 2023), indicating that OXPHOS may also be involved in MMD pathogenesis.However, few studies have examined the relationship between OXPHOS and MMD.Therefore, in this study, we systematically investigated the role of OXPHOS-related genes, signaling pathways, and immune responses in the pathogenesis of MMD using bioinformatics.Four key genes were selected: CSK, NARS2, PTPN6, and SMAD2.
SMAD2 is a signal transducer and transcriptional modulator that mediates multiple signaling pathways (Nomura and Li, 1998).This protein mediates transforming growth factor beta (TGFβ) signal, and thus regulates multiple cellular processes, such as cell proliferation, apoptosis, and differentiation (Kamato and Little, 2020).In our study, we found a significantly downregulation of SMAD2 expression.As for previous studies of SMAD2, Petersen et al. found that SMAD2 knockdown developed an aggressive phenotype in breast cancer cells and enhanced angiogenesis (Petersen et al., 2010).Cao et al. found that activin receptor-like kinase 7 promotes VSMC apoptosis by activating the Smad2/ 3 signaling in diabetic atherosclerosis (Cao et al., 2022).These studies provided evidence that SMAD2 downregulation in patients with MMD might promote angiogenesis, proliferation and phenotype switching of VSMCs, which are the important characteristics of MMD mentioned before (Ma et al., 2022;He et al., 2023c;He et al., 2023d).As for the results in our study, pathway enrichment analysis also showed that SMAD2 was enriched in the Notch, GAP junction, and toll-like receptor signaling pathways, which are also associated with angiogenesis, vascular development, endothelial cell and VSMC growth in previous studies (Beyer and Berthoud, 2018;Shafeghat et al., 2022;Cuervo et al., 2023).In summary, these results generated the evidence that SMAD2 may play an important role in MMD pathogenesis and vascular changes.
CSK is a negative regulator of Src family kinases (SFKs) recruited to the plasma membrane by binding to transmembrane or adapter proteins to suppress signaling via various surface receptors by phosphorylating and maintaining several inactive positive effectors, such as FYN or LCK (D and CJB, 1999;Sun and Ayrapetov, 2023).It plays important roles in cell growth, migration, differentiation, and immune response (Zhu et al., 2023).There were several studies about the CSK on endothelial cells.Alejandro et al. found that CSK activation renders the endothelial barrier more susceptible to TNF-α through p38 activation, which leads to endothelial tight junction loss (Adam et al., 2016).Ulf et al. reported that CSK is involved in endothelial cell proliferation inhibition (Baumeister et al., 2005).Previous studies showed that the CSK has an inhibitory effect on the proliferation of endothelial cells.In this study, we found a downregulated of CSK expression in the MMD group.Accordingly, we inferred that may promote the angiogenesis through stimulating the endothelial cell proliferation, which may cause a intimal thickening of MMD vessels reported in previous study (Takagi et al., 2016;He et al., 2023b).
NARS2 encodes asparaginyl-tRNA synthetase 2, a putative member of the class II family of aminoacyl-tRNA synthetases.This enzyme catalyzes asparagine attachment to mitochondrial tRNA and plays a critical role in mitochondrial protein biosynthesis (Duchene et al., 2009).NARS2 mutations are associated with a combined oxidative phosphorylation deficiency 24 and autosomal recessive deafness-94 (Sissler et al., 2017).In this study, NARS2 was downregulated but there is a lack of research about the relationship between NARS2 and MMD.Due to its association with cellular energy metabolism, we inferred that NARS2 downregulation in MMD may lead to a deficiency in mitochondrial protein synthesis and affect the energy supply within the cell, leading to VSMC and EC dysfunction in MMD vessels.
PTPN6 is a member of the protein tyrosine phosphatase (PTP) family, a group of enzymes that remove phosphate groups from tyrosine residues by hydrolyzing phosphoric acid monoesters.PTPs are signaling molecules that regulate several cellular processes, including cell growth, differentiation, mitotic cycle, and oncogenic transformation (Amin et al., 2021).In this study, PTPN6 was the only upregulated among the four key genes.In previous studies, Clement et al. found that PTPN6 overexpression fully restored PDGF-induced proliferation, migration, and signaling pathways in SMC exposed to high glucose and hypoxia (Clément et al., 2021).As for the role of VSMC in MMD, a previous study reported VSMC over-proliferation in the intima and intimal thickening of pathological blood vessels in MMD (Takekawa et al., 2004).This generates evidence that PTPN6 overexpression in MMD may promote VSMC proliferation and subsequently cause intimal thickening.In addition, PTPN6 expression was significantly positively correlated with the expression of dynamin 2 (DNM2) in this study, a GTP-binding protein subfamily and an important component of the cellular cytoskeleton (Ferguson and De Camilli, 2012).Interestingly, our previous study also found that cytoskeletal remodeling in ECs plays an important role in MMD pathogenesis (He et al., 2023b).This suggests that PTPN6 and DNM2 are involved in MMD by participating in regulating the cellular cytoskeleton.
Our study has limitations.First, the sample size collected in this study was relatively small, which may have reduced the representativeness of the study and weaken the level of evidence when drawing conclusions about causality between identified genes and MMD pathogenesis.Further study needs to be conducted in a larger and ethnically diverse patient population to strengthen the generalizability of the findings.Second, the key genes and their functions in this study are lack of validation through in vitro experiments.Functional experiments are needed to validate the impact of the identified genes on relevant cellular processes in MMD, for example, the cell proliferation assays or migration assays on the gene overexpression or knockdown cell lines.Finally, in order to avoid the omissions of some potential biologically significant genes in the small sample size, this study used the criteria of p < 0.05 and | log2FC | > 1 to identify DEGs, without using adjusted p-values or FDR values.This method may lead to a higher false positive rate, which is also a limitation of this study.

Conclusion
In summary, this study identified four key OXPHOS-related genes associated with MMD: CSK, NARS2, PTPN6, and SMAD2.We systematically investigated the relationship between four key OXPHOS-related genes and the pathogenesis of MMD We also analyzed the profile of immune infiltration and the microenvironment of these genes.Bioinformatic analysis showed that the four key genes are related with angiogenesis, endothelial cell and VSMC proliferation, and VSMC phenotype switching, which may be involved in the mechanisms of MMD.Further studies will be needed to investigate the role of these identified genes in the pathogenesis of MMD in a larger and ethnically diverse patient population, or by in vivo and in vitro experiments.In conclusion, this study provides a comprehensive understanding of the important role of OXPHOS in MMD pathogenesis and will be conducive to the development of new treatment methods for MMD.

FIGURE 1
FIGURE 1Screening of the phosphorylation-associated genes for MMD and functional enrichment.(A) PCA plot for gene expression data of the participants.(B) Volcano plot of DEGs differentially expressed between MMD patients and control participants in the GEO dataset.(C) Venn plot of the overlapping DEGs and phosphorylation-associated genes.Finally, 27 intersecting genes were selected.(D) GO enrichment analysis for the 27 intersecting genes.E, KEGG enrichment analysis for the 27 intersecting genes.MMD, Moyamoya disease; PCA, Principal component analysis; GEO, Gene Expression Omnibus; DEGs, Differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

FIGURE 2
FIGURE 2Identification for the key phosphorylation-associated genes for MMD by machine learning algorithm and functional enrichment.(A, B) Identification for characteristic genes for MMD by LASSO regression algorithm.Finally, 13 intersecting genes were selected.(C) Identification for characteristic genes for MMD by SVM algorithm.Finally, 7 intersecting genes were selected.(D) Venn plot of the overlapping genes identified by LASSO regression and SVM algorithm.Finally, 4 intersecting genes were selected as key phosphorylation-associated genes for MMD.(E) GSEA analysis for gene CSK.(F) GSEA analysis for gene NARS2.(G) GSEA analysis for gene PTPN6.(H) GSEA analysis for gene SMAD2.MMD, Moyamoya disease; LASSO, Least Absolute Shrinkage and Selection Operator; SVM, Support Vector Machine; GSEA, Gene set enrichment analysis.

FIGURE 3
FIGURE 3 Analysis of immune infiltration and microenvironment of the four key phosphorylation-associated genes for MMD.(A) Heatmap for the proportion of 19 kinds of immune cells in each sample.(B) Thermodynamic plot of the relationship between each kind of immune cells.(C) The relationship between the gene CSK (upper of the figure) and NARS2 (lower of the figure), and immune cells.(D) Violin plot for the difference of immune infiltration between MMD and control group.MMD = Moyamoya disease.