Microarray datasets were obtained from the publicly available Gene Expression Omnibus (GEO) database.¹ Eleven datasets were included and analyzed in this study (22–28) (Table 1). The gene expression matrix derived from GSE59867 and GSE62646 was merged and normalized. The ComBat function of the sva package was applied to eliminate batch effects (29). Distribution and variation of the samples before and after batch normalization were evaluated by principal component analysis (PCA). The normalized gene expression matrix was used for subsequent analyses. Each of the remaining data set was processed independently.

DEGs were screened by the Limma package (30). DEGs analyses were performed with a p < 0.05 and |log₂FC| > 0.585 (|FC| > 1.5) signifying statistical significance. Volcano plots for visualization were drawn using the ggplot2 package.

Gene set enrichment analysis was employed to identify enriched gene sets and the leading-edge gene subsets which contributed most to the enriched pathways. GSEA was performed using the clusterProfiler R package (31). Human genes were annotated using org.Hs.eg.db package. Gene Ontology (GO) terms² and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways³ were used as the reference gene sets. Gene sets with an absolute normalized enrichment score (|NES|) > 1, p value <0.05 and adjusted p value (p.adj) < 0.25 were considered to be significantly enriched (32). The Benjamini–Hochberg procedure was used to control the false discovery rate. Visualization was performed using the ggplot2 package.

Gene set variation analysis was performed to quantify pathway activities of the samples using GSVA R package (33). The Hallmark gene sets from MSigDB (34) were selected as the reference gene sets. Limma package was used to define the differentially expressed gene sets. A p value <0.05 was set as the cut-off criterion.

We employed the xCell algorithms to estimate the relative abundance of leukocytes in PBMCs (35) using the online analytical platform.⁴ The Charoentong signatures were selected as the reference gene sets.

Time series analysis was performed using Mfuzz (2.50.0) package (36). Functional enrichment analyses of gene clusters were generated by the online analytical platform METASCAPE.⁵ A p value <0.01 was taken as the threshold for statistical significance.

The consensus clustering was executed using ConsensusClusterPlus package (37). To ensure the classification stability, the analyses was repeated 1,000 times with different initial conditions for modeling fitting to assure internal consistency (parameters: reps = 100, p Item = 0.8, p Feature = 1). Each sample of PBMCs from patients with stable coronary artery disease (SCAD) and MI patients was assigned to either one of two sub-clusters, based on the expression levels of genes belonging to the MSigDB (34) gene sets “GOBP_INFLAMMATORY_RESPONSE” (n = 773). Ward.D and Euclidean distances were used as the clustering algorithm and distance metric, respectively. PBMCs from HF patients were classified into three sub-clusters based on our custom gene sets containing 26 glycolytic genes and 47 oxidative phosphorylation (OXPHOS) genes (Supplementary Data Sheet 1). Ward.D and Spearman distances were used as the clustering algorithm and distance metric, respectively.

Correlation analysis was performed to evaluate correlations between the abundance of and the expression level of pro-inflammatory transcripts in leukocytes. The correlation coefficient was obtained by performing the Spearman’s correlation analysis. Visualization was performed using the ggplot2 package.

Statistical analyses and plotting were performed using R or GraphPad Prism 9 software. Given that the gene expression data do not fit a normal distribution by the Shapiro–Wilk test, two-group comparisons were determined using the Mann–Whitney non-parametric test. Multiple-group comparisons were determined using Kruskal–Wallis test with Dunn’s correction. Comparison of paired samples was performed by using Wilcoxon matched-pairs test. A spectrum of gene expressions was visualized by Box-Whisker plots or violin plots. The Box-Whisker plots denoted median + quartiles (box) and range (whiskers) and the violin plot showed the density distribution of all data points. For all tests applied, statistical differences were determined significant at p < 0.05.

