The blood gene expression profiles of 166 samples which incorporate three phases of MI (D0: acute MI, D30: 30-days post-MI, and Y1: 1-year post-MI) were downloaded with the gene expression omnibus (GEO) under accession number of GSE123342 (Vanhaverbeke et al., 2019). There were 65 D0, 64 D30, and 37 Y1 samples. There were 70,523 probes in Affymetrix Human Transcriptome Array 2.0 corresponding to 30,905 genes. The probes for the same gene were averaged and the data was quantile normalized (Bolstad et al., 2003). We wanted to find the genes with changed expression patterns in post-MI.

Monte Carlo feature selection has been a widely used method for feature selection (Chen et al., 2018a, 2019a,b,d, 2020; Pan et al., 2018, 2019a,b; Wang et al., 2018; Jiang et al., 2019; Li et al., 2019). It was originally developed by Draminski et al. (2008). It randomly constructed many tree classifiers of the sub datasets from the original dataset and assigned the importance to a feature based on how much it participated in the tree classifiers. The java software dmLab¹ with default parameters (Draminski et al., 2008) was used to apply the Monte-Carlo feature selection method.

To be more specific, the original dataset was divided into s subsets of m features (m<<d, where d is the total number of features, i.e., 30,905 genes in this study). Then, for each subset, t trees were constructed. Therefore, a total of s⋅t classification trees were constructed. At last, the relative importance (RI) of each feature was estimated as follows:

where IG(n_g(τ)) was the information gain (IG) of node n_g(τ), (no.inn_g(τ)) was the number of samples in node n_g(τ), (no.inτ) was the number of samples in tree τ, wAcc was the weighted accuracy over all samples, and u and v were two regular factors which were set as default.

After running MCFS, all features can be ranked based on their RI. The higher the RI, the more important a feature was.

With MCFS, all features were ranked. But we still did not know how many genes we should choose. Ideally, we wanted the number of selected genes to be small but their classification performance to be great. To find the balance and the optimal signature, we adopted IFS (Chen et al., 2018b,d; Lei et al., 2018; Li and Huang, 2018; Sieber et al., 2018; Zhang et al., 2018; Wang and Huang, 2019; Yan et al., 2019). During IFS, a serial of feature sets F = [f₁,f₂,…,f_N] were constructed. N ranged from 1 to 1000. For each feature set, we constructed corresponding support vector machine (SVM) classifiers using the R function svm with default parameters in package e1071² and evaluated the performance using leave-one-out cross validation (LOOCV). Therefore, we can get a serial of LOOCV accuracies which corresponded to different feature sets with various numbers of features. With the help of the IFS curve, we can balance the model complexity and classification performance. If the number of features was too small, the performance would be bad. If the number of features was too large, too much noise would be introduced and the performance would decrease. The optimal selection would be achieved when the number of features was small and the accuracy was high.

The biological functions of the optimal MI signature genes were analyzed using hypergeometric enrichment analysis (Shi et al., 2018a, b). The significance of the signature genes onto Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, Gene Ontology (GO) biological process (BP), molecular function (MF), and cell component (CC) were represented with hypergeometric p values.

In this study, we exploited newly published gene expression profiles of patients with MI (Vanhaverbeke et al., 2019). Each patient was represented by 30,905 gene expression features. We integrated expression profiles of all patients into one matrix for quantile normalization followed by applying the MCFS method for ranking analysis. Each feature was assessed by estimating the relative importance (RI) value. After evaluating all features, we generated a feature list F in descending order of RI values of features. The ranked features with RI values were provided in Supplementary Table S1.

According to the feature list obtained by the MCFS algorism, the IFS method was employed to identify optimal feature sets which could train the best performance for SVM. To save computing time, we established the series of feature subsets (F₁, F₂, F₃, …, F₁₀₀₀) based on the top 1 to 1000 genes in F. For each feature set, we established a classifier by SVM algorithm and estimated optimal parameters through Leave-One-Out Cross-Validation (LOOCV). The LOOCV accuracies on multiple feature subsets were shown in Figure 1, from which we can see that the accuracy reached a plateau area when the top 134 features were used for building the classifier. The 134 optimal features were listed in Supplementary Table S2. The confusion matrix of the predicted results using the 134 features was shown in Table 1. It can be seen that all three classifiers had a great performance.

