Identifying communities from multiplex biological networks

Monoplex and multiplex networks

Let X be a (monoplex) network, i.e., an undirected graph. We use X for its incidence matrix too, where the entry X_i,j indicates the presence (X_i,j = 1) or the absence (X_i,j = 0) of an edge between the vertices i and j.

A multiplex network is a collection of undirected graphs (X^(g))_g sharing the same set of vertices (Kivelä et al., 2014). A basic way to apply standard graph approaches on a given multiplex network is to first aggregate all its graph layers into a monoplex network. We considered here three types of aggregation, the intersection X^∩, the union X^∪ and the sum X^Σ of the graphs in (X^(g))_g, respectively defined as: ${\displaystyle X_{i,j}^{\cap }=\underset{g}{\wedge }{X}_{i,j}^{\left(g\right)}\text{,}{X}_{i,j}^{\cup }=\underset{g}{\vee }{X}_{i,j}^{\left(g\right)}\text{and}{X}_{i,j}^{\Sigma }=\sum _g{X}_{i,j}^{\left(g\right)}}$ for all pairs of vertices i and j, where ∧ and ∨ denote the logical conjunction and the (inclusive) disjunction respectively. Note that X^∩ and X^∪ are unweighted graphs whereas X^Σ is a weighted graph.

Modularity

Under the assumption that a network X has a community structure c and by putting c_i for the community of the vertex i, the Newman–Girvan modularity (Newman & Girvan, 2004) is defined as: (1)${\displaystyle \mathcal{Q}\left(X,\mathbf{c}\right)=\frac{1}{2m}\sum _{\genfrac{}{}{0ex}{}{\left\{i,j\right\}}{i\ne j}}\left(X_{i,j}-\frac{k_i{k}_j}{2m}\right){\delta }_{{\mathbf{c}}_i,{\mathbf{c}}_j}}$ (2)${\displaystyle =\sum _a\left[\frac{\sum _{\genfrac{}{}{0ex}{}{\left\{i,j\right\}}{i\ne j,{\mathbf{c}}_j={\mathbf{c}}_j=a}}{X}_{i,j}}{2m}-\frac{\sum _{\genfrac{}{}{0ex}{}{\left\{i,j\right\}}{i\ne j,{\mathbf{c}}_j={\mathbf{c}}_j=a}}{k}_i{k}_j}{{\left(2m\right)}^2}\right]}$ where m is the total number of edges of X, k_i is the number of edges containing the vertex i (i.e., the degree of i) and ${\delta }_{{\mathbf{c}}_i,{\mathbf{c}}_j}=\left\{\begin{array}{cc}1\hfill & \text{if}{\mathbf{c}}_i={\mathbf{c}}_j\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}\right.$.

The modularity was initially conceived as a measure of the strength of a community structure (Newman & Girvan, 2004). In the original article, it is defined as the sum over all the communities a, of the proportions of within-community edges (first term of Eq. (2)) minus the expectation of these proportions in random networks with the same vertex degrees (second term of Eq. (2)).

This modularity has been extended to deal with weighted graphs, by considering the sum of the weights instead of the number of edges (Newman, 2004).

Multiplex-modularity

In order to extend the modularity of Newman & Girvan (2004) to multiplex networks, we measured the strength of a given community structure in all the graphs of a multiplex network (X^(g))_g by considering, for all communities a, the sum of the proportions of within-community edges over all the graphs minus the expectation of this sum under the same model as in Newman & Girvan (2004). Since the expectation of the sum of random variables is the sum of their expectations, we obtained the Eq. (3). Reorganizing this sum shows that it is equal to the sum of the individual modularities of the graphs of the multiplex with regard to the same community partition Eq. (5). The multiplex-modularity of the multiplex network (X^(g))_g is defined as: (3)${\displaystyle {\mathcal{Q}}^M\left({\left(X^{\left(g\right)}\right)}_g,\mathbf{c}\right)=\sum _a\left[\sum _g\frac{\sum _{\genfrac{}{}{0ex}{}{\left\{i,j\right\}}{i\ne j,{\mathbf{c}}_j={\mathbf{c}}_j=a}}{X}_{i,j}^{\left(g\right)}}{2m^g}-\sum _g\frac{\sum _{\genfrac{}{}{0ex}{}{\left\{i,j\right\}}{i\ne j,{\mathbf{c}}_j={\mathbf{c}}_j=a}}{k}_i^g{k}_j^g}{{\left(2m^g\right)}^2}\right]}$ (4)${\displaystyle =\sum _g\frac{1}{2m^g}\sum _{\genfrac{}{}{0ex}{}{\left\{i,j\right\}}{i\ne j}}\left(X_{i,j}^{\left(g\right)}-\frac{k_i^g{k}_j^g}{2m^g}\right){\delta }_{{\mathbf{c}}_i,{\mathbf{c}}_j}}$ (5)${\displaystyle =\sum _g\mathcal{Q}\left(X^{\left(g\right)},\mathbf{c}\right)}$ where m^g is the total number of edges of the graph X^(g) and $k_i^g$ is the degree of the vertex i in the graph X^(g).

