A flowchart of the current study is illustrated in Figure 1 . Dataset GSE66360 based on the Affymetrix platform ( Table 1 ) was downloaded from the Gene Expression Omnibus database (14), consisting of 99 samples, including 49 MI samples and 50 normal samples. The test dataset GSE62646 contains 84 MI samples and 14 patients with coronary artery disease (CAD). The single-cell dataset GSE180678 was obtained from the GEO database, which contained one patient with ischemic cardiomyopathy ( Table 1 ).

Data were preprocessed before performing the DEG analysis. The Affymetrix in R was utilized to improve the database via background calibration, normalization, and log2 transformation (15). We averaged the probe site values if a gene had more than one. According to the annotation files from the platform, the probe IDs were converted to gene symbols, and probes that did not correspond to the symbol of the gene were removed. Finally, the eligible DEGs were screening and electing by the linear models for microarray (LIMMA) data package; P-value < 0.05, and |log2 FC| > 1 (16).

First, we used Seurat (17) for data processing to screen reliable cell subpopulations from single-cell data. The expression of each gene must occur in at least three cells, with at least 200 genes expressed per cell. The percentages of mitochondria and rRNA were calculated by the PercentageFeatureSet function, ensuring 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. The data were then normalized by log normalization and used to identify highly variable genes using the FindVariableFeature function. Principal Component Analysis (PCA) was used to scale down all genes using the ScaleData function. Next, subgroups of cells were determined using the FindNeighbors and FindCluster functions (set resolution = 0.1), and annotations were applied. The R package FindAllMarkers was used to identify markers for each cell subpopulation and the gene expression profile dataset was re-evaluated using the CiberSort (18) algorithm in the context of the expression of marker genes for each cell subpopulation to determine the composition of individual cell subpopulations in different samples in the expression profile. Finally, DEGs were combined with cell subpopulations to analyze the distribution and abundance of DEGs in cell subpopulations.

Differential genes are the basis of our subsequent analysis study, and whether DEGs are associated with MI determines the accuracy of our following study. Therefore, we performed DO enrichment analysis on differential genes to observe the degree of enrichment of DEGs in MI. We used the “cluster Profiler” package (19) and “DOSE” package (20) for DO enrichment analysis to analyze the enrichment degree of differential genes in MI and visualize the results on a bar chart (barplot) and a bubble plot (dotplot).

It has been established that ssGSEA can be used to quantify the correlative scale of immunological cells and functions from DEGs. Thus, ssGSEA was adopted to evaluate the correlation between a single sample gene of immune cells and functions based on the converged biological data (21). Then, the results of immune cells and immune function were visualized in a heat map using the “pheatmap” package (22). Finally, the R package “GGally” was utilized to visualize the correlation between 16 immune cells and 12 immune functions. A P-value < 0.05 was statistically significant.

The divergence in immune scores between normal and MI groups and the difference in immune cells and functions between the cohorts were assessed by data conversion and visualization using ggplot2 and reshape2 statistical packages by R. We further conducted a correlation analysis of significantly expressed genes between immune cells and function

Next, we selected the top five upregulated and downregulated genes among the most significant DEGs to clarify the correlation between MI-significant differential genes and immune cells and function. The top six genes mostly associated with immunity were screened as hub genes. Subsequently, the data were visualized using ggplot2 packages.

The identified hub genes were used to establish a gene-disease diagnostic model to determine the accuracy of the genes mostly associated with immune cells and function for MI diagnosis. The goodness-of-fit between the actual and predictive values and the degree of the fitting was evaluated using a calibration plot and Spiegelhalter’s Z-test. In addition, the receiver operating characteristic curve (ROC), the area under the ROC curve (AUC), and Harrell’s concordance index (C-index) were utilized to appraise the nomogram’s discrimination ability. Finally, we compared the levels of hub gene expression between control participants and MI cases in GSE62646 to validate the accuracy of the model by the same analysis methods as above.

