Single-Cell Sequencing Revealed Pivotal Genes Related to Prognosis of Myocardial Infarction Patients.

Objectives: Myocardial infarction (MI) is a common cardiovascular disease. Histopathology is a main molecular characteristic of MI, but often, differences between various cell subsets have been neglected. Under this premise, MI-related molecular biomarkers were screened using single-cell sequencing.
Methods: This work examined immune cell abundance in normal and MI samples from GSE109048 and determined differences in the activated mast cells and activated CD4 memory T cells, resting mast cells. Weighted gene coexpression network analysis (WGCNA) demonstrated that activated CD4 memory T cells were the most closely related to the turquoise module, and 10 hub genes were screened. Single-cell sequencing data (scRNA-seq) of MI were examined. We used t-distributed stochastic neighbor embedding (t-SNE) for cell clustering.
Results: We obtained 8 cell subpopulations, each of which had different marker genes. 7 out of the 10 hub genes were detected by single-cell sequencing analysis. The expression quantity and proportion of the 7 genes were different in 8 cell clusters.
Conclusion: In general, our study revealed the immune characteristics and determined 7 prognostic markers for MI at the single-cell level, providing a new understanding of the molecular characteristics and mechanism of MI.

1. Introduction

Myocardial infarction (MI) as one of the most common cardiovascular diseases is often caused by myocardial necrosis due to coronary artery occlusion [1]. MI is often clinically diagnosed according to persistent chest pain, elevated cardiac troponin level, regional wall motion abnormalities, ECG recordings, and angiographic detection of coronary thrombus [2]. Compared with healthy people, patients with MI face 4 to 6 times higher chances of death from cardiac arrest [3]. In recent decades, early intervention after MI, especially increased use of myocardial reperfusion therapy, greatly improves MI patients' survival; however, myocardial reperfusion can still induce myocardial damage [4], thereby contributing to high MI-related mortality. New strategies should be developed to limit the MI area and prevent adverse left ventricular remodeling and heart failure after MI [5].

Immune cells have been reported to be involved in the risk of cardiovascular disease. In the process of atherosclerosis, monocyte-derived macrophages and Thelpertype1 cells will be triggered by subendothelial retention of plasma lipoprotein in the arterial wall to form atherosclerotic plaques [6]. McMaster et al. reported that terminal organ dysfunction and damage in hypertension could result from both congenital and adaptive immune system cells [7]. The infiltration of immune cells in the myocardium could adversely affect normal heart function, causing heart failure pathogenesis [8]. The immune system is an important regulator of MI. After MI, the innate immune response could be activated by extensive necrosis of ischemic cardiomyocytes; then a large number of immune cells, such as monocytes and neutrophils, will be recruited to the heart to help clear pathogens and dead tissue and improve healing [9, 10]. Therefore, the characterization of specific immune cell subpopulations may offer a new possibility for MI management.

As a polygenic disease, the occurrence of AMI is the result of the interaction of genetic and environmental factors [11]. At present, several biomarkers have been identified for the clinical diagnosis of AMI-CAD, and cardiac troponin I (cTnI) is considered the gold standard for clinical diagnosis [12]. The limitation of cTnI detection is that the level delay increases 3-4 hours after AMI, and it takes 6-12 hours of continuous measurement to overcome the lack of accuracy of biomarkers [13]. Therefore, the diagnosis efficiency and performance are low and the cost is high. At the same time, cTnI levels increased in patients with heart failure, sepsis, and chronic kidney disease, reducing the specificity of the test [14]. As another basic index of AMI, creatine kinase myocardial band (CK-MB), is also widely used in clinical diagnosis [15]. However, in some diseases, such as congestive heart failure, myocarditis, renal failure, and skeletal muscle tissue injury, serum CK-MB levels are also increased, reducing the diagnostic specificity [16]. Therefore, it is urgent to find new diagnostic biomarkers with high sensitivity, specificity, and accuracy in AMI-CAD. The single-cell RNA-seq (scRNA-seq) technology allows us to have a wide range of molecular definitions and characterizations of immune cell populations [17]. This study analyzed the single-cell sequencing data of 19 MI in GSE180678 to depict immune cell populations and reveal an immune landscape of MI. By comparing markers and enrichment pathways in different cell types, unique components and characteristics of each cell population were determined. In addition, this study also determined new biomarkers for predicting the prognosis of MI.

2. Materials and Methods

2.1. Research Data Acquisition

