Identification of Hub Diagnostic Biomarkers and Candidate Therapeutic Drugs in Heart Failure.

PURPOSE: The objective of this study was to identify the potential regulatory mechanisms, diagnostic biomarkers, and therapeutic drugs for heart failure (HF).
METHODS: Differentially expressed genes (DEGs) between HF and non-failing donors were screened from the GSE57345, GSE5406, and GSE3586 datasets. Database for Annotation Visualization and Integrated Discovery and Metascape were used for Gene Ontology and Kyoto Encyclopedia of Genes and Genomes analyses respectively. The GSE57345 dataset was used for weighted gene co-expression network analysis (WGCNA). The intersecting hub genes from the DEGs and WGCNA were identified and verified with the GSE5406 and GSE3586 datasets. The diagnostic value of the hub genes was calculated through receiver operating characteristic analysis and net reclassification index (NRI). Gene set enrichment analysis (GSEA) was used to filter out the signaling pathways associated with the hub genes. SYBYL 2.1 was used for molecular docking of hub targets and potential HF drugs obtained from the connection map.
RESULTS: Functional annotation of the DEGs showed enrichment of negative regulation of angiogenesis, endoplasmic reticulum stress response, and heart development. PTN, LUM, ISLR, and ASPN were identified as the hub genes of HF. GSEA showed that the key genes were related to the transforming growth factor-β (TGF-β) and Wnt signaling pathways. Sirolimus, LY-294002, and wortmannin have been confirmed as potential drugs for HF.
CONCLUSION: We identified new hub genes and candidate therapeutic drugs for HF, which are potential diagnostic, therapeutic and prognostic targets and warrant further investigation.

Introduction

HF is a clinical syndrome characterized by congestion of the lungs and vena cava, leading to abnormal heart structure or function, which is the final stage of the development of heart disease.1 Approximately 40 million people worldwide suffer from HF, and the incidence rates are steadily rising.2 Despite significant progress in the HF management in recent decades, the treatment options are mainly palliative rather than curative.3 Given the complex pathogenesis of HF, it is essential to elucidate the underlying molecular mechanisms in order to identify potential therapeutic targets and prognostic markers.

Bioinformatics is a high-throughput technique that can screen multiple databases to identify the potential pathological biomarkers of various diseases.4 Weighted gene co-expression network analysis (WGCNA) is a systems biology application that mines the genetic interaction networks to construct highly coordinated gene modules.4 WGCNA has been widely used for detecting disease biomarkers, and elucidating biological mechanisms and drug interactions.5–7 Although biomarkers of HF have been identified, but due to heterogeneity of HF and its complicated pathophysiological manifestations, a single gene cannot accurately predict the characteristics of HF.8,9

In this study, the differentially expressed genes (DEGs) between HF patients and non-failing donors (NFD) were screened from multiple GEO datasets and functionally annotated. The hub genes were then screened through the degree of connectivity in the PPI network, and used to build a co-expression network with WGCNA. The intersecting hub genes between DEGs and WGCNA were identified and validated in other human HF datasets. Gene set enrichment analysis (GSEA) was used to discover the signaling pathways associated with these hub genes. Finally, the potential HF drugs were predicted through the Connectivity Map (cMap) database, and molecular docking between the drug candidates and hub genes was simulated using SYBYL 2.1 software.

We conducted a bioinformatics analysis using DEGS and WGCNA to further investigate the occurrence and development of HF and identify potential therapeutic drugs for biomarkers of HF.

Materials and Methods

Data Collection and Preprocessing

The study design is outlined in Figure 1. We searched the GEO database and included data on human heart tissue samples in this study. The mRNA expression profiles from HF and NFD samples were downloaded from the GSE57345, GSE5406 and GSE3586 datasets of the GEO database (). The subjects included in our study suffered from heart failure with reduced ejection fraction (HFrEF) and heart failure with preserved ejection fraction (HFpEF). Among them were 96 patients with ischemic heart disease, 84 patients with dilated cardiomyopathy (CMP), and 139 non-heart failure samples from GSE57345.10 The GSE5406 dataset included 16 samples without heart failure, 86 patients with dilated cardiomyopathy (CMP), and 108 patients with ischemic heart disease.11 The GSE3586 dataset included 13 patients with dilated cardiomyopathy and 15 patients without heart failure.12 The gene annotation files of GSE57345, GSE5406 and GSE3586 were GPL11532, GPL96, and GPL3050 respectively. GEO2R online tool was used for screening DEGs between HF and NFD with p < 0.05 (calculated by t-test) as the threshold.

