Biomarkers identification for acute myocardial infarction detection via weighted gene co-expression network analysis.

The study aimed to seek potential biomarkers for acute myocardial infarction (AMI) detection and treatment.The dataset GSE48060 was used, consisting of 52 peripheral blood samples (31 AMI samples and 21 normal controls). By limma package, differentially expressed genes (DEGs) between 2 kinds of samples were identified, followed by enrichment analysis, subpathway analysis, protein-protein interaction (PPI) network analysis, and transcription factor network (TFN) analysis. Weighted gene co-expression network analysis was used to further extract key modules relating to AMI, followed by enrichment and TFN analyses. Expression validation was performed via meta-analysis of 2 datasets, GSE22229 and GSE29111.A set of 428 DEGs in AMI were screened out, and the upregulated toll-like receptor (TLR) family genes (TLR1, TLR2, and TLR10) were enriched in wound response, immune response and inflammatory response functions, and downregulated genes (GBP5, CXCL5, GZMA, CCL5, and CCL4) were correlated with immune response. CCL5, GZMA, GZMB, TLR2, and formyl peptide receptor 1 (FPR1) were predicted as crucial nodes in the PPI network. Signal transducer and activator of transcription 1 (STAT1) was the key transcription factor (TF) with multiple targets. The grey module was highly related to AMI. Genes in this module were closely related to regulation of macrophage activation, and spermatogenic leucine zipper 1 (SPZ1) was identified as a TF. Expressions of TLR2 and FPR1 were confirmed via the integrated matrix.Several potential biomarkers for AMI detection were identified, such as GZMB, GBP5, FPR1, TLR2, STAT1, and SPZ1. They might exert their functions via regulation of immune and inflammation responses. Genes in grey module play significant roles in AMI via regulation of macrophage activation.

Introduction

