Three gene expression profiles of DCM with HF (GSE29819, GSE21610, GSE17800) and one single-cell RNA-sequencing dataset of DCM (GSE95140) were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). First, GSE29819 and GSE21610 were merged into the training group for major analysis, containing 20 normal myocardial tissues and 35 myocardial tissues derived from DCM patients with HF (Schwientek et al., 2010; Gaertner et al., 2012). GSE17800 was the test group for validation, including nine normal myocardial tissues and 40 myocardial tissues from DCM patients with HF (Ameling et al., 2013). GSE29819, GSE21610, and GSE17800 were based on the GPL570 sequencing platform. And a total of 21,628 genes and 21,654 genes were obtained in the training group and test group, respectively. The probes were annotated to official gene symbols using the Perl programming. Next, we performed log2 transformation and then normalized the raw count expression data using the function “normalizeBetweenArrays” of the R package “limma”. For the training group, the function “ComBat” of the R package “sva” was then utilized to remove the batch effect of GSE29819 and GSE21610. Moreover, heart specimens from normal individuals and 10 DCM patients were obtained from GSE95140 for further verification and single-cell RNA-seq analysis (Nomura et al., 2018). The sequencing platform of GSE95140 was GPL16791.

To assess the dispersion between the control group and DCM with HF group, PCA was performed and two components were extracted. Then, the R package “limma” was utilized to screen the DEGs between control and DCM with HF samples. The cutoff criteria of DEGs was |log₂Fold Change (log₂FC) | > 1 and adjusted p < 0.05 (Chen et al., 2018). The Volcano plot of all DEGs and the clustering heatmap of the top 50 up-regulated genes and top 50 down-regulated genes were plotted by the R packages “ggplot2” and “pheatmap”, respectively.

To investigate the potential biological function, pathways, and related diseases of up-regulated DEGs, the R packages “clusterProfiler” and “DOSE” were used for Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Disease Ontology (DO) enrichment analyses. The enrichment GO terms, KEGG pathways, and DO terms with q < 0.05 were considered significant. Furthermore, we used the R package “CBNplot” to infer the Bayesian network (BN) on the gene expression profiles and KEGG enrichment analysis, showing the relationships between genes in the enriched pathways (Sato et al., 2022).

The R package “WGCNA” was used to construct the gene co-expression network and identify the modules correlated to DCM with HF (Langfelder and Horvath, 2008). First, a total of 9,379 genes with standard deviation (SD) of expression larger than 0 in the training group were selected for analysis. Second, the function “hclust” was used to cluster the samples to exclude the outliers. Third, a similarity matrix was obtained by calculating each pairwise gene correlation. It was then converted to an adjacency matrix using the soft adjacency function A i j = S i j β , where i and j were any two of 9,379 genes, β is the soft threshold power, S _ij is the similarity matrix, and A _ij is the adjacency matrix. The function “pickSoftThreshold” selected the optimal soft threshold power (β = 1–20) to establish the scale-free topology. The optimal soft threshold power was also graphically validated by the high scale-free topology R ² between log10(k) and log10 (p(k)), where k is the connectivity between genes and p(k) means the probability of connectivity. Fourth, the adjacency matrix was turned into a Topological Overlap Matrix (TOM). The functions “pickSoftThreshold” and “plotDendroAndColors” were conducted to identify the modules and then plot the clustering dendrogram of genes together with dissimilarity on the topological overlap and various module colors. Herein, the parameters “minModuleSize” and “MEDissThres” were set to 60 and 0.3, respectively. Fifth, module-trait correlation analysis was performed to reveal the ME blue module as the most important one related to DCM with HF. Subsequently, the Gene significance (GS) and Module Membership (MM) were calculated by intramodular analysis. Genes with GS > 0.5 and MM > 0.8 were identified as driver genes in the ME blue module for further study (Ai et al., 2020).

