Identification of vital modules and genes associated with heart failure based on weighted gene coexpression network analysis.

AIMS: Heart failure (HF) is a chronic heart disease with a high incidence and mortality. Due to the regulatory complexity of gene coexpression networks, the underlying hub genes regulation in HF remain incompletely appreciated. We aimed to explore potential key modules and genes for HF using weighted gene coexpression network analysis (WGCNA). METHODS AND
RESULTS: The expression profiles by high throughput sequencing of heart tissues samples from HF and non-HF samples were obtained from the Gene Expression Omnibus database. Differentially expressed genes (DEGs) between HF and non-HF samples were firstly identified. Then, a coexpression network was constructed to identify key modules and potential hub genes. The biological functions of potential hub genes were analysed by Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analyses. Finally, a protein-protein interaction (PPI) network was constructed using the STRING online tool. A total of 135 DEGs (133 up-regulated and 2 down-regulated DEGs) between HF and non-HF samples were identified in the GSE135055 and GSE123976 datasets. Moreover, a total of 38 modules were screened based on WGCNA in the GSE135055 dataset, and six potential hub genes (UCK2, ASB1, CCNI, CUX1, IRX6, and STX16) were screened from the key module by setting the gene significance over 0.2 and the module membership over 0.8. Furthermore, 78 potential hub genes were obtained by taking the intersection of the 135 DEGs and all genes in the key module, and enrichment analysis revealed that they were mainly involved in the MAPK and PI3K-AKT signalling pathways. Finally, in a PPI network constructed with the 78 potential hub genes, CUX1 and ASB1 were identified as hub genes in HF because they were also identified as potential hub genes in the WGCNA.
CONCLUSIONS: To the best of our knowledge, our study is the first to employ WGCNA to identify the key module and hub genes for HF. Our study identified a module and two genes that might play important roles in HF, which may provide potential biomarkers for the diagnosis of HF and improve our knowledge of the molecular mechanisms underlying HF.

Introduction

Heart failure (HF) is the most common heart disease worldwide and affects millions of people. In 2016, there were approximately 300 000 deaths associated with HF. HF is an important health problem and causes severe medical and economic burdens worldwide.

Heart failure is a chronic disease and a complication of many other diseases, such as type 2 diabetes mellitus and coronary heart disease. The occurrence and development of HF can be driven by a variety of biological factors, such as myocardial stretch, oxidative stress, inflammatory response, and neuroendocrine system activation. HF involves many pathophysiological mechanisms, accompanied by many molecules entering body fluids and thus becoming biomarkers for diagnosis. Brain natriuretic peptide and N‐terminal proBNP are the two molecular markers considered by the European Society of Cardiology for the clinical diagnosis of HF. In recent years, some new molecular markers have also been used in the diagnosis, classification, and efficacy monitoring of HF, such as the inflammatory mediators IL‐1 and IL‐6, Gal‐3, myocardial remodelling marker matrix metalloproteinases, and various non‐coding RNAs that can stably exist in plasma, such as miR‐133a and miR‐29a. Identifying biomarkers that play key roles in these pathways is of great significance for the early diagnosis of HF, optimizing accurate treatment strategies for patients with HF, and reducing its mortality.

Epigenetic regulation plays an important role in the development of HF. Epigenetic regulation involves RNA, specifically the roles of short and long non‐coding RNAs and endogenous competing RNA regulatory networks, and DNA, including DNA methylation, histone modifications, and chromatin conformational changes. In addition to miRNAs, transcription factors are key regulators of gene expression that enhance or inhibit gene transcription by binding to specific DNA sequences of target gene promoters. If transcription factors and miRNAs jointly regulate a coding gene, and in turn, transcription factors regulate the coding gene of this miRNA, a feed‐forward loop is established between transcription factors, miRNAs, and the coding genes. Due to the regulatory complexity of gene coexpression networks in HF, it is of great importance to explore molecular biomarkers for the diagnosis and treatment of HF.

Weighted gene coexpression network analysis (WGCNA) is a systems biology method for analysing the correlations among genes and the relationships between modules and external sample traits. WGCNA allows the discovery of correlations between gene clusters and clinical traits and between genes and coexpressed modules or clinical traits. In the current study, we hypothesized that gene coexpression networks regulated the hub genes and pathways involved in HF. The results of this study provide potential biomarkers for the diagnosis of HF and improve our knowledge of the molecular mechanisms underlying HF.

Methods

Data collection

Gene expression matrixes of GSE135055 and GSE123976 were downloaded from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo) based on the GPL16791 platform (Illumina HiSeq 2500, ). The GSE135055 dataset consisted of 21 human HF and 9 human non‐HF (NF) left ventricle tissue samples. The GSE123976 dataset consisted of 6 human HF and 3 human NF left ventricle tissue samples. The clinical characteristics of the samples were summarized in Supporting Information, Table . The step‐by‐step flowchart of this research was shown in Figure .