In order to confirm the performance of identified optimal features/genes representing different phases of samples, we performed cluster analysis on expression profiles of 134 optimal genes in 166 samples which incorporate three phases of MI (D0: acute MI, D30: 30-days post-MI, and Y1: 1-year post-MI). We used a heatmap to visualize the expression of such optimal genes among three groups of samples (Figure 2). The cluster tree illustrated that most samples belonging to the same phase can be clustered together and different phases were classified into different branches. In addition, these optimal genes were also classified into three clusters which correspond to high expression in three phases. The largest gene cluster with 90 genes was highly expressed in D0, the cluster with 16 genes had a high expression of D30, and the cluster with 28 genes was highly expressed in Y1.

The expression levels of genes like KLHL8, HCLS1, MOB3A, IL17RA, ETF1, ZFAS1, CRK, MXD1, UBXN2B, FCAR, and EXTL3 decreased in post-MI while the expression levels of genes like DCK and RNU4-7P increased in post-MI. We plotted the boxplots of several representative genes in Figure 3. For example, in Figure 3C, the expression levels of FCAR on D0 was significantly higher than on D30 and the expression levels on D30 was significantly higher than on Y1. There was a consistent post-MI trend of FCAR. These expression patterns may reveal the mechanisms of MI. FCAR is a member of the immunoglobulin superfamily and encodes a receptor for the Fc region of IgA. The cell surface receptors for immunoglobulin, such as the protein of FCAR, can activate many inflammatory processes involved in atherosclerosis and coronary artery disease (Daëron, 1997; Gavasso et al., 2005). The variation in FCAR which causes an amino acid alteration was found to increase the risk of MI and coronary heart disease, indicating the potential functional role of FCAR in the development of cardiovascular disease (Iakoubova et al., 2006, 2008).

We next performed functional enrichment analysis on these 134 optimal features/genes. A hypergeometric distribution test was applied to calculate p value to determine the significantly enriched entries. Firstly, we performed Gene Ontology enrichment analysis on the gene set. In biological progress aspect, the top 3 GO terms were GO: 0044264, GO: 0046903, and GO: 0005976, which correspond to cellular polysaccharide metabolic process, secretion, and polysaccharide metabolic process, respectively (Supplementary Table S3). The top GO term of cellular component was GO: 0005964, corresponding to phosphorylase kinase complex (Supplementary Table S4). The most significantly enriched GO term of molecular function was GO: 0004908, which was annotated to interleukin-1 receptor activity (Supplementary Table S5). Secondly, KEGG enrichment analysis was applied to discover the signaling pathways involved in these optimal genes. In this part, we found the insulin signaling pathway (hsa04910) was the top enriched KEGG pathway (Supplementary Table S6).

To investigate the correlation of optimal genes, we applied gene interaction analysis on 134 features/genes to construct gene interaction networks. Proteins encoded by such classes of genes were input into a STRING database (Szklarczyk et al., 2018), mining interaction relationship. Although part of the genes showed no association with other genes, we found an interaction network consisting of dozens of genes and predicted three hub genes, including IL1R1, TLR2, and TLR4 (Figure 4), which may interact with each other to play a non-negligible role in the progression of MI.

IL1R1, TLR2, and TLR4 showed promising associations with MI. It was reported that the knockout of IL1R1 caused a reduction of leukocyte production after MI, leading to a decreased inflammation with better outcome (Sager et al., 2015). In another mice study, the up-regulated IL1R1 at 7 days post-MI prolonged the inflammation by suppressing neutrophil apoptosis (Iyer et al., 2015).

TLR2 plays a fundamental role in the activation of innate immunity (Binder et al., 2002). There are usually high levels of cytokines that result in inflammation in MI patients; TLR2 served as a key receptor to activate the corresponding pathways (Pagano et al., 2012). The experimental data indicated that circulatory TLR2 is relevant to different manifestations of myocardial I/R injury (Arslan et al., 2010). And the inhibition of TLR2 has beneficial effects on I/R injury in a murine model of MI (Arslan et al., 2009). TLR2 is the key receptor which can induce the inflammation after MI, therefore many MI-related genes show close interactions with TLR2.

