The PhenoClusters workflow described in this study is based on OMIM, Orphanet and HPO data and is outlined in Figure 1. The following sections describe its modules in more detail.

Phenotypes co-occurring in NMDs from Module B and phenotype-function pairs from Module C were investigated in terms of co-mention in the scientific literature, by comparing how many PubMed abstracts mention both terms in a given pair, to how many abstracts mention only one of them. To do so, the NCBI Entrez Programming Utilities API was used as previously described (Díaz-Santiago et al., 2020). In brief, the different textual descriptions for a given i) phenotype as described in HPO, ii) GO term from the biological process vocabulary, iii) KEGG pathways, and iv) Reactome pathways, were retrieved. All these terms were queried separately in PubMed (the complete database) to obtain the lists of PMIDs (PubMed identifiers) of abstracts mentioning each given term. PMID lists were then compared for each pair of terms (phenotype-phenotype, phenotype-GO, phenotype-KEGG or phenotype-Reactome) using Fisher’s exact test to detect the significantly ( P ≤ 0.05 ) co-mentioned pairs. The numbers of significant pairs obtained for each set was then compared to the number of significant pairs obtained from random models, by calculating whether the probability of finding as many significant pairs in the random dataset was at least as high as that detected using real data.

Genes in phenotype-gene associations were further analysed to see if they were more likely to be expressed in neuronal or muscular tissues, both highly relevant for NMDs. To do so, expression location of genes paired with NMD phenotypes was compared to that of genes paired with non-NMD phenotypes. Expression locations were obtained from the normal tissue expression dataset (normal_tissue.tsv.zip file) at The Human Protein Atlas version 20.0 (Uhlén et al., 2015) (http://www.proteinatlas.org). A gene was considered as expressed in neuronal or muscular tissue if it showed medium or high expression in one of the following tissue types: cerebellum, cerebral cortex, hippocampus, caudate neuronal cells, skeletal muscle or heart muscle. The relative proportions of genes from each set that showed expression in neuronal or muscular tissue were compared using Fisher’s exact test.

All the above calculations have been implemented as an automated workflow named PhenoClusters that is based on PhenCo (Díaz-Santiago et al., 2020) and uses additional scripts from (Jabato et al., 2020). The Picasso supercomputer of University of Malaga was used for code implementation and testing. It consists of an OpenSUSE LEAP 12.3 with Slurm queue system and Infiniband network (54/40 Gbps) containing 216 nodes with Intel E5-2670 2.6 GHz cores for a total of 3 456 cores and 22 TB of RAM. The code is available from GitHub at https://github.com/Elenadisa/PhenoClusters. The main requirements are Python 3, Ruby 2.4.1, R 4.0.0 or higher, Bioconductor 3.4 (Huber et al., 2015) or higher, scripts from sys_bio_lab (Jabato et al., 2020), and Anaconda Individual edition (https://docs.anaconda.com/anaconda/install/) as package and environment manager. The workflow was managed using AutoFlow (Seoane et al., 2016).

For co-occurrence analysis, PhenoClusters requires NetAnalyzer (Rojano et al., 2017). Although many co-occurrence methods are available, previous work suggests that the hypergeometric index (HyI) (Fuxman Bass et al., 2013), which can be considered analogous to a contingency table based approach, is the most suitable for co-occurrence studies (Rojano et al., 2017; Bueno et al., 2018). Hence, H y I ≥ 2 was set to calculate associations within the layers of a network and obtain phenotype-phenotype or phenotype-gene associations that were considered significant co-occurrences (Díaz-Santiago et al., 2020).

Functional enrichment was performed using the genes associated with each phenotype. The Bioconductor package clusterProfiler 3.18.0 (Yu et al., 2012) was used to extract the biological process subontology from GO and the KEGG pathways and conduct over-representation analysis for these resources. The Bioconductor package ReactomePA 1.34.0 (Yu and He, 2016) was used to extract Reactome pathway data and perform over-representation analysis for this resource. Both packages were executed with default parameters (pvalueCutoff = 0.05, pAdjustMethod = “BH”, universe, qvalueCutoff = 0.2, minGSSize = 10, maxGSSize = 500) and only the assignments with P ≤ 0.05 after Benjamini-Hochberg multiple testing adjustment (Benjamini and Hochberg, 1995) were considered significant.

Highly interconnected phenotype communities were obtained using the R-CRAN package linkcomm 1.0.13 (Kalinka and Tomancak, 2011) with default parameters for undirected networks. R-CRAN packages required for graphical representations and data management are ggplot2 3.3.2, RColorBrewer 1.1.2, igraph 1.2.6, dplyr 1.0.2 and VennDiagram 1.6.20. The final, user-friendly HTML reports were produced using R markdown packages rmarkdown 2.6, knitr 1.30 and kableExtra 1.3.1.

Details about PhenoClusters execution can be found at https://github.com/Elenadisa/PhenoClusters. As can be seen, diseases from OMIM and Orphanet were considered separately in different scripts. One run usually takes at least 6 h 45 min without the co-mention verification, since literature analysis takes a very long, heterogeneous amount of time depending on the PubMed server overload. The resulting reports are named omim_report.html (Supplementary Material S1) and orphanet_report.html (Supplementary Material S2). To simplify cluster inspection, they are also saved as independent files in omim_clusters_report.html and orphanet_clusters_report.html including the cluster ID, the phenotype graph together with HPO IDs and descriptions, as well as the associated functions and the genes identified within them. Clusters showing functional coherence are given at the top of the file. For the sake of simplicity, “*” is used to tag OMIM clusters, while Orphanet clusters are tagged by “#”, as indicated in Module D.

Diseases in OMIM and Orphanet were classified as NMD or non-NMD (Table 2). All downstream analysis treated the resources separately (Supplementary Material S1, S2). The total numbers of diseases in each resource for which genes and pathological phenotype information was available are shown in Table 2; full details can be found in Supplementary Material S3, S4. It can be seen that relatively few diseases were considered NMDs with known genes and phenotypes in HPO (424 and 125 in OMIM and Orphanet, respectively). Many phenotypes were shared by both NMDs and non-NMDs (row “Common” in Table 2). Fisher’s exact test was then performed to compare the relative occurrence of each phenotype between disease groups to find those relatively more common in NMDs. This resulted in 335 and 216 typical NMD phenotypes for OMIM and Orphanet. Overlap between the two is shown in Figure 2, where Panels B and C show that the proportion of common phenotypes in both sets increases with decreasing cutoff p-values. This reinforces the suitability of the co-occurrence approach to extract typical NMD phenotypes.

The top over-represented phenotypes are shown in Supplementary Material S5. The most significant in both resources is ‘elevated serum creatine kinase’, which reflects the important role of muscle decay in NMDs, given that the normal location of this enzyme is the cytoplasm and mitochondria, not serum (Moghadam-Kia et al., 2016). This result is reassuring with respect to the approach of PhenoClusters, as elevated serum creatine kinase is routinely tested when diagnosing NMDs as a consequence of muscle injury (Thavorntanaburt et al., 2018).

The top phenotypes (Supplementary Material S5) seem to be more indicative of muscular impairment than neuronal dysfunction, pointing to proximal or distal muscular weakness, both very typical in the clinical history of NMD patients (McDonald, 2012; Mukherjee et al., 2019). Other important phenotypes are related to specific diseases, such as Duchenne/Becker muscular dystrophy (‘calf muscle hypertrophy’, ‘Gowers sign’), where muscle hypertrophy can exaggerate postural instability and joint contracture (Kornegay et al., 2012). In contrast, the absence of ‘hypotonia’-related phenotypes in this top ten list is notable. However, it should be made clear that such phenotypes were also significantly over-represented among NMDs, just not included in the top ten (Supplementary Material S3, S4). To conclude, the PhenoClusters approach to obtain statistically over-represented phenotypes for NMDs extracts phenotypes characteristic of these disorders.

Bipartite networks were created (one for OMIM diseases and another for Orphanet), with nodes representing diseases and phenotypes, and edges representing known relations in OMIM or Orphanet respectively. The network was analysed to identify pairs of typical NMD phenotypes that tend co-occur across diseases (Module B in Figure 1). Around a quarter of all phenotype pairs co-occurring in at least one disease showed significant co-occurrence (27.61% in OMIM and 21.20% in Orphanet). Only a small subset (≤1 %) of all phenotype pairs represent significant co-occurrence between NMD phenotypes. However, this value was far lower for the random sets of non-NMD phenotypes, being almost negligible (Table 3), showing that pairs of NMD phenotypes are far more likely to significantly co-occur across multiple diseases than pairs of randomly selected phenotypes. These results point to the presence of comorbidity between NMD phenotypes.

As well as showing that the typical NMD phenotypes tended to co-occur across diseases much more frequently than non-NMD phenotypes, further validation was performed based on co-mention analysis for each phenotype pair within PubMed abstracts. The NMD phenotype pairs were shown to be co-mentioned in the literature much more frequently than equal numbers of pairs of randomly chosen non-NMD phenotypes (Table 3). Fisher’s exact test analysis revealed that 1 173 OMIM pairs and 689 Orphanet pairs (Supplementary Tables S1, S2) were significantly co-mentioned. Conversely, for non-NMDs, the numbers of significant pairs were much lower ( 333.08 ± 34.63 in OMIM and 110.0 ± 13.47 in Orphanet, Supplementary Tables S1, S2, respectively). These differences were also supported by the distribution of their p-values (Supplementary Figures S1, S2) that were clearly more dense at P ≤ 0.05 for NMDs than in non-NMDs. These results demonstrate the validity of the PhenoClusters approach to find pairs of NMD-related phenotypes that tend to co-occur together.

Associations between NMD phenotypes and genes and functions were calculated in order to uncover potential underlying mechanisms leading to the expression of a pathological phenotype. This was done by initially connecting phenotypes with genes using a tripartite network based approach (Module C in Figure 1) by combining all disease-phenotype pairs and all disease-gene pairs in OMIM and Orphanet separately. The networks were then analysed in order to connect genes with phenotypes based on co-occurrence across multiple diseases.

Using a threshold of H y I ≥ 2 , 34 139 and 36 218 phenotype-gene associations were obtained for OMIM and Orphanet, respectively; from these, a total of 1 078 and 611 distinct genes were associated with NMD phenotypes, respectively. There was significant overlap between OMIM and Orphanet in terms of genes associated with NMD phenotypes ( P < 2.2 × 10 − 16 , Fisher’s exact test using all phenotype-associated genes as background). Moreover, a significantly greater proportion of the NMD associated genes (84.41% in OMIM, 86.73% in Orphanet) were shown to be expressed in neuronal or muscle tissue than the non-NMD associated genes (70.32% in OMIM. 72.88% in Orphanet), which is a statistically significant difference ( P = 2.2 × 10 − 16 and P = 6.3 × 10 − 11 , for OMIM and Orphanet, respectively). The same kind of test was performed for the remaining non-relevant tissues, revealing in this case that there was no significant difference between NMD and non-NMD genes after a Fisher’s test ( P > 0.9 for both databases).

The genes associated with each phenotype were then used for functional enrichment analysis using clusterprofiler and ReactomePA to find over-represented ( P ≤ 0.05 ) GO terms (biological process ontology), and KEGG and Reactome pathways (Table 4). As with gene-based analysis, there was significant overlap between OMIM and Orphanet ( P < 2.2 × 10 − 16 , 6.842 × 10 − 8 , and 2.381 × 10 − 14 , respectively, using Fisher’s exact test with all phenotype-associated functions used as background). Supplementary Material S6 shows the overlap between the different resources using Venn diagrams. Therefore, genes associated with NMD phenotypes are significantly expressed in neuronal and muscle tissue, the most relevant for NMDs, supporting the robustness of the phenotype-function pairs obtained by PhenoClusters.

In spite of the robustness of these novel associations between phenotypes and functions, additional validation was performed based on co-mention analysis using PubMed abstracts as a sign of relevance. Term co-occurrence of a given phenotype alongside its corresponding function was looked for within the same abstract (Supplementary Figures S1, S2). Only a subset of the phenotype-function pairs obtained from the network were significantly co-mentioned (“Confirmed” columns in Table 5). However, this subset was several times larger than the number of co-mentioned pairs within the “random” pairs, and this difference was significant for all comparisons ( P < 2.2 × 10 − 16 ) using Fisher’s exact test. Consistent results were found for both OMIM and Orphanet: respectively i) 5% and 4% of phenotype-GO pairs, ii) 10.5% and 9% of KEGG pairs, and iii) 4% and 3% of Reactome pairs were significantly co-mentioned. As such, PhenoClusters (Figure 1) appears to produce robust and reliable associations between NMD phenotypes and functions based on phenotype and gene co-occurrence across diseases. Even though only a small fraction could be confirmed in PubMed abstracts, this was several times more than would be expected by chance. The analysis showed much higher co-mention for the KEGG associations as a proportion of the total pairs. This supports our previous work suggesting that information in KEGG is more reliable for functional studies in certain situations (Luque-Baena et al., 2014).

