Gene expression profiling datasets were obtained from the NCBI Gene Expression Omnibus public database (GEO) (https://www.ncbi.nlm.nih.gov/geo/). Screening was performed in accordance with the following criteria: 1) Tissues originate from skin biopsy on Homo sapiens; 2) At least 10 samples were included; 3) Samples have been treated with no modifying drugs. Finally, the GEO dataset numbered GSE45485, GSE76885, GSE32413, GSE95065, GSE58095, and GSE125362 were selected.

Individual datasets underwent stringent quality control, background correction, log2 transformation, and normalization in the environment of the R software (version 4.2.1). Agilent microarrays (GSE45485, GSE76885, GSE125362, and GSE32413) were normalized using the “limma” package; Illumina microarray (GSE58095) was normalized using the “lumi” package; Affymetrix microarray (GSE95065) was subjected to RMA normalization using the “affy” package, respectively. GSE45485 and GSE76885 were merged, and the batch effects were corrected with the ComBat function of the “sva” package in R. A total of 135 samples (38 HCs and 97 SSc patients) of the merged GSE45485 and GSE76885 were utilized to conduct the WGCNA analysis. And GSE32413, GSE95065, GSE58095, and GSE125362 were utilized for the validation, respectively.

Differential expression analysis of HC and SSc samples was performed using the “limma” package. With |log2 fold change (FC)| > 0.585 and adjusted p < 0.05 as the cutoff threshold, differentially expression genes (DEGs) were detected. To better visualize the results, heatmap and volcano plot of DEGs were generated using the “pheatmap” and “ggplot2” packages.

To investigate the co-expression relationships among the genes and the relationship between the genes and the phenotypes, weighted correlation network analysis (WGCNA) method was applied using the “WGCNA” package in R. After filtering the outlier samples, with an optimum soft threshold was set, the weighted adjacency matrix was transformed into a topological overlap matrix (TOM) to estimate the network connectivity. Then, the co-expression modules were clustered by a dynamic tree-cut approach on TOM-based dissimilarity. Genes with similar patterns were grouped into a module. At last, correlation coefficient analysis of module membership (MM) with gene significance (GS) was implemented.

After identifying the key module that most representing the SSc disease trait, the intra-module connectivity (IMConn) was then calculated to determine the top 30 genes with the highest connectivity within the key module. Besides, the criteria (absolute values of GS > 0.45 and MM > 0.80) was used to screen the genes with biological importance in the key module. The intersection of DEGs, top 30 genes with the highest IMConn, and genes with biological importance in the key module, was taken using the tool on Evenn website (http://www.ehbio.com/test/venn/). The common genes were defined as the final hub genes of SSc.

The functional annotation of DEGs analyzed by Gene Ontology (GO) was reflected in biological processes (BP), cell components (CC), molecular function (MF). Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was performed to analyze related significant pathways. The GO and KEGG pathways were retrieved with a cut-off criterion of p < 0.05 and visualized by the “enrichplot” and “GOplot” packages in R.

The respective functions of each hub gene were revealed by Gene Set Enrichment Analysis (GSEA). After removing HC samples, the left SSc samples were distinguished into two groups, the low expression group and the high expression group, based on median values of hub gene expression levels. Differential expression analysis between the two groups was performed and genes were sorted by logFC from the highest to the lowest. The ridge plots were presented using the “clusterprofiler” and “enrichplot” package in R.

To identify the diagnostic markers among hub genes, LASSO regression was conducted using the “lars” package in R. Next, the sensitivity and specificity of the diagnostic model were evaluated using receiver operating characteristic curves (ROCs) using the “pROC” package in R. The diagnostic performance of the model was assessed by the area under the curve (AUC), and AUC > 0.75 was set as the cut-off value. In general, the gene expression profile of the merged dataset (GSE45485 and GSE76885, n = 135) was used as discovery cohort, and the gene expression profiles of GSE32413 (n = 35), GSE95065 (n = 33), GSE58095 (n = 102), and GSE125362 (n = 12) were used as validation cohorts to verify the ability of diagnostic model. The correlations between diagnostic markers and well-known causative genes (TGFB1, CTGF, COL1A1, COL1A2, IL6, CCL2, VCAM1, and THBS1) in the merged dataset were assessed by Pearson’s correlation test using the “ggplot2” package in R.

Connectivity Map (CMap) is a collection of databases that stores a pool of gene transcription-expression profiles from cultured mammalian cells exposed to active small molecule drugs. Top 30 genes with the highest IMConn were uploaded to the L1000 platform (https://clue.io/) for prediction of potential drugs towards SSc for pharmaceutical development. Compounds with the CMap negative connectivity score of −90 or lower, indicating higher potential anti-SSc effect, were considered to be potential effective drugs. Meanwhile, SwissTargetPrediction online database (http://www.swisstargetprediction.ch/) was utilized to predict the targets of the predicted drugs.

The molecular structures of predicted compounds were obtained from PubChem Compound (https://pubchem.ncbi.nlm.nih.gov/). The 3D coordinates of predicted targets were downloaded from the PDB (http://www.rcsb.org/pdb/home/home.do). AutoDockTools (version 1.5.7), AutoDock Vina (version 1.1.2), PyMOL (version 2.5) and ChemDraw (version 19.0) softwares were utilized for molecular docking studies and model visualization, respectively. The pkCSM (http://structure.bioc.cam.ac.uk/pkcsm) server was utilized to measure ADMET (absorption, distribution, metabolism, excretion, and toxicity) properties of predicted compounds (Pires et al., 2015). Using the canonical smiles strings retrieved from the PubChem database, pharmacokinetic and toxicity properties were calculated.

All statistical analyses were performed using R software. p < 0.05 was considered statistically significant.

The flow chart of the study was summarized in Figure 1. Details of the collected datasets are presented in Table 1. The DEGs were investigated in HC and SSc in the merged microarray dataset (GSE45485 and GSE76885). A total of 86 genes were identified to be differentially expressed between HC and SSc samples, of which 58 genes were upregulated and 28 genes were downregulated. The volcano plot and heatmap of DEGs in each group were presented in Figures 2A, B.

To gain insights into the biological roles of the DEGs, we performed GO categories enrichment analysis. With the criterion of p < 0.05, “extracellular matrix organization”, “extracellular structure organization”, “bone development”, “granulocyte migration”, and “myeloid leukocyte migration” exhibited highly significant enrichment within the BP category. For the CC category, DEGs were significantly enriched in “collagen-containing extracellular matrix”, “collagen timer”, “secretory granule lumen”, “cytoplasmic vesicle lumen”, and “vesicle lumen”. In addition, the MF category contained DEGs significantly enriched in “receptor ligand activity”, “signaling receptor activator activity”, “extracellular matrix structural constituent”, “cytokine activity” and “glycosaminoglycan binding” (Figure 2C). The top enriched KEGG pathways included “cytokine–cytokine receptor interaction”, “viral protein interaction with cytokine and cytokine receptor”, “protein digestion and absorption”, and “biosynthesis of unsaturated fatty acids” (Figure 2D).

WGCNA algorithm was used to identify the co-expressed genes and modules in the merged gene expression datasets of GSE45485 and GSE76885. To construct the scale-free clustering dendrograms, soft threshold power was picked as 5 (Figure 3A). The signed R2 was shown in a log-log linear model for module connectivity analysis is R2 = 0.90, suggesting the successful construction of scale-free correlation (Figure 3B). After merging similar modules, eight modules with high adjacency from the co-expression network were visualized in the dendrograms (Figure 3C). According to the module-trait relationships, the MEblue module was identified as the key module that most associated with SSc disease trait (Cor = 0.64, p < 0.01) (Figures 3D, E). The dendrogram and adjacency heatmap of eigengenes further indicated that the MEblue module was the closest one to reflect the pathogenesis of SSc (Figure 3F). Thus, MEblue module was selected for downstream analysis.

A total of 901 genes were included in the MEblue model. On one hand, with GS > 0.45 and MM > 0.80, 45 genes (top 5% of all genes in the key module) with biological importance in the MEblue module were filtered out. On the other hand, on the basis of the expression values of IMConn, the top 30 highly connected genes in the key module to mine potential key molecules were selected. Venn plot showed the intersection of DEGs, top 30 genes with the highest IMConn, and 45 genes with biological importance (Figure 4A). Thus, 7 overlapping genes (THY1, SULF1, PRSS23, COL5A2, NNMT, SLCO2B1, and TIMP1) were ultimately selected and defined as the final hub genes that might be involved in the pathogenesis of SSc. Detailed information of hub genes was presented in Table 2. Correlation analysis among 7 hub genes revealed that they were highly connected in the expression levels with each other (Figure 4B).

GSEA results revealed the potential biological roles of hub genes (Figure 4C). The ridge plots have shown that they are associated with autoimmune diseases (systemic lupus erythematosus, rheumatoid arthritis, and lipid and atherosclerosis), microorganism infections (tuberculosis, staphylococcus aureus infection, pertussis, and leishmaniasis), inflammatory related pathways (NF-kappa B signaling pathway, chemokine signaling pathway, TNF signaling pathway, AGE-RAGE signaling pathway in diabetic complications), immune responses (phagosome, complement and coagulation cascades, natural killer cell mediated cytotoxicity), and fibrosis process (protein digestion and absorption, ECM-receptor interaction, focal adhesion, osteoclast differentiation).

In the discovery cohort of GSE45485 and GSE76885, LASSO regression algorithm was employed to further screen prognostic SSc-related signature genes. According to the minimum partial likelihood deviance and optimum λ value, THY1 and SULF1 were identified as prognostic markers (Figure 5A). The violin plots revealed that THY1 and SULF1 exhibited higher expression levels in SSc patients than HCs (p < 0.05) (Figure 5B). The AUC values of THY1 and SULF1 in the merged dataset were 0.903, 0.903, respectively, suggesting high diagnostic efficacy, and the combined diagnostic value of THY1 and SULF1was 0.922 (Figure 5B). Then, the relationships between the diagnostic markers and causative genes (TGFB1, CTGF, COL1A1, COL1A2, IL6, CCL2, VCAM1, and THBS1) was verified by Pearson correlation analysis. The expression levels of THY1 and SULF1 were positively associated with these causative genes levels with high relevance (p < 0.05) (Figure 5C). To further verify the diagnostic markers, expression level detection and ROC analysis were conducted in the validation datasets. The expression levels of THY1 and SULF1 were higher in SSc patients than HCs, based on the violin plots of GSE32413, GSE95065, GSE58095 and GSE125362 (p < 0.05) (Figure 5D). Similarly, the AUC values of THY1 and SULF1 were 0.787, 0.963 in GSE32413, respectively; 0.970, 0.996 in GSE95065, respectively; 0.883, 0756 in GSE58095, respectively; and 1.000, 1.000 in GSE125362, respectively (Figure 5D).

Top 30 genes with the highest IMConn were then incorporated into the CMap database for further analysis to screen target drugs related to SSc. According to CMap score, top 10 small-molecule drugs with potential therapeutic effects were ranked as follows: desoxypeganine, clofazimine, BRL-50481, GW-311616, methyllycaconitine, acetyl-geranyl-cysteine, SB-525334, dipyridamole, tomelukast, and warfarin (Table 3). The 3D structure diagrams of these candidate molecule drugs are shown in Supplementary Figure S1. The targets of the predicted drugs analyzed by SwissTargetPrediction online database were also summarized in Table 3. Molecular docking analysis was based on the structure of proteins and structure of drugs. It is believed that the molecular docking binding energy is less than 0, indicating that the ligand and the receptor can spontaneously bind. The results showed that BRL-50481 bounds to PDE4B, PDE4D, and PDE7A with the low binding energy of −4.5, −5.8, and −4.8 kcal/mol, respectively (Figure 6A). Similarly, dipyridamole bounds to PDE2A, PDE5A, PDE10A, SLC29A1 with the low binding energy of −5.5, −4.5, −5.4 and −5.2 kcal/mol, respectively (Figure 6B). SB-525334 bonds to TGFBR1 with the low binding energy of −6.4 kcal/mol (Figure 6C). And, GW-311616 bonds to ELANE with the low binding energy of −5.7 kcal/mol (Figure 6D). These data all indicated the highly stable binding between drugs and proteins.

On the other hand, ADMET properties were predicted to further understand how these compounds behave in a biological system, as shown in Table 4. According to the theory provided by the pkCSM server, the compounds are all well absorbed (Intestinal absorption >30%). Dipyridamole is poorly distributed to the brain (log BB < −1) and unable to penetrate central nervous system (log PS < −3), while SB-525334 can readily cross the blood–brain barrier (log BB > 0.3) and penetrate central nervous system (log PS > −2). Mostly, dipyridamole, GW-311616, and BRL-50481 are not inhibitors of cytochrome P, except that BRL-50481 is an inhibitor of CYP1A2. In terms of toxicity, all compounds are predicted to have no skin sensitization, and no obvious cardiotoxicity. However, most of them may have potential AMES toxicity or hepatotoxicity.