Acute myocardial infarction (AMI) is a major cause for death in the world.[ Although the incidence of myocardial infarction (MI) has a decreased trend between 1997 and 2008 in the United State,[ approximate 155,000 new asymptomatic silent cases occur annually,[ and in developing countries, AMI is still a significant health burden.[ Many advanced approaches have been developed for the management of patients with AMI, such as pharmacological reperfusion therapy and catheter-based interventions. However, AMI remains a major problem worldwide.[

Early detection of AMI contributes to early treatment interventions, which can significantly reduce the mortality rate.[ Many studies have been conducted to seek potential molecular biomarkers for AMI detection. As expected, several genes and proteins such as monocyte-platelet aggregates, heart fatty acid–binding protein, and troponin I have been established as potent markers for AMI diagnosis.[ Moreover, substantial miRNAs have been proposed as crucial markers for management of AMI, such as circulating miR-1, miR-19a, miR-208a (cardiac-specific), and miR-499.[ However, interactions between these molecules and the potentially altered functions and pathways are rarely reported. The pathogenic mechanisms of AMI pathogenesis remain obscure. Suresh et al[ establish a gene expression profiling by microarray, and compare the differentially expression genes (DEGs) between AMI samples and normal samples, and the dysregulated pathways of these DEGs. However, they only focus on the pathways and recurrent events. The regulatory correlations between these genes are not further investigated. Thus, a recent study by Gao et al[ reanalyzes this GSE48060 microarray data, and performs protein–protein interaction (PPI) network analysis and transcription factor network (TFN) analysis after identifying DEGs in AMI samples. Although they predict several crucial genes in AMI progression, such as CCL5, BCL3, and NCOA7, the subpathway analysis and coexpression network analysis are not involved.

Weighted Gene Co-expression Network Analysis (WGCNA) is a useful approach that has been widely used for gene expressions to identify key disease-related modules.[ For instance, using WGCNA, Saris et al[ select two large co-expressed modules that related to amyotrophic lateral sclerosis. By estimating gene expression patterns, Azuaje et al[ identify a WGCN in MI, which is used to determine the potential role of Col5a2 and its transcriptional pattern. Therefore, we reanalyzed Suresh's microarray data, GSE48060, and carried out WGCNA for gene expressions, in addition to enrichment analysis, PPI network analysis, and TFN analysis of the DEGs. Moreover, in the aspects of the DEG screening, we applied more rigorous criteria with fold change (FC) >1.2 and P < 0.05, compared with Gao et al[ with P < 0.01. Furthermore, we performed expression validation of DEGs by using other databases, which could powerfully verify our results from other experimental data. We aimed to reveal molecular underpinnings in AMI development and provide novel biomarkers for the AMI detection and treatment.

Methods

Data resource

The dataset GSE48060,[ which consisted of 52 microarray expression profiles, was downloaded from the Gene Expression Omnibus (GEO, http://www.ncbi. nlm.nih.gov/geo) database. The gene expression data were derived from the peripheral blood tissues of 31 first-time AMI patients within 48 hours after MI with (n = 5) or without (n = 26) recurrent events (AMI samples) and 21 normal cardiac function controls (control samples). The platform of the dataset was GPL570 (HG-U133_Plus_2, Affymetrix Human Genome U133). The samples were recruited from the Mayo Clinic Rochester echocardiography laboratory, and all participants signed informed consent approved by the Mayo Clinic Rochester Institutional Review Board.

Data pretreatment and differentially expressed genes selection

Firstly, the series matrix profile that has been undergone pretreatments such as background correction, quantile normalization, and probe summarization was downloaded. Then, the probes did not correspond to gene symbols were eliminated. When multiple probes corresponded to 1 gene symbol, the expression value of gene was adopted to their mean expression values. Thereafter, the Linear Models for Microarray Analysis (limma, http://www.bioconductor.org/packages/3.0/bioc/html/limma.html) package of R was used to identify DEGs between AMI and control samples.[ The thresholds for DEG screening were the FC > 1.2 (|log2 FC| > 0.263) and P < .05.

Functional and pathway enrichment analyses, protein–protein interaction network analysis and transcription factors identification

Based on the Database for Annotation, Visualization and Integration Discovery (DAVID, http://david.abcc.Ncifcrf.gov/) online tool,[ functional and pathway enrichment analyses for the DEGs were conducted. The significant analysis of the enriched go terms and pathways were conducted by using hypergeometric test. The cut-off value for functional enrichment analysis was P < .05, and that for pathway enrichment analysis was adjusted P < .05.

Combining information in the Search Tool for the Retrieval of Interacting Genes (STRING, http://string-db.org/) database with DEG expressions, we explored potential interactions of these DEGs at protein level.[ PPI score was set as 0.4 to construct a PPI network, which was then visualized by the Cytoscape (http://cytoscape.org/) software.[ In the network, a node represents a protein encoded by corresponding DEG, and the degree of a node is deemed as the number of interactions with other DEGs. A node with high degree serves as a hub gene in the PPI network.

The Human Transcription Regulation Interaction (HTRI, http://www.lbbc.ibb.unesp.br/htri) is an open-access database that contains experimentally validated human transcriptional interactions.[ Based on information in HTRI database, we screened out TFs amongst these DEGs, and identified their downstream targets to build the TF-target regulatory network.

Subpathway analysis

To further elucidate potential pathways of the DEGs, subpathway analysis was carried out for up- and downregulated DEGs, respectively, using the isubpathwayMiner package in R.[ The cut-off value for significant subpathway was P < .05.

Construction of weighted gene co-expression network

WGCNA has been widely applied to establish a scale-free network from gene expression data.[ The WGCNA package in R (http://www.inside-r.org/packages/cran/WGCNA/docs/bicor) was used to construct this network, and to get more genes for WGCNA analysis, the genes with the looser threshold value of P < .05 from results of difference expression analysis between AMI samples and normal samples in GSE48060 were selected. There were 3 steps to construct the WGCNA network: definition of the gene co-expression matrix. In brief, the Pearson correlation coefficient (CC) was calculated for all gene interactions. The correlation between 2 genes, m and n, was deemed as Smn = |cor(m,n)|, based on which the gene co-expressed matrix was built as S = [Smn]; definition of adjacent function. The adjacency coefficient amn was used to determine the correlations of gene interactions, based on the function amn = power (Smn, β)=|Smn|β; determination of parameter in the adjacent function. The weighting coefficient β was determined according to the scale-free principle. Herein, i represents numbers of nodes, and p(i) represents the probability of a node. For WGCNA network construction, it was required that the CC between log (i) and log [p(i)] was at least 0.8.

After the network establishment, the hierarchical average linkage clustering was used to divide different gene co-expression modules.[ Then the correlations between genes in a module and a disease were recognized according to the following methods:P value of a DEG between 2 kinds of samples was calculated; the log P was defined as the gene significance (GS), and the average GS in a module was deemed as module significance (MS). Commonly, a module with a higher MS indicates a greater possibility to be associated with the disease.

Expression validation of differentially expressed genes via other datasets

The expression data of 9 blood tissue samples from healthy individuals in the expression dataset GSE22229 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE22229) and 36 blood tissue samples from patients with AMI in the dataset GSE29111 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE29111) were downloaded from the GEO database. The platforms of the 2 datasets were both GPL570 (HG-U133_Plus_2, Affymetrix Human Genome U133).

CONOR package in R, which is available through the CRAN repository (https://cran.r-project.org/web/packages/),[ was used to perform meta-analysis for the 2 datasets. Then data were normalized using median method, and an integrated expression matrix was obtained as the validation data. Next, limma package was used to select the DEGs between 2 kinds of samples, with the same criteria as aforementioned (|log2 FC| > 0.263 and P < .05). We checked whether DEGs identified in GSE48060 were also differentially expressed in this integrated expression matrix.

Results

Identification of differentially expressed genes between acute myocardial infarction and normal samples

Based on the aforementioned criteria, we obtained a total of 428 DEGs, including 160 upregulated and 268 downregulated genes. Heat map of clustering analysis of gene expression was listed in Figure 1, in which these DEGs could well distinguish the 2 kinds of samples, indicating the DEGs could be used for further analysis.

Figure 1

Heat map of clustering analysis of the gene expressions. X-axis represents the names of each sample in dataset GSE48060, and Y-axis represents the clusters of genes. The color toward red represents high expression values, and color toward green represents lower expression values.

Enrichment and subpathway results

As a result, we found that the upregulated genes were mainly associated with functions such as response to wounding [e.g., toll-like receptor 10 (TLR10), CR1, and TLR2), immune response (e.g., IL1R2, TLR10, and TLR2), and inflammatory response (e.g., TLR10, TLR1, and TLR2). Although the downregulated genes were significantly enriched in immune response (e.g., GBP5, CXCL5, and GZMA), cellular defense response (e.g., KLRC3, CX3CR1, and CCL5), and defense response (e.g., IL15, CCL5, and CCL4) functions (Table 1). In addition, pathways analysis for all the DEGs were significantly associated with IL12-mediated signaling events (e.g., CCL4, GZMA, and GZMB), endogenous TLR signaling (e.g., CD14, TLR1, and TLR2), interferon gamma signaling (e.g., FCGR1B, GBP1, and GBP3), and cytokine–cytokine receptor interaction (e.g., CCL4, CCL5, and CSF2RB) related pathways (Table 2).

Table 1

Significant function of the differentially expressed genes (top 5).

Table 2

Significant pathways of the differentially expressed genes.

Moreover, 32 subpathways in the 8 pathways were identified for upregulated DEGs, such as starch and sucrose metabolism (e.g., HK3, PYGL, and MGAM), pantothenate and CoA biosynthesis (e.g., VNN2 and VNN1), and pentose phosphate pathway (e.g., PGD and TKT). For the downregulated DEGs, only 1 subpathway as retinol metabolism was identified (ALDH1A1) (Table 3).

Table 3

Subpathway analysis results (partly).

Protein–protein interaction network

According to the STRING database, a PPI network of the DEGs was constructed, containing 205 nodes and 590 interactions. In the network, formyl peptide receptor 1 (FPR1) (degree = 34), CCL5 (degree = 32), GZMA (degree = 32), TLR2 (degree = 31), and GZMB (degree = 30) were the highlighted nodes (Fig. 2).

Figure 2

Protein–protein interaction network of the DEGs. The nodes denote the protein product of DEGs, and the lines between 2 proteins denote their correlations. In addition, the grey square-shaped nodes represent downregulated proteins of DEGs, and the white square-shaped nodes represent upregulated proteins of DEGs. The size of each node indicates the number of the interactions with other DEGs, and the boldness of the line is directly proportional to combined score between the proteins. DEG = differentially expressed gene.

Transcription factors and the corresponding targets

Based on HTRI database, a total of 6 TFs that regulated 42 downstream targets were searched. Among them, signal transducer and activator of transcription 1 (STAT1) was the predominant one with multiple targets (Fig. 3).

Figure 3

Transcription factor-target regulatory network. Square shape represents target and rhombus represents transcription factors. Arrows between them predict their correlations. In addition, the grey square shape indicates the downregulated gene, and the white square shape indicates the upregulated gene, whereas the black square shape indicates gene with no differential expression.

Weighted gene co-expression network analysis network and key module selection

A set of 2892 genes (P < .05) from results of difference expression analysis between AMI samples and normal samples in GSE48060 were screened out to establish the WGCNA network. As revealed in the hierarchical clustering tree, there identified a total of 13 modules (Fig. 4), and their correlations with AMI was listed in Table 4. As a result, the grey module has the highest correlation with AMI (CC = 0.88). Therefore, we further analyzed genes in this module, and found that they were significantly related to functions such as multicellular organismal response to stress, regulation of macrophage activation, and spinal cord motor neuron differentiation (Table 5). In addition, spermatogenic leucine zipper 1 (SPZ1) was identified as only 1 TF in the grey module.

Figure 4

Branches of the cluster dendrogram of the most connected genes in 13 modules. A line in the cluster dendrogram represents one gene, and the height in Y-axis indicates the distance between two genes. Different colors represent different modules.

Table 4

Correlations between each module and acute myocardial infaction.

Table 5

Significant pathways of genes in the grey module.

Expression validation of differentially expressed genes by other datasets

As a result, there identified a total of 6249 DEGs from the validation dataset (GSE22229 and GSE29111), including 3409 upregulated and 2840 downregulated genes. Compared with DEGs in the dataset of GSE48060, we found 147 overlapped DEGs, including 53 upregulated and 10 downregulated genes, were both in the GSE48060 dataset and validation dataset. Especially, we focused on the expression of several genes which were predicted as important nodes in the PPI network, such as FPR1, CCL5, GZMA, TLR2, and GZMB. Interestingly, we found that FPR1 and TLR2 both belonged to the overlapped DEGs (Table 6).

Table 6

Overlapped differentially expressed genes in GSE48060 and the integrated matrix.

Discussion

Via a series of bioinformatics methods, we identified a total of 428 DEGs in AMI, which were mainly enriched in IL12-mediated signaling events (e.g., GZMA and GZMB), endogenous TLR signaling (e.g., TLR1 and TLR2) and cytokine–cytokine receptor interaction (e.g., CCL4 and CCL5)-related pathways. Upregulated TLR family genes (TLR1, TLR2, and TLR10) were significantly associated with wound response, immune response and inflammatory response functions, and downregulated genes (GBP5, CXCL5, GZMA, CCL5, and CCL4) were also correlated with immune response. CCL5, GZMA, GZMB, FPR1, and TLR2 were predicted as crucial nodes in the PPI network. STAT1 was the key TF with multiple targets. By WGCNA, we found the grey module was highly related to AMI, and genes in this module were closely related to regulation of macrophage activation. In addition, SPZ1 was identified as a TF in the grey module. Notably, FPR1 and TLR2 as DEGs in AMI samples were validated in other datasets.

The Granzyme 2, Cytotoxic T-Lymphocyte-Associated Serine Esterase 1 (GZMB) encoded protein is an activate enzyme of caspases, which are involved in the signaling pathways of necrosis and inflammation. Reportedly, GZMB is one of the DEGs enriched in the IL12-mediated signaling events pathway, in the comparison of inflammatory with noninflammatory synovial biopsies.[ SOCS3 is a vital protein that regulates the immediate response to cytokines. Through the inhibition of IL12-mediated signaling, SOCS3 plays an important role in Th2 development.[CCL5 is one of the chemokine genes involved in inflammatory process.[GBP5 is an activator of NLRP3 inflammasome assembly that has significant roles in inflammation.[ In our present study, we found that GZMB, GBP5, and CCL5 were linked in the PPI network, and GZMB was enriched in IL12-mediated signaling events pathway, whereas GBP5 and CCL5 in immune response. Inflammation reaction is an important reaction of immune system,[ and immune response is often accompanied by inflammation response during the disease progression. Thus, we speculated that IL12-mediated signaling events pathway might associate with inflammation to cause AMI development, and the above 3 genes might work together in regulating these function and pathways.

The protein encoded by TLR2 is a member of TLR family which has a critical role in activating innate immunity. AMI is reported to have close relationship with increased expression of TLR2.[ During AMI, the high-mobility group box 1 protein can induce myocardial repair, and it can activate NF-κB signaling pathway through the receptors such as TLR2/4/9.[ A study suggests that sTLR2 might involve in the innate immunity response in pathogenesis of heart failure after AMI.[ The inflammatory response is reduced by the activation of TLR, and this may badly influence the reperfusion therapy of MI. Arslan et al[ find that TLR2 regulates myocardial ischemia, and anti-TLR2 antibody causes a pronounced reduction of leukocyte influx and cytokine production. In our study, we discovered TLR family genes (TLR2, TLR10, and TLR1) were crucial DEGs in AMI that enriched in immune response and inflammatory response, suggesting these genes might be closely associated with development of AMI, via involvement of these responses. Moreover, TLR2 was predicted as a hub node in the PPI network. Notably, differential expression of TLR2 in AMI samples was validated by integrating analysis another 2 datasets, and this provided potent evidence that it is an important gene for AMI management and could be used as a therapeutic target.

The FPR1 is one of the FPRs that participant in inflammation reaction. Reportedly, it can protect against myocardial ischemia-reperfusion.[ In adult rat cardiomyocytes, it is found activation of FPR1 contributes to the Ac-ANX-A12–26 cardioprotective actions.[ Notably, FPR1 is also identified as a DEG in human AMI blood tissue, compared with normal blood tissue using microarray data.[ In addition, activation of FPR2, the family member of FPR1, is involved in cardiac repair after MI via mobilization of circulating angiogenic cells.[ Therefore, the deregulated expression of FPRs might account for the occurrence of AMI. In our study, FPR1 was predicted as a crucial node in the PPI network. More importantly, its altered expression was validated by integrating analysis another 2 datasets, which further confirmed our prediction that it was an important gene for AMI prevention.

Activation of the transcription STAT1 is elevated in primary cardiac myocytes under exposure of simulated ischemia. Deficiency of STAT1 exerts a protective function against MI through the control of autophagy.[ Absent in melanoma 2 (AIM2) is a protein involved in inflammatory regulation. Reportedly, AIM2 limits the transcription of proinflammatory cytokine in cardiomyocytes via inhibition of STAT1 phosphorylation.[ In the present study, STAT1 was identified as a critical TF that targeted multiples genes, suggesting its central role in AMI pathology.

Macrophages take part in the progression of various inflammatory-related diseases. The inhibitor of COX-2, whose expression is elevated by macrophage under exposure to proinflammatory stimuli, can increase risk of AMI in patients without cardiovascular risk factors.[ In mice, deletion of matrix metalloproteinase-28 inhibit M2 macrophage activation, and may aggravate left ventricle dysfunction and rupture after MI.[ Our findings by WGCNA indicated that the grey module was highly related to AMI and genes in this module were significantly enriched in regulation of macrophage activation, suggesting genes in this module were highly related to AMI, via regulation of macrophage activation.

SPZ1, a bHLH-zip TF, acts as a downstream gene of mitogen-activated protein kinase (MAPKs) and is important in MAPKs signaling pathway via phosphorylated by MAPK1/ERK2 and MAPK3/ERK1.[ Reportedly, brain MAPK1/ERK2 signaling pathways are activated in AMI rats.[ In our result, SPZ1 was only 1 TF with significant P value in the grey module and we predicted this is a novel TF might be related to AMI due to no direct relationships between its expression alteration and AMI has been reported.

Despite these valuable findings, there remains a limitation that we did not perform experimental validations in vitro. However, we should note that the sample collections are very hard because they are from blood tissues in patients with AMI. The eligible samples from patients with AMI would be collected in our future studies as much as possible.

Our findings provide important values in prediction of molecular events related to AMI as well as potential biomarkers for detection and prevention. However, the results need to be validated by substantial experiments, as well as the performance of those potential biomarkers such as sensitivity, superiority to the current, cost, and convenience should also be analyzed in our future study. In addition, in our future study, prediction of miRNAs that regulate the DEGs in PPI network, as well as submodule analysis and transcriptional regulatory network construction for the genes in the grey module are still needed to potentially provide a more comprehensive insight into AMI.

In conclusion, several novel potential biomarkers for AMI detection were identified, such as GZMB, GBP5, TLR2, FPR1, STAT1, and SPZ1. They might exert their functions via regulation of immune and inflammation responses. The grey module was most related to AMI, and genes in this module play significant roles via regulation of macrophage activation. However, future experimental verification is still needed to confirm those speculations.