The list of diseases available under the Medical Subject Headings (MeSH) terms “neuromuscular,” “musculoskeletal,” and “muscular” diseases was used as a guideline for identifying muscular diseases of interest (Lipscomb, 2000). All available (RNAseq + microarray platform) information from GEO (Barrett et al., 2013) was downloaded and surveyed for maximum coverage of muscle diseases in the MeSH headings identified above. A single platform GPL96 (Affymetrix HG-U133A) offered the highest coverage of muscle diseases surveyed. Choosing studies from one platform alone (GPL96) limited possible noise arising from platform differences. Additionally, studies with non-muscle sample tissue, and less than two-samples/condition were eliminated. Filtering the data for accuracy, and experimental context using our constraints resulted in microarrays from 19 human diseases. In addition to these 19 diseases, samples from “CP” (not available under the MeSH term categories considered above) was also included resulting in 20 diseases for analysis. CP is a movement disorder characterized by contractures of the muscle with its primary insult on the nervous system. Gene expression data for studies with. CEL files (disease and control) are normalized using RMA (Robust Multi-array Average). Studies with series matrix files were downloaded as is. ComBat cross-array normalization is utilized for diseases with more than one associated GSE (e.g., LGMD2A, DMD, and JDM, Table 1), to remove study artifacts (Johnson et al., 2007). Multiple probes were accounted for using the “collapseRows” function of WGCNA library in R (Langfelder and Horvath, 2008). The reduced and processed data sets subsequently included z-transformed expression values of 12,789 genes across 20 diseases.

R (v 3.2.2) (R Core Team, 2015)/Bioconductor (Gentleman et al., 2004) was used for all data processing and pipelines implemented in this analysis.

Cyber-T's (Kayala and Baldi, 2012) regularized t-test performed on z-transformed expression values was used as a measure of a gene's differential activity. The associated T-statistic provided insight into the difference in mean expression of a gene across conditions and was referred to as the differential gene activity (DGA) score for each disease state (Table S0). Partial pearson correlation of DGA scores across diseases, quantified disease similarity based on expression profiles. Use of partial correlation have been shown to be effective in factoring out dependencies such as variation in tissue types and experimental conditions (Suthram et al., 2010). The “pcor” function available through the “ppcor” package in R (Kim, 2015) was utilized to calculate all possible pairwise partial correlations between each pair of samples (here diseases) while eliminating the effect of all other samples. We identified disease clusters using hierarchical (complete-linkage) clustering of partial correlation.

A comprehensive list of genetic factors affecting the 20 muscle diseases from various sources such as OMIM (Amberger et al., 2009), PheGenI (Ramos et al., 2014), ClinVar (Landrum et al., 2014), and DisGeNET (Piñero et al., 2015) (henceforth referred to as the disease-gene list) was downloaded and curated. A hypergeometric model, with a null that disease-genes were randomly drawn from the space of all genes was used to ascertain the statistical overlap between diseases based on the disease-gene list. The function “phyper,” available through base R packages was utilized to calculate the hypergeometric p-values.

We developed a framework of “functional modules” (FM) within muscle that represented significant units required for normal muscle activity (Table 2). Each of the 23 FMs were represented by a group of manually curated list of biomarkers that belonged to a broader functional pathway (family) within muscle (Table S1). This list expanded on an existing framework for muscle functional families (Mukund et al., 2014). Several biomarkers within each FM were multi-functional and were placed in FMs that were most relevant to skeletal muscle. The functional module activity (FMA) score associated with each functional module i in disease k, was calculated as mean DGA score of its component genes. This score reflected the overall state of a module in a particular disease (Table S2). Specifically, a negative FMA values reflected a downregulation of the module genes in disease with respect to controls and the converse for positive FMA values. Significance testing at p < 0.05 identified FMs associated with each disease state.

The human protein-protein interaction network (PPIN) was extracted from the STRING database (v9.1), containing both direct (physical) and indirect (functionally derived) interactions (Franceschini et al., 2013). Limiting the interactions to a combined score cutoff of >0.85 allowed us to account for strong interactions, with sufficient experimental evidence and resulted in a total of 1,48,030 unique interactions among 10,341 proteins. The PPIN was clustered using MCL (Markov Cluster) algorithm a fast, scalable and unsupervised cluster algorithm for networks based on simulation of stochastic flow. The MCL algorithm finds cluster structure in graphs by iteratively computing the probabilities of random walks through the graph (markov matrices) using an alternation of two operators called expansion and inflation. Expansion is the power of a stochastic matrix using the normal matrix product (i.e., matrix squaring) and inflation is the Hadamard power of a matrix (taking powers entrywise), followed by a scaling step. A detailed description and comparison are provided in Enright et al. (2002) and Brohee and Van Helden (2006). Clustering the PPIN using MCL clustering resulted in 1,025 protein modules (filtered for minimum module size >2) with sizes ranging between 3 and 256 genes/proteins. We considered 764/1,025 modules (6,215/8,581 proteins overlapping with our list of 12,789 genes) with at least three genes/module for further analysis (Table S3). Analogous to Suthram et al. (2010) the protein module activity (PMA) score was calculated for each protein module i in a disease k as the mean of DGA scores for its component genes. In the end, we obtained a vector of PMA_ik for each disease, representing the activity level of a given protein module in each disease state (Table S4). Significance testing identified protein modules with a p < 0.05 associated with each disease. We defined the threshold for module expression as upper 50th percentile of |mean PMA| scores from significant modules across diseases.

