The raw data of two eligible microarray datasets (GSE120895 and GSE17800) based on platform GPL570 were downloaded from the GEO database (87 patients with DCM and 16 controls after data merging). Moreover, GSE19303 (73 patients with DCM and 8 controls) was used as the validation set. The Limma package was used to screen DEGs between DCM and controls, followed by data normalization [13]. Bias and variability of the datasets were removed using the Combat function in the sva package. p 0.05 adjusted by the false discovery rate (FDR) was considered significant. The data processing procedure of our research is illustrated in the workflow (Fig. 1).

Fig. 1.

The workflow of our research. DCM, dilated cardiomyopathy; GEO, Gene Expression Omnibus; GO, Gene Ontology; PPI, protein-protein interaction; CMAP, Connectivity Map; KEGG, Kyoto Encyclopedia of Genes and Genomes; GSEA, gene set enrichment analysis; TF, transcription factor; DEG, differentially expressed gene; WGCNA, weighted gene coexpression network analysis; IHC, immunohistochemical; qPCR, quantitative polymerase chain reaction.

A coexpression network of 884 genes with commonly upregulated or downregulated expression in DCM/control was constructed using WGCNA, which is a widely used systems biology approach that helps identify the relationship between genes and disease phenotypes. The main processes were as follows: (1) the hclust function was used for hierarchical clustering analysis; (2) a power of $\beta$ = 7 was selected according to the scale-free topology criterion; (3) genes with coexpression relationships were grouped through gene network connectivity, the minimum module size was set to 30, and each module was assigned a unique color label; (4) the relationship between modules and phenotypes was calculated to identify biologically meaningful modules with Pearson correlation analysis; (5) GO annotation and KEGG pathway enrichment analysis were performed for these functional modules; and (6) a PPI network was constructed.

In WGCNA, module membership (MM) is defined as the correlation of the module eigengene and the gene expression profile, while gene significance (GS) is defined as the correlation between the gene and the trait. Among the modules of interest in our study, GS $>$0.3 and MM $>$0.8 were defined as hub genes in the candidate gene modules.

The least absolute contraction and selection operator (LASSO) method [14], which is suitable for high-dimensional data restoration, was used to select the optimal risk factor prediction characteristics from gene datasets. The A radiomics score (Rad-score) of each gene was calculated by selecting a linear combination of features, which are weighted by their respective coefficients. Then, the variables selected by the LASSO method were used to obtain the final factors for establishing the model. Harrell’s C index and the area under the curve (AUC) were measured to quantify the discrimination performance of the model based on the experimental set and validation set.

GSEA is a computational method used to assess whether a predefined set of genes displays statistical significance and consistency differences between two biological states [15]. The c5.all.v7.4.entrez and c2.cp.kegg.v7.4.entrez datasets in the MsigDB database were used as reference gene sets, and the clusterprofile package was utilized to perform GSEA with integrated gene expression data. p 0.05 was considered significant.

The promoter sequences of target genes were searched using the University of California Santa Cruz (UCSC) and National Center for Biotechnology Information databases (NCBI). Potential binding transcription factors were searched using the JASPAR database and further screened according to transcription direction and related levels. Gene expression profiling interactive analysis (GEPIA) was used to verify the correlations and identify potential transcription factors. Finally, the structural information combined with the interval key structural domain was presented using UniProt.

The CMap database establishes links between genes, compounds and diseases based on similar and opposite gene expression profiles. In this study, the DEGs of the DCM and control groups were grouped according to expression differences. Then, the DEGs were loaded to the “Query” page. In this study, drugs with connectivity scores $>$90 or –90 and p values less than 0.05 were selected as drug candidates. The 2D and 3D structures of the candidate compounds were obtained from PubChem.