To better understand the MI-associated inflammatory processes, we began with exploring alterations of gene expression profiles associated with the MI PBMCs. Microarray datasets of GSE59867 and GSE62646 were obtained from the publicly available GEO database. In the cohort of GSE59867, MI patients at admission (n = 111) who were indicated for percutaneous coronary interventions (PCI) were assigned as the study group. Age-, gender-, and BMI-matched SCAD patients (n = 46) without a history of MI were used as the reference group. The baseline demographic and clinical characteristics of the two groups were included in the original publication (17). In the cohort of GSE62646, MI patients at admission (n = 28) who were indicated for PCI were assigned as the study group. Gender- and BMI-matched SCAD patients (n = 14) without a history of MI were used as the reference group. It was noted that the patients in the SCAD group were older than the patients in the MI group in this cohort (67.8 years of age in the SCAD control group vs. 55.4 years of age in the MI group). Given that age is a non-modifiable risk factor for chronic inflammation associated with age-associated diseases (38), it is possible that the difference in age here might not cause overestimation of the potential association of inflammation with MI. The baseline demographic and clinical characteristics of the two groups were showed in the original publication (23). Given that both GSE59867 cohort and GSE62646 cohort were carried out with similar designs and the PBMCs samples from these two experiments were processed on the same platform (GPL6244), we merged and normalized the two data sets to expand the sample size and to increase the reliability of the analysis. The combined dataset contained 139 PBMCs samples from the MI patients at admission and 60 PBMCs samples from the SCAD patients. Sample distribution without (Figure 1A) and with (Figure 1B) batch-normalization adjustment was visualized using PCA plots. The normalized gene expression matrix was used for subsequent analyses. The alterations of gene expression profiles associated with the MI PBMCs were examined with PBMCs from SCAD patients serving as the control. GSEA revealed that inflammatory response (e.g., cytokine secretion, leukocyte chemotaxis and leukocyte cell–cell adhesion, etc.) and metabolic pathways (e.g., galactose metabolism and glycolysis/gluconeogenesis, etc.) were significantly enriched in the MI group (Figure 1C). Furthermore, the expression patterns of the core leading-edge genes contributed most to glycolysis were similar with that of acute inflammatory response (Figure 1D). Interestingly, several pathways involved in platelet activation, for instance, platelet degranulation, coagulation, and platelet alpha granule lumen, were significantly enriched (Figure 1C). Enrichment of platelet-associated pathways in the MI PBMCs may in part result from increased leukocyte-platelet aggregates and coagulation abnormalities (39, 40). Detailed GSEA enrichment results can be found in Supplementary Data Sheet 2. Moreover, the xCell analysis revealed that in the PBMCs derived from MI patients, the relative abundance of monocytes and neutrophils was significantly increased and the relative abundance of CD4⁺ lymphocytes and CD8⁺ lymphocytes was significantly decreased (Figure 1E), consistent with that reported in the literature (41). This increase in blood neutrophils and monocytes might result from exodus of neutrophils and monocytes from the hematopoietic stem and progenitor cells in the bone marrow (16). It was noted that increased abundance of neutrophils could in part be caused by the presence of low-density granulocytes (LDG) in the MI PBMCs (42). Taken together, the results here indicate that PBMCs from the MI patients at admission are molecularly marked by pro-inflammatory and glycolytic gene signatures.

Consensus clustering was further performed to clarify if augmented inflammation is universally or selectively associated with the MI patients. According to the expression of genes belonging to the MSigDB gene sets “Inflammatory response” (n = 773), each sample was assigned to either one of the two subtypes using consensus clustering (Figure 2A). Based on the expression of selected genes, k = 2 seemed to be an adequate selection with the clustering stability increasing from k = 2–7 (Figures 2B,C). Out of the 139 MI patients, 53 were classified into the subgroup 1 and 86 were classified into the subgroup 2. In contrast, all SCAD patients were classified into the subgroup 1, which was characterized by a hypo-inflammatory phenotype (Figures 2D–G). To further compare the phenotypic differences of PBMCs between the two subgroups, DEGs analyses were performed with a |log₂FC| > 0.585 and a p value of <0.05 signifying statistical significance. GSVA was employed to quantify the hallmark pathway activities of the samples. Detailed GSVA results can be found in Supplementary Data Sheet 3. The results revealed that the expression profiles of MI PBMCs in the subgroup 1 shared more similarity to the profiles in SCAD PBMCs instead of resembling the profiles of the MI PBMCs in the subgroup 2 (Figures 2D–F). Multiple hallmark pro-inflammatory pathway activities were elevated in the MI group, particularly in the subgroup 2, including inflammatory response, IL6/JAK/STAT3 signaling, TNFA signaling via NFκB, Wnt-beta catenin signaling, hypoxia and glycolysis (Figure 2G). These results indicate a portion of the MI patients are classified into the hyper-inflammatory group, while the remaining patients display less-severe inflammation. Furthermore, our analysis here suggests that consensus clustering could be a feasible approach to classify MI patients into the hyper-inflammatory subtype and the hypo-inflammatory subtype.