A background distribution of disease correlations expected at random was generated to assess the significance of observed correlations. Disease and control sample labels were shuffled prior to computing DGA scores and disease similarity (through partial correlations). The whole process was repeated 100 times to create a background distribution of disease correlations. This was utilized to determine a permutation-based p-value (number of the times the permuted statistic exceeded the observed statistic). The background distributions for PMA and FMA scores were similarly generated.

All network visualization was performed using Cytoscape software (Shannon et al., 2003). Enrichment was identified using Gene Ontology's -Biological Process category available via ClueGO-a Cystoscape plugin (Bindea et al., 2009) and DAVID v6.8 (Dennis et al., 2003). Venn diagram was created using the “VennDiagram” package available through R/Bioconductor (Chen and Boutros, 2011).

We utilized the drug gene interaction database (DGIdb) (Wagner et al., 2015) to identify a list of expert-curated proteins/genes that serve as druggable targets. A list of currently approved drugs from the FDA (or at least one jurisdiction) was downloaded from Drugbank (Law et al., 2014), while drugs treating the disease clusters were obtained from Medscape (Medscape, 2017)¹ and UpToDate® (UpToDate, 2017)²

Expression data for diseases affecting muscle (Table 1) was processed (see Methods) resulting in a set of 12,789 genes across the 20 diseases and utilized in our analysis. The change in gene expression for each disease state (with respect to controls) was quantified as the associated T-statistic of z-normalized values (see Methods, Figure 1A, Figure S1) and referred to as differential gene activity (DGA) score. Hierarchical clustering of diseases based on partial correlation of DGA scores resulted in seven disease clusters (Figure 1B, Table 1). Twenty (~10%) of the 190 possible disease pairs were identified as being significant (p < 0.05) via significance testing (see Methods, Figure 1C). Few well-characterized intra-cluster associations marginally missed the significance threshold and were not captured, for instance the association between LGMD2B and DMD, BMD and DMD. We ascertained if the observed significant disease correlations (based on DGA scores) shared known genetic associations. We compiled a list of genes associated with our disease set from publicly available databases (see Methods). Pairs of disease were considered to significantly share disease-genes if the hypergeometric p-value of overlap was < 0.05. We identified that 26/190 possible disease pairs shared a significant genetic basis, and 6/26 interactions overlapped with our 20 significant DGA based associations (Table 3A), one-sided Fisher's exact test p-value of 0.036 (Table 3B).

Protein modules identified within a human-PPIN through modularity detection algorithms (Enright et al., 2002; Brohee and Van Helden, 2006) represented a group of strongly interacting proteins with putative functional associations. We utilized PPIN to examine if there was a common program of regulation underlying different pathologies.

We generated a catalog of 1,025 protein modules by querying a large-scale human PPIN available through STRING (Franceschini et al., 2013) (see Methods) We identified that 764/1,025 modules (containing at least 3 nodes/module), overlapped our list of 12,789 genes, with module sizes varying from 3 to 256 proteins. We defined a protein module activity (PMA) score calculated for each protein module identified per disease as the mean of component protein DGA scores. The common underlying program of regulation or “signature” protein modules were identified utilizing a 2-fold approach-first, all modules with absolute PMA values significantly higher than random (p < 0.05) in more than half the diseases (n > 10) were extracted (35 modules passed this selection criterion). Next, using the threshold for module expression (see Methods), 17/35 modules were identified as the underlying protein signature (Figure 2, Table S5). The complete list of signature modules identified and their top three enrichment terms are provided in Table S6.

Utilizing the human PPIN allowed us to infer a common program of dysregulation underlying muscle diseases. We next sought to elucidate mechanistic similarities between significant disease pairs, in a more skeletal muscle-specific context. We developed a framework of FM that represented significant units required for normal muscle activity (Table 2, see Methods). The set of 23 “FM” captured key biomarkers associated with major processes in skeletal muscle. We calculated a functional module activity (FMA) score that reflected the collective behavior of genes in each functional module, for a given disease state (Table S2). We also computed an associated p-value for each functional module via permutation testing.

