Screening of Secretory Proteins Linking Major Depressive Disorder with Heart Failure Based on Comprehensive Bioinformatics Analysis and Machine Learning

We conducted a two-sample MR analysis to investigate the causal relationship between MDD and HF. The summary data for MDD (exposure, ieu-a-1188) and HF (outcome, ebi-a-GCST009541) were obtained from the IEU OpenGWAS database (https://gwas.mrcieu.ac.uk/↗, accessed on 5 January 2024); the use of which has received informed consent from participants and has been granted ethical approval. These included 59,851 cases and 113,154 controls for MDD and 47,309 cases and 930,014 controls for HF.

Three microarray datasets relating to the HF and control groups, namely GSE19303, GSE1145, and GSE21610, were retrieved from the GEO database (accessible at https://www.ncbi.nlm.nih.gov/geo/↗, accessed on 5 January 2024). The raw expression profile datasets of peripheral blood mononuclear cells (PBMCs) (GSE38206) and brain tissues (GSE80655) including anterior cingulate cortex (AnCg), dorsolateral prefrontal cortex (DLPFC), and nucleus accumbens (nAcc) from MDD patients were also obtained from GEO database as well. The integration of HF-related expression profiles was achieved through the application of the "ComBat" function of the "SVA" package in R (version 4.3.1). This process involved the batch correction of three HF-related datasets, resulting in a final compilation of 76 HF samples and 27 control samples. Detailed descriptive information on datasets and samples is shown in Table S1.

The HF combined dataset and MDD-related datasets underwent background correction, normalization, and gene symbol conversion. Subsequently, the identification of DEGs in HF- and MDD-related datasets was pinpointed through the utilization of the "Limma" and "DESeq2" packages in the R software program. Consequently, the screening of DEGs in HF- and MDD-related datasets was executed based on the criteria of adjusted p-value < 0.05 and |log2 (fold change)| > 0.5. Following this, the expression patterns of the DEGs were visually represented through the utilization of the "ggplot2" package for volcano plots and the "pheatmap" package for heatmaps in the R software program, respectively.

The utilization of WGCNA, recognized as a systematic biological methodology, aimed to disclose patterns of gene associations among various samples and to identify potential biomarker genes or therapeutic targets. The interplay among gene sets and their correlation with phenotypes assumed an integral role in the selection process. Step 1 involves removing genes with a median absolute deviation of 0 from each sample. In Step 2, the "goodSamplesGenes" function is utilized to identify missing values, and samples exceeding a cut height threshold of 20,000 are eliminated as outliers. In Step 3, an optimal soft threshold of 3 is determined using cex1 = 0.85, facilitating the construction of a scale-free co-expression gene network. Following gene clustering, modules are obtained, and similar ones are merged based on MEDissThres = 0.25 criteria. Step 4 entails generating a heat map illustrating module–trait relationships, while Step 5 involves calculating module membership (MM) and gene significance (GS) values and plotting scatter plots depicting MM-GS correlations for each module [11]. The specific codes for WGCNA can be found in Supplementary File S1.

The Human Protein Atlas database (accessible at https://www.proteinatlas.org/↗, accessed on 6 January 2024) provided access to download secretory proteins. The genes encoding secretory proteins, totaling 3947, were acquired from the class of "SPOCTOPUS-predicted secretory proteins" (https://www.proteinatlas.org/search/protein_class%3ASPOCTOPUS+predicted+secreted+proteins↗, accessed on 6 January 2024).

The exploration of interactions between MDD-associated secretory proteins and HF key genes involved the establishment of a PPI network. This network integrated MDD and HF based on data derived from the STRING database (https://www.string-db.org↗, accessed on 30 May 2024). A medium confidence score threshold exceeding 0.4 was applied in this context. Subsequently, the visualization of the PPI network was achieved through the utilization of the Cytoscape software (version 3.9.0). Furthermore, an additional step involved the implementation of the Cytoscape plug-in known as molecular complex detection (MCODE) for the identification of noteworthy subsets in the network. Subsets with scores greater than 10 were singled out as MDD-related pathogenic genes for further analytical procedures.

The exploration of the biological function and mechanistic underpinnings of MDD-related pathogenic genes involved the application of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses. The genes were imported into the Sangerbox database (http://www.sangerbox.com/tool↗, accessed on 30 May 2024) for this purpose. The enrichment significance was determined with a threshold set at p < 0.05. Furthermore, the findings of the functional enrichment analysis were visually represented through the utilization of a lollipop chart and a circos plot.

The cMAP, an online resource (https://clue.io↗, accessed on 30 May 2024), serves as a repository of gene expression profiles that operates by analyzing the impact of gene expression signatures, thereby elucidating connections between genes, diseases, and small-molecule compounds. According to official instructions, the up-regulated genes of cardiac origin in the MDD-related pathological genes were subsequently incorporated into the cMAP online database to unearth prospective small-molecule drugs for HF therapy. Ultimately, the study successfully pinpointed the top 10 compounds characterized by the most elevated enrichment scores.

For the identification of candidate biomarkers and establishment of a diagnostic model for HF, a combination of four different types of ML algorithms, including the random forest (RF), eXtreme Gradient Boosting (XGB), support vector machine (SVM), and generalized linear model (GLM), were employed in this research. We intersected MDD-associated secretory proteins with upregulated DEGs in HF and obtained 16 intersection genes. Then, four ML algorithms were used to perform diagnostic analysis of the intersected genes in the HF combined dataset, obtaining an importance score for each gene, and to rank the importance of each gene with a lollipop chart. In order to get the common important genes recognized by different algorithms, we intersect the top 5 genes scored by each of the four ML algorithms to obtain the hub genes. We used a violin diagram to show the expression pattern of hub genes and used a receiver operating characteristic (ROC) to analyze the diagnostic value of each hub gene for HF. The above four ML algorithms were carried out using the "glmnet", "caret", "random-Forest", "kernlab", and "xgboost" packages in R, and the detailed parameters of the ML algorithms for screening HF biomarkers are shown in. Supplementary File S1

The construction of the nomogram was achieved based on the two hub genes through the application of the "rms" package. The assessment of diagnostic performance for HF entailed the generation of ROC to examine the nomogram. In addition, calibration curves and a decision curve analysis (DCA) enabled the evaluation of the predictive efficiency of the nomogram in HF.

Two independent datasets (GSE5406 and GSE141910) containing HF cases and controls were obtained from the GEO database. The expression patterns of hub genes in external datasets were analyzed by violin diagram, and the diagnostic effect of hub genes and nomogram model was evaluated by ROC analysis.

The quantification of immune cell and stromal cell infiltration from the gene expression profile of HF was performed using the "MCPcounter" package. The box plot, generated using the "ggplot2" package, was employed to present these significant differences. Following that, an examination of the connection between 10 subpopulations of immune cells was demonstrated using the "corrplot" package. Ultimately, Spearman's rank correlation coefficient was employed to present the correlation of the expression of diagnostic biomarkers with the quantity of infiltrated immune cells.

We evaluated the causal relationship between 731 immune cell characteristics and HF based on a two-sample MR analysis. GWAS summary statistics for each immune trait are publicly available from the GWAS catalog (the entry numbers are from GCST0001391 to GCST0002121) [15]. The subsequent analysis process is referred to as MR analysis for MDD on HF.

The establishment of the chronic unpredictable mild stress (CUMS) model adhered to the methodologies delineated in prior research [16]. The CUMS molding time lasted 12 weeks, and the following environmental stress-inducing elements were applied: (1) the confinement of the subjects within small cages equipped with breathing holes; (2) the imposition of continuous overnight illumination; (3) the tilting of the cages to a 40° angle along the vertical axis for a duration of 24 h; (4) the reversal of the light/dark cycle for a duration of 24 h; (5) the introduction of damp bedding, consisting of a mixture of water (200 mL) and sawdust (100 g) for 24 h; and (6) the exposure to noise stimulation at a frequency of 60 Hz for a duration of 1 h. The application of the stressful stimulus occurred daily at 9:00 a.m., with a deliberate randomization of stimulus in order to prevent predictability by the animals. The control group of animals underwent an identical duration of housing in the absence of any form of stressful stimuli. Upon the culmination of the experiment, echocardiography was conducted, followed by the administration of anesthesia to euthanize the animals. Subsequently, serum and heart tissue samples were harvested for further experimental procedures.

The extraction of total RNA from ventricular tissue was carried out utilizing the Trizol reagent (Vazyme, Nanjing, China). Thereafter, reverse transcription was carried out using a Reverse Transcription Kit (Vazyme), and a SYBR Green PCR Kit (Vazyme) was employed for the real-time quantitative PCR (RT-qPCR). Duplicate reactions were carried out to ensure accuracy, and the 2^−ΔΔCt approach was utilized for the quantitation of relative mRNA expression. The primer sequences are provided, including ISLR (F: 5′-CGGCCCCTCATTACTCCATC-3′; R: 5′-AACTCTCAAAGCGGGTCAGG-3′), SFRP4 (F: 5′-TTTGGCAACGTACCTGAGCA-3′; R: 5′-CACTCCTCTGGACCGCTTTT-3′), and 18S rRNA (F: 5′-GAGAAACGGCTACCACATCC-3′; R: 5′-CACCAGACTTGCCCTCCA-3′). Moreover, the serum levels of ISLR and SFRP4 were validated using the designated ELISA kits (Bioswamp, Wuhan, China) in accordance with the guidelines provided by the manufacturer.

To perform the statistical analysis, GraphPad Prism version 9.5.1 (GraphPad Software Inc., San Diego, CA, USA) was utilized in this study. Results were displayed as mean ± SD. The comparisons between the two groups were performed by means of an unpaired Student's t-test. p-values < 0.05 were considered significant and significance for all figures determined as * p < 0.05, ** p < 0.01, *** p < 0.001, and **** p < 0.0001.

The MR and bioinformatics analysis strategy is shown in Figure 1. The IVW analysis showed that MDD was linked to an increased risk of HF (OR, 1.129; 95% CI, 1.052–1211; p < 0.001) (Figure 2A). The effect size of each individual SNP is shown in Figure 2B,C, and the results of the leave-one-out sensitivity analysis (Figure 2D) suggested that the causal associations between MDD and HF were not affected by any individual SNP.

Subsequently, to batch correction, the HF combined dataset, comprising 76 samples from the HF group and 27 samples from the control group, was obtained. As illustrated in Figure S1A–D, the removal of batch effects led to a substantial reduction in differences observed among the three datasets. The DEA between the combined HF and control samples revealed a total of 855 DEGs, with 482 genes found to be upregulated and 373 genes downregulated. The expression pattern of these DEGs in the HF combined dataset is illustrated using a volcano plot and heatmap (Figure 3A,B).

To delve deeper into the identification of HF key genes, a WGCNA was performed to pinpoint the gene modules most closely linked to HF samples. The selection of the soft-thresholding power of 3, based on considerations of scale independence and average connectivity, is highlighted in Figure S2A. This power was applied to generate a total of 11 modules (Figure S2B), and the display of the clustering of module eigengenes can be observed in Figure 3C. Additionally, an investigation into the correlation between HF and 11 gene modules (Figure 3D) was carried out in this study. The analysis unveiled that the green module displayed the most positive correlation with HF (1146 genes, p = 7 × 10⁻⁴). The modules colored in green, yellow, and blue, which exhibited statistical significance, were designated as key modules for subsequent investigations. Moreover, we found a strong association between MM and GS in the green (r = 0.56, p = 1.4 × 10⁻⁹⁵), blue (r = 0.33, p = 3.7 × 10⁻⁵⁹), and yellow modules (r = 0.35, p = 9.9 × 10⁻⁴⁰), respectively (Figure S2C–E). To identify key genes associated with HF, an intersection analysis was performed between DEGs and genes from WGCNA key modules. This process yielded a total of 315 HF key genes for further analysis (Figure 3E).

As illustrated via volcano plots and heatmaps in Figure 4A–H, 874 DEGs were pinpointed in PBMCs of individuals with MDD, and in the regions of AnCg, DLPFC, and nAcc, the DEGs were found to be 499, 683, and 1171, respectively. Given the consideration that the development of HF may be primarily facilitated by the release of secretory proteins in individuals with MDD, the MDD-associated secretory proteins were subsequently acquired by combining the differentially expressed secretory proteins from MDD PBMC (Figure 4I) and brain tissue datasets (Figure 4J–L). This analysis collectively identified 332 MDD-associated secretory proteins.

In order to elucidate the prospective pathogenic genes and mechanistic underpinnings in MDD-related HF, the interplay between MDD-associated secretory proteins and HF key genes was compiled using the STRING database. Through the utilization of the MCODE analysis plugin, subsets with scores greater than 10 were pinpointed, where 78 genes were highlighted as MDD-related pathogenic genes implicated in the development of HF (Figure 5A,B). The analysis of GO terms pertaining to biological process (BP) demonstrated that the MDD-related pathogenic genes exhibited a predominant enrichment in the "chemokine-mediated signaling pathway", "inflammatory response", and "extracellular matrix organization" (Figure 5C). The results of the molecular function (MF) analysis uncovered that the pathogenic genes displayed a strong association with "cytokine activity" and "chemokine activity" (Figure 5C). The findings of the KEGG pathway analysis unveiled robust associations between the MDD-related pathogenic genes and key pathways, namely "cytokine–cytokine receptor interaction" and "PI3K-Akt signaling pathway". Additionally, these genes exhibited enrichment in fundamental factors contributing to HF, including "fluid shear stress and atherosclerosis", "hypertrophic cardiomyopathy", and "dilated cardiomyopathy" (Figure 5D). Furthermore, drug screening analysis of the cMAP database showed that the top 10 compounds with the highest negative scores, namely alpha-linolenic-acid (ALA), nicotine, estrone, CGP-37157, TTNPB, kinetin-riboside, PD-173074, STO-609, methylnorlichexanthone, and vinorelbine, were identified as potential pharmacological agents for the management of MDD-related HF (Figure 5E). The description of the targeted pathways pertaining to these compounds is manifested in Figure 5F.

Given the potential crucial involvement of common differentially expressed secretory proteins in MDD and HF, 16 common genes were identified at the intersection of MDD-associated secretory proteins and up-regulated HF DEGs (Figure 6A). To refine the selection of diagnostic biomarkers, four ML algorithms, namely GLM, SVM, RF, and XGB, were employed to assign rankings to the 16 common genes on the basis of the importance of each gene as a variable indicator (Figure 6B–E). It is worth noting that after the intersection of the top 5 important genes obtained by each of the four algorithms, only two hub genes were obtained, namely immunoglobulin superfamily containing leucine-rich repeat (ISLR) and secreted frizzled-related protein 4 (SFRP4) (Figure 6F). The violin plot shows the expression patterns of ISLR and SFRP4 in the HF combined dataset (Figure 6G,H), with ROC analysis showing that their AUC for HF diagnosis reached 0.831 and 0.881 (Figure 6I), respectively, suggesting their potential as HF biomarkers. To enhance the diagnostic accuracy and predictive capability, a nomogram was constructed incorporating the two hub genes (Figure 6J). The analysis of calibration curves unveiled a close resemblance between the predicted probabilities generated by the established nomogram diagnostic model and those of the ideal model (Figure 6K). Furthermore, the performance of the nomogram was assessed through DCA, revealing that the nomogram model-based decision-making could potentially offer advantages in the diagnosis of MDD-related HF (Figure 6L). In the ROC analysis, the AUC value of the nomogram model reaches 0.904 (Figure 6M).

We obtained two independent HF datasets (GSE5406 and GSE141910) from the GEO database. As shown in Figure 7A–D, ISLR and SFRP4 were significantly up-regulated in the HF group of both datasets. In the GSE5406 dataset, the AUC of SFRP4 and ISLR for HF diagnosis was 0.890 and 0.950, respectively (Figure 7E). In GSE141910, the AUC of SFRP4 and ISLR for HF diagnosis was 0.967 and 0.902, respectively (Figure 7G). In addition, the diagnostic efficiency of the nomogram model in GSE5406 and GSE141910 reached 0.948 and 0.976, respectively, suggesting that the combination of the two biomarkers has excellent diagnostic value.

A close association was revealed between immune processes and the function and pathway analysis of MDD-related pathogenic genes in HF. The characteristics of eight subpopulations of immune cells and two subpopulations of stromal cells were derived using the MCPcounter algorithm. Higher proportions of B lineage, CD8 T cells, T cells, and fibroblasts were observed in HF in contrast to the control group (Figure 8A). Furthermore, a significant positive correlation (r = 0.43, p < 0.05) was found between T cells and B lineage during these 10 cell subpopulations, and fibroblasts demonstrated a positive association with endothelial cells (r = 0.4, p < 0.05) (Figure 8B). Notably, as displayed in Figure 8C, the hub genes, ISLR and SFRP4, showed a significant correlation to immune-related cell accumulation (especially T cells and fibroblasts) in HF. Additionally, the MR analysis showed that 31 immunophenotypes were associated with the development of HF, among which IgD+CD38+ B cells and CD39+CD8+ T cells were risk factors, and activated CD4 regulatory T cells were protective factors (Figure 8D).

A well-recognized animal model of MDD (CUMS model) was established to enhance the validation of the integrated bioinformatics analysis described earlier. The echocardiographic results unveiled a notable reduction in left ventricular ejection fraction (LVEF) in the MDD group when compared to the control group, signifying decreased systolic function of the heart (Figure 9A). Nonetheless, it was observed that certain individuals in the MDD group did not manifest a reduction in LVEF. Consequently, the MDD group was stratified into two subgroups, namely the MDD with reduced LVEF (rLVEF) group and the MDD without rLVEF group on the basis of the mean LVEF value of rats in the MDD group. A notable observation was made wherein no significant difference in LVEF was evident between the MDD without rLVEF group and the control group. In contrast, rats in the MDD with rLVEF group displayed a marked reduction in LVEF compared with the control and MDD without rLVEF groups (Figure 9B). The findings from an RT-qPCR analysis provided confirmation that two hub genes in the MDD with rLVEF group exhibited an upregulation in the expression patterns in contrast to the control and MDD without rLVEF groups (Figure 9C,D). Subsequent correlation analysis unveiled a significant negative correlation between the mRNA expression of ISLR (r = −0.653, p < 0.001) and SFRP4 (r = −0.476, p = 0.008) with LVEF (Figure 9E,F). In addition, the ELISA detection unveiled an elevated trend in the levels of ISLR and SFRP4 in the serum of the MDD with the rLVEF group when compared to non-HF rats (Figure 9G,H). A diagnostic nomogram model, incorporating ISLR and SFRP4, was developed for the prediction of the probability of rLVEF occurrence based on the in vivo experimental data (Figure 9I). The application of the nomogram in decision-making was found to potentially enhance the accuracy of rLVEF prediction, as indicated by the calibration curves and DCA results (Figure 9J,K). The ROC curve analysis unveiled that the AUC for ISLR and SFRP4 in diagnosing rLVEF reached 0.759 and 0.836, respectively. Notably, the AUC of the nomogram model reached 0.903, indicating its significantly advantageous diagnostic performance (Figure 9L,M).