Figure 1

Flowchart of study outline. The flowchart constructed by Draw io online tool ().

Functional Enrichment Analysis of the Overlapping DEGs

The overlapping DEGs were uploaded to the Database for Annotation, Visualization, and Integrated Discovery (DAVID, ) and Metascape () databases for Gene Ontology (GO) function and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses respectively. P<0.05 was considered statistically significant.

PPI Network Analysis

The protein-protein interaction (PPI) network was constructed using the STRING database, and visualized with the Cytoscape software (3.7.2). DEGs with connectivity ≥ 5 in the PPI network were considered to be the hub genes.

WGCNA

The R package WGCNA was used to construct the weight co-expression network of the GSE57345 dataset. The GSE57345 dataset was used in the WGCNA analysis because it contained the largest sample size. The weighing coefficient β was first calculated based on R2> 0.9 of the scale-free real biological network. After determining the adjacency function parameter β, a hierarchical clustering tree of different gene modules was constructed. The Pearson correlation coefficient was then used to transform the correlation matrix into an adjacency matrix and subsequently to a topological overlap matrix (TOM).

Identification of Key Modules and Hub Genes

The correlation among the genes in the aforementioned modules was analyzed and a heatmap was constructed. The module most closely related to the HF state was considered the key module of HF. To verify specific modulus-character associations, the correlation between GS and MM was investigated. Genes with |MM| > 0.8 were subsequently screened out as hub genes in the key module. The intersecting hub genes between the key module and the DEGs were finally defined as the hub genes of HF.

Validation of Hub Genes

We used the GSE57345, GSE5406, and GSE3586 datasets to confirm the identity of the hub genes that may be associated with HF. Then, we used the t-test to determine the significance of the correlation between the expression of hub genes and HF. Receiver operating characteristic (ROC) curves were drawn for the core genes in the three datasets, and the area under the curve (AUC) was calculated. Specificity, sensitivity, and net reclassification index (NRI) were calculated to evaluate the value of the genetic diagnosis.

Functional Analysis of Hub Genes by GSEA

GSEA was performed to clarify the biological functions of the hub genes using the KEGG gene set (c2.cp. kegg. v7.2.) as the default and p < 0.05, as the threshold. The HF samples of the GSE57345 data set were divided into high and low expression groups of each hub gene. The enrichment graph was plotted using the clusterProfiler package of the R language and GSEA function.

Screening of Potential Therapeutic Drugs

The CMap database includes 1309 compounds and expression data of > 7000 human genes. The DEGs intersecting GSE57345, GSE5406, and GSE3585 datasets were uploaded to CMap. The negatively correlated small molecules were screened out using p < 0.0001 and mean < 0.4 as the criteria. ChemBioDrawUltta 17.0 software () was used to draw the 3D structural formulas of potential therapeutic drugs and save them in mol2 format as small molecule compounds. The 3D crystal structures of the core targets were downloaded from the UniProt database (). The Surflex-Dock module of SYBYL2.1 software was used for molecular docking, with total score >4 as the threshold for binding ability. The docking results were visualized using the Pymol software.

Results

Identification and Analysis of DEGs

