The GSE51981, GSE6364, and GSE7307 datasets were downloaded by the R package “GEOquery” (2.60.0) (Davis and Meltzer, 2007) to obtain the expression profile data. Then, the E-MTAB-694 dataset was downloaded through the Array-Express database. The related annotation information including the platforms, the probes, and ID conversion was obtained from the GEO database. When multiple probes corresponded to one gene symbol, the average expression level of multiple probes was used as the expression level of the corresponding gene. ID conversion was conducted with the R package “org.Hs.eg.db” (v3.13.0). Furthermore, the “removeBatchEffect” function in the R package “LIMMA” (v3.48.3) (Ritchie et al., 2015) was used to adjust batch effects, which were evaluated by principal component analysis (PCA).

Differential expression analysis was conducted on 77 EM disease and 71 normal samples of the GSE51981 dataset through the Bayesian analysis of the R package “LIMMA”. The log2FoldChange > 1.5 and p-value < 0.05 were set as the threshold of DEGs. The R package “pheatmap” (v1.0.12) was used to perform clustering analysis of DEGs for the heatmap. To explore the biological significance of these DEGs in the pathogenesis of EMs, GO and KEGG pathway enrichment analyses were performed through the R package “clusterProfiler” (v4.1.3) (Wu et al., 2021) to identify significantly enriched GO terms and significantly enriched KEGG pathways with the threshold of p-value < 0.05.