Figure 1

The flow chart of the study. DEGs, differentially expressed genes; HF, heart failure; NF, non‐HF; PPI, protein–protein interaction; WGCNA, weighted gene coexpression network analysis.

Identification of differentially expressed genes

The differentially expressed genes (DEGs) between HF and NF samples in each dataset were screened using the limma package of R based on a cut‐off value|log2 (fold change)| > 0.05 and an adjusted P < 0.05.

Construction of the weighted gene coexpression network analysis network

Weighted gene coexpression network analysis was performed using the WGCNA package of R to identify potential key modules and potential hub genes. First, all samples were clustered to screen obvious outliers. Second, a soft‐thresholding power β was selected, and a scale‐free network was constructed. Third, hierarchical clustering and the dynamic tree cut function were used to detect modules. Fourth, the correlations between modules and clinical information were analysed to identify key modules, and potential hub genes in the key module were obtained by setting the gene significance over 0.2 and the module membership over 0.8.

Gene ontology and Kyoto Encyclopedia of Genes Genomes pathway enrichment analyses

Potential hub genes were identified by taking the intersection of DEGs and genes in the key module. To explore the biological functions and signalling pathways involved in potential hub genes, gene ontology (GO) analysis and Kyoto Encyclopedia of Genes Genomes (KEGG) pathway analysis were performed using the clusterProfiler package of R. The biological functions of the genes can be divided into three categories: biological process, molecular function, and cellular component. Bar plots and bubble plots were created using the ggplot package of R to visualize the enrichment results. Enriched results with a P value < 0.05 were considered significant.

Construction of the protein–protein interaction network