A total of 583 DEGs were identified between the HF and NFD samples across three GEO datasets (Figure 2A–D). The most significantly enriched GO terms pertaining to biological process (BP) were negative regulation of angiogenesis (p = 1.55E-04), response to endoplasmic reticulum stress (p = 6.63E-04), heart development (p = 0.002463), regulation of ventricular cardiac muscle cell action potential (p = 0.004523), MAPK cascade (p = 0.005025), and blood vessel development (p =0.007107) (Figure 3A and ). The cellular components (CC) terms including extracellular exosome (p =2.72E-08), cytoplasm (p =9.57E-06), actin cytoskeleton p (p =4.94E-05), mitochondrion (p =8.84E-05), nucleoplasm (p =1.03E-04), Golgi apparatus (p =1.21E-04) and lysosome (p =6.84E-04) were significantly enriched (Figure 3B and ). The top enriched molecular function (MF) terms were activating transcription factor binding (p = 0.003852), actin filament binding (p =0.006599), cadherin binding involved in cell-cell adhesion (p = 0.00701), transcription coactivator activity (p =0.018332) and collagen binding (p =0.034108) (Figure 3C and ). In addition, KEGG analysis revealed that the Ras signaling pathway (p =5.15548E-07), Focal adhesion (p =1.11814E-05), Lysosome (p = 2.43906E-05), MAPK signaling pathway (p = 4.23641E-05), PI3K-Akt signaling (p = 4.85396E-05), Protein processing in endoplasmic reticulum (p = 5.21303E-05) and Hippo signaling pathway (p =7.03525E-05) were significantly enriched among the DEGs (Figure 3D and ).

Figure 2

The DEGs between HF and NFD. The volcano plots of DEGs in (A) GSE57345, (B) GSE5406 and (C) GSE3586. (D) Venn diagrams of DEGs in three data sets.

Figure 3

GO and KEGG pathway enrichment analysis. Significantly enriched GO terms for (A) Biological processes, (B) Cellular component, (C) Molecular function. (D) KEGG pathway Molecular function p < 0.05 is considered statistically significant.

Construction of WGCNA Network

A total of 1589 genes were screened from the GSE57345 dataset for WGCNA (p < 0.05, |Log2FoldChange| > 0.5). Sample clustering showed no significant differences in the WGCNA (Figure 4A). At β = 4, the scale-free network fitting index R2 was 0.9, and the average connectivity approached 0, indicating that this value could obtain a scale-free network that met all requirements. Thus, β = 4 was selected to construct a scale-free network (Figure 4B–C). A dynamic shearing algorithm was used to cluster the genes and module divisions.

Figure 4

Determination of soft-threshold power β. (A) Clustering dendrogram of 319 samples. (B) Scale-free topology fit index as a function of the soft-threshold power. The red line indicates that R2 is equal to 0.9. (C) Mean connectivity as a function of the soft-threshold power.

Identification and Visualization of the Key Module

Five gene co-expression modules were finally obtained by calculating the module feature vector of each and merging similar modules (Figure 5A). The genes were clustered in the black, blue, yellow and green-yellow modules, and those that could not be clustered into any module were specified to the gray module. The yellow module with 73 genes showed the strongest correlation with HF (r = 0.77, p = 1e-61) (Figure 5B), as well as with the clinical phenotype as per GS and MM analyses (cor = 0.96, p = 5.5e-41; Figure 5C–G). The genes distributed in the upper right corner were closely related to HF pathogenesis, and are likely the key disease genes. Twenty-one genes in the yellow module were confirmed as hub genes.

Figure 5

Identification of key HF gene modules. (A) Clustering dendrograms of genes and module detecting. (B) Heat map of the correlation between HF modules. (C–G) Correlation of GS and MM in the HF-related module. p < 0.05 is considered statistically significant.

Identification and Verification of Hub Genes

From the 583 overlapping DEGs, 322 hub genes were selected for subsequent analysis (Figure 6A). The key intersecting genes between the 322 hub DEGs and 21 hub genes of the yellow module included PTN, LUM, ISLR and ASPN. The expression levels of these potentially key genes were analyzed in the HF and NFD samples of the GSE57345, GSE5406 and GSE3586 datasets. We further visualized their expression levels in the GSE57345 data set, and found that all four genes were overexpressed in HF samples compared to the NFD group (p < 0.05, Figure 6B). Afterwards, the genes were verified in GSE5406 and GSE3586, which verified higher expression levels in the HF group (p < 0.05, Figure 6C and D).

Figure 6

Analysis of key genes. (A) The Venn diagram of hub genes in the yellow module and hub genes in DEGs. Expression of PTN, LUM, ISLR and ASPN in the (B) GSE57345, (C) GSE5406 and (D) GSE3586 datasets. **p < 0.01 and ***p < 0.001 are considered statistically significant.