Next, we sought to validate the foregoing observations and evaluate changes in inflammatory markers by analyzing GSE123342, GSE48060, GSE141512, GSE19339, GSE66360, GSE109048, and GSE159657 data sets from additional MI cohorts. For each data set, GSEA was performed by ranking the gene sets after differential expression analysis. Detailed GSEA enrichment results can be found in Supplementary Data Sheet 4. Hallmark pathways including acute inflammatory response, cytokine secretion, chemokine production and canonical glycolysis were selected to assess pro-inflammatory activation of the samples. GSEA revealed that gene sets associated with acute inflammatory response, cytokine secretion, chemokine production and canonical glycolysis were enriched in the whole blood (Figures 3A,B), PBMCs (Figure 3C), and thrombus-derived leukocytes (Figure 3D) from the MI patients. Based on the results from the xCell analysis, the increase in the relative abundance of monocytes in the post-MI blood was validated in GSE48060 and GSE141512 data sets and the increase in neutrophils was verified in GSE123342 and GSE48060 data sets (Figure 3E). Genes of interest included those that showed a significant upregulation in over two third of the indicated datasets. The representative leading-edge genes contributed most to the enriched gene sets involved in inflammatory response and glycolysis were listed in Figure 3F. The post-MI peripheral blood was generally marked by an upregulation of pro-inflammatory genes (e.g., CD14, CLEC5A, HP, ICAM1, IL10, IL1B, NLRP3, S100A8, S100A9, TREM1, etc.) and glycolytic genes (e.g., ALDOA, ENO1, GAPDH, GAPDHS, HK1, HK3, PGAM1, PGK1, PKM, TPI1, etc.). In addition, the circulating endothelial cells (CECs; Figure 4A) and platelets (Figure 4B) from the MI patients also manifested a significant enrichment of pro-inflammatory gene sets. Although these specific gene sets did not appear to be enriched in circulating exosomes (Figure 4C), the exosomal level of the classical inflammation markers such as CD14, CLEC5A, IL1B, NLRP3, S100A8, S100A9, and TREM1 was significantly upregulated (Figure 4D). Detailed GSEA enrichment results of these datasets can be found in Supplementary Data Sheet 4. The strength of the relationship between hallmark inflammatory markers and relative abundance of leukocytes was assessed using Spearman’s correlation coefficients. The results revealed that almost all of these pro-inflammatory transcripts exhibited positive correlations with the abundance of monocytes and neutrophils (Figure 4E), indicating monocytes and neutrophils could be key players in the MI-associated acute inflammatory responses. Our results from multiple validation datasets further confirm that MI patients are generally characterized by a phenotype of pro-inflammatory activation in the blood. Furthermore, the results here reveal that the robust post-MI upregulation of the hallmark pro-inflammatory transcripts is a common phenomenon in the peripheral whole blood, PBMCs, thrombus-derived leukocytes, CECs, platelets and exosomes.

The temporal evolution of gene expression profiles of PBMCs from the MI patients