GSE109048 [18] chip data and GSE180678 [19] single-cell sequencing data were downloaded from the NCBI GEO database, and 19 healthy samples and 19 myocardial infarction samples were obtained by retaining the information of healthy samples and myocardial infarction samples. Genes corresponding to the probes were extracted, while the probes corresponding to multiple genes were removed. The average values of multiple probes corresponding to one gene were taken. The flow chart of data analysis in this study is shown (Figure S1A).

2.2. Difference Analysis

For the annotated chip data, the standard deviation of gene expression in all samples was greater than 0.1. The genes in the range of p < 0.05 were screened by the R software package Limma [20].

2.3. Calculation of Immune Cell Abundance

CIBERSORT, a deconvolution algorithm, has 547 genes from 22 human hematopoietic cell phenotypes [21]. The abundance of 22 types of immune cells in healthy samples and myocardial infarction samples was quantified using CIBERSORT, and they were compared by performing the Wilcoxon test.

2.4. Weighted Gene Coexpression Network Analysis (WGCNA)

The WGCNA supports the exploration of module structure in the network and relationship measurement between modules and between genes and modules [22]. The samples were clustered by hierarchical clustering according to the RNA-seq data of 38 samples, and Pearson's correlation analysis was used to measure the distance between each gene. Scale-free topological fitting index (R2) > 0.7 was set as the threshold in WGCNA in the R software package for weight coexpression network development. The weight parameter (β) of the adjacency matrix was selected by the pickSoftThreshold function to generate a weighted adjacency matrix. Based on the dissimilarity of the topological overlap matrix, the gene tree was developed, and there were at least 200 genes in each module. The coexpression modules that meet the conditions (height = 0.25, deepSplit = 2, and minModuleSize = 30) were identified by DynamicTreeCut function [23].

2.5. Functional Enrichment Analysis