The intersection of up-regulated DEGs and driver genes in the ME blue module was identified as feature genes for DCM with HF. Next, the clustering heatmap of the feature gene expression was plotted using the R package “pheatmap”. BPNN, typically composed of processing elements arranged in at least three layers (input, hidden, and output layers), is a multilayer feedforward network trained on an error-back propagation algorithm (Kaundal et al., 2006). Here, BPNN was constructed in the training group based on feature genes and then validated in the test group using the R package “neuralnet”. For BPNN, the input and output layer neurons are equal to the feature gene number and the classification number (control or DCM with HF), respectively. According to common rule-of-thumb methods, the hidden layer neurons were determined and adjusted by 1/3–2/3 of the input layer neurons (Karsoliya, 2012). Accordingly, the neuron numbers of the input layer, hidden layer, and output layer were set as 13, 5, and 2, respectively. The classification accuracies of BNPP in the training and test groups were evaluated by the AUC score of the ROC curve by the R package “pROC”. Additionally, we visualized Pearson’s correlation among the 13 feature genes by the R package “corrplot”.

To narrow down the numbers of feature genes, Lasso was first performed to carry out feature selection. Five machine learning classifiers (RF, GBM, NN, XGBoost, and SVM) were then constructed for further feature selection on the Lasso-selected genes (Han et al., 2014). These models were all constructed through 10-fold cross-validation. The setting paraments of RF, NN, and SVM were 400 decision trees, five neurons in the hidden layer, and radial kernel. And paraments set by default for GBM and XGBoost. By using the R package “glmnet”, Lasso was first applied to reduce gene numbers for screening the feature genes. Then, five machine learning classifiers were established in the training and test groups based on the Lasso-selected feature genes using the R package “caret”. The classification performances of five models in the training and test groups were assessed by AUC scores using the package “pROC”. The residual and feature importance of classifiers were analyzed by the R package “DALEX”. Net Benefit (NB) calculated by Decision Curve Analysis (DCA) and Brier score were also utilized to measure the clinical value and calibration of classifiers, respectively (Fitzgerald et al., 2015; Alba et al., 2017). In general, the classifier with the lowest residual, minimum Brier score, highest AUC, and maximum NB was considered the best predictive classifier for further feature selection. Herein, RF was selected as the best-performing classifier for accurately distinguishing between normal and DCM with HF groups. Using the LASSO-selected genes, the R packages “caret” and “randomForest” were used to build an RF classifier in the training and test groups. The feature importance of RF was investigated by using the Gini gain approach. Feature genes with a mean decrease Gini larger than two were selected as key genes for diagnosing DCM with HF (Tian et al., 2020).