To expand the clinical value of the immune-related gene model, the Enrichr platform (https://maayanlab.cloud/Enrichr/) and the coremine platform (https://www.coremine.com) were used to predict the drug delivery of NR4A2 in the model. Drug predictions are identified through the Drug Signature Database (DSigDB) in Enrichr, a large portal with a large number of different genomic libraries that identify targeted drugs associated with NR4A2. A P-value < 0.05 was used to select drugs with potential clinical application.

Twenty SD rats were divided into an experimental group (n=17) and a control group (n=3). After 7 days of adaptation, the rats underwent ligation of the anterior descending branch of the coronary artery after anesthesia on day 8. An electrocardiogram was performed after ligation to observe whether the ST segment of the rats showed evident elevation, and the ST-segment elevation indicated successful modeling. After the operation, the rats were fed for seven days. The rat’s abdominal aorta was severed on postoperative day 8 to induce massive hemorrhage and death, and the heart tissue was harvested after death. We divided 6 rats into two groups for quantitative real-time PCR (qPCR) and WB (MI=3, control=3). This study was approved by the Ethics Committee of the Affiliated Hospital of Liaoning University of Traditional Chinese Medicine (2022CS(DW)-006-01).

50 mg of tissue from the ischemic portion of the left ventricle of the experimental group and the control group were harvested for qPCR and WB experiments, cut in pre-chilled PBS, and centrifuged for 5 min. After the supernatant was removed, the tissue was homogenized by adding the corresponding lysate volume. Then the tissues were fully lysed at 4°C for 20 min and then centrifuged at 12 000 r/min for 20 min. BCA detected the protein concentration in the supernatant, and the remaining supernatant was mixed with protein loading buffer (5×Loading Buffer) at a ratio of 1:4 and boiled to denature the protein. After separation by protein gel electrophoresis, the proteins were transferred to PVDF membranes and blocked with 5% skim milk. The membranes were incubated with the primary antibody at 4°C at a dilution ratio of 1:5000 for NR4A2 and 1:5000 for GAPDH, incubated overnight, eluted three times with TBST for 10 min each, and 1:5000 HRP-labeled secondary antibodies (goat anti-mouse and goat anti-rabbit) were added, respectively The proteins were incubated for 2 h at room temperature, eluted three times with TBST for 10 min each time, developed by exposure with ECL chemiluminescent solution, photographed by Bio-RAD camera system and analyzed by Image Lab software for relative expression of proteins.

50mg of experimental and control tissues were taken, added to 1mL Trizol, and homogenized on ice. Total RNA was extracted and cDNA synthesized. All primers were designed with Primer5.0 by referring to the sequences provided by Genbank (RNA fragment sequences are shown in Table 2 ). The procedure was conducted as follows: 1st cycle 95°C for 2 min, followed by denaturation at 95°C for 5 s, annealing at 60°C for 10 s, extension at 72°C for 30 s, 40 cycles, and finally terminated at 72°C for 5 min. The Ct value was the number of cycles the fluorescence signal in each reaction tube undergoes when it reaches a set threshold value. The specificity of the reaction products was determined based on the solubility curve and the results of product electrophoresis. The relative ratio of the target gene to GAPDH was used as its expression, and the relative ratio was calculated using the 2-△△Ct method.

R (version 4.1.3) was utilized for statistical analysis. A p-value < 0.05 was set as a significance threshold for screening DEGs, immune infiltration analysis, immune-related gene model, QPCR, WB experimental results data analysis, and drug prediction.

Dataset GSE66360 was selected and underwent differential expression analysis using the “Limma” package in R 4.1.3 software. Three hundred thirty-five DEGs were identified, including 280 up- and 55 downregulated genes (log2FC > 1 and P-value <0.05). These DEGs were visualized in a volcano plot ( Figure 2A ) with downregulated genes in green, upregulated genes in red, and the rest in black. The top 50 upregulated and downregulated genes were selected to create a heatmap ( Figure 2B ).

Cluster analysis categorized scRNA-seq cells into 12 clusters using the “FindNeighbors” and “FindCluster” functions ( Figure 3A ). Marker genes were used to categorize the 12 clusters into 5 types of cells ( Figure 3B ). Next, the “FindAllMarkers” function was used to screen the DEGs and 144 genes intersected with the sequencing results. We found that the DEGs were predominantly distributed in endothelial cells. This implies that endothelial cells may have an important role in MI ( Figure 3C ).

We conducted an enrichment analysis of significant genes and diseases to observe the enrichment degree of DEGs and MI. As shown in Figure 4 , DEGs were significantly enriched in diseases such as arteriosclerosis, bacterial infectious disease, arteriosclerotic cardiovascular disease, atherosclerosis, primary bacterial infectious disease, periodontal disease, lung disease, periodontitis, and myocardial infarction. 29 genes were significantly enriched in MI, suggesting a correlation between these DEGs and MI, highlighting the robust correlation between DEGs and MI. (For more information, please see Supplementary Chart 1 ).

Subsequently, we conducted ssGSEA to appraise immune cell content and function in DEGs. The results showed that the top three immune cells with the highest content of DEGs were CD8+T cells, B cells, and Neutrophils, while the highest immune score was found in cytolytic activity ( Figure 5A ).Then, we analyzed the immune infiltration of immunological cells and functions. The corrplot of immune infiltration divergence illustrated that the levels of 16 immune cells and 12 types of immune functions were high in MI patients.

Analysis of 16 immune cell types revealed that B cells were positively linked to TILs (r=0.738), regulatory T cells (Tregs) (r=0.631), and Natural killer (NK) cells (r=0.557). TILs were positively correlated to B cells, CD8+T cells (r=0.696), and Tregs (r=0.784), while Tregs were positively correlated with TILs, B cells, and Th2 cells (r=0.582). On the other hand, mast cells were negatively correlated with Tregs (r =-0.599), B cells (r=-0.441), TIL (r =-0.624), and CD8+T cells (r=-0.518). Moreover, Dendritic cells (DCs) were negatively correlated to B cells (r =-0.516), Tregs (r=-0.678), and TILs (r =-0.565) ( Figure 5B ). The corrplot of the immune function divergence indicated that parainflammation was positively linked to C-C chemokine receptor (CCR) (r=0.687), inflammation-promotion (r=0.57), and type I interferon (IFN) response (r=0.635), while T cell co-inhibition was positively associated with APC co-inhibition (r=0.507), checkpoints (r=0.68), inflammation-promotion(r=0.747) and T cell co-stimulation (r=0.547), type I IFN response (r=0.581), and type II IFN response (r=0.594) ( Figure 5C ).

After integrating the data of the MI and normal groups from the two datasets, we analyzed the difference in immune cells and immune function. We found that the Cell types that show the most difference in immune cell infiltration between MI and control are macrophages and neutrophils ( Figure 6A ). We also found a significant difference among the seven immune functions, with the most significant difference observed in CCR, Cytolytic activity, Parainflammation, followed by APC co-inhibition, and T cell co-stimulation. This phenomenon suggests that these immune functions may be involved at critical times in the pathophysiology of MI ( Figure 6B ).

The top 5 upregulated and downregulated genes from DEGs were used to analyze the association between immune cells and functions ( Figure 7 ). We found the highest positive correlation between S100A12 vs. CCR, followed by S100A12 vs. Parainflammation and S100A12 vs. Neutrophils. Conversely, the highest negative correlation was found between CCR2 vs. Tfh, followed by GIMAP4 vs. Tfh. The highest correlation between immune cells and function was observed in 6 genes: GIMAP7, GIMAP4, CCR2, NR4A2, CSTA, and S100A12.

We used hub genes to build gene-disease diagnostic models and generated a line chart, calibration map, and AUC curve. The nomogram was utilized to visualize the model, including the 6 hub genes. The fraction of each gene was consistent with the proportion in the nomogram. As shown in the nomogram ( Figure 8A ), GIMAP7, GIMAP4, CCR2, NR4A2, CSTA, and S100A12 were predicting factors of MI. High levels of NR4A2 and S100A12 were positively associated with the diagnosis of MI. Hence, the accuracy of MI diagnosis was predicted based on the total score of the three characteristics. The model’s performance was appraised by C-index, AUC, and calibration plots. The ROC curve of the model yielded an AUC value of 0.978, indicating that the model had good discrimination ability ( Figure 8B ). The calibration plot also indicated the better correction of the model ( Figure 8C ). To validate the accuracy of the prediction model, we performed a validation using a test dataset. The accuracy was evaluated by AUC plots. Interestingly, the AUC value ( Figure 8D ) and calibration plot ( Figure 8E ) was set at 0.9, indicating that the accuracy of the HUB gene is more reliable for the diagnosis of MI. Finally, we further explored the correlation between NR4A2 expression and immune cells ( Figure 8F ). As displayed in Figure 8F , NR4A2 was positively correlated with Tfh, pDCs, Parainflammation, Neutrophils, Macrophages, iDCs, and CCR. and negatively correlated with Tumor Infiltrating Lymphocytes, Th2 cells, T cell co−stimulation, Cytolytic activity, and B cells.

We identified a strong correlation between 180 drugs and MI. Finally, the top five drugs Aldosterone, liothyronine, cholecalciferol, etilefrine, and liothyronine, were selected (adjusted P-value < 0.05) ( Figure 9A ), (For more information, please see Supplementary Chart 2 ). Meanwhile, we found that 11 herbal medicines were associated with NR4A2 ( Figure 9B ); the top 6 herbal medicines were danggui, baizhi, zicao, mudanpi, heshouwu, and youhua.

The transcriptional changes of NR4A2 in the heart tissues of MI rats and controls were detected by qRT-PCR and WB. The results showed that the expression of NR4A2 was increased in both MI groups compared with the control group ( Figures 10A, B ), which was consistent with the in silico analysis results and indicated the potential diagnostic value of NR4A2 for MI.