Development and Validation of a Diagnostic Model Based on Hypoxia-Related Genes in Myocardial Infarction

Introduction

Myocardial infarction (MI) is one of the leading causes of human mortality.¹ MI is characterized by myocardial ischemia and hypoxia caused by coronary artery occlusion, which in turn triggers inflammation and apoptosis.² Extensive studies support that the inflammatory response and myocardial apoptosis are important in the pathogenesis of acute MI.³ When MI occurs, hypoxia leads to injury and necrosis of cardiomyocytes, which release pro-inflammatory chemokines to recruit inflammatory cells to the infarcted area and remove necrotic cardiomyocytes.² Although early treatment with reperfusion and pharmacological approaches significantly reduced mortality after MI, in-hospital mortality remained high.⁴ Therefore, it is desperately warranted to illustrate the underlying pathogenesis of MI.

Hypoxia plays an important role in MI. Myocardial cells undergo ischemia and hypoxia, followed by apoptosis, inflammation of the damaged tissue, and infiltration of immune cells.⁵ Hypoxia was demonstrated to serve as a possible regulator of the activity of epicardial mesothelial cells after MI.⁶ Hypoxia can promote neutrophil survival after acute MI.⁷ Kologrivova et al indicated that cells of the immune system are involved in resolution of inflammation and repair in MI.⁸ Immunoregulatory therapy may have great potential to accelerate cardiac repair and improve ventricular remodeling after MI. Hence, it is important to explore the relationship between hypoxia, immunity and MI for providing new therapeutic methods.

In this study, we identified different MI subtypes based on hypoxia-related genes and investigated the differences among these MI subtypes. A diagnostic model was developed to distinguish different subtypes of MI, which may be conducive to understand the molecular mechanisms underlying MI and potentially provide an approach for personalized therapy of MI patients.

Materials and Methods

Datasets Selection and Data Preprocessing

The gene expression profiles of MI were acquired from the GEO database by searching keywords “MI (myocardial infarction)” and “Homo sapiens”. The datasets that met the following criteria were included in our study: no less than 5 samples and including control samples. After filtering, two datasets (GSE59867 and GSE62646), treated as training sets, and two datasets (GSE123342 and GSE48060), served as validation sets, were included in the study. GSE59867 consisted of 111 patients with MI and 46 controls, and GSE62646 included 28 MI cases and 14 controls. GSE123342 consisted of 67 patients with MI and 22 controls, and GSE48060 included 26 MI cases without recurrent events and 21 controls. All these four datasets examined the blood samples. For the MI group, samples with the earliest sampling time (<48h) available in each dataset were selected for analysis. In addition, 200 hypoxia-related genes (HRGs) in “HALLMARK_HYPOXIA” were retrieved from MSigDB.

Differentially Expressed Genes (DEGs) Analysis

The “limma” package in R was applied to acquire the DEGs between MI and controls with adj_p < 0.05 and |log₂ (Fold Change)| > 0.4. The ggplot2 R package was used to draw volcano plot and heatmap. David was employed to perform GO classification and KEGG pathway enrichment analysis (p < 0.05).

Unsupervised Clustering and DEG Analysis Among Different MI Subtypes

Based on the expression levels of HRGs, the “ConsensusClusterPlus” package was used to classify the patients with MI into different subtypes. With adj_p < 0.05 and |log₂ (Fold Change)| > 0.4, DEGs among different MI subtypes were identified by “limma” package in R. To investigate the difference on biological process among groups, GSVA enrichment analysis was performed using “GSVA” R packages. The gene sets of “c2.cp.kegg.v7.2. symbols” were downloaded from MSigDB database for running GSVA analysis. Differentially expressed pathways were identified using the “limma” package in R, and FDR < 0.05 was considered as statistically significant.

Estimation of Immune Microenvironment (IME) Cell Infiltration

The ssGSEA algorithm was applied to quantify the relative abundance of each cell infiltration in the IME of MI. The gene set for marking each IME infiltration immune cell type was acquired from a previous study.⁹ ImmPort database was applied to acquire immune response-related mRNAs. The enrichment scores calculated by ssGSEA analysis were used to represent the relative abundance of each TME infiltrating cell and activity of specific immune responses in each sample. The Wilcox test was used to compare the differences in immune cell infiltration and immune response activity among different MI subgroups. Furthermore, the expression of HLA-related genes was also compared among different MI subgroups.

