Single cell RNA-seq identifies a FOS/JUN-related monocyte signature associated with clinical response of heart failure patients with mesenchymal stem cell therapy

Sample collection and RNA-sequencing

This study recruited patients at Shanghai East Hospital (Shanghai, China). The research was conducted in accordance with the ethical principles of the Helsinki Declaration [17] and obtained ethical approval from the Ethics Committee of Shanghai East Hospital. Prior to collection of specimens and clinical information, informed consent was obtained from all registered patients. All patients approved the release of personally identifiable information, received verbal and written project information, and signed a written consent form.

Human umbilical cord mesenchymal stem cells (hUC-MSCs) were administered via intravenous injection at a dose of 1×10⁶ cells/kg for heart failure patients [10]. Clinical assessments and evaluations were conducted on patients one week before cell implantation and at 1 month and 6 months after transplantation. These assessments included routine physical examinations, vital sign monitoring, adverse event evaluation, laboratory tests, and functional measurements (Supplementary Tables 1, 2). The above clinical monitoring showed improvement in patients after stem cell therapy. To further investigate the reasons behind the differences in treatment efficacy among patients, single-cell RNA-sequencing was performed on peripheral blood mononuclear cells (PBMCs) collected from four patients at both pre-treatment and one-month post-treatment, resulting in a total of 8 PBMC samples.

Public data collection

The GSE59867 used in present work, including PBMC samples from heart failure patients [18], were collected from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE59867). PBMC samples from patients (n=111) with ST-segment elevation myocardial infarction (STEMI) and the control group comprised patients (n=46) with a stable coronary artery disease (CAD) and without a history of myocardial infarction.

scRNA sequencing for PBMCs from heart failure patients with MSC therapy

The scRNA-seq library was generated using the Chromium Single Cell assay (10× Genomics). The library was sequenced using the NovaSeq 6000 Illumina platform and the raw data were aligned to the human genome (hg38). Quality control and normalization of the scRNA-seq data were performed using the Seurat R package (v4.1.1). The gene expression matrix was normalized to the total UMI counts per cell derived from the filtered cells. We used the FindVariableFeatures function to identify the top 3000 highly variable genes. The FindAllMarkers function in Seurat was executed with default parameters to identify genes exhibiting cluster-specific expression. The cell types of PBMC were annotated based on the WNN algorithm [19]. To determine the composition of samples based on cell types, the cell count for each cell type was calculated from each sample. The aforementioned identical functions were utilized to obtain subclusters of CD14⁺ monocytes for clustering and grouping.

Metacell analysis of scRNA-seq data

Single-cell transcriptomic data are characterized by high background noise, low gene detection rate, and sparse expression matrix. Currently, it is known that functionally similar genes tend to cluster together, thus algorithmically obtaining metacells often better reflect the true expression profile of cells [20, 21]. Therefore, metacell analysis was conducted on single-cell RNA sequencing data of PBMCs before and after MSC therapy for heart failure, thus providing dependable cellular data for subsequent machine learning modeling.

Construction of a machine learning model for therapy-response CD14⁺ monocytes