The multiplex-modularity could easily be adapted to multiplex networks containing weighted edges.

Clustering algorithms for community detection

The network modularity was initially designed to measure the quality of a partition into communities (Newman & Girvan, 2004). Several clustering algorithms were subsequently developed to identify the communities by optimizing the modularity (Blondel et al., 2008; Shiokawa, Fujiwara & Onizuka, 2013). Finding the partition optimizing the modularity is NP-complete (Brandes et al., 2008) and thus requires meta-heuristic approaches when dealing with large graphs such as biological networks. We applied the Louvain algorithm (Blondel et al., 2008). The Louvain algorithm starts from the community structure that separates all vertices. Next, it tries to move each vertex from its community to another, keeping only the changes increasing the modularity, until no change increases the modularity. It then replaces the vertices by the detected communities and iterates the same operations on the new graph obtained, until stability. This algorithm was applied to the networks aggregated as X^∩, X^∪ and X^Σ. Adapting Louvain to the multiplex networks required replacing the optimization of the modularity by the optimization of the multiplex-modularity defined previously.

Resolution parameter

A general issue with modularity-based clustering approaches is that they may fail to detect communities smaller than a certain size, depending on the graph size. This phenomena is known as the resolution limit (Fortunato & Barthélemy, 2007). This issue is overcame by adding a resolution parameter γ to the modularity formula (Reichardt & Bornholdt, 2006). The parametrized modularity${\mathcal{Q}}_{\gamma }$ is defined as: (6)${\displaystyle {\mathcal{Q}}_{\gamma }\left(X,\mathbf{c}\right)=\frac{1}{2m}\sum _{\genfrac{}{}{0ex}{}{\left\{i,j\right\}}{i\ne j}}\left(X_{i,j}-\gamma \frac{k_i{k}_j}{2m}\right){\delta }_{{\mathbf{c}}_i,{\mathbf{c}}_j}.}$

The greater the resolution parameter γ, the smaller the communities maximizing the corresponding modularities. The parametrized multiplex-modularity${\mathcal{Q}}_{\gamma }^M$ of a multiplex network is defined accordingly.

Consensus clustering

Consensus clustering approaches compute a unique community structure from the community structures obtained on each graph independently. We implemented a consensus clustering approach inspired from Lancichinetti & Fortunato (2012). It first computes the community structures of each individual graph with the standard Louvain algorithm. Next, each community structure is transformed into a new graph in which an edge connects two vertices if they belong to the same community. We obtained the same number of new graphs as the initial number of graphs composing the multiplex network. These new graphs form a new multiplex network from which we infer a single community structure by using the sum-aggregation approach and the Louvain algorithm.

Random graph models

Stochastic block models (SBMs, also known as planted partition models) are among the most intuitive models of random graphs with community structures (Holland, Laskey & Leinhardt, 1983). In SBMs, the number of communities is set a priori and each vertex i is assigned to a unique community c_i. The presence of an edge between two vertices i and j is drawn independently conditionally on the respective community assignments of i and j, and with a probability that only depends on these assignments. Note that edges are independent conditionally on their communities, although not independent in an absolute sense. For instance, two edges involving vertices of the same community are more likely to be both present than edges from different communities. Hence, the parameters of the SBMs are the number of vertices, the community structure and the matrix π in which the entry π_ql is the probability of observing an edge between a vertex of the community q and one of the community l.

A natural way to generate multiplex networks (X^(g))_g with a community structure is to draw independent realizations of SBMs sharing the same community structure. This means that all the graphs are drawn independently and, in each graph, all the edges are also drawn independently, conditionally on the community assignments of the vertices. The parameters of the multiplex SBMs are the same as classical SBMs, except that instead of one matrix of edge probabilities, we have a collection (π^(g))_g, the gth matrix corresponding to the gth graph of the multiplex network. To ease the interpretation of the simulations, we parametrized edge probability matrices with only two values: the probability p_I of an internal edge and the probability p_E of an external one. All diagonal entries are equal to p_I and all non-diagonal entries are equal to p_E.

Network incompleteness and missing data