The pairs of significantly co-occurring NMD phenotypes (pale magenta rectangle of Module B in Figure 1) were processed to extract communities of highly interconnected phenotypes (Module D in Figure 1). As summarised in Table 6, the number of communities and average number of phenotypes per community is higher using the modules based on NMD phenotypes compared to non-NMDs. This indicates that the co-occurring NMD phenotypes tend to form groups of related phenotypes.

The phenotype communities were combined with the phenotype-function associations obtained previously (light cyan rectangle of Module C in Figure 1) to evaluate their functional coherence. In total, 40 OMIM communities and 72 Orphanet communities (Table 6) showed coherent function, defined as having shared functional enrichment for at least 70% of their constituent phenotypes; they constituted the “functionally coherent clusters”. The complete list of functionally coherent clusters is provided at the end of the reports given in Supplementary Material S1, S2. In terms of cluster properties, some have high interconnection (*193, *138, *219, #115, #150), others show spoke-hub structure, with a central node connected to a number of other nodes which do not then connect to each other (clusters *58, *1, *3, *4, *91, *26, #144, #2, #3, #21, #60, #22), and others have near linear topology (*146, *85, #97). Compared to their respective random non-NMD models Table 6, OMIM presents a similar number of clusters, in stark contrast with Orphanet, for which 5 times more clusters were found. The coherence threshold must be decreased to 50% to find a substantially different number of coherent clusters in OMIM compared to the random datasets (Supplementary Figure S1), This difference may be due to the OMIM clusters being larger than for Orphanet (16.75 vs. 13.63 phenotypes, respectively) or related to the nature of the diseases considered NMDs in each dataset. These results support the approach presented here in which the information from the different databases is considered separately, with Orphanet arguably providing more confident results with respect to how well the typical NMD phenotypes in these clusters represent realistic and useful groupings of co-occurrent phenotypes that reflect co-morbidity with shared underlying cellular functions.

The functionally coherent clusters found by our approach (Supplementary Material S1, S2) were inspected by clinicians with expertise in paediatric neuromuscular diseases, finding clusters with sets of expected phenotypes for NMDs, and others with unexpected phenotypes. Clusters are grouped in Table 7 based on the area in which the constituent phenotypes manifest, followed by further subdivision by group of NMD. Since they represent direct links between NMDs, their symptoms and their underlying mechanisms, they could be used for clinical diagnosis, by suggesting novel phenotypes to test for, should a patient present other phenotypes within a cluster. They also suggest potentially affected genes, useful for directing genetic analysis, and potentially affected functions. Although this analysis necessarily has a subjective aspect, given that clinicians differ in training and experience, it is remarkable that PhenoClusters produces functionally coherent clusters easily recognised by clinicians. In summary, there are multiple clusters that fit in with our current knowledge about NMDs as shown in Table 7, which reinforces once again the coherence of the results generated by the presented approach.

Other clusters provided by PhenoClusters contained more unusual groupings of phenotypes. Clinicians were surprised by the presence of ‘macroglossia’ (HP:0000158) and ‘arthogryposis’ (HP:0002804) (Figure 3). Macroglossia is usually associated with non-neuromuscular syndromes and in metabolic disorders such as glycogen storage disease type II (Pompe disease), mucopolysaccharidosis, oligosaccharidosis, mucolipidoses, sphingolipidoses and galactosidosis. However, macroglossia appeared as a hub node in the cluster *58 (Figure 3), strongly associated with highly representative NMD phenotypes, such as ‘congenital muscular dystrophy’ and ‘achilles tendon contracture’. This phenotype also appears in clusters *26, and *53, with macroglossia occupying a peripheral position (Figure 3). In macroglossia-containing clusters, 70%–80% of phenotypes show significant enrichment for the “protein O-linked mannosylation” biological process (GO:0035269), suggesting that the clinical observation of macroglossia in patients with NMDs may prompt practitioners to ask for O-linked mannosylation tests. This contrasts with cluster *211 (Supplementary Material S1) that also contains congenital muscular dystrophy and ankle flexion contracture, as well as ‘increased endomysial connective tissue’, but not macroglossia, the functions for this cluster being ‘focal adhesion’ (hsa04510) and ‘ECM proteoglycans’ (R-HSA-3000178) but not “protein O-linked mannosylation”. Hence, cluster *211 seems to be related to processes involved in physically connecting (‘focal adhesion’) cells to the extracellular matrix (‘ECM proteoglycans’). Therefore, functionally related clusters may help to discriminate disorders presenting with the macroglossia phenotype from other muscular dystrophies caused by the disruption of other cell functions. This finding shows the potential clinical utility of our approach to identify phenotypes whose presence alongside a given set of other phenotypes can indicate distinct underlying processes, with implications for diagnosis and treatment.

‘Arthogryposis’ represents another phenotype of potential use for NMD diagnosis (Figure 3). Arthrogryposis multiplex congenita consists of contractures in multiple body areas (Ambegaonkar et al., 2011). Cluster #117 shows a range of clinical phenotypes with arthrogryposis as the central node (Figure 3), while #59 and #88 are cases where this phenotype appears in a peripheral position. Interestingly, cluster #115 (Supplementary Material S2) presents almost all phenotypes observed in cluster #117, except arthogryposis. Hence, it can be inferred that when arthogryposis is present, the pathways associated with the majority of phenotypes in the cluster are more specific to muscle fibre-related processes including contraction dysfunction. However, when arthogryposis is absent, the cluster phenotypes share additional pathways related to cardio-myopathies and the regulation of pseudopodium assembly, indicating a broader aetiology in the genetic systems involved. As such, the presence or absence of arthogryposis alongside other symptoms may help inform NMD diagnosis and treatment, by indicating distinct underlying processes.

The genes underlying the enriched functions were further investigated. It was noted that some of them occur in multiple clusters but are linked to different functions depending on the affected phenotypes. For example, DAG1 is a pivotal component of the dystrophin-glycoprotein complex and its dysfunction is related to many muscular dystrophies, amongst other diseases (Durbeej et al., 1998; Sciandra et al., 2003). It can be considered a pleiotropic gene due to its role in glycosylation, a process which is involved in myriad cellular functions. For NMDs, it seems that glycosylation pathways are affected when DAG1 appears in clusters *1, #29 and #113, in agreement with previous studies (Barresi and Campbell, 2006). However, the involvement of DAG1 in clusters *198, #2, #144, #112, #103, #118, #11, #89 and #63 (Supplementary Material S1, S2), appears to be related to cardiomyopathy (Michele et al., 2009), while it is more focused on ‘ECM proteoglycan components’ in cluster *211. Hence, depending on the cluster in which DAG1 appears, it seems to be involved in different pathways and functions, showing how PhenoClusters facilitates the identification of genes with multiple putative roles in NMDs depending on context.

Other clusters suggest that members of the solute carrier (SLC) family can perform many roles within NMDs. SLCs mediate the transport of a wide range of essential nutrients and metabolites, performing many different functions in cells and tissues (Zhang et al., 2019). SLC proteins were found to be associated with nutrient supply in cluster #93; metabolic transformation in clusters *85, *128, and *218; energy homeostasis in cluster *138; oxidative stress in cluster *205; and neurological regulation in clusters *97 and *166. These findings show how the phenotypic components of the functionally coherent clusters defined using PhenoClusters can point to the disruption of different cell functions in relation to NMDs involving genes from the SLC protein family.

As the use of bioinformatics analysis becomes more routine within biomedicine, it is crucial that published studies are accompanied by workflows that allow the analysis to be recorded and reported in a reproducible manner that can be applied to similar datasets, if applicable. PhenoClusters is based on phenotype co-occurrence, incorporating randomisation using matched non-NMD phenotypes to generate control datasets and including literature co-mention for validation of the approach. It also uses known NMD-causing genes to provide functional support to the co-occurring phenotypes. It produces a HTML report that can be easily interpreted, following the proposed Reproducible Research System (Mesirov, 2010; Piccolo and Frampton, 2016). Although there are many articles in the bioinformatics and systems biology fields presenting modular software tools and packages, in many cases code and data are not provided, despite the importance of reproducible research in health sciences (Harris et al., 2018). As it is essential for the present study, the necessary code to conduct the workflow as described in Figure 1 has been provided in full. Since the analysis starts by downloading diseases from OMIM and Orphanet, it can be adapted to study phenotype/gene/function co-occurrences in other diseases. PhenoClusters is thus an example of the re-analysis of previously published data using a new approach that produces new results and confirms already known facts. As such, it has achieved the goals of reproducible research (Mesirov, 2010; Piccolo and Frampton, 2016) in the health sciences (Harris et al., 2018) and exploiting already published results (Kovalevskaya et al., 2016).

PhenoClusters did not mix data from OMIM and Orphanet based on two main facts. The first one is that OMIM and Orphanet are populated in different ways. For example, the variability of the phenotypic series of OMIM was high, with PS253600 (muscular dystrophy, limb-girdle, autosomal recessive) presenting 28 entries, while PS118220 (Charcot-Marie-Tooth disease) has 70 entries. The second one is that OMIM was searched with somewhat arbitrary keywords, whilst Orphanet NMDs were obtained based on its ontology structure. As a result, there are more OMIM-specific phenotypes and less Orphanet-specific phenotypes (Figure 2), and fewer functionally coherent clusters in OMIM than in Orphanet (Table 6). Besides these discrepancies, both databases rendered a similar ratio of shared phenotype pairs (Table 3) and the same top specific phenotype for NMDs, ‘elevated serum creatine kinase’ (Supplementary Material S5), in agreement with the importance of the creatine kinase test often used in NMD diagnosis (Thavorntanaburt et al., 2018). Other common phenotypes include muscle weakness, contractures, altered gait, functional difficulties, and respiratory issues, in agreement with recognised phenotypes for these diseases (Norwood et al., 2009). However, the interesting phenotypes ‘macroglossia’ and ‘arthogryposis’ (Figure 3) emerged from database-specific group of phenotypes. In any case, despite the mentioned discrepancies, it is striking that, when the phenotypes were associated with genes and functions, there was high overlap (Table 4, Supplementary Material S6). Therefore, we believe it is preferable to use the information from OMIM and Orphanet separately and to compare to contrast the results.

Co-mention validation in Table 5 revealed that even though GO contains many terms to describe biological functions and Reactome contains more pathway descriptions than KEGG, 2-fold more phenotype-KEGG pathway pairs (percentage-wise) were validated by literature co-mention. Moreover, when comparing functional annotation between OMIM and Orphanet, a proportionally larger overlap was found between the two resources for KEGG pathways (Supplementary Material S6). These findings supports our previous finding that KEGG produces better bioinformatic models in genetic algorithms for clinical diagnosis and prognosis (Luque-Baena et al., 2014), as well as for the outcome of diseases (Urda et al., 2018). Hence, KEGG pathways should always be considered in functional analyses.

Due to the etiological heterogeneity of NMDs, each disease is defined by its own set of clinical phenotypes. Phenotypes serve to understand life and disease, but it is not always easy to translate them to molecular mechanisms and vice versa (Yu et al., 2016). Taking advantage of the huge amount of data available in databases, bipartite networks (phenotype-disease, disease-gene) have been used in biomedicine to model factors that influence human diseases, explore their molecular complexity, reveal novel molecular relationships and disease susceptibility genes, uncover the biological significance of disease-associated mutations (Pavlopoulos et al., 2018) and, more recently, discover phenotype/disease clusters that can predict protein-phenotype associations and reveal the underlying mechanisms that link them (Ren et al., 2020). Bipartite gene-disease network analysis was revolutionised with the diseasome and related studies (Goh et al., 2007) which showed that genes associated with similar disease phenotypes have a higher propensity to interact physically with each other, forming distinct disease-specific functional modules, and that disorders tend to form clusters on the basis of similar pathophysiology (Oti and Brunner, 2007). Other studies revealed that phenotypes and gene co-regulation accurately predict unknown disease-gene relations (Deelen et al., 2019) since genes causing a specific disease or disease symptom often have similar molecular functions or are involved in the same biological process or pathway.

With all these facts in mind, the systems biology approach implemented in PhenoClusters lies in bipartite and tripartite networks (Modules B and C in Figure 1) to extract co-occurring phenotypes in NMDs (Supplementary Material S3, S4) that were consistent with i) clinical histories (Norwood et al., 2009) (Table 7), ii) creatine kinase tests (Moghadam-Kia et al., 2016; Thavorntanaburt et al., 2018) (Supplementary Material S5), iii) proximal or distal muscular weakness (McDonald, 2012; Mukherjee et al., 2019), and iv) postural instability and joint contracture (Kornegay et al., 2012), as usually found in the clinical history of NMD patients. The genes significantly associated with NMDs by PhenoClusters were largely expressed in neuronal and/or muscle-related tissues. For those that were not, it is tempting to speculate that this is due to incomplete information in The Human Protein Atlas, however it may also indicate genes involved in development or other regulatory processes that lead to the manifestation of the NMD phenotypes. There is currently much interest in the relationship between tissue/cell-type specificity and disease progression (Hekselman and Yeger-Lotem, 2020). Functional enrichment of the NMD-associated genes was used to provide phenotype-function associations that were then combined with the co-occurrent phenotype communities to produce functionally coherent clusters. Gene-to-function translation reduced the degrees of freedom of the analysis (several genes were required to assign a significantly enriched function to a phenotype), increasing statistical power and allowing the detection of shared functions between phenotypes, even if the underlying genes differ. The comparison of NMD-specific associations with randomised data in Tables 3, 4, 6 demonstrates that the findings were well-founded, as discussed below.