The FM framework was utilized to assess common functional mechanisms underlying 20 significant disease-pairs (with p < 0.05). Table 4 presents a subset of significant disease-pairs, which shared four or more significant FM between them. For instance, JDM and DMD had 15/23 FM overrepresented (p < 0.05). These included modules associated with atrophy, inflammation, ECM and cytoskeleton, all members of the excitation contraction coupling family (FM IDs-3,4,5), members of contraction (7,8), mitochondrial energy metabolism (11,12,13,14,15), inflammation, and fiber type maintenance (19 and 21).

We ascertained if druggable targets were over-represented in the common protein signature modules to support the hypothesis that drugs targeting common targets can treat a variety of muscle diseases. We identified 54/156 proteins in the signature modules belonged to at least one druggable category, with at least one interaction as categorized in the drug gene interaction database (Wagner et al., 2015). We identified 41 of these proteins as targets for 81 approved drugs used in treating a variety of diseases and provided possible avenues for exploration of therapeutic options (Figure 3, top panel; Table S7).

We next hypothesized that protein modules exclusively regulated in each disease cluster might be enriched for drug targets. For instance, we identified protein modules associated with 3 disease clusters (DM/JDM, IBM/PM and DMD/BMD/LGMD) using an approach similar to section Identifying Protein Modules Underlying Muscle Diseases. We identified −20 approved drugs to be targeting proteins regulated in modules unique to the IBM/PM cluster, 26 drugs targeting proteins unique to the DMD/BMD/LGMD cluster and 203 in the DM/JDM cluster (Table S8).

The seven disease clusters identified included several well-characterized clusters such as muscular dystrophies (BMD, DMD, LGMD2I, and LGMD2A), mitochondrial disorders (MELAS and PEO); and a few less characterized clusters such as Sarcopenia, EDMD, LGMD2B and FSHD. Fisher's exact test was performed to capture the significance of overlap between disease pairs identified using DGA scores and those with shared genetic basis (using disease-gene lists). The Fisher's exact p-value of 0.036 implied that the genetic disease similarity was significantly captured by DGA based correlations in this study. A relatively low, albeit, significant p-value of fisher's exact test reflected the possibility that certain disease associations did not pass our significance threshold or were not captured using disease-gene lists. This could also be attributed to certain diseases being better studied than others. For example, genes associated with ALS (895 genes) where far greater than those associated with AQM (26 genes).

Enrichment analysis of 17 signature modules, identified from the protein interaction network indicated an over-representation of modules associated with immune response (e.g., module 14), mitochondrial function (e.g., modules 41, 168, 271, 355, and 340), mitochondrial structure (e.g., 537), metabolism (e.g., 426), calcium homeostasis in muscle (e.g., 334), and the extracellular matrix (e.g., 153, 416) (Figure 2, Table S3).

Early mitochondrial research in muscle disorders have suggested that a widespread occurrence of mitochondrial anomalies did not necessarily imply a primary deficiency in efficacy of mitochondrial function (muscle meeting its energy requirements) (Stadhouders and Sengers, 1987). However, more recent research has repeatedly suggested deficient bioenergetics underlie the pathology of several muscular and neuromuscular diseases in mammalian models (Wallace, 2000, 2013; Pieczenik and Neustadt, 2007; Ramadasan-Nair et al., 2014). Pathology of neuromuscular diseases such as ALS also exhibits mitochondrial dysfunction as a major event in its progression (Dupuis and Loeffler, 2009; Cozzolino and Carrì, 2012). Reduced efficiency in the action of the tricarboxylic acid (TCA) cycle has been also assessed in diseased muscle associated with inflammatory myopathies (Coley et al., 2012), dystrophy (Even et al., 1994), and mitochondrial diseases such as MELAS, PEO (Wallace, 1999). Mounting evidence has suggested that the pathological muscle wasting observed in dystrophies (e.g., DMD) might be due to reduced ATP availability required for maintenance of Ca²⁺ homeostasis and fiber regeneration (Timpani et al., 2015). Bioenergetic pathway enzymes have recently shown to be relevant biomarkers of muscular and neuromuscular disease progression (Santacatterina et al., 2015).