The protein–protein interaction (PPI) network was constructed using the Search Tool for the Retrieval of Interacting Genes online tool (STRING, https://www.string‐db.org) to identify hub genes, and the results were visualized and analysed with Cytoscape.

Results

Identification of differentially expressed genes in the GSE135055 and GSE123976 datasets

FPKM format expression profiles of the GSE135055 and GSE123976 datasets were obtained from the GEO database. A total of 10 316 DEGs were identified between the HF (n = 6) and NF (n = 3) samples in the GSE123976 dataset, including 10 247 up‐regulated and 69 down‐regulated DEGs in HF relative to NF (Figure ). A total of 314 DEGs were identified between the HF (n = 21) and NF (n = 9) samples in the GSE135055 dataset, including 203 up‐regulated and 111 down‐regulated DEGs in HF relative to NF (Figure ). Heatmaps of the top 50 DEGs (ranked by adjusted P value) in the GSE123976 and GSE135055 datasets are shown in Figure , respectively. In addition, a total of 135 overlapping DEGs between the GSE135055 and GSE123976 datasets were identified. Among them, 133 were up‐regulated, and 2 were down‐regulated in HF (Figure ).

Figure 2

Identification of DEGs between the HF and NF groups in the GSE135055 and GSE123976 datasets. (A) Volcano plot showing that a total of 10 316 DEGs between HF (n = 6) and NF (n = 3) samples were identified in the GSE123976 dataset, of which 10 247 were up‐regulated and 69 were down‐regulated in HF. Blue represents down‐regulated genes, red represents up‐regulated genes, and grey represents genes with no significant difference in expression. (B) The top 100 DEGs of the GSE123976 dataset were visualized in a heatmap. (C) Volcano plot showing that a total of 314 DEGs between HF (n = 21) and NF (n = 9) samples were identified in the GSE135055 dataset, of which 203 were up‐regulated and 111 were down‐regulated in HF. Blue represents down‐regulated genes, red represents up‐regulated genes, and grey represents genes with no significant difference in expression. (D) The top 100 DEGs of the GSE135055 dataset were visualized in a heatmap. (E) Venn diagram showing the intersections of up‐regulated DEGs between the GSE123976 and GSE135055 datasets. (F) Venn diagram showing the intersections of down‐regulated DEGs between the GSE123976 and GSE135055 datasets. DEGs, differentially expressed genes; HF, heart failure; NF, non‐HF; PPI, protein–protein interaction; WGCNA, weighted gene coexpression network analysis.

Construction of the gene coexpression network by weighted gene coexpression network analysis

Weighted gene coexpression network analysis was carried out on the expression profile of the GSE135055 dataset. Sample clustering was performed, and the clustering dendrogram revealed no obvious outliers (Figure ); therefore, all samples were included in further analysis. Then, the soft threshold power β was calculated before the construction of the gene coexpression network. According to the results, the network met the scale‐free topology threshold of 0.9 (Figure ), and the mean connectivity was close to 0 (Figure ) when the soft‐thresholding power was 5; therefore, β = 5 was selected to construct a hierarchical clustering tree. Ultimately, 38 modules were identified based on average hierarchical clustering and dynamic tree clipping (Figure ).

Figure 3

Construction of the gene coexpression network for the GSE135055 dataset. (A) Sample clustering of the GSE135055 dataset revealed no obvious outliers. (B) Network topology analysis revealed that the network met the scale‐free topology threshold of 0.9 when β = 5. The x‐axis represents soft‐thresholding power, and the y‐axis represents the scale‐free topology model fit index. (C) Network topology analysis revealed that the mean connectivity was close to 0 when β = 5. The x‐axis represents soft‐thresholding power, and the y‐axis represents the mean connectivity. (D) Clustering dendrogram of DEGs in the GSE135055 dataset based on topological overlap. Each module is given a unique colour and represents a cluster of coexpressed genes.

Identification of key modules associated with heart failure progression

The relationships between each module and all clinical variables, including status (HF vs. NF), age, and sex, in the GSE135055 dataset were investigated. The results showed that among the modules, the yellow module had the strongest association with HF (r = 0.9, P = 9e‐12) (Figure ), indicating that the genes of the yellow module may play important roles in the occurrence and development of HF. Thus, the yellow module was considered as a key module for further analysis.

Figure 4

Identification of the key modules associated with the development of heart failure. Heatmap displaying the correlations and significant differences between gene modules and clinical information. Correlations are displayed in rectangles, and correlation coefficients and P values are displayed in rectangles. Red rectangles denote negative correlations between modules and clinical information, and green rectangles denote positive correlations between modules and clinical information.

Identification of the hub genes in the key module

Pearson's correlation analysis of gene significance vs. module membership in the yellow module showed that they were highly correlated (r = 0.85, P < 1e‐200) (Figure ). To identify the potential hub genes in the yellow module, we selected a threshold of 0.2 for the correlation between genes and clinical variables and a threshold of 0.8 for the correlation between genes and the yellow module. A total of six hub genes were screened: UCK2, ASB1, CCNI, CUX1, IRX6, and STX16. A heatmap of these hub genes was constructed based on the topology overlap (Figure ). The results revealed that STX16 had strong associations with UCK2, IRX6, CCNI and ASB1, and IRX6 was strongly correlated with CCNI.

Figure 5

Hub genes in the yellow module strongly correlated with the progression of heart failure. (A) Scatterplot of gene significance vs. module membership for the progression of heart failure in the yellow module (r = 0.85, P < 1e‐200). (B) A heatmap showing the associations between hub genes based on topology overlap. Dark colours indicate a high degree of topology overlap and strong correlations between the hub genes.

Biological function analysis of potential hub genes

To explore the function of potential hub genes in the yellow module, the intersection of the 135 DEGs identified in the GSE135055 and GSE123976 datasets and the 1877 genes in the key module was obtained. A Venn diagram was constructed to visualize the intersecting genes, and 78 intersecting genes were identified, which were up‐regulated DEGs in the GSE135055 and GSE123976 datasets (Figure ). Subsequently, GO and KEGG pathway enrichment analyses were performed to explore the potential functions of the 78 intersecting genes. The results showed that the 78 intersecting genes were mainly enriched in axonogenesis, regulation of cell shape, Ras protein signal transduction, and regulation of cell morphogenesis in biological process (Figure ) and were mainly located in focal adhesion, microtubules, myofibrils, and synaptic clefts (Figure ). The KEGG pathway results showed that the intersecting genes were significantly enriched in the MAPK signalling pathway, PI3K‐AKT signalling pathway, and ARG‐RAGE signalling pathway in diabetic complications (Figure ), which indicated that these intersecting genes may exert potential roles in the development of HF via these pathways.

Figure 6

Biological function analysis of intersecting genes. (A) Venn diagram showing the intersecting genes between the yellow module and the differentially expressed genes, and 78 differentially expressed up‐regulated genes were screened. (B) Bubble plots of the biological process (BP) enrichment for 78 up‐regulated genes. (C) Bubble plots of the cellular component (CC) enrichment for 78 up‐regulated genes. (D) Bubble plots of Kyoto Encyclopedia of Genes Genomes pathway enrichment for 78 up‐regulated genes.

Identification of hub genes in heart failure

Next, we investigated the correlations between the 78 intersecting genes using STRING and visualized the results using Cytoscape (Figure ). ASB1 and CUX1 were found in both the key module and the PPI network, suggesting that these genes may act as hub genes in the progression of HF.

Figure 7

Identification of hub genes in heart failure using protein–protein interaction analysis. The protein–protein interaction network of the 78 differentially expressed up‐regulated genes. Dots represent genes, and the size of the dot represents the degree of connection with other genes.

Discussion

In the present study, a total of 135 overlapping DEGs were identified in the GSE135055 and GSE123976 datasets, including 133 up‐regulated DEGs and 2 down‐regulated DEGs. We identified a key module and six hub genes in the GSE135055 dataset by constructing a WGCNA network, including UCK2, ASB1, CCNI, CUX1, IRX6, and STX16. Seventy‐eight intersecting genes were identified by taking the intersection of genes in the key module and the 135 DEGs in the HF group compared with the NF group. Enrichment analysis revealed that 78 intersecting genes were mainly involved in the MAPK and PI3K‐AKT signalling pathways. A PPI network was constructed to investigate the PPIs among the 78 genes, and CUX1 and ASB1 were defined as hub genes in the progression of HF.

Weighted gene coexpression network analysis is a method for exploring the relationships among genes, modules, and clinical traits based on a weighted approach. WGCNA involves expression clustering analysis, phenotypic correlation analysis, the analysis of correlations between genes, gene module definition, coexpression network construction, and the identification of connections between modules and phenotypes. Therefore, WGCNA is a reliable and effective tool for exploring pairwise connections among numerous genes. Here, six hub genes (UCK2, ASB1, CCNI, CUX1, IRX6, and STX16) were identified by WGCNA. A previous study discovered the abnormal expression of these genes in advanced‐stage HF patients. CUX1 plays key roles in tumour differentiation and metastasis. , Abnormal CUX1 expression is associated with the Ebstein anomaly, which affects cardiomyocytes and myocardium differentiation. IRX6 was associated with point mutations and small insertions/deletions in a Dutch Brugada syndrome cohort.

In the current study, a PPI network analysis showed that ASB1 was defined as a hub gene and played important roles in PPIs in HF. ASB1 is a member of the ankyrin repeat and SOCS box‐containing (ASB) family of proteins (1‐18). The ASB1 gene, located on chromosome 2q37, encodes a protein with 335 amino acids. ASB1 is associated with several cullin‐associated neddylation proteins that are involved in targeted ubiquitin degradation pathways. Emeny et al. reported that severe anxiety was associated with 48.5% increased methylation at a single CpG site located in the promoter of the gene encoding ASB1. These results indicated that epigenetics regulated the expression of the stress‐responsive ASB1 gene in an anxiety‐related phenotype manner. The regulation of methylation plays an important role in cardiac function. Ortega et al. found that a special ASB1 methylation pattern was linked to left ventricular structure and performance in end‐stage ischaemic cardiomyopathy. The gain of methylation of ASB1 CpG islands is closely related to LV function, dimension, and output. Further study will be needed to elucidate the molecular mechanism of ASB1 in the progression of HF.

In the phosphatidylinositol‐3‐kinase (PI3K)/protein kinase B (PKB/Akt) signal transduction pathway, PI3K is an intracellular phosphatidylinositol kinase. The main effector activated downstream is Akt, also known as serine/threonine protein kinase. In this study, the analyses of the biological functions of the 78 up‐regulated genes in the yellow module showed that these genes were mainly enriched in the MAPK and PI3K‐AKT signalling pathways. Chang et al. reported that tetrahydrobiopterin reverses left ventricular hypertrophy and diastolic dysfunction in spontaneously hypertensive rats through PI3K/Akt. It has been reported that the MAPK signalling pathway is involved in coronary artery disease and that adriamycin induces HF via the MAPK/AMPK pathway. Moreover, cardiac contractility modulation has been shown to alleviate chronic HF in a rabbit model through the PI3K/AKT pathway. Our results supported the role of the MAPK and PI3K‐AKT signalling pathways in the development of HF. However, constructing a WGCNA on HF downloaded from a GO dataset without verification by qPCR is one limitation of the current study. Another limitation is small samples, as well as potential effects of patient race.

Conclusions

One susceptibility module and two hub genes for HF were identified by coexpression network based on bioinformatic analyses, and we found that the module might cause HF through the MAPK and PI3K‐AKT signalling pathways. Our results may provide potential biomarkers for the diagnosis of HF and improve our knowledge of the molecular mechanisms underlying HF.

Conflict of interest

The authors declare that they have no competing interests.

Funding

This working was supported by grants from the National Nature Science Foundation of China (81970342) and Key research and development program of Jiangsu Province (BE2018611).

Table S1. The clinical characteristics of samples in the GSE135055 and GSE123976 datasets.

Click here for additional data file.