Integrating the characteristic genes of macrophage pseudotime analysis in single-cell RNA-seq to construct a prediction model of atherosclerosis

Abstract

Background: Macrophages play an important role in the occurrence and development of atherosclerosis. However, few existing studies have deliberately analyzed the changes in characteristic genes in the process of macrophage phenotype transformation.

Method: Carotid atherosclerotic plaque single-cell RNA (scRNA) sequencing data were analyzed to define the cells involved and determine their transcriptomic characteristics. KEGG enrichment analysis, CIBERSORT, ESTIMATE, support vector machine (SVM), random forest (RF), and weighted correlation network analysis (WGCNA) were applied to bulk sequencing data. All data were downloaded from Gene Expression Omnibus (GEO).

Result: Nine cell clusters were identified. M1 macrophages, M2 macrophages, and M2/M1 macrophages were identified as three clusters within the macrophages. According to pseudotime analysis, M2/M1 macrophages and M2 macrophages can be transformed into M1 macrophages. The ROC curve values of the six genes in the test group were statistically significant (AUC (IL1RN): 0.899, 95% CI: 0.764-0.990; AUC (NRP1): 0.817, 95% CI: 0.620-0.971; AUC (TAGLN): 0.846, 95% CI: 0.678-0.971; AUC (SPARCL1): 0.825, 95% CI: 0.620-0.988; AUC (EMP2): 0.808, 95% CI: 0.630-0.947; AUC (ACTA2): 0.784, 95% CI: 0.591-0.938). The atherosclerosis prediction model showed significant statistical significance in both the train group (AUC: 0.909, 95% CI: 0.842-0.967) and the test group (AUC: 0.812, 95% CI: 0.630-0.966).

Conclusions: IL1RN^High M1, NRP1^High M2, ACTA2^High M2/M1, EMP2^High M1/M1, SPACL1^High M2/M1 and TAGLN^High M2/M1 macrophages play key roles in the occurrence and development of arterial atherosclerosis. These marker genes of macrophage phenotypic transformation can also be used to establish a model to predict the occurrence of atherosclerosis.

Keywords: atherosclerosis; macrophage; phenotypic transformation; pseudotime analysis.

Figure 1

Single-cell transcriptome data: (A, B) The single cells were divided into 45 groups on the basis of their transcriptome data and then ultimately divided into 9 cell populations. (C) Circle chart: The X-axis represents the marker genes that define the cells, and the Y-axis represents the different cell populations.

Figure 2

(A, B) The macrophage population was divided into 13 clusters, which were finally categorized into 3 cell populations (M1, M2, M2/M1). (C) Differential gene analysis among different subtypes of macrophages: The X-axis represents the macrophage subtypes, and the Y-axis represents the differentially expressed genes. (D) Violin plot of macrophage marker genes.

Figure 3

(A–C) Trajectory analysis of M1 macrophages and M2 macrophages. (D) The trajectory analysis of the heatmap of M1 macrophages versus M2 macrophages: The X-axis represents the timeline of trajectory analysis, the left Y-axis represents the KEGG enrichment results, and the right Y-axis represents the differentially expressed genes between the two clusters. (E–G) Pseudotime analysis of genes (MRC1, IL1B, and MARCO): The X-axis represents the cell of trajectory analysis, and the Y-axis represents the relative expression of the gene.

Figure 4

(A–C) Trajectory analysis of M1 macrophages and M2/M1 macrophages. (D) The trajectory analysis of the heatmap of M1 macrophages versus M2/M1 macrophages. The X-axis represents the timeline of trajectory analysis, the left Y-axis represents the KEGG enrichment results, and the right Y-axis represents the differentially expressed genes between the two clusters. (E–G) Pseudotime analysis of genes (MRC1, IL1B, and MARCO). The X-axis represents the cell of trajectory analysis, and the Y-axis represents the relative expression of the gene.

Figure 5

(A–C) Trajectory analysis of M2 macrophages and M2/M1 macrophages. (D) The trajectory analysis of the heatmap of M2 macrophages versus M2/M1 macrophages. The X-axis represents the timeline of trajectory analysis, the left Y-axis represents the KEGG enrichment results, and the right Y-axis represents the differentially expressed genes between the two clusters. (E, F) Pseudotime analysis of genes (MRC1, IL1B, and MARCO). The X-axis represents the cell of trajectory analysis, and the Y-axis represents the relative expression of the gene.

Figure 6

(A) The leftmost color block represents the module, and the rightmost color bar represents the correlation range. In the heatmap in the middle part, the darker the color is, the higher the correlation. Red indicates a positive correlation, and blue indicates a negative correlation. The numbers in each cell indicate relevance and significance. The X-axis represents the sample type. (B) A scatterplot of gene significance (GS) for treat vs. module membership in the turquoise module. There is a highly significant correlation between GS and MM in the module. (C) The left circle represents the disease-characteristic genes screened using the WGCNA method, and the right circle represents the characteristic genes that change most clearly over time between macrophage subtypes. The intersection of the two circles represents the intersecting genes. (D, E) Boxplot of the residual and reserve cumulative distribution of the residual. (F) The ROC curve shows the difference between SVM and RF. (G, H) RF analysis results and screening for important genes.

Figure 7

(A) A nomogram was created to represent the disease model. It uses the X-axis to display the expression of a single gene, as well as the score scale of a single gene, the total score scale of all genes, and the disease incidence scale. Meanwhile, the Y-axis shows individual genes, points, total points, and risk of disease. (B) A graph was used to plot the predicted event rate (Predicted Probability) on the abscissa and the observed actual event rate (Actual Rate) on the ordinate, ranging from 0 to 1. This can be interpreted as the event rate in percentage. The diagonal dashed line serves as the reference line, representing the scenario where the predicted value equals the actual value. (C) The DCA graph employs the threshold probability (ThresholdProbability) on the abscissa and the net profit rate after subtracting the disadvantages on the vertical axis. (D) A graph was used to represent the high-risk threshold and benefit rate on the abscissa, and the number of high risks on the ordinate.

Figure 8

(A) Heatmap of intersecting genes in the training group. (B) Heatmap of intersecting genes in the test group. (C–H) ROC curves of disease signature genes in the test group. (I, J) The area under the curve of the AS prediction model in the train group and test group.