Ca²⁺ homeostasis in muscle largely determines its contraction and relaxation properties. This is tightly regulated by the Ca²⁺ signaling apparatus within muscle comprising of the ryanodine receptors, sarcoplasmic endoplasmic reticulum calcium pumps (SERCA), troponin complex, calsequestrin; in addition to Ca²⁺ binding proteins such as parvalbumin, sarcolipin, phospholamban and calpains. We observed a strong dysregulation of several of these proteins- ATP2A1 (SERCA pump), sarcolipin (SLN, which inhibits SERCA) and calsequestrin [CASQ, restrains Ca²⁺ to the sarcoplasmic reticulum (SR)], in ALS and DMD (Wang et al., 2012; Mukund and Subramaniam, 2015). Likewise, regulation of ASPH (regulator of ryanodine receptors) and SLN have also been observed in muscle from diseases such as CP (Smith et al., 2011). Figure 3 provides a representative sample of 4/23 signature modules and their enrichment.

Although existing research on several muscle diseases (such as ALS, DMD, BMD, and CP) has shown varying extents of mitochondrial dysfunction and calcium dysregulation in their pathomechanism, our approach points to widespread, statistically non-random dysregulation of mitochondrial function and calcium homeostasis associated with most muscle diseases including relatively less characterized diseases such as AQM and CFS. Further, the absence of modules associated with structural sarcomeric proteins (myosins, z–disc proteins, dystroglycan) at our significance threshold emphasizes the vital role of muscle bioenergetics, calcium signaling and homeostasis pathways in the pathogenesis of diseases affecting muscle.

The functional module framework utilized, allowed us to capture common functional mechanisms underlying disease pairs. Table 4 represented a subset of significant disease-pairs, which shared four or more significant FM between them (e.g., JDM and DMD having 15/23 FM overrepresented; p < 0.05). DMD and JDM represent myopathies, where the primary insult is on the skeletal muscle however, JDM is a systemic autoimmune vasculopathy characterized by weakness of proximal muscles and skin rashes with its histopathology showing evidence for necrosis, fiber size variation, and a muscle degeneration/regeneration phenotype (Peloro et al., 2001). JDM shares many pathologic similarities with muscle of children affected by DMD. A comparison of the expression profiles of children with DMD and JDM have revealed similarities in gene cascades involving muscular atrophy, deficits in mitochondrial metabolism and contraction, along with upregulation of extracellular matrix and cytoskeletal cascades (Tezak et al., 2002) consistent with the functional overlap observed in our study (see Results).

Finding relevant FM consistent with the current understanding of similarities between JDM and DMD further justified the efficacy of the adopted approach in identifying FM affected in more than one diseases state, in a context specific manner. To further elucidate disease associations much less explored, we focused on two diseases ALS and CP and their overlapping FM and associated FMA scores.

Current advancements in understanding muscular/ neuromuscular disease pathophysiology have allowed for drastic improvements in drug therapy, however, several of the diseases discussed here are as yet untreatable with high rates of morbidity and mortality with limited therapeutic options. Gene therapy and precision medicine are yet to be realized in their full potential for several of the diseases considered here (Nightingale et al., 2016; Dalakas, 2017). Given this shortage of drug/therapeutic availability for muscle diseases, we aimed to identify if drug-repurposing opportunities could be inferred from our quantitative framework across the 20 diseases. We identified 81 approved drugs to be targeting proteins contained within the signature modules (Table S7) warranting further analysis.

Further, several of the available treatments for muscle diseases currently only offer symptomatic relief, for instance, no specific therapeutic treatments exist for dystrophies such as DMD, BMD LGMD, or EDMD, with patients requiring aggressive supportive care. These patients are only often treated for associated conditions of the heart and lung. Likewise for neurological diseases which are multi-symptomatic such as ALS and CP, patients are mostly provided symptomatic relief with a multitude of agents such as antiparkinsonian, anticonvulsant, antidopaminergic, antispasticity, anti-sialorrhea, or antidepressants. This prompted us to further explore if there were drugs uniquely shared by disease clusters that could provide opportunities for drug repositioning within disease clusters.

Our results from three associated disease clusters (DM/JDM, IBM/PM and DMD/BMD/LGMD clusters) supported this hypothesis (Figure 5, Table S8), for example, several drugs such as diltiazem, cyclophosphamide, cyclosporine—prescribed for management of DM and JDM were identified associated within the DM/JDM cluster (“UpToDate” 2017)³. We also observed in the cluster presence of Ruxolitinib- a JAK inhibitor developed to treat neoplastic diseases but suggested as investigational therapy for DM (Hornung et al., 2014). Interestingly, Sirolimus a FDA approved drug for prophylaxis against organ rejection, is currently suggested for symptomatic relief in patients with DM (Nadiminti and Arbiser, 2005). Sirolimus appears to also target proteins associated in the DMD/BMD/LGMD clusters suggesting possible therapeutic opportunities in different muscle disease categories for Sirolimus (and possibly other shared drugs identified in Figure 5). These results further emphasized opportunities and a need for exploring repurposing of therapeutics in diseases affecting the muscle.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.