Weighted Gene Co-Expression Network Analysis (WGCNA)

To construct a scale-free gene co-expression network, the top 25% genes were subjected to the subsequent analysis with WGCNA package in R. First, the samples were clustered with “hclust” function to eliminate the samples with obvious outliers. Then, the “pickSoftThreshold” function is used to select an appropriate soft threshold power regulator to build the scale-free topology. The height is set to 0.90. Then, the adjacency matrix is calculated according to the kernel value, and the adjacency matrix is transformed into topological overlap matrix and the corresponding dissimilarity matrix. Genes with similar expression patterns were grouped together, and modules were divided according to the “cutreeDynamic” function with default parameters. The minimal gene module size was set to 100 to obtain appropriate modules, and the threshold to merge similar modules was set to 0.25.

Identification of Hub Genes

The gene significance (GS) and module membership (MM) were calculated to select candidate hub genes. We defined the genes with |GS| > 0.3 and |MM| > 0.5 in hub modules as MI subtype-associated genes. The intersection of MI subtype-associated genes and DEGs between MI and controls was defined as hub genes. The GO and KEGG enrichment analyses were conducted for hub genes with David (p < 0.05).

Drug Prediction

We intersected the DEGs between MI group and control group with the DEGs among different MI subtypes. Then, potential candidate drugs were predicted based on intersection genes using the DGIbd database (https://dgidb.org) to offer novel insights for disease diagnosis, treatment and research.

Construction of Diagnostic Model for MI Subtype

The “pROC” package in R was used to draw ROC curve for each hub gene in training sets and validation sets. Then, the genes with area under the ROC curve (AUC) >0.85 in all subtypes were subjected to the subsequent analysis. The LASSO regression using “glmnet” package in R was utilized to develop a diagnostic model for MI subtype. ROC analysis was performed to determine the accuracy of the diagnostic model.

Real-Time qPCR (RT-qPCR)

A total of 28 blood samples were obtained from 13 patients with MI and 15 normal controls for RT-qPCR. Every individual provided the written informed consent. Relative gene expression was analyzed by 2^−ΔΔCT method. GAPDH and ACTB were used as internal reference. The study was approved by the ethics committee of the university.

Results

Differentially Expressed Analysis

Compared with normal controls, 389 DEGs, including 191 up-regulated mRNAs and 198 down-regulated mRNAs, were identified in MI with adj_p < 0.05 and |log₂ (Fold Change)| > 0.4 (Figure 1A and B). These 389 DEGs were enriched in protein binding, transmembrane signaling receptor activity, RNA processing, signal transduction and plasma membrane (Figure 1C). KEGG analysis revealed several significantly enriched pathways, including cytokine–cytokine receptor interaction, human T-cell leukemia virus 1 infection and complement and coagulation cascades (Figure 1D).

Figure 1 Differentially expressed gene (DEG) analysis. (A) The volcano plot of DEGs between MI and control group. (B) The heatmap of top 25 up- and down-regulated DEGs between MI and control group. (C) GO enrichment analysis of DEGs between MI and control group. (D) The KEGG analysis of DEGs between MI and control group.

Identification of Two Different Subtypes in MI

Consensus clustering revealed that the k = 2 was identified with optimal clustering stability (Figure 2A and B). Then, 139 patients in training sets were clustered into two subtypes, namely, cluster A (n = 82) and cluster B (n = 57) (Figure 2C). Consistent with analysis in training sets, the validation in validation sets displayed the similar results (Figure 2D and E). Then, 93 patients in validation sets were clustered into two subtypes, including cluster A (n = 30) and cluster B (n = 63) (Figure 2F). Consensus clustering revealed significant differences in molecular features between the two MI subtype. A total of 211 DEGs, including 180 up- and 31 down-regulated mRNAs, between two subtypes were identified in training sets (Figure 3A and B). GSVA analysis indicated that several pathways, including ECM receptor interaction, phenylalanine, metabolism, and complement and coagulation cascades, were activated in cluster A, and several pathways, including mismatch repair, DNA replication, aminoacyl tRNA, and biosynthesis, were activated in cluster B (Figure 3C).

Figure 2 Consensus clusters by HRGs in training sets (A–C) and validation sets (D–F). (A and D) Consensus clustering cumulative distribution function (CDF) for k = 2 to 5. (B and E) Relative change in area under the CDF curve for k = 2 to 5. (C and F) Consensus clustering matrix for k = 2.

Figure 3 Differentially expression analysis between two different subtypes in training sets. (A) The volcano plot of DEGs between cluster A and cluster B. (B) The heatmap of top 25 up- and down-regulated DEGs between cluster A and cluster B. (C) GSVA enrichment analysis showing the activation states of biological pathways between cluster A and cluster B.

The Difference of Immune Characteristics Between Two Subtypes

The results of ssGSEA analysis indicated that cluster A group has a higher level of infiltrated macrophage, activated dendritic cell, MDSC, eosinophil, mast cell, and monocyte, while activated CD4 T cell, activated B cell, and activated CD8 T cell are enriched in cluster B (Figure 4A). Meanwhile, more immune response activity was increased in cluster A, while significantly higher immunoscore in cluster B was detected (Figure 4B and C). Similarly, the expression of several HLA-related genes, such as HLA-A, HLA-C, HLA-DRB1 and HLA-DRB, were increased in cluster A (Figure 4D). As for the hypoxia and EMT status, we found that the enrichment scores of hypoxia and EMT2 were significantly higher in cluster A, while EMT1 and EMT3 were statistically nonsignificant between the two groups (Figure 4E). The same analysis for infiltrating immunocytes, activity of specific immune responses, immunoscore, HLA gene expression, hypoxia and EMT status were performed in validation sets. Similar trends were observed in validation sets (Figure S1A–E).

Figure 4 Characteristics of IME cell infiltration between cluster A and cluster B group in training sets. (A) The abundance of each IME infiltrating cell in cluster A and cluster B group. (B) Difference activity of specific immune responses between cluster A and cluster B group. (C) Difference in immune score between cluster A and cluster B group. (D) Difference in the HLA-related gene expression between cluster A and cluster B group. (E) Differences in EMT pathways and hypoxia condition between cluster A and cluster B group. *Indicates p < 0.05; **Indicates p < 0.01; ***Indicates p < 0.001; ****Indicates p < 0.0001; ns indicates p > 0.05.

WGCNA and Identification of Hub Genes

A co-expression network analysis was performed using the WGCNA R package. The power β = 13 was screened to ensure a scale-free network (Figure 5A). A total of 8 modules were determined, including a red module, green module, pink module, turquoise module, grey module, yellow module, brown module, and black module (Figure 5B). The correlation between module eigengenes and MI subtypes was analyzed. Compared with the other modules, the turquoise module and brown module had higher correlations with MI (Figure 5C). The results indicated that the turquoise module, including 880 genes, was positively associated with cluster B (r = 0.62, p = 3e-16) and the brown module, including 414 genes was negatively associated with cluster B (r = −0.58, p = 8e-14) (Figure 5D and E). Thus, the turquoise module and brown module were selected for further analysis. With |GS| > 0.3 and |MM| > 0.5, 1085 genes were identified as MI subtype-associated genes. Then, 211 hub genes were obtained by overlapping MI subtype-associated genes and DEGs between two subtypes. GO analysis indicated that these 211 hub genes were enriched in transcription factor activity, sequence-specific DNA binding, RNA processing, regulation of cholesterol efflux, and positive regulation of transcription from RNA polymerase II promoter (Figure S2A). KEGG analysis revealed that these 211 hub genes were enriched in Herpes simplex virus 1 infection and NF-kappa B signaling pathway (Figure S2B and C).

Figure 5 Hub module selection. (A) Determination of soft thresholding power in the WGCNA. (B) The cluster dendrogram of module eigengenes. (C) The module trait relationships. (D) A scatter plot of gene significance for MI versus the module membership in the turquoise module. (E) A scatter plot of gene significance for MI versus the module membership in the brown module.

Drug Prediction

The potential drugs for MI were analyzed with the DGIbd database by overlapping the DEGs between the MI group and the control group with the DEGs between two subtypes. The results indicated that 11 genes were observed as targets of the 60 predicted drugs, such as testosterone, pinocembrin and cilostazol (Figure S3).

Construction of Diagnostic Model for MI Subtype

A total of 9 genes with AUC > 0.85 in both two MI subtypes were obtained. Then, 5 genes, ANKRD36, HLTF, KIF3A, OXCT1 and VPS13A, were further identified as potential diagnostic markers (Figure 6A and B). Based on these 5 genes, we established logistic regression model for the diagnosis of MI subtypes. The ROC analysis revealed that the AUC of the diagnostic model is 0.975 in training sets, and the AUC is 0.920 in validation sets, proving that the diagnostic model has high credibility (Figure 6C and D).

Figure 6 Construction of diagnostic model. (A) LASSO coefficient profiles of 9 hub genes. (B) Selection of the optimal parameter (lambda) in the LASSO model. (C) ROC curves showing the predictive efficiency of diagnostic model in training sets. (D) ROC curves showing the predictive efficiency of diagnostic model in validation sets.

Correlation of Candidate Genes and Immune

The expression levels of 5 hub genes in training sets and validation sets were displayed in Figure S4A and B. The results exhibited that compared with control group, these 5 genes were significantly lowly expressed in cluster A, VPS13A, OXCT1 and HLTF were significantly decreased in cluster B in training sets; while compared with cluster B, these 5 gene expression levels were decreased in cluster A. In the validation sets, a similar situation exists. Compared with the control group, they were also significantly lower expressed in cluster A, but no significant difference was found in cluster B. Compared with cluster B, these 5 genes were also significantly lower expressed in cluster A. Correlation analysis of 5 hub genes was performed, suggesting that there was a strong positive correlation between the five genes (Figure S5). VPS13A had the strongest correlation with OXCT1 (0.88). The correlation between ANKRD36 and HLTF was the smallest, but also 0.65. Furthermore, correlation analysis of 5 hub genes and immune cells and immune score was performed (Figure 7). The results indicated that all five genes had strong positive correlations with immune scores and various immune cells, such as Activated B cell, Activated CD4 T cell, Activated CD8 T cell, Immature B cell, and Type 2 T helper cell, and strong negative correlations with other immune cells, such as Monocyte, Macrophage, Activated dendritic cell, Regulatory T cell, Neutrophil, Mast cell, Plasmacytoid dendritic cell, MDSC, and CD56bright natural killer cell. Among them, VPS13A had the strongest positive correlation with Activated CD4 T cell and the strongest negative correlation with Monocyte.

Figure 7 Correlation of candidate genes and immune cells and immune score. *Indicates p < 0.05; **Indicates p < 0.01; ***Indicates p < 0.001.

RT-qPCR Validation

The expression of ANKRD36, HLTF, KIF3A, OXCT1 and VPS13A were validated by RT-qPCR. As shown in Figure S3A, compared with control group, the expression of ANKRD36, HLTF, KIF3A, OXCT1 and VPS13A were reduced in cluster A and cluster B. Consistent with the results in Figure S3A, all five genes were down-regulated in the MI group compared with the control group (Figure 8).

Figure 8 Expression validations of ANKRD36, HLTF, KIF3A, OXCT1 and VPS13A by RT-qPCR. (A) ANKRD36, (B) HLTF, (C) KIF3A, (D) OXCT1, (E) VPS13A. *p < 0.05.

Discussion

The underlying pathological mechanism of MI remains unclear, and needs to be clarified urgently. In this study, we identified two different subtypes of MI, and the two subtypes showed different expression patterns and hypoxic-immune status. The diagnostic model, which can well distinguish the two subtypes of MI, was developed and validated.

Consensus clustering divided MI samples into two subtypes, each with its own unique immune characteristics. The results of ssGSEA analysis showed that more immune cell infiltration and more active immune responses were detected in cluster A, while pronouncedly higher immunoscore in cluster B was detected. In addition, the expression levels of several HLA-related genes were increased in cluster A. Hypoxia is a central pathophysiological component in MI, which can promote the activation of epithelial-to-mesenchymal transition (EMT) regulatory genes.⁶ The epicardium is a potential source of cardiac progenitors to support reparative angiogenesis after MI through EMT.¹⁰ Hence, we compared the differences of hypoxia and EMT status between two subtypes and found that the enrichment scores of hypoxia and EMT2 were significantly higher in cluster A. All the above results indicated that compared with cluster B, cluster A was characterized by high hypoxia status and low immunity status.

As for potential drugs, 11 genes were observed as targets of the 60 predicted drugs, such as testosterone, pinocembrin and cilostazol. Previous studies had explored the relationship between testosterone and cardiovascular disease. A review suggested that testosterone may perform a critical function in the treatment of angina, heart failure, and myocardial ischemia.¹¹ Low endogenous bioavailable testosterone levels had been associated with higher rates of cardiovascular-related mortality.¹² A Mendelian randomisation study indicated that endogenous testosterone was positively associated with thromboembolism, heart failure, and MI in men.¹³ Lungkaphin et al reported that pinocembrin can reduce cardiac arrhythmia and infarct size in acute myocardial ischemia/reperfusion rats.¹⁴ Similarly, Zhang et al indicated that pinocembrin can decrease cardiac infarction area.¹⁵ Ye et al suggested that pinocembrin can ameliorate atrial electrical remodeling, attenuate autonomic dysfunction, and suppress MI-induced inflammatory responses in a MI rat model.¹⁶ Tanaka et al suggested that cilostazol might be used for temporary adjunctive treatment of ST-segment elevation MI.¹⁷ Lee et al reported that cilostazol has protective effects on myocardial function and myocardial remodeling in MI with restrained stress rat model.¹⁸ The above findings may suggest a potential novel strategy for MI treatment.

A total of five genes, including ANKRD36, HLTF, KIF3A, OXCT1 and VPS13A, were included in the diagnostic model, and showed good diagnostic effect in distinguishing the two subtypes. Yan et al found that Ankyrin repeat domain 36 (ANKRD36) was significantly lower expressed in hypertension.¹⁹ An integrated genomic analysis indicated that ANKRD36 was a novel biomarker of early chronic myeloid leukemia progression.²⁰ Helicase-like transcription factor (HLTF) encodes a member of the SWI/SNF family, which has helicase and ATPase activities. Liu et al indicated that HLTF is associated with the migration and invasion of colorectal cancer.²¹ Kinesin family member 3A (KIF3A) overexpression was related to tumor recurrence, lymph node metastasis and chemoresistance.²² The gene, 3-oxoacid CoA-transferase 1 (OXCT1) has been associated with tumor proliferation and metastasis. A review showed that the elevation of OXCT1 expression was detected in different categories of cancer cells and considered OXCT1 as a novel target for cancer therapy.²³ A bioinformatics analysis on pan-cancer suggested that VPS13A is associated with the tumor development.²⁴ Elevated VPS13A might contribute to increased risk of stroke in current smokers.²⁵ In addition, the VPS13A and OXCT1 are the most correlated hub genes, indicating that they may function together in MI. Our findings indicated that the model that included five genes had a better diagnostic effect than each individual gene.

The results of correlation analysis of hub genes and immune cells showed that VPS13A had the strongest positive correlation with activated CD4 T cell and strongest negative correlation with monocyte, indicating that the expression of VPS13A was associated with activated CD4 T cell and monocyte. Inflammation has been reported to be involved in the initiation and progression of numerous diseases, including MI.²⁶ Regulatory T cells (Tregs) can prevent the accumulation of proinflammatory cells and the secretion of proinflammatory cytokines, and Tregs cell activity might be pivotal for the maintenance of cardiovascular homeostasis.²⁷ In addition, Zhao et al indicated that activated CD4+ T cells promote cardiac remodeling following MI.²⁸ Monocytes play an integral role in the innate immune system in human and participate in inflammation.²⁹ Previous study suggested that monocytes are involved in both initiation and resolution of inflammation and are critical mediators of the inflammatory response following MI and ischemia-reperfusion injury.³⁰

Conclusion

In conclusion, we divided MI samples into two subtypes, cluster A and cluster B, based on the expression levels of HRGs. Among them, cluster A was characterized by high hypoxia status and low immunity status. In addition, some drugs (testosterone, pinomicin, cilostazol, etc.) have been found to be effective in the treatment of MI but may produce different results in the two different subtypes. Finally, a diagnostic model, including ANKRD36, HLTF, KIF3A, OXCT1 and VPS13A, was established, which can distinguish the two subtypes well.