Unique and shared transcriptomic signatures underlying localized scleroderma pathogenesis identified using interpretable machine learning

Abstract

Using transcriptomic profiling at single-cell resolution, we investigated cell-intrinsic and cell-extrinsic signatures associated with pathogenesis and inflammation-driven fibrosis in both adult and pediatric patients with localized scleroderma (LS). We performed single-cell RNA-Seq on adult and pediatric patients with LS and healthy controls. We then analyzed the single-cell RNA-Seq data using an interpretable factor analysis machine learning framework, significant latent factor interaction discovery and exploration (SLIDE), which moves beyond predictive biomarkers to infer latent factors underlying LS pathophysiology. SLIDE is a recently developed latent factor regression-based framework that comes with rigorous statistical guarantees regarding identifiability of the latent factors, corresponding inference, and FDR control. We found distinct differences in the characteristics and complexity in the molecular signatures between adult and pediatric LS. SLIDE identified cell type-specific determinants of LS associated with age and severity and revealed insights into signaling mechanisms shared between LS and systemic sclerosis (SSc), as well as differences in onset of the disease in the pediatric compared with adult population. Our analyses recapitulate known drivers of LS pathology and identify cellular signaling modules that stratify LS subtypes and define a shared signaling axis with SSc.

Keywords: Autoimmune diseases; Autoimmunity; Bioinformatics; Immunology.

Figure 1. scRNA-Seq defines cell populations in the skin of patients with LS compared with age- and sex-matched healthy controls.

(A) The 3 major steps of the pipeline are illustrated: tissue collection, dissociation, and library preparation; single-cell transcription; and data analysis. (B) The uniform manifold approximation and projection (UMAP) of the 122,809 cells derived from the skin of 27 LS samples and 17 healthy samples is displayed, with the final clustering into 15 cell type subclusters. (C) The marker genes specific to each cell cluster that led to the final cellular annotations are displayed. (D) The proportions of LS and healthy that comprise each of the 15 subclusters are illustrated, with overall higher proportions of inflammatory cells (i.e., T cells and B cells) and stromal cells (keratinocytes) in LS.

Figure 2. Cell-specific expression modules differentiate patients with LS from healthy controls.

(A) SLIDE model performance measured by area under the receiver operating characteristics curve (AUC) using k-fold cross-validation. Significance assessed using permutation testing (negative control). Asterisks indicate Wilcoxon P < 0.0001. (B) Levels of significant latent factor 1 in LS and healthy patients. Asterisks indicate Wilcoxon P < 0.001. (C) Correlation network representation of significant latent factor 1. Purple and green edges indicate positive and negative correlations, respectively. Triangles indicate genes with higher expression in LS, and squares indicate genes with higher expression in control. (D) Bubble plots showing expression of genes in significant latent factor 1. Bubble size indicates average expression for within each group. (E) Levels of significant latent factor 2 in LS and healthy patients. Asterisks indicate Wilcoxon P < 0.01. (F) Correlation network representation of significant latent factor 2. (G) Bubble plots showing expression of genes in significant latent factor 2. (H) AUCs for classifying each group from healthy controls. Box plots show the interquartile range, median (line), and minimum and maximum (whiskers).

Figure 3. Cell-specific expression modules in adult LS resemble expression in SSc.

(A) SLIDE model performance measured by AUC using k-fold cross-validation. Significance assessed using permutation testing (negative control). Asterisks indicate Wilcoxon P < 0.0001. (B) Latent factors trained to classify adult patients with LS can classify patients with SSc in cross-prediction. Cross-prediction AUCs for classifying each group from healthy controls. (C) Levels of significant latent factor 4 in adult patients with LS and healthy controls. Asterisks indicate Wilcoxon P < 0.001. (D) Correlation network representation of significant latent factor 4. Purple and green edges indicate positive and negative correlations, respectively. Triangles indicate genes with higher expression in LS, and squares indicate genes with higher expression in control. (E) Levels of significant latent factor 5 in adult patients with LS and healthy controls. (F) Correlation network representation of significant latent factor 5. Box plots show the interquartile range, median (line), and minimum and maximum (whiskers).

Figure 4. Adult and pediatric transcriptomic signatures in LS.

(A) SLIDE model performance measured by AUC using k-fold cross-validation. Significance assessed using permutation testing (negative control). Asterisks indicate Wilcoxon P < 0.0001. (B) Levels of significant latent factor 1 in adult LS and pediatric LS. (C) Correlation network representation of significant latent factor 1. Purple and green edges indicate positive and negative correlations, respectively. Triangles indicate genes with higher expression in adult LS, and squares indicate genes with higher expression in pediatric LS. (D) Levels of significant latent factor 2 in adult LS and pediatric LS. (E) Correlation network representation of significant latent factor 2. Box plots show the interquartile range, median (line), and minimum and maximum (whiskers).

Figure 5. Transcriptomic signatures associated with mLoSSI score in LS.

(A) SLIDE model performance measured by correlation with mLoSSI using k-fold cross-validation. Asterisks indicate Wilcoxon P < 0.001. (B) Latent factors trained to predict mLoSSI score in LS (top) can cross-predict MRSS in SSc (bottom). Plot shows true and predicted scores with P values measured by Spearman correlation. (C) Spearman correlation between latent factor 1 levels and mLoSSI score. (D) Correlation network representation of latent factor 1. Purple and green edges indicate positive and negative correlations, respectively. Triangles indicate genes associated with higher mLoSSI, and squares indicate genes associated with lower mLoSSI. (E) Spearman correlation between latent factor 2 levels and mLoSSI score. (F) Correlation network representation of latent factor 2. Box plots show the interquartile range, median (line), and minimum and maximum (whiskers).

Figure 6. Spatial transcriptomics images of a representative LS from the dataset.

(A) The H&E stain of the representative LS tissue shown with areas of interest blown out. Marker genes (MUCL1 for eccrine, KRT10 for keratinocytes, PECAM1 for endothelial cells, and CD163 for macrophages) were used to look for latent factors identified in the corresponding cell types. (B) Colocalization and intensity of signals from all latent factors around sebaceous glands/follicles with inflammation pointing toward possible interaction. (C) KRT14 and ACTG1 upregulated in eccrine cells and (D) KRT5 and AQP3 in keratinocytes. These were colocalized in the same areas. (E) SOX4 and KLF2 were upregulated in endothelial cells, and (F) HNRNPA2B1 was upregulated in macrophages and colocalized in the same area.