The module showing the strongest correlation with activated CD4 memory T cells was screened by correlation analysis, and the Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) with the WebGestaltR [24] (http://www.webgestalt.org) tool were applied to explore the genes in this module. The statistical cutoff value was P < 0.05. Pathway and process of enrichment were shown in the form of a bubble diagram after visualization.

2.6. Reduction and Clustering of Single-Cell Data

The single-cell data were filtered according to the threshold that each gene was expressed in at least 3 cells and at least 250 genes were expressed in each cell. The ratio of mitochondria to rRNA was calculated by the PercentageFeatureSet function, and the cells with more than 30% mitochondrial genes but fewer than 500 expressed genes were filtered. The NormalizeData function was used for data normalization. All genes were scaled by the ScaleData function, and cells were classified by PCA. When dim = 30, the cells were clustered according to the FindNeighbors and FindClusters functions, and the FindAllMarkers function in the Seurat R software package was applied to count the difference multiple at logfc = 0.5. The minimum expression ratio of differential genes Minpct = 0.25 was the threshold to select marker genes from each cell subpopulation.

2.7. Cell Types in Different Subgroups

Human cell marker genes were downloaded from the CellMarker [25] official website (http://biocc.hrbmu.edu.cn/CellMarker/). The corresponding tissues of the myocardium, heart, and blood were selected, and the cell types of different subsets were determined by the enricher function of the clusterProfiler package [26].

2.8. Prediction of Cell Development Trajectory

For different subsets of cells, the monocle was employed to predict trajectories of pseudotime development [27] by using reversed graph embedding to show multiple fate decisions under full unsupervision. Different branches in developmental trajectory were closely observed to identify cells on the same branch, which we defined as in a shared state. Cell differentiation degree was measured by Pseudotime.

2.9. ReactomeGSA Pathway Analysis

ReactomeGSA, a tool for multispecies, multigroup approach analysis, reveals key biological mechanisms through integrating large “genomics” data sets [28]. We used the reactome database to analyze single-cell data. In this study, data of different species were mapped to a common path space. The enrichment scores of different subpopulations in all pathways were calculated, and the top 10 pathways with the greatest differences in enrichment scores were compared through drawing heat maps.

2.10. Statistical Analysis

All bioinformatics analysis of the scRNA-seq data was carried out in the R software package. P value was adjusted according to the Benjamini-Hochberg program. Statistical significance was defined when the value was lower than 0.05.

3. Results

3.1. Screening of Modules Related to Immune Cells by WGCNA

The abundance difference of 22 types of immune cells from MI and normal samples was analyzed. Significant differences in the abundance of resting mast cells, activated CD4 memory T cells, and activated mast cells between the two groups were detected (Figure 1(a)). The soft threshold of 4 was decided for ensuring a scale-free nature of the network (Figure 1(b)). Genes in GSE109048 were clustered into 13 modules (Figure 1(c)). The relationship between the modules and the three types of cells was analyzed. Interestingly, the Brown module showed the most significant correlation with resting mast cells, and the turquoise module was strongly associated with activated CD4 memory T cells. All 13 modules showed no significant correlation with activated mast cells (Figure 1(d)).

Figure 1

Screening of modules associated with immune cells. (a) Wilcoxon's test compared the abundance of 22 types of immune cells between normal and myocardial infarction samples. (b) Average connectivity of selected β and different β values when R2 reaches 0.7 for the first time. (c) Clustering dendrogram of comethylation modules from GSE109048. (d) Correlation analysis of 13 modules with resting mast cells, activated mast cells, and activated CD4 memory T cells.

3.2. Function of the Turquoise Module and Identification of Hub Genes

GO and KEGG analyses were performed on 5200 genes in the turquoise module that was associated with activated CD4 memory T cells. The results showed that the cells were annotated to 117 molecular functions (MFs), 251 cellular components (CCs), 301 biological processes (BPs), and 24 KEGG pathways. These GO items and KEGG pathway were related to RNA activity and translation during cytogenesis (Figures 2(a) and 2(b)). 341 genes were screened from the turquoise module based on the absolute value of correlation between phenotypic traits (GS) and gene > 0.2 and relationship between gene expression profile (MM) and modular eigengene > 0.4 (Figure 2(c)). Among them, only 10 genes were found to be significantly different between normal samples and myocardial infarction samples. To understand the interaction between the 10 hub genes, a gene-gene interaction network was developed (Figure 2(d)). Among the 10 hub genes, PCDHA6 was expressed higher in normal samples than that in samples of myocardial infarction, and the other 9 genes were significantly expressed higher in patients with MI (Figure 2(e)). Pearson's correlation analysis showed a negative correlation between PCDHA6 and 9 hub genes, and only a significant negative correlation of PCDHA6 with RINT1 (Figure 2(f)).

Figure 2

The function of the turquoise module and the interaction of hub gene. (a) The GO analysis of genes in the turquoise module includes three aspects: molecular function, cellular component, and biological process. (b) The KEGG pathway enrichment analysis on the turquoise module genes. (c) Scatter plot of turquoise-GS-MM distribution. (d) The turquoise module gene-gene interaction network diagram. (e) The differential expression of 10 hub genes in normal and myocardial infarction samples was analyzed. (f) Pearson's correlation analysis among 10 hub genes.

3.3. Single-Cell Data Analysis Revealed the Immune Landscape of Human Myocardial Infarction

We obtained single-cell data to characterize the immune landscape of human MI. After initial quality control, 1084 cells and 15726 genes were included in our analysis (Figure S1BC). According to the FindVariableFeatures function, the distribution of hypervariable genes and nonhypervariable genes was analyzed, and the first 20 hypervariable genes were shown (Figure S2). We used t-distributed stochastic neighbor embedding (t-SNE) for cell clustering into 8 subgroups (Figure 3(a)). Cell types of 8 subpopulations were annotated, and subpopulations 0-7 were hematopoietic stem cell, CD1C+_B dendritic cell, MYH11+ CNN1+ smooth muscle cell, CDH5+ HEXIM1+ smooth muscle cell, natural killer cell, CD1C−CD141− dendritic cell, endothelial cell, and cardiomyocyte, respectively (Figure 3(b)). The marker gene of each subgroup was identified (Table S1), and the expression of the five most prominent marker genes of each subgroup in different subgroups was analyzed by a thermograph. It was found that each subgroup had different gene expression characteristics. Markers including TAGLN, MYL9, ACTA2, TPM2, RGS5, HIGD18B, NDUFA4L2, COX4I2, and GJA4 were shared by subpopulation 2 and subpopulation 3 (Figure 4). In addition, subpopulation 2 specifically expressed MYH11 and CNN1, while subpopulation 3 specifically expressed CDH5 and HEXIM1 (Figure 3(c)). The developmental trajectory of different cell subsets was predicted by the monocle. The developmental trajectories of all cell subsets were distributed in three branches. We found that CD1C+_B dendritic cell and MYH11+ CNN1+ smooth muscle cell were mainly distributed in branch 1, CDH5+ HEXIM1+ smooth muscle cell, hematopoietic stem cell, and endothelial cell were mainly distributed in branch 2, and cardiomyocytes were distributed in both branches. In addition, natural killer cells and CD1C-CD141-dendritic cells appeared in the same branch (Figure 3(d)).

Figure 3

scRNA-seq data analysis of human myocardial infarction. (a) The cell samples were divided into 8 subgroups by the t-SNE algorithm. (b) Eight cell subpopulations were annotated according to the marker. (c) The violin picture showed the expression patterns of specific genes in two kinds of smooth muscle cells. (d) The cell development trajectories of 8 subpopulations were predicted by the monocle. Each point represents a cell and is marked with pseudotime (right) and the cell state (left).

Figure 4

Expression heat maps of the first five marker genes in each cell subpopulation in different cell clusters.

3.4. Expression Distribution of Hub Genes and Chemokines in 8 Cell Clusters

There were 10 hub genes related to activated CD4 memory T cells in the model of turquoise. Of the 10 hub genes, 7 were detected in single-cell sequencing analysis. Figure 5(a) showed the distribution of these 7 genes in 8 cell clusters, and the proportion of their expression in subgroup 0 (hematopoietic stem cell) and subgroup 1 (CD1C+_B dendritic cell) was not high. The proportion of high LDOC1 expression was observed only in subgroups 2 (MYH11+ CNN1+ smooth muscle cell) and 7 (cardiomyocyte). The expression rates of PCBD2, PGBD4, and RINT1 in all cell clusters were low. The high proportion and expression level of RMDN1, SH3BP5, and SLC43A3 in are shown subgroup 5 (natural killer cell). The expression proportion and expression level of SLC43A3 in subgroup 2 (CDH5+ HEXIM1+ smooth muscle cell) were the highest among all 7 hub genes.

Figure 5

Expression distribution of Hub gene and chemokines in 8 cell subpopulations. (a) Bubble map of 7 hub genes expressed in 8 cell subpopulations. (b) The expression level of chemokines and their receptors in each cell cluster. The dot color demonstrates the average intensity of expression, and the dot size shows the quantity of cells expressing the gene.

In addition, 41 chemokine and 18 chemokine receptor genes were extracted from the study of Ru et al. [29] to their expressions in 8 cell clusters. The results showed that CCL3, CXCL3, CXCL8, and CXCL16 were highly expressed in CD1C-CD141- dendritic cells and that CCL4 and CCL5 were highly expressed in natural killer cells. Moreover, the CXCR4 receptor showed a high expression in CD1C-CD141- dendritic cells (Figure 5(b)).

3.5. Enrichment Pathway Analysis of Different Cell Subpopulations

We also calculated the enrichment scores (ES) of different subpopulations in all pathways and found the 10 pathways with the greatest difference in ES. Analysis on the enrichment pathway activity of each cell subpopulation (Figure 6(a)) showed a limited difference between the two smooth muscle cell enrichment pathways, but only in the resolvin conjugates and biosynthesis of protectin in tissue regeneration (PCTR and RCTR) pathway, the ES of MYH11+ CNN1+ smooth muscle cell was higher than that of CDH5+ HEXIM1+ smooth muscle cell. A great similarity was observed between endothelial cell and hematopoietic stem cell abundance. The enrichment pathways of cardiomyocyte and CD1C B dendritic cells were similar, but the main difference between the two was that the score of COX reactions enriched by cardiomyocytes was higher, and the score of FGFR1c and Klotho ligand binding and activation enriched by CD1C B dendritic cell was higher. The enrichment pathways of CD1C-CD141- dendritic cells and natural killer cells were also similar, and the difference between the two was largely in classical antibody-mediated complement activation and metabolism of vasopressin. Moreover, cardiomyocytes, endothelial cells, and CD1C B dendritic cells showed higher ES in the metabolism of serotonin, choline catabolism, and PCTR and RCTR pathways, indicating that these three kinds of cells were important contributors to anabolism (Figure 6(b)).

Figure 6

Characteristics of enrichment pathways in different cell subpopulations. (a) Enrichment activity of each cell subpopulation in the first 10 pathways. (b) Specific enrichment scores of 8 kinds of cell subpopulations on the 10 pathways with the greatest difference.

4. Discussion

Great efforts have been devoted to the understanding of different types of cardiovascular disease and to identifying the presence of high-risk cardiovascular disease. Nevertheless, we still need highly sensitive biomarkers to make the diagnosis of MI early and accurately, which at the same time can reduce the exposure of patients to extensive testing or radiation. But so far, biomarkers of this kind are rare [30]. Here, we first divided the genes from MI samples into 13 modules by the WGCNA. In view of differences in the abundance of activated mast cells, resting mast cells, and activated CD4 memory T cells between normal samples and MI samples, the turquoise module, which was the most strongly related to activated CD4 memory T cells, was detected. In the turquoise module, there were 10 genes differentially expressed between MI and normal samples. It is preliminarily predicted that these 10 genes were potential markers of MI.

Previous studies have mostly identified disease markers from tissues of different states, such as primary and normal tissues, with high and low immune infiltration, but in this way, the differences between various cell subsets could be easily neglected and may miss important biomarkers [31]. The application of scRNA-seq has enabled the use of single-cell precision to treatment and diagnosis of cardiovascular diseases. scRNA-seq allows analysis of individual cell phenotypes to discover cell subpopulations that may result in cardiovascular pathogenesis [32]. In this study, we distinguished differentially expressed genes according to different tissue states, analyzed gene expression of 1084 cells in MI patients at the single-cell level, and finally determined 8 cell clusters. A single cell constantly goes through a dynamic process and responds to various environmental stimuli, which is particularly reflected in the molecular characteristics of the cell [33]. Molecular marker analysis of subpopulations confirmed that each subpopulation had different gene expression characteristics and regulated the activity of multiple pathways, including immune-related pathways, and synthetic metabolism-related pathways.

Of the 10 hub genes identified by the WGCNA in the turquoise module, 7 were detected by single-cell sequencing analysis. The expression quantity and proportion of them were different in 8 cell subpopulations. It is reported that after MI, dead cells and damaged matrix release danger signals, activate complement cascades, and stimulate Toll-like receptor/interleukin-1 signal, resulting in nuclear factor-κ B system activation and induction of chemokine, cytokines, and adhesion molecules [34]. Our analysis on the expression of chemokine and its receptor genes in MI samples showed that chemokines CCL3, CXCL3, CXCL8, and CXCL16 were highly expressed in CD1C-CD141- dendritic cells and that CCL4 and CCL5 were highly expressed in natural killer cells. These chemotactic factors are involved in the progress of MI [35-38]. Therefore, dendritic cells and natural killer cells may have a critical function in MI-induced chemokines. In addition, we also observed that seven genes were significantly overexpressed in AMI samples (Figure S3A), which was more conducive to clinical detection. It is worth mentioning that SH3BP5 was significantly positively correlated with naïve B cells and activated CD4 memory T cells, and SLC43A3 was significantly positively correlated with gamma delta T cells and CD4 memory resting T cells, LDOC1 and PCBD2 are significantly positively correlated with activated CD4 memory T cells, and PGBD4 is significantly positively correlated with gamma delta T cells (Figure S3B). These results show that these genes are fully involved in the process of immune regulation. The interaction between these seven genes is analyzed by using the STRING database, so it can be observed that there is no direct interaction between these genes. This suggests that these genes may be independently involved in different biological processes (Figure S3C). Literature mining shows that LDOC1 is associated with many diseases, such as reproductive-related tumors [39, 40]. RINT1 is a new moderately permeable cancer susceptibility gene. It is common in breast cancer and possibly associated with Lynch's syndrome. The loss of RINT1 affects the fate of PDAC cells by interfering with the SUMOylation pathway, and in acute liver toxicity [41, 42]. In acute hepatotoxicity, mitochondrial SH3BP5, as the target of JNK, maintains its activation by promoting the production of reactive oxygen species [43]. The overexpression of Slc43a3 reduces FA uptake in differentiated OP9 cells and leads to reduced lipid droplet accumulation [44]. These results show that the seven genes play an important role in a variety of complex diseases.

The current work also contained certain limitations that should be equally noted. Specifically, the sample size of the current work was small and came from a public database. Secondly, we lack a systematic correlation analysis investigating the association of hub gene expression with clinical features, and also, the results were not verified in solid MI tissues.

5. Conclusions

In summary, this study described the immune landscape of MI and identified the marker genes related to the prognosis of MI at the level of a single cell, offering references for understanding the influencing factors as well as characterization of MI patients.