There is an absence of gold-standard datasets for phenotypic relationships, particularly for rare diseases. As such, Ren et al. (2020) applied systems biology methods to obtain phenotype/disease clusters to seed machine learning models to predict protein-phenotype associations without any objective validation, assuming that they were relevant since authors can find biological support for some of the clusters. Here, abstract co-mention was proposed as a source of external validation for phenotype-phenotype and phenotype-function pairing approach. There were far fewer co-mentions for randomly paired terms than for specific pairs (Table 5), even though only a small fraction of phenotype-phenotype and phenotype-function pairs ( < 10   % in the second) showed significant co-mention in the literature. This is likely due to a combination of the following reasons: 1) term relationship existed, although not yet described in literature; 2) the co-mention appeared in the main body, but not in abstracts, of articles; 3) phenotypes were co-mentioned with genes rather than functions; and 4) abstract co-mention existed, but using slightly different terms or natural language synonyms that escape the search. In any case, these issues also apply to the randomised/non-NMD associations. In conclusion, although the abstract co-mention is not perfect, the fact that the amount of significant literature co-mention was much higher for NMD-specific pairs than controls (Table 5) shows it to provide suitable validation of the described approach and that real relationships were found.

Many PhenoClusters clusters are in clear agreement with current medical knowledge about the relationships between phenotype comorbidity, NMDs, their pathophysiology and genetic causes (Table 7). This serves to further support the bioinformatic approach presented here. Regarding differential diagnosis and treatment, the clusters for which the presence/absence of a specific phenotype (Figure 3) could change the nature of the associated functions are particularly interesting: these have potential applications for differential diagnosis and treatment selection. As such, ‘macroglossia’ can be used to identify O-mannosylation alterations in a given NMD dysfunction (Spence et al., 2002; Martin, 2007; Moore and Winder, 2010; Dobson et al., 2013; Goody et al., 2015), while ‘arthogryposis’ (Ambegaonkar et al., 2011) may facilitate NMD diagnosis in other situations, as its presence may indicates that the disease is more specific to muscle fibres and contraction dysfunction.