TLR4 regulates the cytokines after cardiac damage (Arslan et al., 2010). Activation of TLR4 was related to myocytic inflammatory reaction in MI patients 14 days after onset, suggesting that TLR4 signaling plays a role in the progress after MI (Satoh et al., 2006).

Using the feature selection, 134 genes were extracted and exhibited an excellent performance in our prediction model of SVM, suggesting that these genes may participate in the progression of MI. Here, we took some of the selected genes as examples to give a detailed discussion to validate the relevance of a given gene in distinguishing different pathological phases of MI. Through a literature review, several experimental evidences or analysis results have been found to confirm the reliability of our prediction.

Given that the selected 134 genes were deemed as important features for the classification of different phases of MI, we performed GO and KEGG functional enrichment analysis to explore the key biological processes or pathways during the progress of disease. As shown in Supplementary Tables S3–S6, we analyzed the enriched GO terms and KEGG pathways which showed statistical significance. A detailed discussion was given about the linkage between certain functional sets and MI.

Based on the enrichment results of 134 selected genes, we found some GO biological process terms with high scores turned out to be involved in the polysaccharide metabolic process, including GO: 0044264 and GO: 0005976. As early as 1965, scientists have noticed the important role of glucose load in MI (Cohen and Shafrir, 1965). Recent studies reported that certain polysaccharide compounds can affect myocardial injury via regulating the inflammation response (Li et al., 2011; Lim et al., 2016). As demonstrated by experiments on rat, the polysaccharide extract from Momordica charantia down-regulated the expression of NF-kappaB and ameliorated oxidative stress and inflammation, which caused a cardioprotective effect against MI (Raish, 2017). Polysaccharide metabolism plays an important role during the progression of MI, so the biological processes related to polysaccharide metabolism are meaningful and can be used to indicate the progression of disease based on its specific pattern.

Apart from GO terms that belong to biological processes, we found these 134 genes are also enriched in a cellular components term GO: 0005964 with the highest probability. GO: 0005964 refers to phosphorylase kinase complex. For cardiomyocytes, the storage of glycogen is important during the emergency situation. Increasing Ca²⁺ concentration in cytosol can induce glycogenolysis by the activation of phosphorylase kinase, which can alleviate myocardial damage during MI or cardiac surgery (Raish, 2017). In fact, some phosphorylases have been applied in the diagnosis of myocardial ischemic injury and infarction since the serum level of phosphorylase showed a signature with the diseases (Krause et al., 1996). It is reasonable for the MI-related genes to be enriched in such GO term that would mean the phosphorylase play a crucial role during the progression of MI.

The most enriched GO terms of molecular function turned out to be interleukin-1 (IL-1)-related functions including GO: 0004908 and GO: 0019966, which represent IL-1 receptor activity and IL-1 binding, respectively. An interleukin-1 receptor gene ST2 was increased in the serum after MI, suggesting that this gene may participate in innate immunity during myocardial injury (Weinberg et al., 2002). What’s more, ST2 was reported to be able to predict the clinical outcome in AMI due to its role in cardiac pathophysiology (Shimpo et al., 2004). Many publications have observed the elevated serum level of IL-1 receptor in patients with AMI (Shibata et al., 1997; Balbay et al., 2001). These findings proved the important role of IL-1 in the progression of MI, and confirmed the relation between selected genes and MI.

The KEGG pathways enrichment analysis provided various pathway results. Among these, the highest enriched pathway turned out to be hsa04910, which is an insulin signaling pathway. Increased insulin can promote the metabolism of glucose to maintain the balance of blood glucose. The connection between abnormal insulin signaling and heart disease has already been reported, in that diabetes mellitus significantly increased the risk of ischemic heart disease (Miettinen et al., 1998). Insulin can protect cardiomyocytes from apoptosis through activating downstream pathways such as PI3K and Akt (Yao et al., 2014). It was reported that impaired insulin signaling will cause the dysfunction of mitochondria after MI due to the reduced glucose transport and oxygen content (Sena et al., 2009). Thus, the insulin signaling pathway is important during the progression of MI and influences the pathological degree of disease.