We further integrated the data sets of GSE59867 (n = 436) and GSE62646 (n = 98) and incorporated all samples from four time points to investigate the temporal evolution of gene expression profiles of PBMCs from the MI patients. Dataset of GSE59867 contained a total of 436 PBMCs samples, including 46 samples from the SCAD patients, 111 samples from the MI patients at admission (on the 1st day of MI), 101 samples from the MI patients at discharge (after 4–6 days of MI), 95 samples from the patients 1 month after MI, and 83 samples from the patients 6 months after MI (22). The number of patients included decreased from one time point to the next because some patients failed to follow-up. Dataset of GSE62646 contained 98 PBMCs samples, including 14 samples from the SCAD patients, 28 samples from the MI patients at admission (on the 1st day of MI), 28 samples from the MI patients at discharge (after 4–6 days of MI) and 28 samples from the patients 6 months after MI (23). No patients from the original cohort of GSE62646 had been lost to follow-up. We merged the two datasets and normalized the matrix using sva package to control for variability between batches. The distribution of samples without (Figure 5A) and with (Figure 5B) batch-normalization adjustment was visualized by PCA plots. The normalized gene expression matrix was used for subsequent analyses. The combined dataset contained 534 PBMCs samples, including 60 samples from the SCAD patients, 139 samples from the MI patients at admission, 129 samples from the MI patients at discharge, 95 samples from the patients 1 month after MI and 111 samples from the patients 6 months after MI. GSVA score was computed per sample to quantify the pathway activities of each patient. Detailed GSVA results can be found in Supplementary Data Sheet 5. Next, we selected several hallmark pathways involved in inflammatory response and compared their GSVA score among different groups. As shown in Figure 5C, the GSVA score of inflammatory response and TNFA signaling via NFκB were increased at admission and returned to the baseline at discharge, while those of IL6/JAK/STAT3 signaling and glycolysis returned to the baseline 1 month after MI. The GSVA score of coagulation remained elevated until 1 month after MI, while GSVA score of OXPHOS was below the baseline level 6 months after MI. Next, the relative abundance of immune cells generated by xCell algorithm was compared among groups. As shown in Figure 5D, increased abundance of monocytes and neutrophils and decreased abundance of CD4⁺ and CD8⁺ T lymphocytes were observed in PBMCs from the MI patients at admission. The relative abundance of neutrophils returned to the baseline at discharge, while the level of monocytes seemed to be above the baseline even 6 months after MI. The relative abundance of B lymphocytes did not differ among groups. In addition, we visualized the dynamic expression patterns of the post-MI PBMCs using GSEA method. The MI groups from different time points were compared separately to the SCAD group. Detailed GSEA results can be found in Supplementary Data Sheet 6. As shown in Figure 6, the signals of most of the hallmark pathways related to inflammation and metabolic alterations peaked at admission and declined thereafter. In contrast, the signals of natural killer (NK) cell-mediated immunity were lower at admission and returned to the baseline over time. Interestingly, we noted a significant downregulation of OXPHOS in PBMCs 6 months after MI, which was consistent with the GSVA results shown in Figure 5C. Thus, both GSEA and GSVA algorithms confirmed that the signals of OXPHOS in PBMCs 6 months after MI was downregulated. To further characterize the temporal dynamics of the post-MI PBMCs, we used Mfuzz method to cluster the DEGs according to their expression patterns throughout the time course. A total of 94 DEGs with 1.5-fold enrichment (p value<0.05 and |log₂FC| > 0.585) in at least one time point were taken into Mfuzz analysis. Two clusters of DEGs were identified by Mfuzz (Figures 7A,B). Cluster 1 contained 45 genes whereas cluster 2 contained 49 genes. Detailed Mfuzz results can be found in Supplementary Data Sheet 7. The expression level of genes in cluster 1 (e.g., SOCS3, HP, ECRP, APQ9, and FAM20) was increased on the first day of MI and decreased thereafter, while the expression in cluster 2 (e.g., SNORA16A, SNORD45B, CLC, KLRC2, and SONRA20) were decreased on the first day of MI and increased thereafter. The dynamic expression patterns of the representative leading-edge genes contributed most to the enriched GSEA terms were depicted using boxplots (Figure 7C). These results suggest a possibly dynamic expression pattern of genes involved in inflammatory response, metabolic reprogramming and platelet activation in PBMCs during a 6-month observation period.

In addition, as shown in Figure 6, the signals of platelet activation peaked at discharge rather than at admission. To our best knowledge, this observation has not been previously reported. To avoid error possibly caused by consolidation of datasets, GSEA analysis for the dataset GSE59867 and GSE62646 was performed separately. As shown in Supplementary Figures 1A,B, the signals of hallmark platelet activation pathways were higher at discharge than admission. It has been demonstrated that the number of the circulating reticulated platelets (RPs) is elevated in MI patients who have received PCI and anti-platelet therapies and remained relatively stable during the first few months after the onset (43). Given that RPs are newly-formed immature and hyperactive platelets containing significantly more RNA compared to mature platelets (MPs) (44), we further tested the possibility that RPs might contribute to higher signals of platelet activation. First, we compared the transcriptomic profiles of RPs and MPs (45). The results showed that multiple pathways were enriched in RPs, including coagulation, platelet activation, platelet aggregation and platelet degranulation (Supplementary Figure 1C). In particular, transcripts involved in platelet reactivity, clotting and cell–cell interactions (e.g., DMTN, SELP, ITGA2B, ITGB3, GP6, TREM1, TUBB1, and C6orf25) were upregulated in the RNA-rich RPs (Supplementary Figure 1D), indicating that RPs have a prothrombotic transcriptomic profile Moreover, the expression levels of these hallmark RP genes in MI PBMCs were higher at discharge than admission (Supplementary Figure 1E). Thus, these results suggest that higher platelet activation could be in part explained by possibly enhanced platelet turnover and increased RPs in MI patients at discharge.