ROC Curve Analysis of Key Genes

The potential diagnostic value of PTN, LUM, ISLR and ASPN was ascertained by plotting the ROC curve based on the expression data in GSE57345, GSE5406 and GSE3586. As shown in Figure 7A–C, the AUC of all genes in all datasets exceeded 0.9, except for 0.785 calculated for ISLR in GSE3586. We uploaded AUC and Standard Error data into MedCalc software, and applied the Z test to compare the expression of PTN, LUM, ISLR and APSN between the datasets. The results showed that there was no statistical difference (p >0.05). We used the NRI to analyze differences in the expression levels of the four predicted HF hub genes in the GSE57345, GSE5406 and GSE3586 datasets. The PTN and ISLR of the GSE5406 dataset showed significant differences in predicting HF (NRI [95% CI]: 0.3228 [0.0361–0.6095], p: 0.027); the prediction effects of the PTN and ISLR genes were significantly different, NRI [95% CI]: 0.3905 [0.0073–0.7736], p: 0.045). The prediction effects of ISLR and LUM gene in the GSE5406 dataset were significantly different (NRI [95% CI]: −0.5077 [−0.9361–0.0793], p: 0.020). Subsequently, ROC analysis was performed to determine the diagnostic value of the four key genes, and the results suggested that these four hub genes can diagnose HF with high sensitivity and specificity (Table 1).

Figure 7

The ROC curve of hub genes (A) GSE57345. (B) GSE5406. (C) GSE3586. The x-axis shows specificity, and the y-axis shows sensitivity.

Table 1

ROC Curve Analysis of Hub Genes

Data set	Gene	Specificity (%)	Sensitivity (%)	Youden Index	AUC	Standard Error	95% CI
GSE57345	PTN	87.5	81.9	0.6942	0.912	0.0171	0.875–0.941
	LUM	83.8	93.2	0.7704	0.917	0.0177	0.881–0.945
	ISLR	83.1	88.7	0.7179	0.921	0.0155	0.886–0.949
	ASPN	80.1	91.0	0.7111	0.912	0.0182	0.875–0.941
GSE5406	PTN	91.8	93.7	0.8550	0.939	0.0303	0.898–0.967
	LUM	90.7	87.5	0.7822	0.960	0.0187	0.924–0.982
	ISLR	93.8	93.7	0.8756	0.940	0.0438	0.899–0.968
	ASPN	98.5	87.5	0.8595	0.952	0.0324	0.914–0.977
GSE3586	PTN	73.3	100	0.7333	0.918	0.0508	0.751–0.988
	LUM	86.7	100	0.8667	0.938	0.057	0.779–0.994
	ISLR	73.3	84.6	0.5795	0.785	0.0909	0.589–0.916
	ASPN	86.7	100	0.8667	0.938	0.0529	0.779–0.994

Abbreviations: AUC, Area under the ROC curve; 95% CI, 95% Confidence interval.

GSEA of Hub Genes

GSEA of PTN, LUM, ISLR and ASPN revealed direct involvement in the pathogenesis of HF. As shown in Figure 8A–D, all genes were enriched in arrhythmogenic right ventricular cardiomyopathy (ARVC), dilated cardiomyopathy, ecm receptor interaction, focal adhesion, gap junction, hypertrophic cardiomyopathy (HCM), regulation of actin cytoskeleton and TGF-β signaling pathway. In addition, PTN, LUM and ASPN were enriched in the WNT signal pathway, LUM in the calcium signaling pathway, and ASPN is likely involved in seleno-amino acid metabolism.

Figure 8

Results of GSEA. (A) PTN. (B) LUM. (C) ISLR. (D) ASPN.

Identification of Potential Therapeutic Drugs

There were 264 upregulated and 499 down-regulated genes intersecting across the three datasets. The genes were uploaded to the cMap database to filter out negatively related small molecule compounds (p < 0.0001 and mean < 0.4). Sirolimus, LY-294002, and wortmannin were identified as potential drugs of HF (Figure 9A–C). Molecular docking showed that sirolimus had good affinity for PTN, ISLR, LUM, and ASPN, and wortmannin for PTN and LUM (Figure 9D). The molecular docking diagrams of potential compounds and hub targets are shown in Figure 9E–J.

