Juvenile dermatomyositis (JDM) is a rare chronic childhood-onset autoimmune disease characterized by inflammatory infiltration in small vessels and tissues within skin and muscle. The incidence of JDM is 2–4 per million per year in the United States (Feldman et al., 2008), with female:male ratios ranging from 1.5:1 to 5:1 (Lindsley et al., 1995). The major manifestations of JDM patients consist of symmetrical proximal muscle weakness, skin rashes, and internal organs involvement (Crowe et al., 1982). Up to 30% of JDM may present with calcifications, one of the prognostic factors of long-term disability (Arabshahi et al., 2012; Li & Zhou, 2019; Ravelli et al., 2010). Adults with JDM in childhood are susceptible to premature cardiovascular damage (Gitiaux et al., 2016).
Pathological state and treatment have been reported to affect growth and puberty in the active phase of JDM (Nordal et al., 2019). Ongoing disease activity, irreversible damage, and aggressive immunosuppressive therapy remain major challenges for long-term outcomes and quality of life in JDM patients (Hoeltzel et al., 2014). The etiology of JDM remains ill-defined although genetic and environmental factors are suspected to be involved in its pathogenesis. It has been reported that JDM patients had higher incidence of Epstein-Barr virus infection (Zheng et al., 2019), and the prominent type 1 interferon (IFN) signature was shown to affect the vasculature JDM (De Paepe, 2017; Greenberg, 2010). Adaptive and innate immune mechanisms involving IFN-associated molecules appear to mediate endothelial tubule-reticular formations and peri-fascicular atrophy.
Weighted gene co-expression network analysis (WGCNA) algorithm is a powerful bioinformatic method that mines practical information from gene expression profiles by constructing of gene modules, thereby interpreting the biological significance of a gene (Langfelder & Horvath, 2008). WGCNA has been widely used in various diseases (Zhao et al., 2010), including malignancies, cardiovascular diseases and autoimmune diseases, where it has provided useful information for understanding pathological process and for discovery of diagnostic and prognostic biomarkers. Nevertheless, WGCNA has never been applied to JDM.
Therefore, we used WGCNA for the first time to analyze pathological state and gene expression data in JDM muscular samples to explore and validate hub genes associated with JDM, as well as to predict small-molecule compounds to treat JDM with promising perspectives.
Materials & Methods
Data collection and differentially expressed genes screening
The flowchart of the study is shown in Fig. S1. Microarray profiles of JDM were retrieved from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) of the National Center for Biotechnology Information using the search terms of “juvenile dermatomyositis” restricted in the title. The datasets enrolled in this study must contain musclular specimens with three biological replicates at least. The “affy” package in R environment (version 3.6.1) was used to quantile normalize the expression within each dataset (Sasik, Calvo & Corbeil, 2002). The corresponding platforms were applied to annotate each probe according to Entrez ID, and the average expression value was calculated if several probes corresponded to the same Entrez ID (Table S1). The “limma” R package was performed for identifying differentially expressed genes (DEGs) between JDM samples and normal samples under cut-off criteria of false discovery rate (FDR) <0.05 and —log2fold change— ≥ 1.
Co-expression network construction
The variance of each gene expression value was calculated and the genes with variance ranked in the top 25% were selected for the construction of WGCNA (Langfelder & Horvath, 2008). The “WGNCA” package was used to construct the co-expression network. In detail, the function goodSamplesgenes was used to include the qualified genes and samples, followed by choosing an appropriate soft-thresholding power to construct the weighted adjacency matrix by the function pickSoftThreshold. The adjacency matrix was transformed into the topological matrix (TOM), and TOM-based dissimilarity (1-TOM) measure was used to cluster the genes using the flashClust function. Genes in the same module were highly interconnected. Then, phenotype (clinic traits) was imput into the co-expression network, and the following parameters were calculated: module eigengene (ME), gene significance (GS), and module membership (MM). ME represents the significant component in the principal component analysis for each gene module, and MM refers to the connectivity between genes and modules. GS was representative of correlation strength between gene expression and clinical traits, which was calculated by log10 transformation of the P-value (GS = lg P) in the linear regression. Key modules were considered based on the criteria that the correlation coefficient ≥ 0.80 and P-value <0.05.
Functional enrichment analysis
All genes in key modules were uploaded to the g:Profiler online (Reimand et al., 2007) database to perform Gene Ontology (GO) functional annotation (Ashburner et al., 2000) and the Kyoto encyclopedia of genes and genomes (KEGG) enrichment pathway analysis (Kanehisa & Goto, 2000). GO functional analysis consists of biological process (BP), cellular component (CC), and molecular function (MF). Analysis results were extracted under the condition of adjusted P- value <0.05. The top five terms were visualized if there were more than five terms.
Selection and validation of hub genes
Genes with high correlation in candidate modules were defined as candidate hub genes. High connectivity was considered when the connectivity ranked in the top 2%. Candidate hub gene met the absolute values of MM >0.80 and GS >0.20. After identifying hub genes highly associated with clinical traits, the search tool for the retrieval of interacting genes (STRING) database was used to construct a protein-protein interaction (PPI) network for the candidate hub genes, and molecular complex detection (MCODE, a plugin in Cytoscape) was used to further select the real hub (Shannon et al., 2003; Szklarczyk et al., 2015). Genes with MCODE score ≥ 0 in the PPI network were selected as the final hub genes. A separate dataset (GSE11971) was used to validate the differential expression of the final hub genes.
Related small-molecule compounds screening
Connectivity map (CMap) database (https://portals.broadinstitute.org/cmap) was used to screen out small molecule compounds based on the real hub genes associated with JDM, because most compounds in this database are the United States Food and Drug Administration-approved drugs (Lamb et al., 2006). First, real hub genes were divided into upregulated and downregulated groups. Next, these probe sets were used to query the CMap database based on the platform of the Affymetrix Human Genome U133 Plus 2.0 Array (http://www.affymetrix.com/analysis/netaffx/index.affx). Finally, enrichment scores representing similarity were calculated, ranging from −1 to 1. Small molecules generated from up-regulated genes suggested therapeutic goals, while down-regulated genes predicted inhibitors of therapy for the disease. Potential compounds were selected based on connectivity score, P-value and correlation.
Two-tailed Student’s t-test was applied to the significance of differences between groups, and P-value less than 0.05 was considered as statistically significant. Statistical analyses were performed using Graphpad Prism 8.0.
Data collection and differentially expressed genes
We employed two datasets on JDM muscular expression profiles. Dataset GSE3307 was used as the training set (Bakay et al., 2006). The original study enrolled 39 muscular biopsy samples, including 21 JDM patients and 18 healthy controls (HC). Dataset GSE11971, including nineteen JDM patients and four normal controls, was used as the validating set (Chen et al., 2008). The gene expression profiles of all tissue samples were analyzed based on the platform of the Affymetrix Human Genome U133 Plus 2.0 Array. A total of 2,834 differentially expressed genes between JDM and HC were identified, including 1,888 down-regulated genes and 946 up-regulated genes. The DEGs are listed in Table S2.
Construction of a weighted co-expression network and identification of key modules
5103 genes whose variance ranked in the top 25% with 21 JDM samples and 18 control samples in GSE3307 were used for WGCNA construction. The “WGCNA” R package was used for expression matrix of GSE3307, and soft-thresholding power β value equal to 10 was selected to ensure a scale-free network with scale-free R2 equal to 0.90 (Figs. S2A–S2B) (Langfelder & Horvath, 2008). A total of 13 modules were returned by WGCNA analysis (Figs. 1A–1B).
The interaction relationship of 12 modules was analyzed using network heatmap plots (Fig. 1C). The division of all modules was highly independent from our analysis. The module eigengene dendrogram showed that 12 modules were divided into two clusters, and the adjacency heatmap of eigengene showed a similar result (Fig. 2A). Based on the criteria that correlation coefficient ≥ 0.80, P value <0.05, blue, lightgreen and midnightblue modules were identified as key modules for further analysis (Fig. 2B). Therefore, we selected the blue, lightgreen and midnightblue modules for subsequent analysis, to identify the relevance between key modules and the pathological state of JDM with substantial biological significance (Figs. 2C–2E).
Functional and pathway enrichment analysis
GO and KEGG pathway enrichment was performed for all genes in the key modules to mine the biological functions associated with JDM. Biological process of GO analysis showed genes in the blue module were associated with generation of SRP-dependent cotranslational protein targeting to membrane, cotranslational protein targeting to membrane, protein targeting to ER, nuclear-transcribed mRNA catabolic process and establishment of protein localization to endoplasmic reticulum; and that in the lightgreen module was relevant to response to type I interferon, type I interferon signaling pathway, cellular response to type I interferon, defense response to virus and response to virus.The top five pathways related to the midnightblue module were cellular response to chemical stimulus, extracellular structure organization, extracellular matrix organization, response to organic substance and cell motility (Fig. 3A). Pathway enrichment results of MF and CC in three key modules are presented in Figs. 3B–3C. The results of the KEGG pathway enrichment analysis in three modules are shown in Fig. 3D.
Identification of hub genes
Based on the criteria that MM >0.80 and GS >0.20, a total of 45 DEGs with the high connectivity in key modules were screened as candidate hub genes. Then, a PPI network was constructed for candidate hub genes using Cytoscape, consisting of 42 nodes and 80 edges according to STRING database (Fig. 4). We conducted molecular complex detection (MCODE) (a plugin in Cytoscape) analysis for 45 candidate hub genes, and 28 genes (blue = 15, lightgreen = 11, midnightblue = 1) were considered hub genes according to the criteria of MCODE score ≥ 0. Table 1 shows 28 hub genes in the three modules.
All hub genes were validated using JDM data from another GEO database (GSE11971). Because of the differences in microarray probes used in two data sets, boxplots were used to show the validation results for the final 22 hub genes (Fig. S3). We found that seven genes, SP110, SAMHD1, IFIT5, PLSCR1, IFI16, MX2 and CLIC1, were significantly upregulated in JDM compared to HC, while thirteen genes, COX5B, COX6A2, COX7C, NDUFA4, NDUFB4, MDH2, ATP5O, ATP5B, RPL21, TPI1, SLC25A3, VDAC1 and EIF4B, were significantly downregulated in JDM in comparison of HC. Figure 5 summarizes the cross-talk pathways involved in the pathogenesis of JDM by hub genes and literature (Miller et al., 2018; Thompson, Piguet & Choy, 2017).
Related small-molecule compounds screening
The CMap database was used for small molecule drugs screening based on 20 real hub genes associated with JDM. Based on the criterion that the number of instances exceeds five and P-value less than 0.05, twelve small-molecule compounds were identified (Table 2). Among these compounds, acacetin, helveticoside, lanatoside C, deferoxamine, famprofazone, tanespimycin and LY-294002 may perturb the development of JDM, while betonicine, felodipine, valproic acid, and sirolimus might provide potentially therapeutic goals for JDM.
|Cmap name and cell line||Mean score||Number||Enrichment||P-value||Specificity||Percent non-null|
|trichostatin A - HL60||0.283||34||0.623||0||0.0798||52|
|trichostatin A - MCF7||0.172||92||0.399||0||0.673||56|
|LY-294002 - PC3||−0.243||12||−0.494||0.00325||0.2249||58|
|valproic acid - HL60||0.254||14||0.417||0.01034||0.1812||50|
|sirolimus - HL60||0.305||10||0.534||0.00359||0.0347||70|
|vorinostat - MCF7||0.22||7||0.5||0.0378||0.7655||71|
In this study, we used WGCNA to construct a co-expression network, detect key gene modules and identify hub genes in JDM for the first time. Our research provides some potential biomarkers or molecular targets for JDM through the Cmap database. We found that three modules highly correlated with JDM. The expression of 28 genes in these three modules showed significant changes in patients with JDM compared to control individuals in the training period, and 20 genes were validated as the real hub genes in the GSE11971 dataset, including the downregulation of NADH dehydrogenase, ATP synthase and cytochrome c oxidase and upregulation of IFN-stimulated genes. However, few of them were identified as biomarkers or crucial genes in JDM yet.
Functional enrichment analysis indicated that type I interferon signaling and various virus infection pathways were strengthened in JDM compared to HC, which is consistent with findings of previous studies (Moneta, Marafon & Marasco, 2019; Piper et al., 2018). IFIT5, IFI16 and MX2, interferon-stimulated genes, both nuclear transcriptional factors, were found to be upregulated in other autoimmune diseases but not in JDM (Wang et al., 2019b; Zhang & Xu, 2019). In the present study, we found that interferon-stimulated genes were significantly upregulated in JDM patients, as well as the so-called “interferon signature”, demonstrating a possible mechanism that viral mimics or other stimuli may play a crucial role in the pathogenesis of JDM. Viral mimics are thought to participate in the pathogenesis of JDM (Musumeci & Castrogiovanni, 2018) and other autoimmune diseases (Christen et al., 2004; Sellami et al., 2019), consistent with the notion that JDM patients have higher rates of viral infections (Tansley, McHugh & Wedderburn, 2013; Zheng et al., 2019). This may suggest that the prevention of certain viral infections would decrease the incidence of autoimmunity by inhibiting self-antigenic mimics. NADH dehydrogenase (NDUFA4 and NDUFB4), ATP synthase (ATP5O and ATP5B) and cytochrome c oxidase (COX) family (COX5B, COX6A2 and COX7C) are crucial molecules involved in the oxidative phosphorylation in mitochondrial metabolism, and the decreased levels of these molecules suggested a crucial role of impaired mitochondrial phosphorylation and lower oxidative capacity in the pathogenesis of JDM, accounting for the extremity weakness in JDM patients.
Hypoxia caused by suppressed oxidative phosphorylation induces changes in reactive oxygen species (ROS) generation, whereby severe hypoxia in skeletal muscle results in elevated H2O2 generation. ROS accumulation produced by mitochondrial dysfunctions, in turn, drives type I interferon responses and muscle inflammation, and may thereby self-sustain the disease process (Wang et al., 2019a). Similar to other autoimmune diseases, high-dose glucocorticoids, used alone or in combination with immunosuppressive agents are routine treatment for JDM patients wheras some refractory patients may develop functional limitations. It has been suggested that refractory JDM patients, have lower maximal oxygen uptake (Drinkard et al., 2003; Hicks et al., 2002) than do healthy children and with children with juvenile dermatomyositis in remission (Takken et al., 2008), suggesting that mitochrondrial dysfuction may contribute to the severity of JDM. Current concepts on the therapy of muscle weakness in JDM focus on induction of partial recovery and exposure to serious adverse events (including muscular toxicity). Our data suggest a novel therapeutic perspective for JDM by protecting mitochondria from dysfunction.
Bioinformatics combined human and material resources to develop more efficient tools with lower error rates (Irizarry et al., 2003). WGCNA is an efficient approach to construct co-expressed modules and hub genes in several diseases. Previous studies using microarray expression profiles from adult-onset DM patients showed that IFN-stimulated genes were upregulated (i.e., MX2, STAT1 and OAS3), suggesting that the IFN signature overlapped the pathogenesis both in adult and juvenile DM. Nevertheless, mechanisms linked to hypoxia are less prevalent in adult-onset DM, suggesting mitochondrial dysfunctions contribute more to juvenile-onset DM rather than adult-onset DM.
We used the CMap database to predict several kinds of small-molecule compounds with promising capacity as therapeutic goals or inhibitors on treatment for JDM. No evidence has demonstrated the direct association between these compounds and JDM, while they hinted indirect link to JDM, according to the literatue. Among these compounds, acacetin, helveticoside, lanatoside C, deferoxamine, famprofazone, tanespimycin and LY-294002 showed negative enrichment scores and thus may have the potential to perturb the development of JDM, while betonicine, felodipine, valproic acid, and sirolimus showed positive enrichment scores and might provide potentially therapeutic goals for JDM. Acacetin, an inhibitor of lipopolysaccharide-induced inflammation, can promote the expansion of Treg cells and supress the differentiation of Th17 cells in a dose-dependent manner in collagen-induced arthritis (Liu et al., 2018). Helveticoside can regulate metabolism and signaling processes as a biologically active component, but little is known in inflammatory reactions (Kim, Lee & Kim, 2015). The iron chelator deferoxamine was shown to reduce mitochondrial oxidative stress in a transient cerebral ischemia model as well as the release of pro-inflammatory molecules including matrix metalloproteinase-9 and hypoxia inducible factor-1 (Im et al., 2012). LY294002, a kind of PI3K inhibitor, has potential against experimental autoimmune myocarditis (Liu et al., 2016). The heat-shock protein 90 inhibitor tanespimycin has been shown to inhibit cutaneous inflammation in experimental epidermolysis bullosa acquisita (Tukaj et al., 2017) and other experimental autoimmune models (Dello Russo et al., 2006). Felodipine, commonly used to treat hypertension and angin, has been evidenced to inhibit oxidative stress and inflammation in endothelial cells, which is consistent with our results (Qi et al., 2017). Valproic acid is a histone deacetylase inhibitor (HDACI), can suppress the inflammatory responses mediated by cytokines, oxidative stress molecules (ROS, NO), activating receptors (NK, T γ δ, and cytotoxic lymphocytes), perforin, granzyme, costimulatory molecules, and autoantibodies (Soria-Castro et al., 2019). Sirolimus can restore immune balance in rheumatoid arthritis patients by expanding the pool of circulating Treg cells (Niu et al., 2019). Our results based on the CMap database might provide hints as to future therapy for JDM; nevertheless, studies in vitro and in vivo are necessary.
This study has some limitations. First, this is retrospective, with all data in this study being retrieved from a public database. A multicenter, prospective study is needed to evaluate the significance of these hub genes in terms of long-term outcomes and possible applications of molecular drugs for therapy. Second, experiments in vivo and in vitro are necessary to interpret potential mechanisms of real hub genes and small-molecule compounds for future clinical translation. Third, clinical traits cannot correlate with gene modules when performing WGCNA because of lack of clinical trait data in these GEO datasets.
Based on weighted gene co-expression analysis, three key modules and 20 real key genes associated with the pathological state of JDM were identified, suggesting pivotal roles of mitochondrial dysfunction and the interferon signature in JDM. This analysis provides several candidate small-molecule compounds for use as targeted therapy of JDM.
Clustering dendrogram of all samples and determination of soft-thresholding power
(A) Sample clustering was conducted to detect outliers. All samples are located in the clusters and pass the cutoff thresholds. (B) Analysis of the scale-free fit index for various soft-thresholding powers (β). (C) Analysis of the mean connectivity for various soft-thresholding powers.
Boxplots of 22 hub genes
(A) ATP5O, ATP5B, NDUFV2, COX7C, COX5B, NDUFB4 and COX6A2; (B) NDUFA4, VDAC1, SLC25A3, MDH2, TPI1, RPL21, MRPS7 and MX2; (C) SAMHD1, EIF4B, IFIT5, SP110, IFI16, PLSCR1, and CLIC1.