To gain additional insights into the molecular signatures of PBMCs associated with the development of post-MI HF, we analyzed GSE59867 data set and extracted samples from the patients with or without clinical HF 6 months after MI. Samples from the SCAD patients (n = 46) in GSE59867 data set served as controls. According to the original study, the post-MI patients who volunteered for a visit 6 months after the onset were divided into the HF group (n = 9) and the non-HF group (n = 8) (22). Patients in the HF group were characterized by decreased LVEF% (39.3 ± 8.4%) and elevated level of NT-pro-BNP (918.3 ± 848.5 pg./ml) 6 months after the onset of MI. The sample size of the HF group at admission, at discharge, 1 month and 6 months after MI was 9, 9, 8, 8, respectively. The number of patients included may decreased from one timepoint to the next because of the missing data in the original study. Patients in the non-HF group did not develop HF 6 months after MI. The NT-pro-BNP levels and LVEF values in the non-HF group 6 months after the onset of MI were 62 ± 14.1 pg/ml and 66.8 ± 1.9%, respectively. The sample size of the HF group at admission, at discharge, 1 month and 6 months after MI was 8, 6, 8, 8, respectively. The number of patients at discharge decreased because of the missing data in the original study. The baseline demographic and clinical characteristics of the HF group and the non-HF group were showed in the original publication (22). There was no significant difference between the two groups in the following parameters: gender, age, BMI, hypertension, diabetes, hypercholesterolemia, previous MI, and smoking (22). All patients in the HF group and the non-HF group received standard medications, including aspirin, beta blockers, ACE inhibitors, statins, etc. Firstly, the HF group and non-HF group were analyzed against the SCAD group using GSEA method. Detailed GSEA results at various time points can be found in Supplementary Data Sheet 8. Based on the aforementioned results, acute inflammatory response, coagulation, canonical glycolysis and OXPHOS were chosen as representative gene sets. As shown in Figure 8, MI patients at admission were marked by a significant enrichment of genes in inflammatory response, coagulation, canonical glycolysis and OXPHOS in their PBMCs, irrespective of the development of HF. Gene sets of acute inflammatory response were consistently enriched until 1 month after MI in the HF group, whereas this enrichment was not detected at discharge in the non-HF group. Meanwhile, coagulation gene set was positively enriched in the HF group and negatively enriched in the non-HF group 1 month after MI. Next, we visually compared the difference in PBMCs gene expression profiles between the HF group and the non-HF group. Detailed GSEA results at various time points can be found in Supplementary Data Sheet 9. As shown in Figure 9A, compared to the non-HF group, the HF group had 142 DEGs at admission (61 upregulated and 81 downregulated), 42 DEGs at discharge (23 upregulated and 19 downregulated), 45 DEGs 1 month after MI (25 upregulated and 20 downregulated) and 14 DEGs 6 months after MI (1 upregulated and 13 downregulated). GSEA revealed that acute inflammatory response and coagulation gene sets were significantly enriched in the HF group at admission, at discharge and 1 month after MI (Table 2), indicating that post-MI development of HF is associated with elevated and prolonged inflammation and coagulation abnormalities. Furthermore, OXPHOS gene set was negatively enriched in the HF group 6 months after MI, suggesting potential mitochondrial dysfunction of PBMCs in the post-MI HF patients. The core genes contributed most to the enriched gene sets were further assessed by GSEA leading edge analysis. The results were visualized using a heatmap (Figure 9B). The temporal expression patterns of representative pro-inflammatory transcripts (e.g., CD163, CXCR2, ICAM1, IL1B, NLRP3, PTGS2, S100A9, and TNF) in the PBMCs from the post-MI HF and non-HF groups were depicted using line plots (Figure 9C). Compared to the non-HF group, all of these genes exhibited a trend of upregulation at the early time points in the HF group, although not all differences were statistically significant. No differences were observed between the two groups 6 months after MI. Taken together, these results suggest that PBMCs from the post-MI HF patients are characterized by excessive activation of inflammation and coagulation pathways, indicating that prolonged elevation in inflammatory response and coagulation in PBMCs may predict the development of post-MI HF.