The STRING (v11.5) (https://string-db.org/cgi/input.pl) (Szklarczyk et al., 2021) has been widely applied to construct a protein–protein interaction (PPI) network. Based on those DEGs, the “Multiple proteins” option was selected. In the PPI network, the minimum required interaction score was set as “high confidence (0.700)”. Then, the cytoHubba (Chin et al., 2014) was employed to identify hub genes. The eccentricity algorithm was selected and 15 top-ranked genes were chosen as hub genes. Finally, the hub gene network was visualized with Cytoscape (v3.9.0) (Demchak et al., 2014).

The R package “randomForest” (v4.6.14) (Liaw et al., 2014) was used to construct a random forest model to screen DEGs. The number of random seeds and decision trees was set as 1–5,000 and 3,000 in the random forest classifier originally, respectively. Finally, the number of random seeds and decision trees was set as 4,543 and 219, respectively, which represented higher accuracy of the constructed model and stable model error. The Gini coefficient method was used to obtain the dimensional importance value of all variables from the constructed random forest model. Those DEGs with an importance value greater than 4 were screened as important genes of EMs for subsequent model construction and verification. The R package “pheatmap” was used to perform clustering analysis of the screened important genes for the heatmap in this dataset.

The GSE6364 dataset downloaded through the R package “GEOquery” was selected as the training set for the construction of the artificial neural network model. After the data normalization, the R package “neuralnet” (v1.44.2) (Fritsch and Guenther, 2016) was used to construct an artificial neural network model of those important variables. The number of hidden neuron layers should be two-thirds of the number of the input layer plus two-thirds of the number of the output layer. Therefore, six hidden layers were set as the model parameter to construct a classification model of EMs through the predicted gene weight information. The R packages “pROC” (v1.18.0) (Robin et al., 2011) and “ggplot2” (v3.3.5) (Gómez-Rubio, 2017) were used to calculate the verification results of AUC classification performance and draw the ROC curve. Another two datasets E-MTAB-694 and GSE7307 were used to verify the accuracy of the constructed neural network model for the diagnosis of EMs. The R package “pROC” was used to draw ROC curves for each dataset, and the AUC value was calculated to verify the classification efficiency. Meanwhile, the sensitivity and specificity in distinguishing the disease samples from normal samples were calculated.

The R package “GEOquery” was used to download the GEO dataset GSE51981 (77 EM disease samples and 71 normal samples) and obtain detailed information. We used the “removeBatchEffect” function in the R package “LIMMA” to adjust batch effects and then conducted principal component analysis (PCA) analysis to evaluate the performance of batch effect adjustment. PCA results ( Supplementary Figure S1) indicated that the disease samples were mixed with the normal samples, which suggested the challenge of diagnosing. We also used the R package “LIMMA” to perform differential expression analysis for the dataset GSE51981 through the Bayesian test. We finally identified 2,267 significantly upregulated and 285 significantly downregulated expressed genes between the disease samples and the normal samples with the threshold of fold change values of >1.5 and p < 0.05. The detailed information of all DEGs is listed in Supplementary Table S1. The results of these DEGs and the heatmap of these DEGs are visualized in Figures 2A and 2B, respectively.

To explore the biological significance of these DEGs in the pathogenesis of EMs, we performed GO and KEGG pathway enrichment analyses through the R package ‘clusterProfiler’. GO terms were classified into three categories: biological process (BP), cellular component (CC), and molecular function (MF). The top five GO terms of genes with significantly upregulated and downregulated expression levels were visualized in Figures 3A,B. The GO enrichment analysis results indicated that these significantly upregulated expressed genes were mainly involved in the transmembrane transporter activity, ATPase activity, metallopeptidase activity, aldehyde dehydrogenase NADP⁺ activity, and lipid transporter activity (Supplementary Table S2), while these significantly downregulated expressed genes were mainly involved in the flavin adenine dinucleotide binding, acyl-CoA dehydrogenase activity, phosphatidylcholine transporter activity, extracellular matrix structural constituent, and ATPase-coupled intramembrane lipid transporter activity (Supplementary Table S3). For KEGG pathway enrichment analysis, the results indicated that these upregulated expressed genes were significantly associated with the cAMP signaling pathway, adrenergic signaling in cardiomyocytes, aldosterone synthesis and secretion, ABC transporters, and salivary secretion (Supplementary Table S4), while these downregulated expressed genes were significantly associated with fatty acid degradation and metabolism; valine, leucine, and isoleucine degradation; lysosome; the PPAR signaling pathway; and the Hippo signaling pathway (Supplementary Table S5). Furthermore, we constructed a PPI network through the STRING database. The hub genes selected from the PPI network are shown in Supplementary Figure S2. According to the eccentricity scores, we identified 15 hub genes from the network, which had highest confidence scores.

To screen DEGs, we put these DEGs into the random forest classifier and set the number of random seeds to 4,543. By referring to the relationship between the model error and the number of decision trees (Figure 4A), we selected 219 trees as the parameter of the random forest model, which represented a stable error in the model. In the modeling process, we used the Gini coefficient method to measure the importance of all variables according to decreased mean square error and model accuracy (Figure 4B). Finally, we selected seven DEGs (AHI1, DMXL2, NAA16, CCDC22, CISD3, COMT, and EIF3E) with a mean decrease of Gini index greater than 4 as important variables for subsequent analysis. Interestingly, all these DEGs were included in the 15 hub genes identified from the constructed PPI network. Among these variables, AHI1 was the most important, with the mean decrease of the Gini index being much higher than other variables (Supplementary Table S6). A small number of variables meant a small out-of-band error, which represented a high accuracy of the constructed random forest model. Based on these seven variables, we performed the k-means clustering of the dataset. The results suggested that these seven genes could be used to distinguish the disease sample from the normal samples (Figure 4C). Furthermore, AHI1, DMXL2, and NAA16 genes were clustered as a group with low expression in the normal sample and high expression in the disease sample. On the contrary, CCDC22, CISD3, COMT, and EIF3E were clustered as another group with high expression in the normal sample and low expression in the disease sample.

Based on the R package ‘neuralnet’, we use the GSE6364 dataset (21 disease samples and 21 normal samples) as the training set to construct the artificial neural network model. Firstly, we performed the preprocessing and normalization of this dataset. According to the output results of the neural network model (Figure 5A), it is illuminated that the entire training was performed in 11,684 steps. Among the output results, the predicted weights of each hidden neuron layer were −3.97906, 1.04457, 2.76611, −2.00181, −11.84206, and −0.90829 (Supplementary Table S7). Next, we drew the ROC curve to evaluate the predicted performance; the AUC values of AHI1, COMT, DMXL2, CISD3, NAA16, EIF3E, and CCDC22 were 0.7150, 0.7809, 0.6927, 0.7266, 0.7217, 0.7093, and 0.7050, respectively (Figure 5B). The larger the AUC value of each DEG is, the higher the credibility of the constructed diagnostic model will be. We also used another two datasets E-MTAB-694 (18 disease samples and 17 normal samples) and GSE7307 (18 disease samples and 23 normal samples) to verify the accuracy of the constructed neural network model. In the E-MTAB-694 dataset (Figure 5C), the AUC values of the seven DEGs were 0.8226, 0.6623, 0.6836, 0.6625, 0.8367, 0.8471, and 0.8617. In the verification results of the GSE7307 dataset (Figure 5D), the AUC values of the seven DEGs were 0.7464, 0.6484, 0.7020, 0.6300, 0.9075, 0.8295, and 0.8327. In general, we constructed a novel diagnostic model of EMs and verified its diagnostic efficacy through the constructed artificial neural network in two public datasets.