The advantages of using functions rather than genes only is particularly evident in the findings related to pleiotropic genes and gene families involved in multiple functions. DAG1 is involved in many diseases because it participates in glycosylation, whose malfunction affects many cellular processes and can produce several NMDs (Durbeej et al., 1998; Sciandra et al., 2003). However, its specific involvement in the proteoglycan components of the ECM seems to be related to cardiomyopathy. Additionally, despite the clear relevance to health and disease of the 190 different SLCs found mutated in human diseases, it is considered an understudied family (César-Razquin et al., 2015). Functionally related clusters from PhenoClusters reveal which particular role can be putatively assigned to the present members of the SLC family. Taken together, these results encourage future research lines directed to systematically exploit co-occurrence and clusters of co-occurrent phenotypes and functions in other diseases, as well as further study of a number of presented clusters to obtain more details on the functional implications.

It has been shown that functionally coherent clusters have potential utility in NMDs in terms of better understanding clinical presentation in these diseases and obtaining clues as to the underlying molecular genetic mechanisms. The former has important implications for diagnosis and patient classification; the latter will help researchers better understand these diseases and search for potential new therapeutic targets.

Their potential for diagnosis is clear: given a patient that presents a number of phenotypes corresponding to one of the clusters, one can make inferences about other clinical phenotypes that the patient may also suffer and that should be tested for. This will help obtain a better clinical profile for the patient, aiding diagnosis. Moreover, genes and functional terms associated with the same cluster can be of use for guiding genetic diagnosis and indicate further pathological examination. For example, multiple clusters in the peripheral neuropathy group in Table 7 contain phenotypes associated with myelination and axon ensheathment related genes (Zhou and Notterpek, 2016; Kamil et al., 2019). As such, given a patient showing phenotypes belonging to one of these clusters, but for whom full diagnosis has not yet been achieved, the clinician could refer the patient for the appropriate assays, such as MRI, to look for demyelination and guide genetic studies to focus on related genes.

A further clinical application is related to differential diagnosis based on the phenotypic profile presented by a patient. The presence/absence of phenotypes ‘macroglossia’ and ‘arthogryposis’ (Figure 3) has been thoroughly explained as examples of phenotypes indicating differing underlying mechanisms that lead to disease manifestation. However, many other examples undoubtedly exist among the clusters presented in Supplementary Material S1, S2.

Nevertheless, there are several limitations in our bioinformatic approach, largely related to the external databases from which the information is acquired. OMIM and Orphanet are both incomplete resources: there are undoubtedly multiple diseases suffered by people around the world whose description and genetic basis have not made their way into these databases. This means that we are likely to miss co-occurrent phenotypes as the diseases in which they co-occur have not yet been added. Similarly, the resources used for functional enrichment analysis, used to ascribe function to the NMD phenotypes in this study, are also incomplete. We do not yet know the function of all genes in the genome, and this has effects on the enrichment procedures. As these resources improve, the potential of our workflow to find clusters of co-occurrent phenotypes and fully characterise their underlying basis will no doubt improve with them. However, regarding the information of clusters, this approach can only extract functions from genes known to be involved in NMDs, and cannot hypothesise about other genes even if they are working in the same pathway. And as already mentioned above, there is a limitation concerning literature co-mention validation (Table 5), as a higher proportion of confirmed pairs it would be desirable.

In conclusion, our approach and the resultant functionally coherent clusters of NMD phenotypes (Supplementary Material S1, S2) can 1) relate phenotype co-occurrence across NMDs to the underlying genes and mechanisms involved in the NMD (or other diseases); 2) provide clinicians with hints about clinical tests to produce a more reliable diagnosis based on the presence or absence of some specific phenotypes that have not yet been reported by the patient or the clinician, in the comorbidity context of the patient; 3) give researchers clues to perform new experiments to discover the underlying biological mechanisms of a disease; and 4) help in selecting optimal treatment. Therefore, PhenoClusters can be considered a new tool for more accurate diagnosis and an advance towards personalised medicine for NMDs.