Metabolic reprogramming is mechanistically implicated in the inflammatory activation of immune cells (46). We thus further zoomed in on glycolysis and OXPHOS signals in PBMCs. GSEA revealed a significant upregulation of glycolysis and OXPHOS signals in PBMCs from the MI patients at admission (Figure 10A), consistent with the results shown in Figure 6. However, a significant downregulation of OXPHOS signals was noted in the patients 6 months after MI compared to the SCAD patients (Figure 10B), which was in agreement with the results shown in Figures 5B, 6. It is worth noting that similar observations were made in the HF patients 6 months after MI compared to those without HF (Figure 9C). The top five leading-edge genes contributed most to the canonical glycolysis and OXPHOS gene sets were marked on the GSEA plots (Figures 10A–C). The expression level of glycolytic genes such as HK3, PGAM1, PKM, ALDOA, and HK2 was increased at admission and decreased thereafter in both HF and non-HF groups (Figure 10D). Compared to the non-HF group, these glycolytic genes exhibited a trend of upregulation at the early time points in the HF group, although not all differences reached statistical significance. The expression level of OXPHOS genes such as COX6C, NDUFV2, NDUFA1, UQCRB, and COX7B was below the baseline after MI (Figure 10E). No significant differences in the expression level of these genes were noted between the HF group and non-HF group. However, although statistically insignificant, lower expression of OXPHOS genes COX7A2, NDUFB5, NDUFC2, NDUFB6, and UQCR10 was noted in the HF group 6 months after MI (Figure 10E). The statistical insignificance was possibly due to small sample sizes and the heterogeneity of the HF patients. Nevertheless, these results suggest a possibility that PBMCs from the post-MI HF patients are characterized by prolonged metabolic abnormalities.

To further clarify if the downregulation of OXPHOS genes in PBMCs is universally or selectively associated with the HF patients, and to determine whether metabolic reprogramming is correlated with inflammatory response in PBMCs from the HF patients, GSE77343 data set containing 197 PBMCs samples from chronic HF patients was analyzed. Consensus clustering was performed to stratify the HF patients based on their metabolic characteristics. After developing a custom gene set containing 26 glycolytic genes and 47 OXPHOS genes (Supplementary Data Sheet 1), a group of robustly co-expressed glycolytic (n = 21) and OXPHOS (n = 29) genes identified by consensus clustering was subjected to molecular subtyping (Figure 11A). Based on the median relative expression levels of the co-expressed glycolytic and OXPHOS genes, PBMCs from the HF patients were divided into four subgroups based on their metabolic profiles, namely, subgroup 1 (S1: low OXPHOS and low glycolysis), subgroup 2 (S2: low OXPHOS and high glycolysis), subgroup 3 (S3: high OXPHOS and low glycolysis), and subgroup 4 (S4: high OXPHOS and high glycolysis; Figure 11B). Heatmap was generated to depict the expression patterns of the glycolytic and OXPHOS genes (Figure 11C). Subgroup 2 defined the largest group of cases (75/197; 38.1%), followed by subgroup 3 (73/197; 37.1%), subgroup 4 (30/197; 15.2%), and subgroup 1 (19/197; 9.6%). Since it has been reported that under aerobic conditions, the energy metabolism of PBMCs primarily relies on mitochondrial OXPHOX rather than cytoplasmic glycolysis (47), the S3 subgroup most closely resembled the metabolic profiles of the normal PBMCs. Thus, further comparisons were made against the S3 subgroup. To answer whether the metabolic profiles of PBMCs was associated with inflammation, blood routine analysis and GSEA were carried out for each group. According to the results detailed in the original study, the S1 subgroup had a lower proportion of monocytes and higher neutrophil-to-monocyte ratio (NMR) than the S3 subgroup; the S2 subgroup had a higher proportion of neutrophils, lower proportion of lymphocytes, higher neutrophil-to-lymphocyte ratio (NLR) and higher monocyte-to-lymphocyte ratio (MLR) than the S3 subgroup; the S4 subgroup had a higher proportion of monocytes and higher MLR than the S3 subgroup (Figure 11D). Neither age, nor the proportions of basophils and eosinophils differed among the four subgroups based on their PBMCs metabolic profiles. GSEA was performed to further validate the heterogeneity of PBMCs in different subgroups. Detailed GSEA results can be found in Supplementary Data Sheet 10. As shown in Figure 11E, canonical glycolysis gene set was positively enriched in the S2 and S4 subgroups, whereas OXPHOS gene set was negatively enriched in the S1 and S2 subgroups. The S2 and S4 subgroups, which were characterized by a highly glycolytic phenotype, exhibited a positive enrichment of gene sets involved in inflammatory response and phagocytosis. However, the S1 and S2 subgroups, which displayed a low OXPHOS phenotype, displayed a negative enrichment of gene sets implicated in mitochondrial function. Thus, the results here demonstrate the diversity and heterogeneity of PBMCs from HF patients.