The study was approved by the Medical Ethics Committee of Nanjing Medical University, and a total of 16 patients were recruited to provide blood samples, and informed consent was obtained from patients before participation. Expression levels of 4 hub genes were verified in 8 DCM blood samples and 8 normal blood samples (Supplementary Table 1). Total RNA was extracted using TRIzol reagent according to the manufacturer’s instructions. Total RNA was reverse transcribed into cDNA by PrimeScript RT Master Mix (TaKaRa, Japan) after measuring the corresponding concentration. Then qRT-PCR was performed using Power SYBR Green PCR Master Mix (No. A25742; Thermo Fisher Scientific, Waltham, MA, USA). Finally, the relative expression levels of miRNAs and target genes were calculated according to the 2 -$\mathrm{\Delta }$$\mathrm{\Delta }$Ct method. GAPDH levels were used to normalize mRNA expression levels. The procedure for PCR is shown below: In brief, samples were incubated at 95 °C for 5 min, followed by 40 cycles at 95 °C for 10 s and finally at 60 °C for 20 s.The sequences of the primers are shown below: 5′-CAACGAATTTGGCTACAGCA-3′ and 5′-AGGGGAGATTCAGTGTGGTG-3′ for GAPDH, 5′-GGACATCGACACCATCTCCC-3′ and 5′-TAGTTCCGCCAGTAGGACCG-3′ for ARL2; 5′-TCCCTTTGCTATTCACCGTCA-3′ and 5′-CACCCGACATTCCCAGCATC-3′ for MLF2; 5′-GCCGGATGGCAGAATGGAG-3′ and 5′-GGAGGGAGACTAGGACGATGG-3′ for GYPC; 5′-ATGAGTGCGGAAGTGAAGGTG-3′ and 5′-GCTCTCTAACATTGGGCATGTC-3′ for COPS7A.

The Human Protein Atlas was used to validate the immunohistochemistry of potential target genes. This database facilitates the systematic study of transcriptome and gene pathology expression of coding genes in different tissue types. Staining of core gene proteins in human myocardial tissue based on immunohistochemical techniques. In addition, tissue types and staining levels were retrieved from the database to analyze the quality of the data to interpret the results [16, 17].

We retrieved the 3D structure of the receptor from the Uniprot and RCSB protein databases, and the corresponding simplified molecular-input line-entry system (SMILES) for the receptor ligand was obtained from the PubChem database. After the ligand energy was adjusted to a minimum by using Chem 3D software 18.0 (PerkinElmer, Waltham, MA, USA), the crystals were imported into Pymol 2.4.0 software (Schrödinger, L. & DeLano, W., CA, USA) for dehydration, hydrogenation and ligand separation, followed by AutoDockTools-1.5.7 [18] to construct docking grid boxes for each target. Docking is done with Autodock Vina 1.1.2 (The Scripps Research Institute, San Diego, CA, USA) [19].

To identify DEGs associated with DCM, we screened the integrated normalized data and obtained 12,748 genes after adjusting for batch effects. As shown in the volcano map (Fig. 2), compared to controls, there were 75 DEGs in the DCM sample, of which 39 were upregulated and 36 were downregulated. The 20 most upregulated and 20 downregulated genes in DCM are visualized in the heatmap.

Fig. 2.

The result of DEGs identification. (a) Principal component analysis (PCA) before batch effects removement. (b) PCA after batch effects removement. (c) A heatmap of 20 most up-regulated and 20 most down-regulated genes. (d) Volcano plot of the DCM-Control. Abbreviations: DCM, dilated cardiomyopathy; logFC, log2 fold-change.

