The single-cell dataset GSE180678 was obtained from the GEO (https://www.ncbi.nlm.nih.gov/gds) database, which contained a total of 1 sample from patients with ischemic cardiomyopathy. The gene expression profile dataset GSE182923 containing 19 healthy tissues and 19 patients with MI was also downloaded (14). For the single-cell dataset GSE180678, we first downloaded the raw data for quality control and data filtering using the R package Seurat (15). For the expression profile dataset, we first obtained the annotation information of the probes, mapped the probes to genes, removed multiple matches, used the median value as the gene expression when multiple probes matched to a gene, and finally obtained the gene expression profile. In addition, MI-related genes were acquired from the DisGeNET database (http://www.disgenet.org/) and the Comparative Toxicogenomics Database (http://ctdbase.org/). We selected the set of genes common to both databases as highly reliable MI-associated genes.

To obtain reliable cell subpopulations, we first filtered the single-cell data using the R package Seurat (15) for data processing, setting each gene to be expressed in a minimum of 3 cells, and each cell to express at least 250 genes. The percentage of mitochondria and rRNA was calculated by the PercentageFeatureSet function and ensured that each cell expressed more than 500 genes and <4,000 genes, with <30% mitochondrial content and at least 100 unique molecular identifiers (UMIs) per cell. Then, the data were normalized by log-normalization, and the FindVariableFeatures function was used to find highly variable genes. All the genes were scaled using the ScaleData function, and principal component analysis (PCA) downscaling was performed. Finally, the cells were clustered using the FindNeighbors and FindClusters functions (set Resolution = 0.1) to obtain cell subgroups, and the cells were annotated by annotation.

The marker genes for each cell subpopulation were identified using the R package FindAllMarkers function, and the gene expression profile dataset was re-evaluated using the CIBERSORT (16) algorithm based on the expression of marker genes for each cell subpopulation as a background to determine the components of individual cell subpopulation in different samples in the expression profile.

Seven genes related to drug activity were first obtained from the database ChEMBL database (https://www.ebi.ac.uk/chembl/) (17) and further intersected with MI-related genes. Then, these genes were mapped to the String (https://cn.string-db.org/) database (18), and a confidence score of 0.4 was set to obtain their interactions, and after filtering the nodes with fewer edges in the network, the potential set of key genes with a high degree of interaction was finally obtained.

To elucidate the multiple mechanisms of action of these drugs on MI at the systemic level, GO BP, cellular component (CC), molecular function (MF), and KEGG pathway enrichment analyses were performed using WebGestaltR (19) on the potential key gene sets. We set p < 0.05 as the significance threshold, and visualized the most significant top 10 pathways and GO Term using the R the software package ggplot2.

We performed a module-based network analysis using Cytoscape (20). Specifically, we obtained gene interaction data from the String database, set a confidence score of 0.4, and used Cytoscape to visualize and analyze the topological properties of the network. We further used the MCODE plugin (21) to find tightly connected protein clusters in the target network to obtain network modules.

Based on the gene expression profile dataset, 10 immune cell infiltration analyses were performed for each sample using the R package MCPcounter (22) to obtain 10 immune infiltrating cell scores for each patient. The Pearson correlation coefficients and significance of gene expression with immune cells were further calculated by the R package Hmisc package and visualized using the R package ComplexHeatmap (23).

For key genes, because drugs targeted by these genes may have a greater impact on the developmental process of MI, we obtained drug target datasets from the Drugbank database (https://go.drugbank.com/). A protein–protein interaction network (PPI) was based on the String database by setting a threshold score of 400 to create a drug-PPI network, in which we calculated the proximity of drugs and MI disease. Here, we were able to calculate the drug-PPI network with a given S (the set of MI-related genes), D (the degree of nodes in the set of B MI-related genes in the PPI), T (the set of drug target genes), and the distance d(s, t) as the shortest path between node s and node t (where s ∈ S, the MI-related genes; t ∈ T, the drug target genes). The formula is shown below:

Here, ω referred to the weight of the target gene. If the target gene was a gene in the set of MI-associated genes, it was calculated by ω = –ln (D+1 ), otherwise ω = 0.

A simulated reference distance distribution corresponding to the drug was generated. Briefly, we randomly selected a set of protein nodes in the network as stimulated drug targets with the same number of nodes as the target size (denoted as R). Then, the distance d(S,R) between these simulated drug targets (representing the simulated drug) and the associated genes was determined, and the simulated reference distribution was developed after 10,000 random repetitions. Meanwhile, we used μ_{d(S, R)} and σ_{d(S, R)} with references to the mean and standard deviation (SD) corresponding to the actual observed distance and converted it to a normalized score, which was proximity degree z:

We used the random data obtained from the reference to perform multiple hypothesis testing, and selected drugs with close distances and false discovery rate (FDR) <0.001 as the set of drug candidates related to the genes obtained from the analysis.

We used Autodock Vina software for molecular docking. First, AutoDockTools 1.5.6 (https://autodock.scripps.edu/) was used to process receptor proteins and small molecule ligands, such as adding polar hydrogen and adding charge. In molecular docking, the coordinates of the grid in each direction of XYZ were −10, 18.5, and 13.5, respectively. The length of the grid in each direction of XYZ was 20 Å. The Lamarckian algorithm was used to identify the most binding mode of the ligand molecule. The exhaustiveness was set to 8, the maximum number of conformations output was set to 10, and the maximum energy range allowed was set to 3 kcal/mol. Pymol was used for the processing of the resultant plots. In addition, molecular dynamics simulations of 100 ns were performed using the Gromacs2019 package (https://manual.gromacs.org) to assess the binding stability of the receptor-ligand complex. The CHARMm36 force field was used in the molecular dynamics simulations. The str files of the ligands were obtained using the CHARMM General Force Field (CGenFF) program. The system was solventized in a dodecahedral box in a TIP3P water molecule, and sodium and chloride ions were added to neutralize the charge of the system at a concentration of 0.154 M. The energy minimization of the solventized system was performed using the steepest descent algorithm with a cutoff of 5,000 steps. The bond length of the covalent bonds was limited using the LINCS algorithm. The PME algorithm was introduced to calculate overall electrostatic interactions. Subsequently, NVT and NPT simulations were performed for 100 ps, at constant temperature (300 K) and pressure (1 bar), with the confined atoms of the compound equilibrating the system at their initial coordinates. The product MD ran for 100 ns with a time step of 2 fs. The root mean square deviation (RMSD) values of the ligands were calculated using the Gromacs built-in tool.

Single-cell filtration analysis and the PercentageFeatureSet function yielded 1,124 cells. In Supplementary Figure S2A, UMI and the amount of mRNA were significantly correlated. After quality control, the distribution of mRNA/UMI/mitochondrial content/rRNA content of the samples was uniform (Supplementary Figures S2B,C). A “ScaleData” function was performed to scale all genes extracted from the scRNA-seq dataset GSE180678 and performed PCA dimensionality reduction to find anchor points (Supplementary Figure S2D). Clustering analysis of 1,124 cells base on “FindNeighbors” and “FindClusters” function were classified into 7 clusters (Figure 1A). These 7 clusters were labeled as 5 cell types by marker genes (Figure 1B). An overview of 1,124 cells from the G1 phase, G2/M phase, and S phase of samples is shown in Figure 1C. Next, the “FindAllMarkers” function was conducted to screen marker genes of 5 cell types, the top 5 marker genes in 5 cell types are listed in Figure 1D. The G1 phase, G2/M phase, S phase, and cell numbers were characterized by 5 cell types (Figure 1E). KEGG enrichment analysis showed that maker genes of 5 cell types were enriched in different various pathways, indicting the five cell subtypes obtained based on scRNA-seq analysis showed significant heterogeneity (Figure 1F).

To further screen out the subgroups of differences in patients with AMI, we calculated the abundance of the 5 identified cell types in AMI (19 cases) and healthy tissues (19 cases) from the GSE109048 dataset by the CIBERSORT method (Figure 2A). The results showed that hematopoietic stem cell frequency was higher in AMI than in healthy tissues (Figure 2B).

Through overlap analysis of genes associated with 7 drugs activity, genes associated with MI in DisGeNET, and genes associated with MI in Comparative Toxicogenomics Database, we found a total of 54 common genes (Figure 3A and Supplementary Table S1). Moreover, 1 gene of 54 genes belongs to the marker gene of hematopoietic stem cell types (Figure 3B). The 1 gene, MGLL, single-cell expression level was shown in Figure 3C. PPI analysis based on the String database showed there were 52 common genes used for the next analysis (Figure 3D). According to the analysis of network topological properties, it can be observed that the degree distribution mainly presents the dark rate distribution, which is in line with the characteristics of a biological network (Supplementary Figure S3A). The closeness, betweenness, and eigenvector distributions show that only a few nodes have high closeness, betweenness, and eigenvector (Supplementary Figures S3B–D).

To explore the functional annotation of those genes, we performed KEGG and GO enrichment analyses. For the GO functional annotations of genes, 795 terms were enriched in BP (p < 0.05) and the top 10 annotation results are shown in Figure 4A. The 2 terms enriched in CC are shown in Figure 4B. In total, 114 terms were enriched in MF, and the top 10 terms enriched in MF are shown in Figure 4C. The results of the KEGG pathway enrichment analysis showed a total of 30 pathways were enriched, and the top 10 annotations are shown in Figure 4D.

Cytoscape was used to analyze the network based on modules, and the mature MCODE algorithm was used to find the closely connected protein groups in the target network. The results showed that 5 clusters were obtained by the MCODE algorithm (Figure 5A). KEGG enrichment analysis of genes in 5 clusters showed those genes involved include immune and inflammatory pathways (Figure 5B). In addition, the correlation analysis between 5 clusters and hematopoietic stem cell type using the “Hmisc package” indicated that three genes (ATM, CARM1, and CASP8) in 5 clusters were negatively correlated with hematopoietic stem cell type, while genes, CASP3 and PPARG, were positively correlated with hematopoietic stem cell type (Figure 6A). Next, correlation analysis of 5 hub genes (ATM, CARM1, CASP8, CASP3, and PPARG) and immune cell score calculated by the MCPcounter package showed that ATM, CARM1, CASP8 were positively, while CASP3 and PPARG were negatively related to T cells and B cells (Figure 6B). In addition, we used ssGSEA to analyze the enrichment scores of each patient in different KEGG pathways, and calculated the correlation between the expression of these five genes and the KEGG pathway. We observed ATM, CARM1, CASP8, CASP3, ACUTE_MYELOID_LEUKEMIA, and ADIPOCYTOKINE_SIGNALING_PATHWAY_. There was a significant positive correlation between pathway, PPARG and ACUTE_MYELOID_LEUKEMIA, and ADIPOCYTOKINE_SIGNALING_PATHWAY was significantly negatively correlated (Supplementary Figure S3E). The interaction relationship between these five genes was obtained by using a String database. It can be observed that most of them have a direct regulatory relationship (Supplementary Figure S3F). Co-expression analysis shows that ATM, CARM1, and CASP8 are highly positively correlated, and PPARG is highly negatively correlated with ATM, CARM1, and CASP8 (Supplementary Figure S3G).

Based on the above analysis, ATM, CARM1, CASP8, CASP3, and PPARG were possible key genes associated with AMI. Drugs targeting these genes may have a greater impact on the occurrence and progression of AMI. According to previously described methods, distance, the density distribution of drugs to related gene sets (Figure 7) identified 5 drugs that may be correlated with PPARG (Table 1).

To confirm their potential as molecular drugs for treating patients with AMI, molecular docking of these components to PPARG was achieved using Autodock Vina software. Of the 4 compounds, DB05490 had the lowest docking score (−9.9 kcal/mol, Table 2), which had a hydrogen-bonding interaction with SER289 (Figures 8A,B). The molecular dynamics simulation of 100 ns showed that the RMSD of compound DB05490 was relatively stable in general, except for a significant increase in the first 5 ns (Figure 8C), to some extent, reflecting that compound DB05490 could stably bind to the active site of PPARG, and then play an inhibitory role on PPARG.