Figure 9

The potential therapeutic drugs of HF. (A–C) The 2D structure of Sirolimus, LY-294002 and Wortmannin. (D) Heat map of the docking score between potential drugs and hub targets. The intensity of red color indicates binding ability. (E–J) Molecular docking diagram of certain core compounds and hub targets.

Discussion

In this study, we combined WGCNA and DEGs to screen for genes associated with HF and found that the expression levels of the PTN, ISLR, LUM, and ASPN genes were all upregulated in HF. Further analysis using the ROC curve showed that these four genes may be potential biomarkers of HF. At present, PTN and ISLR have not been reported to be associated with HF, but there is evidence that they may be potentially associated with HF. Single-gene GSEA showed that the hub genes of HF are related to arrhythmic right ventricular cardiomyopathy, dilated cardiomyopathy and hypertrophic cardiomyopathy, which is consistent with reports demonstrating that these cardiovascular disorders precede the final HF stage.13–15 Single-gene GSEA showed that the hub genes of HF are related to arrhythmic right.

GO analysis of the DEGs showed that endoplasmic reticulum stress (ERS) is closely related to HF pathogenesis, which is consistent with previous reports.16–18 The risk factors of HF can induce ERS in myocardial cells, which culminates in apoptosis and cardiovascular dysfunction. In addition, the DEGs were enriched in ferroptosis, MAPK signaling pathway, PI3K-Akt signaling pathway, and the Hippo signaling pathway, all of which are involved in HF. Liu et al19 detected a high level of ferroptosis in the cardiomyocytes of a rat model of pressure overload-induced HF. Exogenous expression of ferritin FTH1 and GPX4 and reduction in ROS levels through NOX4 knockdown inhibited ferroptosis in cardiomyocytes and improved cardiac function. Fang and Koleini et al20,21 found that doxorubicin induced the accumulation of oxidized phospholipids in undifferentiated cardiomyocytes and up-regulated heme oxygenase1 (HMOX1), resulting in heme degradation, free iron overload, and ferroptosis, which eventually leads to HF. The loss of myeloid differentiation protein 1 (MD1) activates ROS and exacerbates autophagy induced by the MAPK signaling pathway. Therefore, the MD1-ROS-MAPK axis is a novel therapeutic target for HF that can preserve the ejection fraction.22 Mao Liu et al23 showed that paeoniflorin reduced myocardial fibrosis and improved cardiac function in rats with chronic HF by regulating the p38/MAPK signaling pathway. Apelin-13 can slow down oxidative stress by inhibiting the PI3K/Akt signaling pathway in the rat HF model and ameliorate angiotensin II–induced cardiac insufficiency, impaired cardiac hemodynamics, and fibroblast fibrosis.24 Hou and Li et al25,26 showed that YAP/TAZ can initiate the transcription of connective tissue growth factor by interacting with the TEAD domain family, increase the expression of extracellular matrix genes, promote cardiac remodeling and fibrosis, and thus delay the progression of HF. Leach et al27 found that knocking out the SALV gene increased the number of left ventricular myocardial cells in mice with myocardial infarction, which reduced ventricular fibrosis and increased the number of new capillaries around the injured myocardium, indicating that the Hippo signaling pathway can enhance heart function.