It is reasonable to consider that certain genes with similar expression patterns may perform their functions as a global network and participate in similar biological processes. The 883 DEGs with absolute deviation were used for WGCNA. A total of six gene modules were obtained, which are represented by branches of the clustering tree, and different colors are shown in Fig. 3. In addition, modular enrichment analysis showed that primary enrichment in biological process (BP) corresponding to each module included mitotic DNA damage checkpoint, neutrophil degranulation, protein-containing complex disassembly and muscle system process. Primary enrichments in cell component (CC) involved integral components of endoplasmic reticulum membrane, tertiary granule and mitochondrial inner membrane. Primary enrichments in molecule function (MF) consisted of E-box binding and structural constituent of ribosome. Primary enrichments in the KEGG pathway were chemical carcinogenesis - reactive oxygen species and diabetic cardiomyopathy (Fig. 4). The PPI network of DEGs was constructed with STRING and visualized with Cytoscape (3.8.2). As shown in the Supplementary Fig. 1, the PPI network contained nodes and edges in a circular distribution. The nine most important modules were identified by the Cytoscape plugin MCODE.

Fig. 3.

Construction of gene co-expression modules. (a,b) Analysis of network topology for various soft-thresholding powers. (c) Hierarchical cluster dendrogram of DCM-relate genes based on one dissimilarity measure. (d) Hierarchical cluster heatmap of the adjacencies in the eigengene network. (e) Topology overlay heat map: Both rows and columns indicate individual genes, and dark yellow and red indicate a high degree of topological overlap. (f) Module-phenotype associations. Each row corresponds to a module eigengene, and each column corresponds to a clinical feature. Abbreviations: DCM, dilated cardiomyopathy; BMI, body mass index; EF, ejection fraction; Lvidd, left ventricular internal diameter at end-diastole.

Fig. 4.

GO annotation and KEGG pathway enrichment analysis. (a) Biological process of module. (b) Cell component of module. (c) Molecule function of module. (d) KEGG pathway of module. Abbreviations: GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Gene and Genomes.

Correlations between the modules of the DCM and control groups and the phenotypes of the six clinics (group, age, sex, body mass index, ejection fraction, left ventricular internal diameter at end-diastole) were calculated, and their corresponding p values were analyzed for magnitude. Except for the gray module, the blue module (r = –0.14) had the strongest negative correlation, while the turquoise module (r = 0.16) had the strongest positive correlation (Fig. 3).

Based on GS and MM values (Supplementary Fig. 2), 12 genes in the yellow module were identified as hub genes RPS15, OXA1L, MLF2, GYPC, GNB2L1, EMC4, COPS7A, BANF1, ATP5G2, ARL2, AP2S1 and AK2. Through LASSO regression, four variables, GYPC (glycophorin C), MLF2 (myeloid leukemia factor 2), COPS7A (COP9 signalosome subunit 7A) and ARL2 (ADP ribosylation factor like GTPase 2), were ultimately identified and used to construct the model. ROC curve analysis was performed on the regression model to predict DCM in the training set with an area under the curve (AUC) of 0.853. To further test the diagnostic efficacy, we used the GSE19303 datasets as the validation set, and the AUC was 0.89, indicating that the model gene had a high diagnostic value. The C-index for the prediction model of the cohort was 0.90, and it was 0.85 through bootstrapping validation, demonstrating the model’s good discrimination. In addition, we found that there were significantly upregulated in the expression levels of GYPC (p 0.0001), COPS7A (p = 0.011), ARL2 (p 0.0001) and MLF2 (p 0.0001) in validation set. These findings suggest that these four genes are highly associated with the occurrence of DCM and could be used as further detectable biomarkers. Moreover, We compared the mRNA expression levels of the 4 hub genes in DCM tissues and paired normal tissues. As shown in Fig. 5, compared with their expression in paired normal samples, the expression of ARL2 (p = 0.012), MLF2 (p 0.0001), GYPC (p 0.0001) and COPS7A (p = 0.026) was significantly increased in DCM samples In addition, we observed hub gene protein expression changes in myocardial tissue by immunohistochemistry of tissue sections from the Human Protein Atlas database (Fig. 6).

Fig. 5.

Screening and verification of key genes. (a,b) LASSO model. (c,d) ROC curves for training set and validation set. (e,f) Expression of hub genes in the validation set and qPCR. Abbreviations: LASSO, least absolute shrinkage and selection operator; ROC, receiver operator characteristic; qPCR, quantitative polymerase chain reaction.

