Revealing the roles of glycosphingolipid metabolism pathway in the development of keloid: a conjoint analysis of single-cell and machine learning

Abstract

Keloid is a pathological scar formed by abnormal wound healing, characterized by the persistence of local inflammation and excessive collagen deposition, where the intensity of inflammation is positively correlated with the size of the scar formation. The pathophysiological mechanisms underlying keloid formation are unclear, and keloid remains a therapeutic challenge in clinical practice. This study is the first to investigate the role of glycosphingolipid (GSL) metabolism pathway in the development of keloid. Single cell sequencing and microarray data were applied to systematically analyze and screen the glycosphingolipid metabolism related genes using differential gene analysis and machine learning algorithms (random forest and support vector machine), and a set of genes, including ARSA,GBA2,SUMF2,GLTP,GALC and HEXB, were finally identified, for which keloid diagnostic model was constructed and immune infiltration profiles were analyzed, demonstrating that this set of genes could serve as a new therapeutic target for keloid. Further unsupervised clustering was performed by using expression profiles of glycosphingolipid metabolism genes to discover keloid subgroups, immune cells, inflammatory factor differences and the main pathways of enrichment between different subgroups were calculated. The single-cell resolution transcriptome landscape concentrated on fibroblasts. By calculating the activity of the GSL metabolism pathway for each fibroblast, we investigated the activity changes of GSL metabolism pathway in fibroblasts using pseudotime trajectory analysis and found that the increased activity of the GSL metabolism pathway was associated with fibroblast differentiation. Subsequent analysis of the cellular communication network revealed the existence of a fibroblast-centered communication regulatory network in keloids and that the activity of the GSL metabolism pathway in fibroblasts has an impact on cellular communication. This contributes to the further understanding of the pathogenesis of keloids. Overall, we provide new insights into the pathophysiological mechanisms of keloids, and our results may provide new ideas for the diagnosis and treatment of keloids.

Keywords: glycosphingolipid; immune; keloid; machine learning; single cell.

Figure 1

The study workflow.

Figure 2

Combining different datasets. (A, B) Boxplots of mRNA expression distribution before and after removing batch effects. (C, D) PCA plots before and after removing batch effects.

Figure 3

Diagnostic model of keloid was constructed and evaluated. (A) Heatmap for differential analysis of GSLERGs between keloid and normal samples. (B) Random forest algorithm screening for gene importance ranking. (C) SVM-REF algorithm screening for genes. (D) Venn plot showed the intersection genes of the top 10 of RF, SVM-REF and DEGs. (E) ROC curves under AUC values in the diagnostic model built using the 6 GSLERG. (F) Bootstrap resampling algorithm to validate the model. (G) Keloid prediction by nomogram. (H) Calibration curve to evaluate the nomogram. (I) Decision curves to assess the predictive performance of the model.

Figure 4

Correlation between candidate genes and immune cell infiltration. (A) Co-expression patterns of 6 GSMRG in all samples based on the Pearson correlation analysis. (B, C) Heatmap of correlation between 6 GSMRG and immune cell infiltration in Cibersort and EPIC algorithm. (D) Heatmap of correlation between GSMRG and inflammatory factors. (E–F) Top 5 correlation plots in Cibersort and EPIC algorithm. *p < 0.05, **p < 0.01, ***p < 0.001.

Figure 5

Unsupervised clustering analysis in keloid. (A) The empirical cumulative distribution function (CDF) plots revealed the consensus distributions for each k. (B) The area change under CDF curve when k=2-10. (C) The circular manhattan (CM) plot exhibited the clusters at k = 2. (D) The bar plot showed the score of each subtype for the number of clusters k from 2 to 10. (E) Heat map showing the distribution of GSLMRG expression between different clusters. (F–H) Box plots of the distribution of immune cell infiltration and inflammatory factors expression in different clusters.

Figure 6

Functional enrichment analysis. (A, C) GO enrichment analysis of DEGs between fibroblasts of different subtypes. (B, D) KEGG enrichment analysis of DEGs between fibroblasts of different subtypes. (E) GSEA enrichment analysis of biological functions between two clusters. (F, G) Functional categories of differentially expressed genes between fibroblasts of different subtypes.

Figure 7

Cell populations and marker genes in keloid and normal skin. (A) After standard quality control of all cells from three keloids and three normal tissues, 43,910 cells were included in the analysis. (B) The number of genes detected was significantly correlated with the sequencing depth, with a Pearson correlation coefficient of 0.91; the same number of mitochondria was detected at different sequencing depths. (C) The cell clusters visualized by the dimensional reduction of t-distributed stochastic neighbor embedding (t-SNE). (D) Heatmap showing the top 5 genes per cell cluster after differential analysis to obtain marker genes. (E) Dot plot showed annotation of cell clusters by known markers. (F) tSNE plot presented cell type annotation for each cluster. (G) Proportions of distinct cell types for different samples. (H) Heatmap showed the top 5 marker genes between cell types.

Figure 8

Progression of fibroblast cell profiles revealed by pseudotime analysis. (A) GSL metabolic pathway scores of fibroblasts in keloid and NS. (B) tSNE plots showed the expression of marker genes in fibroblasts. (C) Fibroblasts were clustered again by downscaling and shown by tSNE plots. (D) tSNE plots of the GSL metabolic pathway scores of individual fibroblasts. (E) The GSL metabolic pathway scores of the fibroblast subpopulations that were downscaled again. (F) Heatmap showing the expression changes of genes in GSL metabolic pathway with fibroblast differentiation. (G–J) Trajectory differentiation maps according to cell differentiation status, cell development time coloring, cell cluster and tissue type. ****p < 0.0001.

Figure 9

(A) Integrated cell-cell communication networks drawn by number and weight of interactions. (B) The heatmap of outgoing/incoming interaction strength for 10 cell types. (C) The dot plot of outgoing and incoming interaction signal pathways for fibroblasts of two subtypes.