Employing machine learning and leveraging single-cell RNA sequencing to identify crucial regulatory genes may potentially assist in screening for subgroups more suitable for stem cell therapy. The Least Absolute Shrinkage and Selection Operator (LASSO) analysis was utilized to obtain the core genes, and seven machine learning algorithms were employed, resulting in appropriate algorithms and feature genes to better characterize the state of therapy-response CD14⁺ monocytes. The machine learning algorithms were conducted using the mlr3verse package (version 0.2.7) (https://CRAN.R-project.org/package=mlr3verse).

Gene set enrichment analysis (GSEA) and function enrichment analyses

Gene set enrichment analysis (GSEA) can score sorted gene lists (usually based on fold-change) and calculate if a certain gene set is significantly enriched by permutation test [22]. We used GSEA to investigate the potential role of response-related monocytes in the process of MSCs treating heart failure. In parallel, the Gene Ontology and KEGG enrichment analyses were also performed using clusterProfiler package [23].

Trajectory analysis for CD14⁺ monocytes

The Monocle 2 algorithm can elucidate the changing trends of gene expression as cellular states transition, even revealing hidden patterns of variation [24]. We employed pseudo-time analysis to investigate the process of patient response to MSC treatment before and after, and depicted the expression profiles of 7 key genes that change over time.

Single-cell metabolism analysis

The scMetabolism package was used to calculate the single-cell metabolic activity of cell subpopulations. In particular, the VISION algorithm was adopted and the Kyoto Encyclopedia of Genes and Genomes (KEGG) metabolic gene set was employed for analysis [25]. The DotPlot.metabolism function is used for visualization.

SCENIC analysis for CD14⁺ monocytes

Single-cell regulatory network inference and clustering (SCENIC) is a GRN (gene regulatory network) algorithm specifically developed for single-cell data [26]. Its innovation stems from the incorporation of a gene co-expression network, deduced through transcription factor motif sequence validation statistical methods, which are capable of identifying highly reliable, TF-driven GRNs. Using SCENIC analysis, we further emphasized the importance of FOS.

Statistical analysis

Statistically significant difference is calculated in R. P-values of less than 0.05, 0.01 and 0.001 were considered statistically significant.

Data and code availability

The data supporting the research findings can be found in the main text and supplementary manuscript. Other data could be obtained by contacting the corresponding author.

Consent for publication

The authors agree with the terms of conditions for publication.

Characteristics of PBMCs before and after MSC therapy for HF patients

To better understand the circulating system changes of heart failure patients before and after stem cell therapy, we conducted scRNA-seq for bettering depicting the expression pattern of peripheral blood mononuclear cells (PBMCs) of HF patients. As showed in Figure 1A and Supplementary Tables 1, 2, we collected four patients who have received MSC therapy (Supplementary Tables 1, 2). As a result, there were 81226 single cells obtained after stringent criteria (Supplementary Figures 1, 2). As shown in Figure 1A, 1B, the detected single cells were annotated into different cell types using “Weighted nearest neighbor” analysis and the referenced UMAP plot from the multimodal PBMC reference was adopted for better visualizing the scRNA-seq data [19]. Furthermore, we tried to explore the expression characteristics of different cell types, especially changes of cell proportion and cell count. Intriguingly, the proportion of CD14⁺ monocytes was down-regulated after MSC therapy, no matter in the response group or non-response group (Figure 1C, 1D). Moreover, the total cell number of CD14⁺ monocytes (17048 cells) was more than other cell types (Figure 1E). Then, the metabolism analysis of CD14⁺ monocytes was performed and results showed that many metabolism pathways varied among four groups (Supplementary Figure 3), In brief, CD14⁺ monocytes might be affected most with MSCs therapy.

Figure 1. scRNA-seq reveals the different cell types and distribution among patients with or without MSCs treatment. (A) Cell annotation using the azimuth package (https://app.azimuth.hubmapconsortium.org/app/human-pbmc). And the referenced UMAP was leveraged to better visualized the results. (B) Cell distribution for all four groups. (C) Cell proportion of all cell types in four groups. (D) Bar plot shows the cell proportion of all cell types in four groups. (E) Cell distribution analysis indicates the largest cell proportion of CD14⁺ monocytes.

Nonnegative matrix factorization (NMF) analysis showed the heterogeneity of response and non-response CD14⁺ monocytes

To compare the different of CD14⁺ monocytes before and after MSC therapy, we utilized three methods to obtain conserved markers of CD14⁺ monocytes in different groups: (i) method 1: the “findallmarker” function in Seurat package [19] was used, which would compare the interested cell cluster to other all cells using Wilcox algorithm; (ii) method 2: the “findmarker” function Seurat package was adopted, and then we compared the interested cell cluster one by one using Wilcox algorithm. At last, the intersected marker genes were considered specific to the interested cell cluster; (iii) method 3: Metacell algorithm was applied and we obtained the “de-noised” metacells [20, 21] and thereafter we could use the DESeq2 package [27] for differential expression analysis. With these three methods, we could obtain the highly conserved markers for the four group of CD14⁺ monocytes.

Subsequently, with the highly conserved markers, we run the NMF analysis for understanding the heterogeneity and potential functions of CD14⁺ monocytes (Figure 2B, 2C). Of interest, the response and non-response groups were obviously separated using the UMAP plot (Figure 2B), indicating the difference of CD14⁺ monocytes in response and non-response group. However, the CD14⁺ monocytes before and after MSC therapies in the response group were quite similar while relatively distinct in the non-response group (Figure 2C), demonstrating that CD14⁺ monocytes might be capable of serving as important indicators of whether and which HF patients might be suitable for MSC therapy or could be benefit from MSC therapy.

Figure 2. Subcluster analysis of the CD14⁺ monocytes. (A) Identification of highly-conserved DEGs/markers for CD14⁺ monocytes in the four groups, respectively. Three approaches for exploring DEGs were adopted. DEGs: differentially expressed genes. (B) NMF analysis for the reduction analysis of CD14⁺ monocytes based on all highly-conserved DEGs/markers. NMF: Nonnegative Matrix Factorization. (C) Distribution of all CD14⁺ monocytes from four groups. (D) Cell proportion of CD14⁺ monocytes was visualized, comprising 11 subclusters. (E) Cell distribution of all 11 identified subclusters. (F) Cell distribution of all 11 identified subclusters in the four groups. The CD14⁺ monocytes from response and non-response patients were quite different. (G) GSEA pipeline for showing the enriched terms for CD14⁺ monocytes from the “before response” group.

For better investigating the characteristics of CD14⁺ monocytes in the “before_response” group, we utilized the subcluster analysis based on the NMF dimension reduction. As a result, 11 subclusters of CD14⁺ monocytes were obtained (Figure 2D–2F). Then, we used GO (Supplementary Figure 4), KEGG and Reactome (Supplementary Figure 5) databases for function enrichment analyses. Interestingly, we confirmed the various function of different subclusters. For instance, the cluster 0 was enriched in “IL−17 signaling pathway”, “Neutrophil degranulation” and “cellular response to chemical stress”, while cluster 1 was enriched in “Viral myocarditis”, “Interferon Signaling” and “positive regulation of leukocyte cell− cell adhesion”. Then, we also found that cluster 0 (2786 cells) was predominantly enriched in the “before_response” group and would be down-regulated after MSC treatment (Figure 2F). Thus, we identified cluster 0 as “therapy-response CD14⁺ monocytes”. Further, we leveraged the GSEA analysis for running the enrichment analysis of HALLMARK datasets (downloaded from http://www.gsea-msigdb.org). Of interest, the potential function of response CD14⁺ monocytes were mainly related to “COMPLEMENT”, “INFLAMMATORY RESPONSE”, “INTERFERON_GAMMA_RESPONSE” and so on, which were highly related to the therapeutic mechanism of MSCs (Figure 2G).

Machine learning algorithms were leveraged to construct the predictive model for therapy-response CD14⁺ monocytes and identified potential indicators

To identify the therapy-response CD14⁺ monocytes in single cell levels, we used the metacell algorithm [20, 21] and divided CD14⁺ monocytes into training cohort and test cohort. Using LASSO algorithm, we had filtered the 77 conserved markers and only 20 features were left (Figure 3A, 3B), including IFI6, PTAFR, S100A12, FOXN2, MAP3K1, CAST, F13A1, HLA-DPB1, AKIRIN2, ANXA1, DUSP6, EPSTI1, FOS, SOCS3, RFX2, TPM4, ZFP36, PPDPF, MX2 and RPL39. Then, we used the “Bootstrap” algorithm for internal validation for 1000 times and the performance was reliable (Figure 3C). To avoid the bias effect of different algorithms, we further used another seven machine algorithms for constructing the classifier for therapy-response CD14⁺ monocytes (Figure 3D). Intriguingly, the performance of all seven machine algorithms was quite robust, especially SVM algorithm (Figure 3E). Then, we used the test cohort for further validation, similar result was obtained that the performance of the constructed classifier was good (AUC: 0.98) (Figure 3F).

Figure 3. Machine learning algorithms were leveraged to construct a predictive model for therapy-response CD14⁺ monocytes and explored potential indicators. (A) LASSO algorithm was adopted for filtering optimal feature genes. Only the 21 highly-conserved makers/DEGs for therapy-response CD14⁺ monocytes were adopted. (B) Coefficients of identified feature genes from LASSO algorithm was shown. (C) ROC analysis for the LASSO model. LASSO: Least absolute shrinkage and selection operator. (D) Machine learning algorithm for constructing a predictive model for therapy-response CD14⁺ monocytes. (E) ROC analysis for the machine learning model. (F) ROC analysis for the machine learning model in test cohort. The AUC value was 0.98. AUC: area under curve. (G) Intersection of genes identified by the LASSO algorithm with the highly conserved markers of the therapy-response CD14⁺ monocytes. Seven overlapped genes were delineated. (H) Protein-protein interaction network for identified genes. The network prediction was based upon an online web-server: GeneMANIA (http://www.genemania.org).

Then we also intersected the machine learning-identified feature genes and the conserved up-regulated genes in the “before_response” group (Figure 3G). Result showed that PTAFR, FOXN2, AKIRIN2, ANXA1, DUSP6, FOS and RFX2 were obtained (Figure 3G and Supplementary Figure 6). Then, we used the GeneMANIA for constructing real-time multiple association network integration [28] (Figure 3H). Of note, the FOS/JUN-related transcription factors were enriched.

Trajectory analysis of therapy-response CD14⁺ monocytes

Then, we used the “addmodule” function in Seurat package for calculating the signature score of the seven therapy-response-specific and up-regulated genes (Figure 4A–4C). Of note, the signature score was higher in CD14⁺ monocytes from the “before_response” group (Figure 4A) and was also higher in cluster 0 (therapy-response CD14⁺ monocytes) (Figure 4B, 4C). Then, we performed the pseudotime analysis using monocle2 algorithm [24]. Intriguingly, we could find that the cluster 0 was mainly distributed in the early and middle timepoint of the pseudotime plot (Figure 4D, 4E). Then, we detected the expression pattern of the seven genes, most of which were consistent with the pseudotime distribution (Figure 4F). Considering that the middle stage character and potential FOS/JUN-related function of therapy-response CD14⁺ Monocytes, which indicated that the therapy-response CD14⁺ Monocytes might be much more flexible and its state would be shifted, it triggered us to understand whether its status might be similar to the published monocyte signature (MS) reported by Reyes et al. [29], whose research has found four monocyte signatures in the PBMC of sepsis patients and considered MS1 as a unique disease-related CD14⁺ monocytes state that is expanded in sepsis patients and has validated its power in distinguishing sepsis patients from healthy controls [29]. Interestingly, we found that the therapy-response CD14⁺ monocytes also had higher MS1 score than other CD14⁺ monocytes (Figure 4G) as well as positively correlated with the MS1 score (Figure 4H). In summary, the state of therapy-response CD14⁺ monocytes might also be different from other CD14⁺ monocytes, and might be an indicator of those HF patients with could potentially benefit from MSC therapy.

Figure 4. Trajectory analysis of CD14⁺ monocytes and the pseudotime changes of the feature genes identified by machine learning algorithms. (A) Violin plot displays the response-related CD14⁺ monocyte signature in pseudo-bulk levels of CD14⁺ monocytes. (B) Violin plot displays the response-related CD14⁺ monocyte signature across all identified subclusters of CD14⁺ monocytes. “AddModule” function in Seurat package was leveraged to estimate the signature score. (C) Signature score of response-related CD14⁺ monocytes was mainly activated in cluster 0. (D) Trajectory analysis of CD14⁺ monocytes. Monocle2 package was adopted. (E) Pseudotime distribution of other and response-related CD14⁺ monocytes. (F) Expression changes of the feature genes following pseudotime. (G) Violin plot shows the monocyte signature scores. MS: monocyte signature. (H) Correlation analysis between the response signature and the other four published MS phenotypes.

SCENIC analysis revealed potential regulators for therapy-response CD14⁺ monocytes

Then, we tried to explore the transcriptional regulatory network of CD14⁺ monocytes. As a result, the activities of several TFs were different between therapy-response CD14⁺ monocytes and other CD14⁺ monocytes (Figure 5A). The top TFs were showed in Figure 5B, 5C. Data showed that therapy-response CD14⁺ monocytes were mainly regulated by JUN, FOS, CEBPD while other CD14⁺ monocytes were predominantly controlled by IRF7, STAT1 and BCL3 (Figure 5C). More than this, therapy-response CD14⁺ monocytes could be aggregated in the UMAP plot based on regulon activity scores calculated by SCENIC (Figure 5D). Noteworthily, the FOS regulons (FOS_extend_(37g) and FOS_(12g)) and FOS mRNA expression level were all highly up-regulated in therapy-response CD14⁺ monocytes (Figure 5E, 5F), underpinning the crucial roles of FOS in regulating the functions of therapy-response CD14⁺ monocyte signature. Taking it a step further, we downloaded published sequencing data (GSE59867) of PBMC related to post-MI HF patients and result showed significant correlations between detected TFs and the seven genes (Figure 5G). To further understand the role of FOS, we also predicted the potential drugs for both capable of activating FOS and also treating heart diseases (Supplementary Table 3). Patients who have adopted these drugs for HF therapy might be much suitable for MSC therapy.

Figure 5. Transcription factor analysis for CD14⁺ monocytes. (A) SCENIC pipeline for identifying potential regulatory TFs for therapy-response CD14⁺ monocytes and heat map was used for visualization. SCENIC: Single-Cell rEgulatory Network Inference and Clustering; TF: Transcription factor. (B) Violin plot for showing the identified TFs between other and therapy-response CD14⁺ monocytes. (C) Ranking of the potential TFs for other and therapy-response CD14⁺ monocytes. (D) Cell distribution based on the identified regulons/TFs. (E) Activities of FOS regulons in UMAP plots. Two regulons of FOS were considered significant. (F) Expression levels of FOS shown by UMAP plot. (G) Correlation analyses among identified TFs and hub feature genes.

NMF subtype analyses reveals monocyte signatures were capable of classifying post-MI HF patients into three groups

Considering the pivotal roles of monocytes in post-MI heart failure, we tried to further discover whether the conserved markers we identified were powerful enough to distinguish HF patients from healthy patients. Herein, we also used the GSE59867 data, which contained PBMC samples from stable coronary diseases to heart failure. Using NMF analysis [30] (Supplementary Figure 7), we had classified stable CAD patients and post-MI patients into three subgroups based on the 77 conserved markers of CD14⁺ monocytes (Figure 6A, 6B). Interestingly, FOS was highly expressed in subgroup 3 (Figure 6B). Then, we used CIBERSOFT algorithm in IOBR packages [31] to deconvolute the bulk sequencing data for understanding the immune filtration characters (Figure 6C) response CD14⁺ monocyte signature and the signature score was higher in HF patients compared to stable CAD patients (Figure 6D). Correlation analysis reveals that CIBERSOFT-inferred monocytes were positively correlated with all seven therapy-response-related genes (Figure 6E) and monocytes were also highly expressed in subgroup 3 (Figure 6C, 6F) compared to subgroup 1 and 2. Consistently, the therapy-response CD14⁺ monocyte signature was also highly expressed in subgroup 3 (Figure 6G). Finally, the distribution of patients was showed (Figure 6H). In brief, CD14⁺ monocyte participated in the development of post-MI HF patients and the therapy-response CD14⁺ monocyte signature might be considered as typical changes of subgroup 3.

Figure 6. NMF subtype analyses uncovers that therapy-response monocyte signatures were capable of partitioning post-MI HF patients into three groups. (A) Consensus clustering for the GSE59867 data based on the hub genes, which contained PBMC samples from stable coronary diseases to heart failure. (B) Heatmap shows the differences among the three patient groups. (C) Heat map displays the infiltration levels of immune cells among the three patient groups. CIBERSOFT algorithm in IOBR packages was adopted. (D) Boxplot shows the signature score of therapy-response monocytes in bulk RNA-seq levels between control and heart failure groups. (E) Correlation analysis of the hub genes and monocytes in bulk RNA-seq levels. The monocytes were highlighted. (F) Infiltration levels of deconvoluted immune cells in the NMF groups. (G) Box plot reveals the signature score of therapy-response monocytes among the three patient groups. (H) Distribution of different clinical information in the NMF groups. NMF: Nonnegative Matrix Factorization.