We simulated missing data by setting an arbitrary probability for a vertex to be absent in a graph, and using it to randomly remove some vertices. The community detection algorithm was then applied to the incomplete multiplex networks. In practice, whenever a vertex is missing in a graph layer, it is considered as not interacting in this graph.

Adjusted Rand index

The comparison between community partitions is performed by computing the adjusted Rand index (Santos & Embrechts, 2009), a widely used measure of the similarity between two partitions. The closer two community partitions, the greater their adjusted Rand index.

Annotation enrichment

The actual biological communities are unknown, but the accuracy of the detected communities can be assessed by verifying their consistency with known biological processes. We performed exact Fisher tests (Fisher, 1954) to search for enrichment in Gene Ontology Biological Process (GOBP) terms (Ashburner et al., 2000) associated with the genes/proteins of each community. The Bonferroni correction for multiple testing was applied on the Fisher p-values by taking into account both the number of tested ontology terms and communities. Similarly, the communities were screened for enrichment in disease genes extracted from the CTD database (Davis et al., 2014) and filtered to select the curated gene-disease associations (marker/mechanisms/therapeutic).

The intersection and union aggregations, and consensus clustering, fail to recover the community structure

We first observed that the identification of communities from the intersection of the graphs composing the multiplex network eventually leads to detect less accurate community structures (Fig. 1). An exception occurs for the simulation of dense graphs, for which we observed an improvement with the intersection of multiplex networks composed of 2 and 3 graph layers (Fig. 1C-1). In all others cases, the accuracy of the communities detected by the intersection aggregation decreases with the number of graphs composing the multiplex networks. Indeed, as the probability of observing an edge between two given vertices in all the graphs decreases with the number of graphs, the intersection tends to be empty.

A similar behavior is expected for the identification of communities from the union of the graphs composing the multiplex network, as the union aggregated network tends to be complete. This is observed for the simulation of dense and mixed-densities multiplex networks (Figs. 1C-1,2 and 1B-1,2). However, the adjusted Rand index increases with the number of graph layers on sparse simulations (Fig. 1A-1, 2). This unexpected behavior comes from the fact that these simulated networks are so sparse that their unions do not saturate, given the number of graphs considered here. In summary, both the union and intersection aggregations ultimately lead to lose the community structure.

The performance of the consensus clustering increases slowly when considering more than 2 graph layers, showing that a general community structure may eventually emerge from that detected on the individual graphs (Figs. 1A-1, 1B-1 and 1C-1). However, the accuracy of the consensus clustering never outperformed nor got close to the Multiplex-modularity accuracy in our simulations.

The multiplex-modularity better recovers communities in heterogeneous density and missing data contexts

In our simulations, the accuracy of the communities detected by the sum aggregation and the multiplex-modularity approach is consistently better than the accuracy of the union and intersection aggregations, or the consensus clustering. Importantly, the accuracy of the communities detected by the sum aggregation and the multiplex-modularity approach increases with the number of graph layers considered. Although the community structure seems undetectable from a single graph, it is almost perfectly inferred by the multiplex-modularity approach when a sufficient number of graphs is reached. For instance, communities are recovered with Rand Indexes over 0.9 with 3 graphs for sparse multiplex networks (Fig. 1A-1), 4 graphs for mixed multiplex networks (Fig. 1B-1) or 9 graphs for dense graphs and sparse multiplex networks with missing data (Figs. 1C-1 and 1A-2).

The sum aggregation and multiplex-modularity approaches show very similar performances on sparse and dense homogeneous multiplex networks (Figs. 1A and 1C). However, these two approaches do differ for mixed-densities multiplex networks and when missing data are simulated to account for network incompleteness (Fig. 1, Methods, Section ‘Network incompleteness and missing data’). In these cases, the accuracy of the multiplex-modularity approach is better than the sum aggregation. As expected, in a missing data context, the community structures are more difficult to detect. All the approaches based on the aggregation of layers composing a multiplex network are sensitive to missing data. Nonetheless, the multiplex-modularity approach still identifies adequate communities when the number of graphs is sufficiently large, and is particularly more accurate than the sum aggregation approach. Additionally, we compared our multiplex-modularity approach to GenLouvain, a general approach to detect communities from multiplex networks incorporating inter-slice connections (Mucha et al., 2010). In most simulated cases, GenLouvain has a behavior similar to the the sum aggregation approach (Fig. S1).

The identification of communities from 4 biological networks and the issue of the resolution parameter

We constructed 4 biological networks from physical and functional interaction data: a protein–protein interaction network, a network of pathways, a network of protein complex, and a mRNA co-expression network (Methods, Section ‘Biological multiplex networks’).