PTN is a highly conserved proto-oncogene closely related to tumor angiogenesis and metastasis.28 It is highly expressed in various malignant tumors, such as breast cancer, prostate cancer and rectal cancer,29,30 and promotes the proliferation, mitosis, differentiation, and migration of vascular endothelial cells.31 Overexpression of PTN gene can promote bone formation, whereas PTN gene knockout mice have dysfunctional bone growth and remodeling.32 LUM is a member of the SLRP family of leucine-rich proteins that are secreted by the extracellular matrix. It is widely distributed in various tissues, and shows aberrant expression levels in pancreatic cancer, colorectal cancer, breast cancer and cervical cancer.33 LUM has both oncogenic and tumor-suppressive functions depending on the cancer type. For instance, LUM facilitated the metastasis of colon cancer cells by reconstructing the actin cytoskeleton, but inhibited the adhesion of osteosarcoma cells via the TGF-β2 signaling pathway.34,35 ISLR is a conserved immune-related protein that is mainly expressed in stromal cells.36 Xu et al37 showed that ISLR can inhibit Hippo signal transduction during intestinal regeneration and tumorigenesis and activate YAP factor in epithelial cells. Knocking out ISLR in mouse stromal cells significantly affected intestinal regeneration and inhibited colorectal tumorigenesis. Zhang and Hara et al38,39 further showed that the ISLR can promote muscle regeneration and improve myocardial tissue repair. However, little is known regarding the correlation between ISLR and HF. ASPN is an extracellular matrix protein and a member of the leucine-rich small proteoglycan family.40 Sasaki et al41 showed that ASPN protected gastric tumor cells against oxidative stress by up-regulating HIF1α and reducing the levels of mitochondrial ROS. It also increased the expression of CD44 to accelerate the migration and infiltration of gastric cancer cells. However, other reports indicate an anti-tumorigenic role of ASPN in breast cancer.42,43 Studies also show that ASPN is up-regulated during aortic stenosis or coronary artery ligation in ischemic cardiomyopathy patients and animal models.44 ASPN may also increase the apoptosis and fibrosis of H9C2 cardiomyocytes.45 However, the exact role of ASPN in HF pathogenesis needs further investigation.

Lu A et al46 found that Wnt3a binds to FZD and LRP5/6 receptors, thereby activating the classic Wnt-Dvl-β-catenin signaling pathway and promoting myocardial hypertrophy. Wnt signaling can inhibit Na+ channels by directly or indirectly inhibiting the expression of Scn5a. Thus, blocking these intracellular cascades is a rational therapeutic strategy against HF.46 He et al47 found that Wnt3a and Wnt5a ligand were up-regulated in a mouse model of cardiac hypertrophy, underscoring the role of the Wnt signaling pathway in HF pathogenesis. TGF-β is an important factor regulating myocardial fibrosis, which gradually worsens during HF and alters cardiac function from the compensatory phase to the decompensated phase.48 Kakhi et al49 found that sirolimus, an mTOR inhibitor, reversed new HF after kidney transplantation in mammals. Gao et al50 also showed that rapamycin (sirolimus) can reduce cardiomyocyte apoptosis and promote autophagy by regulating mTOR and ERS, thus preventing myocardial damage caused by chronic HF. LY294002 and wortmannin are protein kinase inhibitors that block the PI3K signaling pathway. Melatonin alleviates cardiac hypertrophy by inhibiting the Akt/mTOR pathway and reducing Atg5-dependent autophagy, which can be reversed by LY294002.51 Studies show that apelin may reduce the myocardial damage caused by acute HF by regulating the APJ/Akt/ERS signaling pathway. However, wortmannin and LY294002 can reverse the cardioprotective effects of apelin.52 The PI3K-Akt signaling pathway was also enriched among the HF-related DEGs, indicating a vital mechanistic role in its pathological progression. Molecular docking showed that sirolimus and wortmannin had a high affinity to the hub targets, LY-294002 bound weakly and may therefore have other targets. Nevertheless, all three drugs could be potentially effective for treating HF.

There are several limitations in this study. First, the data used in this study was obtained from the GEO database, which lacks clinical, in vivo, and in vitro experimental research certifications for pivotal genes and HF-associated genes. Second, the datasets used in this study were relatively small, and a larger sample size is needed to verify our results. However, our findings provide new insights into the underlying molecular mechanisms of HF, along with potential diagnostic biomarkers and candidate therapeutic drugs, which will help provide new clues for HF research, diagnosis and treatment, and target selection.

Conclusion

PTN, LUM, ISLR, and ASPN are overexpressed in HF patients compared to NFD, and are mainly related to the TGF-β and Wnt signaling pathways. Sirolimus, LY-294002, and wortmannin are potential drug candidates for HF treatment. The in silico data will need to be verified by functional and clinical studies.