The DCM RNA expression datasets were collected from the online GEO database (http://www.ncbi.nlm.nih.gov/geo/). The keywords “dilated cardiomyopathy”, “Homo sapiens”, and “expression profiling by array” were used on the initial search, and 184 DCM related studies were found. Then the following criteria were used to further screen datasets: 1) The study includes DCM case vs healthy control; 2) tissue samples obtained from left ventricle; 3) sample size was bigger than 10. Three datasets qualified for the above criteria and performed on the same platform were combined for analysis. Another two datasets derived from other platforms were used as validating datasets. The processed data of GSE1145 (platform: GPL570, including 11samples of control and 12 samples of DCM), GSE21610 (Schwientek et al., 2010) (platform: GPL570, including eight samples of control and 21 samples of DCM), GSE29819 (Gaertner et al., 2012) (platform: GPL570, including six samples of control and seven samples of DCM), GSE5406 (Hannenhalli et al., 2006) (platform: GPL96, including 16 samples of control and 86 samples of DCM) GSE57338 (Liu et al., 2015) (platform: GPL11532, including 132 samples of control and 82 samples of DCM) were downloaded as expression matrix with R package ‘GEOquery’.(Davis and Meltzer, 2007). The mRNA expression profiles of controls and targeted patients were extracted and were performed log2 transformation before further analysis (Only if they have not be log2 transformed).

The ‘limma’ (Ritchie et al., 2015) and ‘sva’ (Parker et al., 2014) packages in R software (R version 4.0.3) were used to correct intra- and inter-batch effect. Prcomp function was used to perform principal component analysis (PCA). Next, probes annotations were performed. Probes annotated to >1 gene were removed. For multiple probes annotated to the same gene, the first one appeared retained. Finally, the ‘limma’ package was also used for DEGs identification between different groups with cut-off values of adjusted p-value < 0.05 and |fold change| ≥ 1.5. In addition, DEGs were also identified by applying robust rank aggregation (RRA) algorithm (Kolde et al., 2012) with the same criteria. Venn diagrams (http://bioinformatics.psb.ugent.be/webtools/Venn/) were used to summarize the overlapping DEGs between ‘sva’ and RRA algorithms. Boxplot function, ‘ggplot2’ and ‘pheatmap’ packages were then used to plot gene expression boxplot, volcano plot and heatmap.

To further analyze the functions of DEGs, the R package ‘clusterProfiler’ (Yu et al., 2012) was used to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. The cut-off values for GO and KEGG were set as p < 0.05. The ‘enrichplot’ package was used to draw dot plots for the results of functional enrichment analysis.

xCell (Aran et al., 2017) is a new gene signature-based method to estimate the content of immune and stromal cells. It was validated using extensive in-silico simulations and cytometry immunophenotyping. XCell package was applied to normalized merged data to portray the cellular heterogeneity landscape of left ventricular expression profiles. We compared the cell distribution differences between the two groups through t-test, and the cutoff value was set as p < 0.05. The results were significantly different between the two groups, categorized based on their traits into three categories: “lymphoid and myeloid cells”, “stem cells” and “stromal cells and others.” These were visualized by using ‘ggplot2’ and ‘ggpubr’ packages. Correlation analyses of immune cell subtypes and hub genes were performed with ‘psych’ and ‘corrplot’ packages. Pearson correlation coefficient was used to assess the strength of correlation. Hub genes whose absolute correlation coefficient with immune cells >0.5 and p-value < 0.05 are selected to further study.

Search Tool for the Retrieval of Interacting Genes (STRING, https://string-db.org) is a webtool that provide validated and predicted information of PPIs (Szklarczyk et al., 2019). The list of DEGs was uploaded into STRING website to detect significant protein interactions with minimum interaction score >0.7. The network was then exported and visualized by Cytoscape 3.7.1 software (Cline et al., 2007). The CytoHubba plugin (Chin et al., 2014) was used to calculate hub genes with a high degree. The results were directly visualized by Cytoscape. Additionally, the MCODE plugin (Bader and Hogue, 2003) was used to identify highly interconnected clusters with the cutoff parameters set as follows: degree cutoff = 2, node score cutoff = 0.2, k-Core = 2, Max. Depth = 100. The results were further screened with criteria set as MCODE score >4 and nodes number >5. Gene ontology biological process enrichment analysis was performed on the significant modules.

Hub genes with correlation coefficient >0.5 and p < 0.05 with immune cells were selected. The expression profiles of the selected biomarkers were visualized in boxplot and validated in another two datasets GSE5406 and GSE57338. Then, a diagnosis model combining the selected biomarkers was constructed in the merged dataset by logistic regression using ‘glm’ function and verified in GSE5406 and GSE57338. Receiver-operating characteristic (ROC) curves were used to assess the discrimination ability of the key genes and the diagnosis model.

Expression profiles of healthy controls and DCM patients from GSE1145, GSE21610 and GSE21869 were combined into a merged dataset. The merged dataset contains 25 healthy controls and 40 DCM patients. As shown in Figures 2A,C, the merged dataset had strong batch effect. After normalization, it was effectively removed (Figures 2B,D). Then in the differential expression analysis, 1,005 DEGs were recognized in the merged dataset, which includes 385 down-regulated genes and 620 up-regulated genes by integration method of ‘sva’ (Figure 2E). 179 DEGs including 101 up-regulated genes and 78 down-regulated genes were identified by integration method of RRA. Venn diagrams depict DEGs across two different integration methods (Supplementary Figure S1). DEGs obtained by ‘sva’ method were used in the following analyses. All of the DEGs were displayed in heatmap (Figure 2F).

To assess the functions of DEGs, enrichment analyses of GO and KEGG were performed. The categories of GO analysis include biological process (BP), cellular component (CC) and molecular function (MF). The leading 10 enriched terms of each GO categories and KEGG with p value <0.05 were visualized in Figure 3.

The mainly enriched BP terms included epithelial cell proliferation and its regulation, extracellular matrix and structure organization, cell chemotaxis, leukocyte migration and regulation, and inflammatory response (Figure 3A). The mainly enriched MF terms were extracellular matrix structural constituent, glycosaminoglycan, collagen and fibronectin binding, G protein-coupled purinergic nucleotide receptor activity (Figure 3B). The mainly enriched CC terms contained collagen-containing extracellular matrix, endoplasmic reticulum lumen, external side of plasma membrane, myofibril, sarcomere, basement membrane, and collagen trimer (Figure 3C). In whole, GO results indicated that the gene function of DEGs mainly associated with both extracellular structure reorganization and fibrosis, immune, and inflammatory abnormalities.

In KEGG pathway enrichment analysis, DEGs mainly enriched in PI3K-Akt signaling pathway, Cytokine-cytokine receptor interaction, AGE-RAGE signaling pathway in diabetic complications, Phagosome, viral protein interaction with cytokine and cytokine receptor, and HIF-1 signaling pathway (Figure 3D).

xCell was used to estimate the cell composition heterogeneity of left ventricle between DCM and controls. As shown in Figure 4, 18 cell types were significantly changed in DCM cardiac tissue compared to control samples, among which the scores of CD8⁺ T-cells, cDC, adipocytes, fibroblasts, and smooth muscle in DCM were significantly increased, while the scores of macrophages M1, Th1 cells were significantly decreased.

Moreover, the correlation among immune cells were calculated by using the Pearson’s correlation coefficients. As shown in Figure 5, iDC has the highest positive correlation with DC (Pearson’s coefficient = 0.89). The correlation between macrophages and macrophages M1 was the second strongest positive (Pearson’s coefficient = 0.78). Additionally, aDC, DC, CD8⁺ Tcm, macrophages, monocytes, and macrophages M1 had strong correlation coefficient with most of the remaining immune cells.

As shown in Figure 6A, there were 633 edges among 350 proteins in the PPI network. The top 10 genes with the highest degree in the above network were selected as the hub genes. These are STAT3, IL6, CCL2, PIK3R1, ESR1, CCL5, IL17A, TLR2, BUB1B, and MYC (Figure 6B). The expression of the top 10 hub genes across all datasets was displayed in supplementary material. Then, the correlation coefficients between 10 hub genes and significantly changed immune cells were calculated. As shown in Figure 6C, CD8⁺ T-cells and macrophages M1 had significant correlation with most of the hub genes. Moreover, with the screening rules as |Pearson’s coefficient| > 0.5 and p < 0.05, we obtained three immune-related hub genes with potential diagnostic value: TLR2, CCL2 and CCL5.

Additionally, the densely connected modules were identified by using MCODE plug-in. Nine modules were obtained from the PPI network (Supplementary Figure S2). The most enriched GO BP terms of each module were listed in Table 1. Module_1 was enriched in negative regulation of transcription from RNA polymerase II promoter. Module_2 was enriched in sister chromatid cohesion. Enriched BP terms of Module_3 and four were associated with metabolism. Module_5 was enriched in proteolysis. Module_6 and eight were enriched in extracellular matrix organization. Module_9 was enriched in immune response.

As it is shown in Figure 7A, CCL2 and TLR2 were significantly downregulated, while CCL5 was significantly upregulated in the merged dataset and verified in another two datasets GSE5406 and GSE57338 (Figures 7B,C). Next, ROC curve analysis was performed to verify the diagnostic value of the selected biomarkers. In the merged dataset, the ROC curves of the three biomarkers revealed high diagnostic value for DCM (Figure 8A). The AUC (area under the curve) scores of CCL2, CCL5, TLR2 were 0.963, 0.803, 0.766 respectively. In verification dataset GSE5406, the AUC scores of CCL2, CCL5, TLR2 were 0.733, 0.821, 0.794, respectively (Figure 8C). In verification dataset GSE57338, the AUC scores of CCL2, CCL5, TLR2 were 0.738, 0.831, 0.836, respectively (Figure 8E). Finally, we used the three biomarkers to construct a diagnosis model by logistic regression and visualized in ROC curves. The AUC score of the diagnosis model was 0.981, 0.867 and 0.946 in the merged dataset, GSE5406 and GSE57338, respectively (Figures 8B,D,F).