The GSE61145, GSE60993, and GSE34198 gene expression datasets related to human MI, which were developed based on the GPL6106, GPL6884, and GPL6102 platforms, respectively, were downloaded from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) database. The GSE61145 dataset contained data on 14 blood samples from patients with MI and 10 samples from normal controls. The GSE60993 dataset included data on a total of 24 samples (7 and 17 blood samples from normal controls and patients with MI, respectively). Finally, the GSE34198 dataset contained 97 samples (48 and 49 blood samples from normal controls and patients with MI, respectively).

The raw data were downloaded and the probes were annotated into gene symbols based on platform annotation information. Because a single gene could correspond to several probes (multiple values), the average gene expression values were calculated for each gene. Afterward, log2 conversion was performed to transform the gene expression data from a skewed distribution to an approximately normal distribution. The data were then normalized using the limma package (MetaDE) (http://www.bioconductor.org/packages/2.9/bioc/html/limma.html) in R language.

DEGs were screened by using the MetaDE package (Wang et al., 2012) in R based on the GSE61145 and GSE60993 datasets. The raw data were downloaded and the probes were annotated into gene symbols based on platform annotation information. The average gene expression values were calculated for each gene. Afterward, log2 conversion was performed to transform the gene expression data from a skewed distribution to an approximately normal distribution. The data were then normalized using the limma package in R language. The heterogeneity of gene expression data based on different platforms was analyzed using the MetaDE.ES method (Kim et al., 2022), with tau² = 0 and Qpval >0.05. Differential expression analysis of genes with homogeneous expression was then performed between the disease and control groups, with an FDR (false discovery rate) of < 0.05 defined as the threshold value.

Weighted gene co-expression network analysis (WGCNA) (MetaDE) is a typical system biology algorithm used to construct gene co-expression networks based on high-throughput mRNA expression data. The genes and modules related to MI in this study were analyzed for DEGs based on the WGCNA algorithm (Langfelder and Horvath, 2008). The correlation coefficient between gene expression was calculated using the function S m n = | c o r ( m , n ) | . Then, the coefficient was then weighted by the exponential adjacency function a m n = p o w e r ( S m n , β ) . According to the principle of scale-free networks, the weight coefficient β was determined for the adjacency function. To measure the dissimilarity between nodes, the correlation matrix (S_mn) was transformed into an adjacency matrix (a_mn). A hierarchical clustering tree was then constructed based on the dissimilarity coefficients between genes, with different branches of the clustering tree representing different gene modules. Finally, t-tests were used to analyze the correlations between network modules and disease states.

The co-expressed modules that were closely associated with disease state were analyzed and the module genes were collected to construct the co-expression network. Genes related to disease were subjected to gene ontology (GO) enrichment analysis using the clusterProfiler package (Yu et al., 2012) in R based on the hypergeometric distribution algorithm. P < 0.05 was defined as the threshold value. The formula for the hypergeometric distribution algorithm is shown in Eq. 1. p = 1 − ∑ i = 0 H − 1 ( M i ) ( N − M K − H ) ( N K ) , (1) where N represents the number of genes with GO functional annotations; K represents the number of DEGs among N genes, and M represents the number of genes that are annotated with a particular GO function.

To classify the samples, an optimal classification hyperplane must be selected from numerous options to maximize the distance δ between the sample set and the classification hyperplane. When ε = |wx _i + b| = 1, the distance between the two types of sample points is 2 ((|wx _i + b|)/‖w‖) =(2/‖w‖). The goal is to construct an optimal classification hyperplane under the constraint of Eq. 2 to maximize (2/‖w‖) and minimize (‖w‖2/2). { w x i + b ≥ 1 , y i = 1 w x i + b ≤ − 1 , y i = − 1 i = 1,2 ⋯ , l . (2)

Most classification issues can be treated as nonlinear separable problems, and quadratic programming problems can be modified as follows by introducing the slack variable ξi in Eq. 3. { min 1 2 ‖ w ‖ 2 + C ∑ ξ i , ξ i ≥ 0 constraint   condition :   y i ( ( w x i + b ) ) ≥ 1 − ξ i , i = 1,2 ⋯ , l (3) where ξ _i is the slack variable and C is the penalty coefficient.

With GSE61145 as the training dataset and all genes of interest as classification factors, the SVM model was established using the e1071 package (MetaDE) in Rto distinguish the disease and control samples. The classification factors were added individually until all of them had been added to the SVM classifier. The classification accuracy of the SVM classifier was then calculated and the genes that affected classification accuracy were removed. The SVM model was then validated in the GSE60993 and GSE34198 gene expression datasets.

A total of 1,231 DEGs were identified by the MetaDE package (Langfelder and Horvath, 2008; Wang et al., 2012; Yu et al., 2012; Meyer, 2013). The top 10 DEGs are listed in Table 1, including GZMK (granzyme K), HLA-DQA (histocompatibility complex, class II, DQ alpha), and EOMES (eomesodermin). First, the heterogeneity of gene expression data based on different platforms was analyzed using the MetaDE.ES method, with tau2 = 0 and Qpval >0.05. Then, the differential expression analysis of genes with homogeneous expression was conducted between the disease and control groups, with an FDR (false discovery rate) of <0.05 defined as the threshold value. A total of 1,231 DEGs were identified. The top 10 DEGs with the smallest p-values in the gene difference analysis between the disease and control groups were selected; that is, the genes with the largest difference between disease and control groups. The present study analyzed the co-expressed modules that were closely associated with the disease state and identified the module genes to construct a co-expression network. Genes related to disease were subjected to gene ontology (GO) enrichment analysis. KEGG pathway enrichment analysis was not performed.

To satisfy the precondition of scale-free network distribution, we selected a power of 18 as the adjacency parameter. The results of the consistency analysis showed a high correlation between the GSE61145 and GSE60993 datasets (correlation coefficient = 0.86, p-value < 1e-200). Additionally, GSE61145 was used as a training set to identify disease-associated modules (Figure 1A). Module partitioning for the GSE60993 dataset (Figure 1B) showed high consistency with the GSE61145 dataset. We then calculated the correlation coefficient between module and disease state (normal and MI samples) for the GSE61145 (Figure 2A) and GSE60993 (Figure 2B) datasets, respectively (Table 1). According to the correlation coefficients, the top three modules (black, pink, and red) were identified.

The correlation coefficients between genes in the top three modules and disease state were calculated, which revealed 98 genes with correlation coefficients >0.5. These included 30 genes (11 up-regulated and 19 down-regulated) in the black module, 19 genes (9 up-regulated and 10 down-regulated) in the pink module, and 49 genes (22 up-regulated and 27 down-regulated) in the red module. The gene co-expressed networks of the 98 genes were then constructed (Figure 3). GO analysis showed the enrichment of 10 GO terms among the genes in the black module (Table 2) and 15 GO terms among the genes in the red module (Table 2). The GO terms enriched in the black module included negative regulation of cell proliferation (p-value = 0.009704), regulation of cell proliferation (p-value = 0.014724), and positive regulation of macromolecule metabolic process (p-value = 0.019608). The GO terms closely related to the genes in the red module mainly included positive regulation of I-kappaB kinase/NF-kappaB cascade (p-value = 0.024598), regulation of I-kappaB kinase/NF-kappaB cascade (p-value = 0.029497), and positive regulation of signal transduction (p-value = 0.037345). No GO terms were significantly enriched among the genes in the pink module.

Based on the SVM classification model, we removed genes that could not distinguish between the disease and control samples. Finally, seven genes were obtained: ACOX1 (Acyl CoA oxidase 1), ADCK2 (aarF domain containing kinase 2), AFF3 (AF4/FMR2 family member 3), BCL6 (B-cell lymphoma 6), CEACAM8 (Carcinoembryonic antigen-related cell adhesion molecule 8), CUGBP2 (CUG triplet repeat-binding protein 2) and GPX7 (glutathione peroxidase 7). The SVM classification model of these seven genes could distinguish all samples in the GSE61145 dataset. The scatterplot of the GSE61145 dataset is shown in Figure 4A. The GSE60993 and GSE34198 datasets were then used as validation datasets to confirm the SVM classification model. As shown in Figure 4B, the SVM classification model correctly distinguished 23 (17 disease and 6 normal samples) of 24 samples in the GSE60993 dataset. Additionally, the scatterplot of the GSE34198 dataset indicated that the SVM classification model correctly distinguished 90 (48 disease and 42 normal samples) of the 97 samples (Figure 4C). The efficiency receiver operating characteristic (ROC) curves of the SVM classification model are shown in Figure 5 and the efficiency parameters of each dataset are listed in Table 3.