These networks are sparse, heterogeneous in densities and incomplete. Their densities range from 0.001 to 0.023, with an average of 0.01. They are also prone to missing data as their sizes range from 2,528 to 12,110 nodes (Table S1). Missing data is indeed an important feature in biological networks since insufficient depth and coverage of the interaction space leads to incomplete interaction datasets (Venkatesan et al., 2009; Braun et al., 2009).

Altogether, the 4 biological networks contain 17,003 nodes and 1,371,739 edges, among which 1,338,086 appear in 1 network only (Table S1). The Louvain algorithm adapted to the multiplex-modularity was applied to identify communities from this biological multiplex network. The classical Louvain algorithm was applied to the 4 networks individually, and to their union, and sum aggregations. Considering their weak performance in the simulations, the intersection-aggregation and consensus clustering approaches were not taken into account in the following analyses.

As some of the biological networks considered here are very large, we used the resolution parameter γ to reduce the sizes of the detected communities (Methods, Section ‘Resolution parameter’). We tested the clustering of the 4 individual networks, sum and union aggregations, and multiplex-modularity approach with γ parameter values ranging from 1 to 15. We next elected the parameter value γ = 5 that gives module average sizes lower than 50 proteins for the sum, union and multiplex partitions (Fig. S3). The number of modules obtained from each individual network ranges from 85 communities for the Pathway partition to 1,320 communities for the Co-expression network partition (Fig. S2). Using this value of resolution parameter, more than 800 communities are obtained from the sum and union aggregation networks, and 350 from the multiplex-modularity approach. The communities obtained with the multiplex-modularity approach display a median size of 17 proteins per community as opposed to a median sizes of 2 proteins per community for the sum or union aggregations (Fig. S4). The multiplex-modularity communities are hence more balanced and easier to interpret biologically.

The multiplex-modularity clustering identifies communities closer to the individual partitions

We computed the adjusted Rand Indexes between the community partitions obtained from the 4 individual networks and from their sum and union aggregations, and the multiplex-modularity approach (Methods, Section ‘Adjusted Rand Index’). We observed that, on average, the multiplex-modularity partition is closer to the individual network partitions than the aggregations (Fig. 2).

The multiplex-modularity identifies communities enriched in gene ontology terms

To assess the consistency of the communities, we tested their enrichment in Gene Ontology Biological Process (GOBP) terms (Ashburner et al., 2000, Methods, Section ‘Annotation enrichment’). The communities detected by the multiplex-modularity approach are more often associated to at least one significant (q-value <0.05) GOBP term than the sum or union aggregation communities (Fig. 3). For γ = 5, for instance, 38% of the multiplex-modularity communities are annotated, compared to 14% and 10% of the sum and union communities. Overall, these results suggest that the multiplex-modularity approach is extracting relevant biological communities.

Clustering of proteins involved in the Coffin-Siris syndrome and other syndromes with overlapping clinical features

As an illustration, we screened the modules obtained with the multiplex-modularity approach (γ = 5) for enrichment in disease-related proteins (Methods, Section ‘Annotation enrichment’). Twenty-nine communities are significantly associated to at least one disease. Among those, one module contains the 6 gene products known to be associated to the Coffins-Siris syndrome (q-value =5 × 10⁻¹⁰, Fig. 4, Methods, Section ‘Annotation enrichment’). This syndrome is a rare congenital disorder characterized by specific cranofacial and digital features, microcephaly and intellectual disability, among other phenotypes (Kosho & Okamoto, 2014; Vergano & Deardorff, 2014).

The module is enriched in proteins involved in transcription regulation, and participating in different cellular processes such as chromatin remodeling, transcription from RNA-polII promoter and DNA repair. The 6 Coffin-Siris syndrome proteins are components of the SWI/SNF chromatin remodeling complex (Tsurusaki et al., 2012). Interestingly, the module also contains other gene products implicated in other syndromes, which clinical features overlap with the Coffin-Siris syndrome. MCPH1 is involved in primary autosomal recessive microcephaly (Venkatesh & Suresh, 2014), BAZ1B is involved in the William-Beuren syndrome that also presents neuropsychological deficits (Fusco et al., 2014; Xiao et al., 2009) and ACTG1 mutations are involved in the Baraitser-Winter syndrome, another developmental disorder associated to intellectual disability (Rivière et al., 2012). It is to note that all these syndrome-related proteins are clustered only when the multiplex-approach is used. Indeed, when the different network layers are aggregated, only 4 out of the 6 Coffin-Siris syndrome proteins belong to the same module.