Fig. 6.

Immunohistochemistry of the hub genes based on the Human Protein Atlas database (HPAD). (a) Protein levels of ARL2 in myocardial tissue. (b) Protein levels of COPS7A in myocardial tissue. (c) Protein levels of MLF2 in myocardial tissue. (d) Protein levels of GYPC in myocardial tissue.

After validating the JASPAR preselected transcription factors by gene expression correlation (p 0.05 and R $>$ 0.5), the following results were obtained that the transcription factor upstream of ARL2 was CCCTC-binding factor (CTCF), GYPC was zinc finger protein 384 (ZNF384), COPS7A was zinc finger and BTB domain containing 7A (ZBTB7A) and MLF2 was Sp2 transcription factor (SP2) (Fig. 7).

Fig. 7.

Transcription factor validation and binding threshold. (a–d) Scatter plots of gene expression correlations corresponding to the most likely transcription factors bound by the hub genes, respectively. (e–h) Binding sites for transcription factors.

To find drugs for DCM, we searched the Cmap database and subsequently obtained seven drug candidates, namely, trihexyphenidyl, meclofenamic acid, daunorubicin, simvastatin, doxorubicin, dirithromycin and spiperone (Table 1). The 2D (Supplementary Fig. 3) and 3D structures of these drugs were provided by PubChem, and their corresponding active human validated targets were identified. The spatial structure of the gene and the spatial structure of the drug can be used for further molecular docking simulations to seek possible mechanisms of occurrence.

Molecular docking simulations were used to delve into the possible therapeutic mechanisms of these drugs. The binding energy between two counterparts was calculated to predict their affinity. Binding energies below 0 indicated that the two molecules bind spontaneously, with smaller binding energies leading to a more stable conformation. The 3D structure of MLF2 was not available in the Protein Data Bank and ligands corresponding to GYPC and COPS7A were not suitable for molecular docking simulation. The binding energies of ARL2 with Trihexyphenidyl, Meclofenamic acid, Daunorubicin, Simvastatin, Doxorubicin, Dirithromycin and Spiperone were –6.0, –6.7, –7.2, –7.6, –7.4, –6.7 and –5.0 kcal/mol, respectively. Fig. 8 detailed the local structure of molecular docking. It has been illustrated that Doxorubicin acts as an induction target for DCM as positive controls [20]. From the point of view of binding energy this suggested that Simvastatin may be a target for the treatment of DCM (the binding energy of the protein is lower than that of the positive control).

Fig. 8.

The results of the molecular docking simulations. (a) ARL2 with Trihexyphenidyl. (b) ARL2 with Meclofenamic acid. (c) ARL2 with Daunorubicin. (d) ARL2 with Simvastatin. (e) ARL2 with Doxorubicin. (f) ARL2 with Dirithromycin. (g) ARL2 with Spiperone.

Gene set enrichment analysis showed that, compared to control samples, the DCM group was significantly enriched in cell component (CC) such as pseudopodium, lumenal side of membrane, respiratory chain complex IV, MHC protein complex and high-density lipoprotein particle; in molecule function (MF) such as chemokine receptor binding, chemokine activity, glutathione binding, oligopeptide binding, MHC class II protein complex binding and MHC protein complex binding; in signaling pathways such as Chemical carcinogenesis - DNA adducts, Drug metabolism - cytochrome P450, Glyoxylate and dicarboxylate metabolism and Asthma; and in reactome such as Chondroitin sulfate biosynthesis, CS/DS degradation, Defective B3GAT3 causes JDSSDHD and Defective B3GALT6 causes EDSP2 and SEMDJL1 (Fig. 9 and Table 2).

Fig. 9.

GSEA. (a) Cell component of gene set.(b) Molecule function of gene set. (c) KEGG pathway of gene set. (d) Reactome of gene set. Abbreviations: GSEA, gene set enrichment analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes.