The detection and analysis of community structure in networks has received considerable attention in recent years [1-3]. Community structure based clustering is the division of a network into groups of nodes with relatively dense connections within the group and sparser connections to other groups in the network. In many networks, finding the underlying groups provides practical and important information. Groups within social networks for example, may correspond to social units; in biological networks it helps to reveal substructures with common functionality or association with diseases [4-6]. Many algorithms now exist to solve the problem of breaking up the network into its communities, but most struggle to scale with current biological datasets.

We focused on a sub-set of community detection algorithms, with the goal of applying these to protein interaction networks. Nodes in these networks represent biomolecules, such as genes or proteins, and the edges the structural or protein interactions connecting them. We implemented in C++ a suite of three algorithms, optimised for the scalable analysis of large interaction networks based on the Modularity measure containing: Geodesic and Random Walk edge Betweenness, and Spectral Modularity; we also include a Cytoscape [7] App (version 3.0) for the Spectral Modularity algorithm. We tested the suite’s performance on synthetic networks, on the scale of 1000 proteins to represent realistic biomolecular complexity on multi-core workstations. We also tested the suite’s clustering performance on two previously published datasets [8,9], and investigated the enrichment of the clusters obtained for functional annotations and common neurological and neurodevelopmental diseases/disorders. Evaluation of enrichment was performed using a combination of Hypergeometric and permutation tests. We evaluated the robustness of clusters found with our suite’s boot-strap facility, and against commonly used Modularity based clustering algorithms available [10,11].

Synthetic networks were generated using the [12] benchmark and modelled on the order of 1000 proteins and 5000 interactions as described in Table 1 (Networks 2, 3 and 4, see also Supplementary Table S1). Three previously studied biological datasets were included. The MASC complex, representing a protein complex surrounding the mammalian NMDA receptor [8] consists of 101 proteins and 246 interactions (Network 1 in Table 1, see also Supplementary Table S2). The second, from a list of 1461 proteins obtained from a study of the PostSynaptic Density (PSD) in the human brain [9]. Protein-protein interactions were obtained by mining publicly available databases: HIPPIE [13], BioGRID [14], IntAct [15] and performing an InterologWalk over different species using Bio::Homology::InterologWalk [16]. This second network (Network 5 in Table 1, see also Supplementary Table S3) consists of 1312 proteins and 8031 protein interactions. The third is the Human interactome network BioPlex [17]. This network (Network 6 in Table 1, see also Supplementary Table S7) contains 7668 proteins and 23744 protein interactions, found using high-thoughput affinitypurification mass spectrometry in human embryonic kidney (HEK) 293T cells.

Table 1: Network.

Modularity (Q) measures the quality of a particular network division into communities. It measures the number of edges found within the communities relative to the expected number of edges within the communities, if these edges had been placed at random. The Modularity based algorithms we implemented [18,19] are largely designed for single workstation deployment which is not ideal for bootstrap analysis of large networks. Further we found no readily available implementation of the Random Walk algorithm [18]:

(1)

In eqn-1 the degree of the i^th vertex is given by . Where Aij is the interaction matrix, m the total number of edges in the network and the final term the Kronecker delta function, which equals 1 if the i^thand j^thnodes fall into the same community and 0 otherwise. The maximum value of Modularity is 1, with typical values for real networks ranging from 0.3 to 0.7. We can optimise clustering in the Geodesic and Random Walk algorithms by searching for the global maximum.

The third algorithm, Spectral Modularity, makes use of the spectral properties of the network [19]. The spectral properties are obtained using the Modularity matrix (B) using eqn-2:

(2)

The network is expressed in terms of its eigenvectors and eigenvalues and recursively partition into two communities; positive values of the eigenvector giving a set of nodes belonging to one community, and the negative values to another. The partition can be expressed as a column vector (S) of the nodes with values ± 1. Partitioning occurs if the maximum positive eigenvalue is greater than the tolerance (10-5) for the current partition, and if it results in a positive contribution to the Modularity. The change in Modularity is calculated from this leading eigenvector and the generalised Modularity matrix (B(g)) for the split is given by eqn-3:

(3)

where the generalised Modularity matrix for a given split is expressed as,

(4)

Given an initial separation using the Spectral method, it is possible to maximise for the change in Modularity using a fine-tuning step. The first stage here is to find the node which, when moved from one community to the other, gives the maximum change in Modularity. This node’s community is then fixed and we repeat the process until all nodes have been moved. The whole process is repeated from this new state until the change in the Modularity, between the new and old state, is less than the predefined tolerance. Our implementations were tested on the widely used benchmark, Zachary’s “karate club” network, which represents observed social patterns between members in a university sports club. Our implementations reproduced the results as reported in [18,19].

The suite is designed to run on computing clusters and includes a boot-strap facility, allowing a random sub-sample of the data to be selected; currently set to 80% of the network’s node size. Sampling from multiple bootstrap runs allows the robustness of each algorithm applied to the data to be investigated. The package clusterCons [20] has been used in conjunction with the suite, to build a consensus matrix from which to test the robustness of discovered communities, and proteins found inside each community.

Our implementations were found comparable in terms of speed and scalability with other implementations, when a direct comparison could be made, as shown in Figure 1a. However, performance of the sequential Geodesic and Random Walk implementations was found inadequate over our datasets, as shown in Table 2. To improve their performance required re-programming, by parallelising the implementations using OpenMP. This allows the implementations to make use of the multiple processors on common a workstation, and still remain portable to run on computing clusters. The OpenMP implementations of the Geodesic and Random Walk algorithms parallelise the betweenness score calculations of the edges. Both implementations were tested and found to perform optimally for multi-cores as discussed in Section Speed and Scalability.

Figure 1: CPU-time verses network size for the Geodesic and Random Walk edge Betweenness algorithms.
CPU-time verses number of network edges for Geodesic (A) and Random Walk (B) implementations as given in Table 2. Showing R’s igraph package [10] implementation (crosses), single core (circles), 4 cores (squares), 8 cores (up-triangles) and 12 cores (down-triangles). Figures show both implementations scale to large network size using multiple cores. Figure (B) shows the sequential Random walk implementation unscalable at largest network size, and a predicted linear time (dotted line) given.

Table 2: Single core CPU-times.

A quantitative and qualitative discussion of our suite’s clustering performance applied to the MASC, human PSD and BioPlex datasets is given in Sections MASC network, Human PSD network and BioPlex network. For this we made use of the functional annotation studies [8,21] and the topOnto package (https://github.com/statbio/topOnto), to perform the disease ontology enrichment analysis. Annotation enrichment of the clusters was carried out using a combination of the Hypergeometric test and permutation study as described in Section Clustering Performance. We find our implementations reproduce results found in the MASC network study, whilst our Spectral Modularity algorithm out performs others when applied to the human PSD network. The increase in diverse functional clusters, found in the human PSD network, with our Spectral Modularity algorithm is a result of our fine-tuning step as discussed in Section Comparison with Existing Cluster Algorithms. We further tested and compared the clustering performance of our suite’s Spectral Modularity algorithm on the larger BioPlex network [17].

Speed and scalability

We checked the suite’s performance against single and multi-core instances (2.4 GHz Intel Xeon “Westmere” with 12-cores) using gcc version 4.4.7. Table 2 details the CPU time in minutes for each algorithm with respect to varying sized synthetic networks as given in Table 1. Figure 1 shows the execution times for the sequential and OpenMP implementations of the Geodesic and Random Walk algorithms (see Supplementary Table S1 for details). Our Geodesic edge Betweenness algorithm was tested against the current known R implementation [10] with comparable performance (Figure 1a). For the Random Walk algorithm on Network 5, the execution time using a single core was over 14 days (20160 minutes), a limit beyond which the algorithm was considered unscalable. The predicted time of 26 days (37889 minutes) is shown by the dotted line in Figure 1b. For Network 4, the Random Walk algorithm took approximately 5 days (6827 minutes) to complete. This execution time dropped to 23 hours (1360 minutes) using the OpenMP implementation with 12 cores. As illustrated in Figure 1 there is no significant performance difference for networks smaller than 600 nodes (or less than 4000 edges). For networks larger than this size the performance difference is significant.

Clustering performance

The performance of each algorithm on the MASC and human PSD networks was studied. Tables 3 and 4 summarise the results with clustering consistency tested using the Normalised Mutual Information (NMI) and Adjusted Rand Index (ARI) metrics to the MASC and human PSD networks respectively. The original MASC study provides an analysis of the organisation and underlying functionality of the modularised MASC complex, clustered using the Random Walk algorithm and shown in the Figure 3 of [8]. We replicated this using our Random Walk algorithm in Figures 2a, and with our Spectral Modularity algorithm as shown in Figure 2b, visualised using the package Visone [22]. We generally find one or two large clusters (coloured in blue) containing molecules commonly associated with signal processing. Feeding into these are the primary signal reception clusters (coloured in red), formed around ionotropic and metabotropic receptors. Also feeding into signal processing clusters are several intermediate clusters (coloured in orange) regulate overlapping sets of pathways, and numerous clusters (coloured in green) are enriched for common downstream effector pathways.

Table 3: Algorithm characteristics applied to MASC complex.

Table 4: Algorithm characteristics applied to human PSD network.

Figure 2: Clustering of the MASC complex.
Clustering of the MASC PPI dataset [8] using the Random Walk edge Betweenness (A) and Spectral Modularity (B) algorithms. Protein communities have been colour coded in accordance with the clustering found in Figure 3 of [8]. Here red coloured communities (1 and 2) indicate proteins involved in input, blue (community 3, and 4 in (B)) information processing, orange communities (4 and 5 in (A), 5 and 6 in (B)) for sets of pathway outputs, and the remaining green communities for specific individual effector responses. The Spectral Modularity algorithm (B) captures many of the clustering results found in [8], and shown in (A). Key difference includes the separation of the central processing unit (community 3 in (A)) into two (communities 3 and 4 in (B)), diversifying communication between communities.

Figure 3: Clustering of the human PSD network.
The human PSD network is shown as a furball (A) with 291 Schizophrenia related genes, as given in Table 5, highlighted in magenta. Clustering from commonly available algorithms (Random Walk, Geodesic, Walktrap and MCL algorithms) tends towards many single (or two gene) communities, the MCL (B) algorithm is shown. Algorithms also tend towards the other extreme, generating a small number of communities of large size. These include the FC, SG, GLay and LEC (D) algorithms. Clustering results from our implementation of the Spectral Modularity (C) algorithm is shown to lie between these two extremes. Further details for each algorithms clustering characteristics can be found in Table 8.

For each cluster in the MASC and human PSD network found, we further tested the significance of enrichment for function and disease, using the Hypergeometric distribution formula for sampling without replacement:

(5)

Where in eqn-5 N is the total number of genes in the network; Cn the number of genes in the community; F the total number of functional annotated genes in the network, and f the number of functional annotated genes per community. P-values, ≤ 10^-2, were tested for their strength of significance (sig), by recording the percentage of P-values found from every community/annotation combination, lower than or equal to the observed P-value, when 1000 random permutations of the annotation labels were made. P-values found with a strength of significance <5% where consider statistically significant. P-values lower than the more stringent Bonferroni correction at the 0.05 significant levels is highlighted throughout the enrichment tables.

Performance of the suite’s boot-strap facility was also tested on the MASC and human PSD datasets. We ran each implementation on a distributed computing facility provided by ECDF (Edinburgh Compute and Data Facility, U of Edinburgh. 2013, www.ecdf.ed.ac.uk) 500 times, randomly selecting 80% of the network node size each time. The package clusterCons [20] was used in conjunction with the suite, to build a consensus matrix from which to test the robustness of discovered communities, and proteins found inside the communities. Community and protein robustness values range from 0, indicating no confidence in existing to 1, indicating absolute confidence in the cluster existing.

MASC network

We first discuss the clustering performance of the MASC network. Supplementary Table S4 gives the community enrichment values for the three algorithms relative to the functional family, sub-family, disease and phenotype annotation as found in the original MASC study [8]. The community enrichment values confirm the expected similarity between the Geodesic and Random Walk algorithms, and agree with the enrichment results reported in [8]. For example, community 2 in Supplementary Table S4 for the Geodesic algorithm and Random Walk algorithm appears to be enriched with proteins annotated as involved in G-protein coupling and in Depression. Similarly community 4, is enriched for Erk1/2 MAP kinases (Kinases P=4.7 × 10^-4 in both algorithms respectively) and also Depression (P=4.2 × 10^-2 in both respectively). These findings are supported by the high similarity between gene clusters generated by the Geodesic and Random Walk algorithms as shown in Table 3.

The Spectral Modularity algorithm also captures the important enrichment results found in both Geodesic and Random Walk algorithms. Enrichment for Kinase and Depression functional families can be found in community 5 (see Supplementary Table S4) for the Spectral Modularity algorithm, and correspondingly communities 4 in the Geodesic and Random Walk algorithms. As shown in Figure 2 this community encapsulates the well-studied MAPK/ERK signalling pathway [23]. There also exists a noticeable structural difference when clustering with the Spectral Modularity, compared to the Geodesic and Random Walk algorithms. The large community 3, observed in Figure 2a, was found to assimilate signals from various sources, including ionotropic and metabotropic signals from communities 1 and 2 [8]. This community, and the multiple signals it processes, separates as shown in Figure 2b into communities 3 and 4. With SH2 motif proteins (Src,Grb2) and PI3K/AKT pathway proteins (PI3-K,Akt2) bound to community 3, and Cytoskeletal and the well-studied Ser/Thr kinase Pka protein into community 4.

We further investigated the phenotype enrichment results with updated gene-disease annotation data collected from the OMIM [24], GeneRIF [25] and Ensembl variation [26] databases using the topOnto package, as shown in Supplementary Table S4. The annotation data was standardised using MetaMap [27] and NCBO Annotator [28,29] to recognise terms found in the Human Disease Ontology (HDO) [30]. Recognised disease ontology terms were then associated with gene identifiers and stored locally. Disease term enrichment, for a given dataset, could then be calculated using the Topology-based Elimination Fisher method [31] found in the topGO package (http:// topgo.bioinf.mpi-inf.mpg.de/), together with the standardised genedisease annotation data and the full HDO tree. As shown in Table 5, this approach allows us to compare the gene-disease association for the MASC dataset using the combined OMIM/GeneRIF/Ensembl variation data for common neurological and neurodevelopmental diseases and disorders: Schizophrenia, Alzheimer’s disease, Autistic disorder and Bipolar disorder. We further tested the enrichment results against a Neurological and Immunological clipped HDO tree following results found in [32], and a more stringent Neurological-specific tree, by clipping the full HDO tree for the “disease of mental health” and its children terms. From Table 5, we see the MASC dataset remains enriched for Schizophrenia and Alzheimer’s disease (P=1.5 × 10^-9 and P=3.1 × 10^-7) respectively. The following community enrichment values applied to the MASC network, use the combined OMIM/GeneRIF/ GeneRIF gene-disease annotation data and the neuro-specific ontology tree.

Table 5: Top disease enrichment values for MASC and the human PSD datasets.

Supplementary Table S4 shows that community 1 using the Geodesic, Random Walk and Spectral algorithms is enriched for Schizophrenia (P=5.5 × 10^-3, P=1.4 × 10^-3 and P=4.0 × 10^-3 respectively). Evidence linking Alzheimer’s disease with crosstalk between the ERK/ MAPK and PI3K/AKT signalling pathways through activation of AKT proteins (Akt1) has been observed [33]. Using the Spectral algorithm, we find further evidence supporting this observation, with enrichment for Alzheimer’s disease (P=4.3 × 10^-2) in the MAPK/ERK signalling pathway cluster, community 5 in Figure 2(b).

Typical community robustness values for the MASC network in our studies were found to range from 0.4 to 1.0, and from 0.26 to 1.0 for protein robustness. We used Cytoscape [7] to visualise the robustness information for the Geodesic (A), Random Walk (B) and Spectral Modularity (C) algorithms, as shown in Supplementary Figure S1. The robustness of each community is shown by the opacity of the community’s colour; the more opaque a community’s colour, the more confidence we have in it existing. The robustness of proteins inside the community is shown by the nodes size; bigger the node size, the more confidence we have in the protein existing within the community. In Supplementary Figure S1 we see the inter-connected nature of the enriched communities 1 and 2, and MAPK/ERK signalling pathway community 4 (community 5 for the Spectral Modularity algorithm). The large community 3 (communities 3 and 4 for the Spectral Modularity algorithm) is less robust. It is influenced by changes in the network structure, supporting its central role connecting communities in the network.

Human PSD network

Scalability of the suite’s performance to extract functional and disease enrichment was tested on the larger human PSD network. In a similar approach to that carried out on the MASC complex, we made use of the topOnto package to extract gene-disease annotation data from the OMIM, GeneRIF and Ensembl variation databases. The original study of the human PSD dataset was shown to exhibit a high density of neural disease enrichment [9]; reporting 269 diseases associated with 199 genes obtained using the OMIM database. A comparison of frequency of hits between the top 5 diseases in the original study and using the OMIM database with the topOnto package shows comparable results, see Supplementary Table S5. We extended the study to find enrichment for our diseases of interest in the PSD dataset, using the more populated OMIM/GeneRIF/Ensembl variation gene-disease annotation data and full HDO tree as shown in Table 5. We found the human PSD dataset significantly enriched for Schizophrenia and Alzheimer’s disease (P=1.0 × 10^-30 and P=1.3 × 10^-23 respectively), this remained true using the our clipped HDO ontology’s: P=6.6 × 10^-24 and P=2.7 × 10^-5 for Schizophrenia and Alzheimer’s disease respectively using the Neurological clipped tree.

For each algorithm we tested the disease enrichment per community for the top diseases shown in Table 5 (see Supplementary Table S6 for details). The community enrichment values for Schizophrenia, the most significantly enriched disease found in the human PSD network, are presented in Table 6. The P-values were calculated using eqn-5 and tested for their strength of significance. Table 4 compares the clustering performance on human PSD network using each algorithm. Where the Modularity values in Table 4 show evidence of greater clustering of the data using the Spectral Modularity algorithm (0.36), relative to the Geodesic (0.27) and Random Walk (0.18) algorithms.

Table 6: Community enrichment for Schizophrenia in the human PSD network.

Communities in Table 6, with size greater than 2, were further studied for possible enrichment of Schizophrenia related synaptic functional groups. Use was made of the 18 synaptic functional gene groups, covering 1026 pre- and post-synaptic genes, significantly associated with Schizophrenia [21]. It should be noted that the 18 functional gene groups are only used as classifications of functional type, and not group-truth clusters to test against. Three of the synaptic groups were found significantly associated with increased risk of Schizophrenia: Intracellular signal transduction (P=2.0 × 10^-4), Excitability (P=9.0 × 10^-4) and Cell Adhesion and Trans-synaptic (CAT) signalling (P=2.4 × 10^-3). As shown in Table 7 (see also Supplementary Table S6), we found evidence supporting communities enriched for CAT signalling in each algorithm: community 106 (P=1.5 × 10^-4) with the Geodesic algorithm, community 21 (P=1.3 × 10^-2) with the Random Walk algorithm, community 25 (P=2.9 × 10^-5) and community 17 (P=1.4 × 10^-8) with the Spectral Modularity algorithm. The Spectral algorithm further supported enrichment of Intracellular signal transduction in community 14 (P=1.3 × 10^-2), Structural plasticity in community 7 (P=2.8 × 10^-5) and Protein cluster in community 22 (P=3.0 × 10^-9). Community 22 shows postsynaptic adhesion molecules (including the DLG’s: DLG1/2/3/4) binding to specific presynaptic membrane proteins (CASK,LIN7A/B/C), providing evidence of a transynaptic link between the PSD and presynaptic active zone [34,35]. Our Spectral Modularity algorithm therefore (also the Geodesic algorithm through community 21 Table 7), captures tentative evidence for trans-synaptic signalling molecules involved in Schizophrenia clustered together solely driven by data.

Table 7: Enrichment values for synaptic functional groups associated with Schizophrenia for the human PSD network.

Robustness of the communities found with the Spectral algorithm applied on the human PSD network was further investigated using our boot-strap procedure applied to the MASC complex. Typical values for community robustness range from 0.13 to 0.95 and from 0.04 to 0.98 for protein robustness as shown in Supplementary Figure S2. The robustness of each community is shown by the opacity of the community’s colour, while robustness of proteins inside the community is shown by the node’s size. Schizophrenia related genes are highlighted in magenta. Communities corresponding to those shown in Table 7 are highlighted along with the community robustness values. Edges have been made opaque to allow the networks clusters to be visualised.

Comparison with existing cluster algorithms

Enrichment results for the MASC and human PSD networks were also compared against popular Modularity based clustering algorithms available in R’s igraph package [10] and Cytoscape [11], shown in Supplementary Tables S4 and S6 respectfully. Algorithms tested include implementations of the agglomerative fast-greedy community algorithms (FC and GLay) [36,37], the agglomerative Random Walk algorithm (Walktrap) [38], the coupled Potts model and simulated annealing algorithm (SG) [39,40], the Markov clustering algorithm (MCL) [41], and implementation of the Spectral Modularity algorithm (LEC) [19] from R’s igraph package [10].

The parameters used in the FC, LEC, SG and Walktrap algorithms have been chosen to maximising the Modularity; further details of each algorithms parameter sets can be found in Supplementary Table S4 and S6. We used these algorithms with their optimised settings to test against the clustering results from our implementations. We find functional cluster similarities exist between each algorithm when applied to the smaller MASC network. However, when applied to the larger human PSD network, evidence for structurally relevant communities is restricted to the Spectral Modularity algorithm with fine-tuning step. Evidence for this is summarised in Section Human PSD network, Table 8 and shown in Figure 3. The Random Walk, Geodesic, Walktrap and MCL algorithms were found to generate skewed distributions of community size; one or two large communities with many single (or two gene) communities. The skewed distribution is a known cause of the hierarchical approach to clustering and the distribution of the hub (highest degree) nodes inside the network [42,43]. The FC, LEC, SG and GLay algorithms are observed to produce a small number of large communities. The agglomerative fastgreedy community algorithms, FC and GLay, are known to produce a small number of large communities [44,45]; the unbalanced nature of the algorithms merging step, results in a few large communities growing fast, by merging in many smaller communities. The SG and LEC algorithms are divisive in nature. Clustering results of the SG algorithm was found to be sensitive to the parameter gamma, which is set to 1.0 to maximise the Modularity. However increasing gamma to 5, which favouring more edges between communities than edges inside (lowering the Modularity score), we find the SG algorithm can also lead to improved functional enrichment results. We also note differences in the clustering results between our Spectral Modularity and the LEC algorithm.

Table 8: Algorithm cluster characteristics applied to the human PSD network.

There are differences in implementation between the two Spectral Modularity algorithms. We provide here a descriptive and quantitative account for why, in our opinion, our Spectral Modularity implementation delivers functionally more relevant clusters when applied to the proteomic datasets. We do not claim to cover all differences here. Both implementations use the same divisive hierarchical method to clustering data [19]. We therefore tested whether our fine-tuning step gave closer agreement to the LEC results or not. If the two implementations were identical, with the fine-tuning step in operation or not, we would expect a NMI value of 1. We note we reproduced the Modularity values quoted in [19] for the Zachary’s ``karate club’’ network, with and without the fine-tuning step (see Supplementary Table S6). We further find agreement with the LEC algorithm results with our fine-tuning step switched off (NMI value of 0.72 compared to 0.31, see Supplementary Table S6) when applied to the human PSD network. This would suggest it is the addition of a fine-tuning step, which has a positive effect on the clustering of proteomic data and quality of functional clustering results, as highlighted in Section Human PSD network and Supplementary Tables 7 and 8. Differences in speed between the two implementations can be explained in the different eigenvector solver used. We make use of Numerical Recipes [46], which computes the full eigenvalue spectrum. The LEC implementation makes use of ARPACK [47], which calculates only a few of the leading eigenvalues, and is commonly found to perform faster than Numerical Recipes (Performance of ARPACK, Eigen, and Numerical Recipes, https://github.com/meznom/arpaca/tree/ performance/).

BioPlex network

To test how our findings translate to other large proteome datasets we looked at the BioPlex network [17]. This network is made from data obtained from a large-scale proteomic analysis of 2594 affinity purified baits that identify 7668 proteins from HEK cells. Therefore, we used a similar approach to that carried out on the MASC complex and human PSD network and extracted gene-disease annotation data for the BioPlex network as shown in Table 5 (see also Supplementary Table S8) and tested for enrichment. Unexpectedly we found evidence that the BioPlex network is enriched for proteins implicated for Alzheimer’s disease (P=3.7 × 10^-7), this remained true using our Neurological and Immunological clipped HDO tree (P=3.7 × 10^-7) and Neurological clipped tree (P=1.3 × 10^-5). We note that this result supports in part a possible neuronal origin of the HEK293 cell line [48,49].

The clustering performance of our Spectral Modularity algorithm was tested against the clustering found in the BioPlex study [17]. The BioPlex network was first clustered using clique percolation, before further sub-divided using the fast-greedy community algorithm into 345 communities. A priori we would not expect similarity between the clique-based BioPlex clustering and our Modularity-based Spectral algorithm, as both starts from different definition of what a community is. However, if we compare the BioPlex clustering to our Spectral Modularity algorithm as shown in Table 9, we find our Spectral Modularity algorithm shows evidence for cluster similarities (NMI 0.7 and ARI 0.38) with those found in the BioPlex study. The coupled Potts model and simulated annealing algorithm (SG) [39,40] also shows evidence for cluster similarity, but as we discuss in Section Comparison with existing cluster algorithms this algorithm is sensitive to tuning its parameter gamma.

Table 9: Algorithm characteristics applied to BioPlex network.

We further tested the communities found using our Spectral Modularity algorithm, and those in the BioPlex study, for enrichment of our diseases of interest (see Supplementary Table S8). As shown in Table 10, we found evidence supporting communities statistically enriched for Alzheimer’s disease using our Spectral Modularity algorithm. It was noted that by clustering the BioPlex network using the Spectral Modularity algorithm, we started to reach the known Modularity resolution limit [50]. However, we found many (110) of the BioPlex communities (greater than 2 in size) contained within our own, and by re-running our Spectral Modularity algorithm over these communities, could recover the sub-structure contained within them, including core components of the original BioPlex clusters. An example of this is shown in Supplementary Table S8, where BioPlex’s cluster 12 has been extracted from within our original community 6.

Table 10: Community enrichment for Alzheimer’s disease in the BioPlex network.

We present a suite of C++ implemented Modularity based community detection algorithms. The implementations have been optimised for speed using OpenMP and tested on varying sizes of synthetic networks. The suite is designed to run on clusters and includes a boot-strap facility, allowing a random sub-sample of the data to be selected; currently set to 80% of the network node size. Sampling from multiple bootstrap runs allows the robustness of each algorithm applied to the data to be investigated. The package clusterCons [20] has been used in conjunction with the suite, to build a consensus matrix from which to test the robustness of discovered communities, and proteins found inside each community. The package topOnto has been used in conjunction with this suite, to study the disease enrichment values of clusters obtained from published studies. We find our Spectral Modularity algorithm out performs the Geodesic and Random Walk edge Betweenness algorithms in terms of speed as shown in Table 2. It also out performs commonly available clustering algorithms in unveiling functional enrichment in proteomic datasets, as shown in Table 8 and in Supplementary Tables S4 to S6. Further investigation reveals the fine-tuning step, in conjunction with the Spectral Modularity method, is of key importance.

The datasets supporting the results of the article are included within the article (and its additional files). A C++ version of the suite is available to download at SourceForge following: http://sourceforge.net/projects/ cdmsuite/. The C++ version requires a standard gcc compiler; we tested against version 4.4.7, and has no external dependencies. A Cytoscape App for the Spectral Modularity method is also available. The App requires Cytoscape version 3.0.0 or higher, and was tested against version 3.2.1.

None declared

Implementing OpenMP versions of Geodesic and Random Walk edge Betweenness algorithms. Implementation of Spectral Modularity algorithm.

We wish to thank Theologos Kotsos, Oksana Sorokina, Andrew Pocklington, Giuseppe Gallone and Vathsala Achar with help in testing these algorithms and with general discussions.

The research leading to these results received funding from the European Union Seventh Framework Programme under grant agreement nos. HEALTH-F2-2009-241498 (“EUROSPIN”), HEALTH-F2-2009-242167 (“SynSys”), and European FET flagship project “Human Brain Project” (Subproject 1, Strategic Human Brain Data, WP1.3: T1.3.1 and T1.3.2).