Then, within the training and test groups, we analyzed the expression distribution and diagnostic accuracy of key genes selected by RF. By using the R package “ggpubr”, the boxplots were plotted and demonstrated the differential expression distributions of key genes between control and DCM with HF groups. The diagnostic accuracies of key genes for DCM with HF were evaluated by the AUC score of the ROC curve using the R packages “pROC”. Furthermore, we established a logistic regression model for evaluating the diagnostic accuracy of the key genes’ combination using the R package “glmnet”. Furthermore, the nomogram was used to illustrate and visualize the logistic regression model by the R package “rms”. The calibration curve and NB calculated by DCA were utilized to evaluate the prediction effect of the nomogram. The gene with an AUC score ranging from 0.85 to 1.0 was considered to have excellent specificity and sensitivity to distinguish between the control group and DCM with HF group. The Human Protein Atlas (HPA) database (https://www.proteinatlas.org/) was used to explore the extracellular locations and tissue-specific expression levels of the key genes. Moreover, Drug Signatures Database (DSigDB) (https://maayanlab.cloud/Enrichr/) (Kuleshov et al., 2016) was applied to predict potential therapeutic drugs target for three key genes, with adjusted p < 0.05 as the filtering threshold.

To further confirm the differential expression of key genes in DCM with HF, we performed single-cell RNA-seq analysis on GSE95140, mainly including cell type identification and Pseudotime-ordered analysis. First, we performed quality control and data cleaning using the R package “dplyr”. The screening criteria were set: the number of genes in each cell ranges from 200 to 7,500, and the percentage of mitochondrial genes in each cell is lower than 20%. Second, the function “NormalizeData” of the R package “seruat” was used to perform log2 transformation on the raw count data. Next, the top 2000 hypervariable genes were selected by the function “FindVariableFeatures”. Z-score standardization was conducted using the function “ScaleData”. Third, PCA was carried out on the hypervariable genes in the control and DCM groups. Next, two PCs (control group, all p < 0.01) and four PCs (DCM group, all p < 0.01) were selected. Then, two clusters (control group) and four clusters (DCM group) were obtained and projected onto UMAP plots, respectively. Fourth, the marker genes of cardiomyocyte (NNT2, MYH7), fibroblast (DCN, COL1A1), and endothelial cell (CDH5, VWF) obtained from CellMarker 2.0 database (http://yikedaxue.slwshop.cn/) (Zhang et al., 2019) were used to identify cardiomyocytes in different states. We then explored the expression distribution of key genes and NPPB in the control and DCM groups. NPPB, an important paralog of NPPA with up-regulation in DCM with HF (Barth et al., 2006), was selected as an indicator to reveal whether the key genes showed co-expression and dynamic expression patterns of NPPB. Fifth, the R package “monocle” was used to conduct Pseudotime-ordered analysis using the “DDRTree” method on the DCM group. Using the function “plot_cell_trajectory”, the differentiation trajectory of cardiomyocytes was visualized in the tree diagram. Additionally, we analyzed and visualized the dynamic expression patterns of key genes and NPPB during the cardiomyocyte transitions across the Pseudotime using the function “plot_genes_in_pseudotime”.

Considering the limitation of performing KEGG on small gene sets, GSEA was utilized to explore the potential pathways of key genes. We conducted GSEA on the training group using the R package “clusterProfiler”. The KEGG gene sets annotated for GSEA (c2.cp.kegg.v7.5.1.symbols) were downloaded from the MSigDB database (http://www.gsea-msigdb.org/gsea/msigdb/index.jsp). First, the samples were divided into low- and high-expression groups according to the median expression of NPPA, OMD, and PRELP, respectively. Then, we ranked 21,628 genes in descending order by calculating the log₂FC for each gene. Second, the GSEA was performed on these ranked genes using the KEGG gene sets. The enrichment KEGG pathways with p < 0.05 and |Normalized Enrichment Score (NES)| > 1 were considered significant (Zheng et al., 2022). NES was the primary statistic to examine the enrichment results, the value indicates the relationships between enrichment pathways and key genes. Herein, the intersection of enrichment KEGG pathways of each key gene was determined as the potential pathway in DCM with HF. Subsequently, we carried out Pearson’s correlation analysis between each key gene and the potential pathway-related genes. |R| > 0.3 and p < 0.05 was the cutoff criteria. The intersection of pathway genes significantly correlated with each key gene was considered important and used for further analysis. Next, the Pathview online tool (https://pathview.uncc.edu/home) (Luo et al., 2017) was used to perform Z-score normalization on the expressions of these intersection genes between control and DCM with HF groups and projected their relative expression values onto the KEGG pathway.

Statistical analysis was conducted using R programming language (version 4.2.1). Wilcoxon test was performed to analyze the differential expression of key genes in DCM with HF group versus control group. The diagnostic ability of each key gene was evaluated by the AUC score. P < 0.05 was considered statistically significant.

The gene expression matrices of the training group were obtained after data cleaning, normalization, batch effect removal, and merging, including 20 normal myocardial samples (control group) and 35 myocardial samples derived from DCM patients with HF (DCM with HF group). The PCA score plot showed the relatively consistent clustering of two groups from two components (Figure 1A), which can be used in subsequent analysis due to between-group differences. A total of 218 DEGs (135 up-regulated genes and 83 down-regulated genes) were screened under the threshold (|log₂FC| > 1 and adjusted p < 0.05). The Volcano plot of DEGs was then plotted using the log₂FC and -log10 (adjusted p) (Figure 1B). The top 50 up-regulated genes and the top 50 down-regulated genes were shown in the heatmap (Figure 1C).

The GO, KEGG, and DO enrichment analyses were conducted to annotate the potential functions, pathways, and diseases of 135 up-regulated DEGs, respectively. As shown in Figure 2A, the enrichment GO terms of biological process (BP) included muscle cell migration, striated muscle adaptation, and skeletal muscle adaptation; the enrichment GO terms of cellular component (CC) included collagen−containing extracellular matrix, myofilament, and striated muscle thin filament; the enrichment GO terms of molecular function (MF) included integrin binding, cytokine activity, and extracellular matrix structural constituent. DO enrichment analysis showed that the DEGs were associated with many heart diseases, such as arteriosclerotic cardiovascular disease, coronary artery disease, and myocardial infarction (Figure 2B). Additionally, the enriched KEGG pathways mainly contained muscle contraction, post-translational protein phosphorylation, and regulation of Insulin-like Growth Factor (IGF) transport and uptake by Insulin-like Growth Factor Binding Proteins (IGFBPs) (Figure 2C).The BN plot further demonstrated the enriched gene interactions in muscle contraction (Figure 2D). Furthermore, NPPA acted as the core gene with higher expression and more edges, indicating its crucial role in muscle contraction.

WGCNA was utilized to construct the gene co-expression network and identify the trait-correlated modules in the training group. First, the similarity matrix was calculated and then transformed into the adjacent matrix via the optimal soft threshold power (β = 2) (Figures 3A,B). Besides, the inverse relationship of p(k) and k was observed (Figure 3C), and the scale-free topology R ² of log10(k) and log10 (p(k)) was 0.87 when β = 2 (Figure 3D), which illustrated β = 2 was suitable for constructing a scale-free topology. Second, the TOM was converted from the adjacent matrix and turned into the dissTOM (dissTOM = 1-TOM) for hierarchical clustering. Dynamic tree shearing clustered the similarity genes based on the topological overlap and then divided them into various modules (Figure 4A). Third, we related the modules to clinical trait and further obtained the driver genes in the strongest positive correlation module. As shown in Figure 4B, we ultimately identified the ME blue module strongly positively related to DCM with HF (R = 0.77, p < 0.001), containing 2,411 genes. Fourth, the significant GS distribution across the seven modules was observed (Figure 4C). Additionally, the intramodular analysis demonstrated the close relationship between MM and GS of the ME blue module (R = 0.81, p < 0.001; Figure 4D), thus relieving the module genes highly associated with DCM with HF. Finally, a total of 15 genes with GS > 0.5 and MM > 0.8 (C14orf132, CFH, COL8A1, CTGF, ETNPPL, FIBIN, FRZB, ITGBL1, LTBP2, MFAP4, NPPA, NRK, OMD, PRELP, SFRP4) were identified as driver genes in the ME blue module.

We selected the intersection of 135 up-regulated DEGs and 15 driver genes in the ME blue module as the feature genes of DCM with HF, with a total of 13 genes (OMD, PRELP, NPPA, LTBP2, C14orf132, FIBIN, SFRP4, FRZB, MFAP4, COL8A1, CFH, CTGF, ITGBL1) (Figure 5A). These 13 genes were previously proved to participate in muscle contraction, striated muscle adaptation, and extracellular matrix structural constituent (Figures 2A,B), indicating their potential crucial roles in DCM with HF. Figure 5B showed the differential expressions of feature genes between control and DCM with HF groups. Subsequently, a BPNN with 13 neurons in the input layer (13 feature genes), five neurons in the hidden layer, and two neurons in the output layer (classification numbers: control or DCM with HF) was constructed in the training group and then validated in the test group (Figure 5C). The BPNN showed prediction accuracy with a high AUC score of 0.988 in the training group (Figure 5D). However, the predictive power of BNPP decreased significantly in the test group compared to the training group, merely achieving an AUC score of 0.745 (Figure 5E). These results suggested that the collinearity between the screened 13 feature genes resulted in the constructed BPNN over-fitting for general prediction (Supplementary Figure S2). Thus, feature selection was crucial for further obtaining the key genes for the diagnosis of DCM with HF.

To narrow down the range of feature genes, Lasso within the training group was first performed. As shown in Figure 6A, Lasso retained five genes from 13 feature genes under the optimal log(λ) as 0.0156. Through the 10-fold cross-validation, the minimum binomial deviation with optimal log(λ) was observed, indicating the five genes were optimally selected by Lasso (Figure 6B). Second, five machine learning models (RF, GBM, NN, XGBoost, and SVM) were established on the five LASSO-selected genes and validated their classification performances based on a series of indexes (residual, Brier score, AUC score, NB) in the training and test groups. Strikingly, the RF model showed robust classification ability, with the lowest residual (Figure 6C) and Brier score (Supplementary Table S1), highest AUC achieving 1 (Supplementary Figures S3A,C), and NB under threshold probability (Supplementary Figures S3B,D) in the training and test groups. Thus, the RF classifier was considered the most suitable performance classifier for further feature selection. Moreover, the importance of feature genes selected from five classifiers was illustrated by Root Mean Square Error (RMSE) loss (Figure 6D). Third, we constructed the RF model with 400 decision trees (Figure 6E), and three features with mean decrease Gini coefficients larger than 2 (OMD, PRELP, NPPA) were selected as key genes for diagnosing DCM with HF (Figure 6F).

To test whether NPPA, OMD, and PRELP could serve as the diagnostic biomarkers of DCM with HF, their mRNA expression levels and diagnostic abilities were then analyzed. As shown in Figures 7A,B, the over-expression levels of NPPA, OMD, and PRELP were observed in DCM with HF in the training group and test group. Furthermore, the AUC scores of NPPA (AUC = 0.967), OMD (AUC = 0.984), PRELP (AUC = 0.974), and their combination (AUC = 0.986) for diagnosing DCM with HF were all larger than 0.9 in the training group (Figure 7C). The AUC scores of NPPA (AUC = 0.863), OMD (AUC = 0.722), and PRELP (AUC = 0.816) decreased in the test group (Figure 7D). However, the AUC score of the combination of NPPA, OMD, and PRELP markedly achieved 0.922 within the test group. These results suggested NPPA, OMD, and PRELP be promising diagnostic biomarkers for DCM with HF. To measure the combined diagnostic accuracies of NPPA, OMD, and PRELP, a logistic regression model was constructed and displayed in the nomogram (Supplementary Figure S4). Interestingly, the combination of NPPA, OMD, and PRELP improved the diagnostic effect within the training and test groups. We also observed the heart muscle-specific expressions of NPPA, OMD, and PRELP, and their extracellular locations were all predicted to be secreted in the HPA database (Supplementary Figure S5), which provided strong evidence that NPPA, OMD, and PRELP served promising biomarkers of DCM with HF.

Herein, we used DSigDB to predict the potential therapeutic drugs associated with NPPA, OMD, and PRELP. A total of 17 target drugs related to NPPA were finally predicted, but not associated with OMD and PRELP. And the specific information on the top 10 predicted drugs was shown in Supplementary Table S2. Among these drugs, the antihypertensive efficacy of alprostadil, labetalol, felodipine, and irbesartan has been reported (Johnson and Fugman, 1983; MacCarthy and Bloomfield, 1983; Cattaneo et al., 2003). Additionally, previous studies have shown that irbesartan, furosemide, and Bonuten can treat HF (Johnson and Fugman, 1983; Swedberg et al., 1987; Thompson et al., 1999). Several pieces of evidence have supported that these predictive drugs targeting NPPA could cure some cardiovascular disorders such as hypertension and HF, thus indicating their value in the treatment of DCM with HF.

To investigate the expression patterns of NPPA, OMD, and PRELP in DCM patients, we conducted single-cell RNA-seq analysis on GSE95140. The quality control, data cleaning, and PC selection were shown in Supplementary Figure S6. The entire cell profiles of the control group and DCM group were classified into 2 clusters (Figure 8A) and 3 clusters (Figure 8C) and projected onto UMAP plots, respectively. We identified the cell type using specific gene markers (NNT2 and MYH7 for cardiomyocyte; DCN and COL1A1 for fibroblast; CDH5 and VWF for endothelial cell) in the control and DCM groups (Supplementary Figure S7). We annotated one cardiomyocyte type and three distinct cardiomyocyte types in the control group and DCM group, respectively (Figures 8A,C). Then, we analyzed the expression of NPPA, OMD, PRELP, and NPPB in cardiomyocytes of control and DCM groups. As shown in Figures 9B,D, higher expression levels of NPPA, OMD, PRELP, and NPPB were observed in most cardiomyocytes in DCM compared to the control group. Interestingly, we noticed that type 2 of the DCM group showed high expression of NPPA and NPPB, suggesting that cardiomyocytes with higher degrees of HF in these types. We also performed Pseudotime-ordered analysis on the DCM group to explore the expression changes of NPPA, OMD, PRELP, and NPPB during the cardiomyocyte transitions along the Pseudotime. Figure 9E showed the differentiation trajectory path of cardiomyocytes ordered by different types. Figure 9F further displayed the stable expression patterns of NPPA, OMD, PRELP, and NPPB in diverse cardiomyocytes by Pseudotime. Notably, the dynamic expression changes of NPPA and NPPB were similar, indicating an expression linkage between NPPA and NPPB in DCM. Overall, these results showed NPPA, OMD, and PRELP were enriched in cardiomyocytes of DCM and displayed stable expression patterns, which fluctuated slightly with differential cell states.

To gain insight into the potential pathways of NPPA, OMD, and PRELP, GSEA was performed for each gene. As shown in Figure 9A, the top three positively significant KEGG pathways for NPPA were “citrate cycle (TCA cycle)”, “basal cell carcinoma”, and “TGF-β signaling pathway”. The top two positively significant KEGG pathways for OMD were enriched (Figure 9B), including “basal cell carcinoma” and “TGF-β signaling pathway”. The top three positively significant KEGG pathways for PRELP shown in Figure 9C were “citrate cycle (TCA cycle)”, “basal cell carcinoma”, and “TGF-β signaling pathway”. Interestingly, the terms “basal cell carcinoma” and “TGF-β signaling pathway” were enriched in the KEGG pathways of NPPA, OMD, and PRELP, and acted as the positive pathways. The TGF-β signaling pathway was considered the most likely mechanism involved by NPPA, OMD, and PRELP. We then analyzed Pearson’s correlation of NPPA, OMD, and PRELP and the TGF-β signaling pathway gene sets (Supplementary Table S3), respectively. And 16 NPPA-correlated, 20 OMD-correlated, and 26 PRELP-correlated genes were obtained (Supplementary Table S4). The intersections of these genes were INHBA, ID4, TGFB3, SMAD7, and SMAD9 (Figure 9D). Then, the expression levels of INHBA, ID4, TGFB3, SMAD7, and SMAD9 between control and DCM with HF groups were normalized and visualized in the TGF-β signaling pathway (Figure 9E). Additionally, we observed that ID4, SMAD7, and SMAD9 were up-regulated in the TGF-β signaling pathway. While INHBA and TGFB3 were down-regulated. Accordingly, we inferred that the TGF-β signaling pathway positively associated with NPPA, OMD, and PRELP might play a crucial role in the process of DCM with HF.