Multiscale Embedded Gene Co-expression Network Analysis — TEI-tagged view

Won-Min SongpersName, Bin ZhangpersName (2015date). “Multiscale Embedded Gene Co-expression Network Analysis.” PLOS Computational Biologytitle·j.

Object of study — MEGENAterm·mth, a planar-embedding co-expression networkterm·key framework, tested on simulated and TCGAterm data

Research question — how to build co-expression networksterm·key that avoid PMFGterm·mth's cost, false positives and clustering limits

Contribution — MEGENAterm·mth outperforms established clustering methods and co-expression networkterm·key construction approaches

Full text — tagged

The reference list is not part of the text we hold for this paper. Our copy was extracted from the publisher's page and the bibliography did not survive that extraction, so the article below ends before its references. Nothing has been removed by us.

Previously, a graph filtering technique called Planar Maximally Filtered Graphterm·mth (PMFGterm·mth) has been applied to many real-world data sets such as financial stock prices and gene expression to extract meaningful and relevant interactions. However, PMFGterm·mth is not suitable for large-scale genomic data due to several drawbacks, such as the high computation complexity O(|V| 3 ), the presence of false-positives due to the maximal planarity constraint, and the inadequacy of the clustering framework. Here, we developed a new co-expression network analysisterm·mth framework called Multiscale Embedded Gene Co-expression Network Analysisterm·mth (MEGENAterm·mth) by: i) introducing quality control of co-expression similarities, ii) parallelizing embedded network construction, and iii) developing a novel clustering technique to identify multi-scale clustering structures in Planar Filtered Networksterm·mth (PFNsterm·mth). We applied MEGENAterm·mth to a series of simulated data and the gene expression data in breast carcinoma and lung adenocarcinoma from The Cancer Genome Atlasterm (TCGAterm). MEGENAterm·mth showed improved performance over well-established clustering methods and co-expression networkterm·key construction approaches. MEGENAterm·mth revealed not only meaningful multi-scale organizations of co-expressed gene clusters but also novel targets in breast carcinoma and lung adenocarcinoma.

Introduction Often, complex diseases involve multiple intertwined signaling circuitries. Cancer is an excellent example with a number of biological machineries activated in tumor pathogenesis including proliferation, angiogenesis, avoidance of cell death, evasion of tumor suppressing mechanisms, immortality, invasion etc[ 1 ]. The complexity of cancer further manifests via “tumor microenvironment”, a concept that incorporates interactions between not only the tumor cells, but also normal cells that contribute to the expression of the cancer hallmarks[ 2 ].

In many cases, networks of these intertwined signaling cascades, such as protein-protein interaction networks and metabolic networks are highly heterogeneous[ 3 – 5 ]. Particularly, these networks share certain characteristics such as the scale-free propertyterm·dom (the degree distribution follows a power law), small world effectterm·dom (diameter of network scales with logarithm/double-logarithm of the number of nodes)[ 3 , 5 ], assortativityterm·dom (preference for a network’s nodes to attach to others that are similar in some ways, i.e., high degree nodes tend to attach to high/low degree nodes)[ 6 ], and community structuresterm·dom [ 7 , 8 ]. These observations suggest that the biological networks may follow the similar evolutionary dynamics, and thus network analysis approaches from other domains are very helpful for understanding biological networks[ 4 ].

These organizational principles are reflected in transcriptional control of cells: highly modular and yet diverse functional patterns emerge by means of “co-expression” [ 9 , 10 ]. Co-expressed gene clusters represent coherent unique functional pathways not only in normal conditions [ 10 , 11 ], but also in disease states[ 9 , 12 – 14 ]. These “guilt-by-association” approaches were further extended to encapsulate gene-gene interactions by regarding genes as nodes and interactions as links, known as “co-expression network analysisterm·mth”. These methods first evaluate the association strength between each gene pair by a similarity score (e.g., Pearson’s correlation coefficient) or statistical significance of the association, then identify co-expressed clusters or communities in the context of network topology[ 9 , 15 , 16 ].

However, the existing techniques to construct co-expression networksterm·key suffer from a number of drawbacks. For instance, some popular co-expression networksterm·key such as those from Weighted Gene Co-expression Network Analysisterm·mth (WGCNAterm·mth) enforce the connectivity to exhibit a power-law distribution[ 16 ], unweighted networks by hard thresholds contain a large number of false positive interactions[ 17 ], k-nearest-neighbor networks require the number of neighbors to connect by subjective criteria such as connectedness[ 15 ], and partial correlation based co-expression networksterm·key require at least O(|V| 3 ) computational complexity[ 18 ], limiting the practical applications to |V| < 10 4 . These are further complicated by clustering analysis to identify modular organization of these networks. Some widely used clustering methods such as k-means and spectral clustering require predefined number of clusters[ 19 ].

More importantly, many of network-theoretic clustering methods are incapable of different levels of aggregations of clusters co-existing within a single network. There are several factors accounting for this particular drawback. Firstly, Newman’s modularityterm·key measure suffers the inherent resolution limitterm·dom that fails to differentiate certain configurations of obvious clusters[ 20 ]. Secondly, they are often restricted to identify a single partition of a network by optimizing for the modularityterm·key, thus overlook multiscale organization of complex networks where coarse-grained and compact clusters co-exist[ 8 ].

In order to account for these shortcomings, we adopted a network embedding paradigm on a topological sphere. In other words, a co-expression networkterm·key is embedded on a spherical surface such that one link does not cross the others. Planar Maximally Filtered Graphterm·mth (PMFGterm·mth) was developed to extract most relevant information from similarity matrices based on topological sphere, and has been applied mostly in financial domain[ 21 ]. PMFGterm·mth becomes an ideal platform to construct co-expression networksterm·key due to the following attractive features: i) the preservation of hierarchy by retaining Minimum Spanning Treeterm·mth (MSTterm·mth) as a subgraph, ii) the correspondence between a coherent cluster (if any) and a connected subnetwork, iii) the abundance of 3- and 4-cliques and exhibition of rich clustering structures[ 21 ], and iv) the possession of a wide spectrum of fundamental network characteristics in embedded networks such as transitions between scale-freeterm·key to exponential degree distributions, and large-world to semi-ultra-small world[ 22 , 23 ]. Applications to financial data have revealed that characteristic features of complex systems such as emergence of bubbles[ 24 , 25 ], aggregation of similar firms in same sectors[ 24 ], highly connected hubs and hierarchical organizations[ 26 , 27 ]. Furthermore, an embedded network inference framework called “Directed Bubble Hierarchical Treeterm·mth” (DBHTterm·mth)[ 27 ] was developed to infer meaningful clustering and hierarchical structures in PMFGsterm·mth from gene expression, financial, and simulated data [ 27 ].

However, the existing PMFGterm·mth embedding technique cannot efficiently handle large-scale genomic data. Firstly, pair-wise similarities are noisy and redundant, yielding high false-positive rates in identifying gene-gene interactions[ 17 ]. Enforcing maximal planarity inevitably introduces a significant number of these redundant links in a filtered network and may obscure the underlying “true” interactions. Secondly, the computation complexity for testing planarity is too high ( O (| V | γ ), 2 ≤ γ ≤ 3) for large scale network analysis. Thirdly, clustering analysis in PMFGterm·mth via DBHTterm·mth framework is not optimal. DBHTterm·mth framework is based on inference of the patterns between separating triangles in PMFGterm·mth, and requires that every node belongs to at least one triangle. Noting that gene-gene interactions do not necessarily form triangles, DBHTterm·mth framework may not assign these genes to appropriate clusters. Lastly, a rigorous and formal definition of multiscale organization in these networks has been ignored. Although hierarchical structures have been exploited via agglomerative hierarchical clustering within bubbles and bubble clusters, they are inherently limited by aforementioned drawbacks of the bubble topology, and require a more rigorous algorithm to extract the full information encoded in embedded networks.

Here we developed a new network construction and analysis framework named Multiscale Embedded Gene Co-expression Network Analysisterm·mth (MEGENAterm·mth) to resolve the aforementioned issues with PMFGterm·mth and DBHTterm·mth, and more broadly with the existing co-expression network analysisterm·mth methods. In the rest of the paper, we will briefly overview MEGENAterm·mth, and then perform a comprehensive performance comparison of MEGENAterm·mth and the established network construction and clustering analysis approaches using a series of simulated data as well as the real-world large-scale gene expression data. Finally, we will address the advantages of MEGENAterm·mth and highlight some novel insights derived by MEGENAterm·mth.

Results In this section, we first give an overview of MEGENAterm·mth and address its algorithmic scalability. Then, we perform a comprehensive evaluation of MEGENA-derived Planar Filtered Networksterm·mth (PFNsterm·mth) and MEGENA-derived multi-scale clusters. We evaluate interactions captured by PFNsterm·mth using a series of simulated and real-world gene expression data in comparison with False Discovery Rate (FDRterm·mth) based networks (FDRNterm·mth). We further compare the functional relevance of MEGENA-derived multi-scale clusters and those identified by several established clustering analysis approaches. Simulated expression data as well as the gene expression data in breast carcinoma (BRCA) and lung adenocarcinoma (LUAD) from The Cancer Genome Atlasterm (TCGAterm) were used. A detailed description about the data acquisition and preprocessing can be found in S1 Text . We proceed to highlight some novel insights revealed by multi-scale clusters and then conclude the paper by summarizing the key contributions and pointing out future work.

Overview of MEGENAterm·mth MEGENAterm·mth consists of four major steps: 1) Fast Planar Filtered Network constructionterm·mth (FPFNCterm·mth) by introducing parallelization, early termination and prior quality control; 2) Multiscale Clustering Analysisterm·mth (MCAterm·mth) by introducing compactness of modular structures characterized by a resolution parameter; 3) Multiscale Hub Analysisterm·mth (MHAterm·mth) to identify highly connected hubs of each cluster at each scale and 4) Cluster-Trait Association Analysisterm·mth (CTAterm·mth) to explore the relevance of cluster to clinical outcomes. Fig 1 shows the overall analysis flow of MEGENAterm·mth. Below we give a brief description of FPFNCterm·mth, MCAterm·mth and MHAterm·mth. The details about these steps are presented in Methods .

Flow chart of MEGENAterm·mth. A) Fast planar filtered network constructionterm·mth. Significant interactions are first identified and then embedded on topological surface via a parallelized screening procedureterm·mth described in the text. On the right, a toy example is illustrated to show construction of PFNterm·mth from a thresholded network by FDRterm·mth (top left), and gradual construction of PFNterm·mth with number of included links and screened pairs shown on the top of each. B) Multi-scale clustering: Beginning from connected components of the initial PFNterm·mth as the parent clusters, clustering is performed for each parent cluster and compactness of the sub-clusters are evaluated. These steps are described in the dotted box. The clustering is performed iteratively until there remains no further parent clusters meaningful to split. C) Downstream analyses: Multiscale Hub Analysisterm·mth (MHAterm·mth) is performed to detect significant hubs of individual clusters and across α, characterizing different scales of organizations in PFNterm·mth. Then, clusters are ranked by associations to clinical traits including enrichment of differentially expressed gene (DEG) signatures, and correlations to survival end-point etc.

FPFNCterm·mth constructs PFNterm·mth by mostly following the network embedding rationale from the PMFGterm·mth algorithm. All pairs of genes are first ranked via a similarity measure quantifying respective interaction strengths and then iteratively tested for planarity to grow the embedded network that favors inclusion of pairs with larger similarities[ 21 ]. To make the PFNterm·mth construction scalable for whole genome co-expression network analysisterm·mth, two techniques were developed. Firstly, insignificant interactions are removed before the network embedding step by controlling the False Discovery Rate (FDRterm·mth) of similarity for each gene pair. However, such a filtering may not be necessary since we will show in the subsequent section of Evaluation of PFNsterm·mth that PFNsterm·mth are very robust with respect to different FDRterm·mth thresholds. Secondly, a parallelized screening procedureterm·mth (PCPterm·mth) is developed to extract a subset of gene pairs that are more likely to be embedded. Such procedures enable FPFNCterm·mth to efficiently and effectively construct embedded co-expression networksterm·key by capturing significant interactions at the whole genome level.

PFNterm·mth constructed through FPFNCterm·mth is then input to MCAterm·mth to identify multiscale clusters. MCAterm·mth incorporates three distinct criteria to identify locally coherent clusters while maintaining a globally optimal partition. First, shortest path distancesterm·mth (SPDterm·mth)[ 28 ] are utilized to optimize within-cluster compactness. Second, local path indexterm·mth (LPIterm·mth) is used to optimize local clustering structure. Third, overall modularityterm·key (Qterm·mth)[ 29 ] is employed to identify optimal partition. Specifically, MCAterm·mth adopts a hierarchical divisive approach to dissect complex interactions in PFNterm·mth into coherent interactomes across different resolution scales by iterating two steps, k-split and compactness evaluation. k-split identifies the clusters that lead to an optimal partition of a parent network via optimization of SPDterm·mth, LPIterm·mth and Qterm·mth. In the step of compactness evaluation, individual clusters from k-split are compared to the parent network via a measure of network compactness defined below, υ = S P D ¯ log ( | V | ) α (1) where, V is the set of nodes in the network, S P D ¯ is the average of shortest path distancesterm·mth of all node pairs, and α is the resolution parameter. Given that the denominator log(| V |) α is the hallmark of the small-worldterm·key property represented by the scaling relation S P D ^ ∼ log ( | V | ) when α = 1, ν measures the coherence of a network’s topology. Therefore, a smaller α identifies more compact clusters. For a given cluster (network), MCAterm·mth searches through a range of α values for a resolution scale that leads to more compact clusters than the parent cluster (network). These clusters are further split by k-split until no more compact clusters can be identified. Each split represents a finer picture of modular structure of the given PFNterm·mth. The output of MCAterm·mth is a hierarchy of clusters at various levels defined by α.

Finally, MHAterm·mth and CTAterm·mth constitute the downstream analyses in MEGENAterm·mth. MHAterm·mth first identifies significant hubs within each cluster with respect to an established random model of planar networks[ 23 , 29 – 31 ]. The nodes that are hubs at multiple scales are called multiscale hubsterm·dom. CTAterm·mth evaluates the relevance of individual clusters to clinical outcomes through principal component and correlation analyses.

Scalability of MEGENAterm·mth: Effectiveness of PCPterm·mth PCPterm·mth is a key technique developed to speed up PFNterm·mth construction to overcome the worst case O(|V| 3 ) complexity of the existing serial PMFGterm·mth algorithm (See Methods for detailed discussion) [ 21 , 27 ]. In conjunction with correlation screening, PCP-mediated FPFNCterm·mth dramatically increases its efficiency in construction whole genome co-expression networkterm·key. In order to verify this, we compared PCP-mediated network embedding and the existing serial PMFGterm·mth using the TCGAterm gene expression data that involve over 20,000 genes. We compared the acceptance rateterm·var of pairs filtered by PCPterm·mth in MEGENAterm·mth and that of non-filtered pairs by PMFGterm·mth. The acceptance rateterm·var, defined as |E|/|E| max , where |E| is the number of edges embedded in a PFNterm·mth, and |E| max = 3(|V| - 2) is the maximal number of edges embeddable in a planar network by Euler relation [ 32 ].

As shown in Fig 2 , the acceptance rateterm·var by the serial PMFGterm·mth algorithm quickly decreases close to 0% as the number of links in PFNterm·mth reaches the maximal number of links. The finding indicates that PMFGterm·mth performs exponentially increasing number of computations to embed more edges as the number of links in PFNterm·mth saturates towards the maximal number. On the contrary, PCPterm·mth remedies the problem by dramatically boosting the acceptance rateterm·var close to 100% as the number of links in PFNterm·mth increases. These results demonstrate the effectiveness of PCPterm·mth in reducing the overall computation time by leveraging parallel computation capability, and scalability of FPFNCterm·mth for whole-genome co-expression networkterm·key.

Comparison of acceptance ratesterm·var of correlation pairs into PFNterm·mth links. A,B) Results from PFNterm·mth construction from TCGAterm lung squamous cell carcinoma (LUSC) data including 20523 genes. 57562 links out of maximal possible link number of 61563 are embedded. The left panel (A) shows the acceptance ratesterm·var without PCPterm·mth (denoted as “serial”, and colored as blue), and after performing PCPterm·mth (denoted as “PCPterm·mth”, and colored as red), as a function of number of links already embedded on the PFNterm·mth, normalized by the maximum possible number of embedded links. The right panel (B) shows the ratio of acceptance ratesterm·var after PCPterm·mth to the acceptance ratesterm·var without PCPterm·mth is plotted as a function number of links already embedded on the PFNterm·mth, normalized by the maximum possible number of embedded links. C,D) Results from TCGAterm thyroid carcinoma (THCA) data including 16639 genes. 44802 out of maximal possible link number of 49911 are embedded. The right and left panel show the same plots as described in the case of LUSC.

Evaluation of PFNsterm·mth We evaluated the performance of PFNsterm·mth from multiple aspects. We first evaluated capacity of PFNsterm·mth in capturing underlying regulatory interactions by comparing to golden standard networks using simulated datasets from DREAM challenge[ 33 ]. We then compared the network neighborhoods of a number of genes in PFNsterm·mth with their actual targets derived from the perturbation experiments. Furthermore, we compared the global topological properties of inferred PFNsterm·mth with the established hallmark signatures of complex networks.

Evaluation by simulated data To evaluate the accuracy of PFNterm·mth in capturing the underlying gene regulatory network, we leveraged GeneNetWeaver[ 34 ] to generate 10 time series datasets for each of the golden standard networks from DREAM challenge[ 33 ]. The generation of the simulated data was detailed in S1 Text . We systematically compared MEGENA-derived PFNsterm·mth with those based on the state-of-art methods including ARACNE[ 35 ] and Random Forest (RF)[ 36 ] across various FDRterm·mth thresholds and similarity/dissimilarity measures. Specifically, we tested Pearson’s correlation coefficient (denoted as PCC ), Mutual Information (denoted as MI ), and Euclidean distance (denoted as euclid ), and applied FDRterm·mth thresholds (0.01, 0.05, 0.1, 0.2, 0.4, 0.6, 0.8 and 1) on each similarity/dissimilarity measure to filter out insignificant interactions. FDRterm·mth thresholded interactions were then used to construct PFNsterm·mth and ARACNE networks. Note that the three measures are all applicable to PFNterm·mth, but ARACNE uses only MI as it assumes Data Processing Inequality (DPI) imposed by MI[ 35 ]. As RF is not dependent on any similarity/dissimilarity measure, FDRterm·mth of the interactions determined by RF is estimated by permuting input data.

For each inferred network from a simulated data, we compared the weighted shortest path distancesterm·mth for all pairs of the nodes in the network and those in the underlying gold standard network. The node pairs connected in a gold standard network were treated as a positive class and the pairs not connected were taken as a negative class. We then calculated Area Under Curve in Receiver Operating Characteristic curve (AUC-ROC) for shortest path lengths in each inferred network.

As shown in Fig 3A , the PFNsterm·mth from PCC and MI consistently outperform the RF and ARACNE networks across various FDRterm·mth thresholds. Table 1 shows the best average AUC-ROC scores, indicating that PFNsterm·mth from PCC and MI across various FDRterm·mth thresholds show consistently the best performance except for InSilicoSize100-Yeast2 data set where PFNsterm·mth are only slightly outperformed by RF based networks at an FDRterm·mth threshold of 1. At FDRterm·mth thresholds of 0.2 or less are practically used in almost all cases, PFNsterm·mth from PCC and MI show the best overall performance.

Validation of PFNsterm·mth in comparison to various network inference methods. A. Comparisons of AUC of ROC for weighted shortest path distancesterm·mth of inferred networks from simulated data from various golden standard networks (labeled on the top), in comparison to ARACNE and RF. Different combinations with Pearson’s correlation coefficient (Pearson), mutual information (MI) and Euclidean distance (Euclid) were tested. B-C. Comparison of BRCA TF knock down signatures on BRCA PFNterm·mth (red) and FDRNterm·mth (green) neighborhoods of the target TFs, inferred from MI. The strips on the top of each plot shows expression fold changes (1.3 and 1.5 respectively) to derive these signatures. B shows FDRterm·mth corrected FET p-values against the number of significantly enriched signatures. C shows enrichment fold change cut-off against the number of significantly enriched signatures. D-E. Comparisons of BRCA TF knock down signatures on inferred networks from PCC. D and E correspond to FDRterm·mth corrected FET p-values and enrichment fold changes, similarly to B and C.

Table of best average AUC-ROC across various FDRterm·mth thresholds. Each column represents the combination of network inference method and similarity/dissimilarity measure tested, and each row represents gold standard networks from which time series were generated. The best performing methods are highlighted by bold font.

To assess the stability of the performance of each method, we calculated Coefficient of Variation (CV) of the average AUC-ROC scores across different FDRterm·mth thresholds. PFNsterm·mth from different similarity/dissimilarity measures have the most stable performance among the tested networks across different FDRterm·mth thresholds ranging from 0 to 1 ( S1 Fig ) . The performance of PFNsterm·mth peaks at FDRterm·mth thresholds around 0.01 and/or 0.05 in most cases. However, PFNsterm·mth from Euclidean distance have relatively poor performance in general.

Evaluation by disease-specific data We further evaluated the accuracy of PFNsterm·mth in detecting the true interactions in a disease dataset, we constructed PFNterm·mth from BRCA gene expression data set from TCGAterm (hence, BRCA PFNterm·mth). The details about the data acquisition and preprocessing can be found in S1 Text . For comparison with the BRCA PFNterm·mth, we also constructed FDRterm·mth thresholded network (FDRNterm·mth) with FDRterm·mth < 0.05 for PCC and MI. These co-expression networksterm·key were then tested for the enrichment of the siRNA knockdown signatures of key transcription factors (TFs) of breast carcinoma in MCF7 cells [ 37 ]. Of 78 TFs in the knockdown experiments, 32 TFs remained after data processing (see Data Acquisition and Processing in S1 Text ), and appeared in the BRCA PFNterm·mth and their siRNA signatures were used for testing the networks. Namely, the 32 TFs are: BCL2 , BRCA1 , BRCA2 , CCL5 , CCNA2 , CCNB1 , CDC20 , CDC25A , CDC25B , CDKN2A , CEBPB , CEBPD , CENPE , CENPF , CHEK1 , E2F1 , E2F5 , ERBB2 , ESR1 , FOS , FOXC1 , GATA3 , HIF1A , HOXB7 , ID1 , MYBL2 , MYC , PAX3 , SKP2 , STAT1 , TOP2A , and WT1 . Two differentially expressed gene (DEG) signatures for each experiment were identified based on T-test p value < 0.05 and fold changes ≥ 1.3 or 1.5. Two different fold change cut-offs were chosen to give a more comprehensive performance.

For each co-expression networkterm·key, we tested the enrichment of a TF’s DEG signature in its l -layer network neighborhoods ( l = 1,…, l max , where l max is the largest shortest distance between the given TF and the other nodes in the network). The optimal layer l optimal was chosen based on FET p-value that shows the most significant enrichment of the DEG signature. The final FET p-values were adjusted for multiple testing. In the case of MI, Fig 3B and 3C show the number of signatures enriched in the BRCA PFNterm·mth and FDRNterm·mth at various thresholds for the corrected FET p-values and enrichment fold changes, respectively. Fig 3D and 3E show the results for the PCC based coexpression networks. Clearly, these experiment-derived gene signatures are far more significantly enriched in the PFNsterm·mth than in the FDRNsterm·mth, implying that the PFNsterm·mth constructed by MEGENAterm·mth can better capture true gene regulatory relationships using either PCC or MI.

Evaluation by topological characteristics Degree distributions and diameters of PFNsterm·mth were computed since these are the key network topological characteristics representing the hallmark features such as scale-freeterm·key and small-worldterm·key networks [ 28 ]. We present the PFNsterm·mth constructed from PCC in this section. Fig 4 shows the global PCC-based BRCA PFNterm·mth and its clusters identified at a scale defined by α = 1.3, and S2 Fig shows the global PCC-based LUAD and its clusters at α = 1. We reserve “k” to denote the node degree only in this section, but use k to denote the number of clusters in the other sections.

The global BRCA PFNterm·mth. Different node colors represent different clusters identified at a scale of α = 1.3. Node size and label size are proportional to node degree.

As shown in Fig 5A , the BRCA PFNterm·mth is scale-freeterm·key by following a typical power-law degree distribution with exponent γ ≤ 3, consistent with the frequently observed range of exponents, i.e., 2 ≤ γ ≤ 3, in real-world complex networks [ 28 ]. However, the LUAD PFNterm·mth does not exhibit the characteristics of scale-freeterm·key degree distribution across all k though the distribution between 3 ≤ k ≤ 50 is scalefreeterm·key ( Fig 5B ). The decaying tail at k ≥ 50 in the LUAD PFNterm·mth shows the characteristics of exponential distributions. The diameters of the BRCA and LUAD PFNsterm·mth are ~11.3, which is in accordance to the hallmark feature of small-worldterm·key networks with diameters around log(|V|).

Degree distributions of the BRCA PFNterm·mth (A) and the LUAD PFNterm·mth (B). The x-axis is the logarithm of degree k and the y-axis is the logarithm of inverse cumulative degree distribution, P(k’ > k). Red straight line is fitted distribution for P(k^'>k)~k^(γ+1), where γ is the estimated exponent of the underlying degree distribution. Respective γ value is displayed at the top.

In summary, these PFNsterm·mth possess the hallmark features of complex networks. The qualitative difference between the degree distributions of the two PFNsterm·mth demonstrate a broad range of network characteristics from scale-freeterm·key to exponential distributions[ 23 ]. These results support the breadth of topological diversity in embedded networks such as PFNsterm·mth. Several studies of statistical mechanics of embedded networks on topological sphere showed that those networks can possess either exponential degree distribution [ 31 ], or scale-freeterm·key degree distribution with various exponents[ 23 ] in thermodynamic limits, depending on the underlying evolutionary dynamics.

Evaluation of Multicale Clusters We further compared the clusters derived from MEGENAterm·mth and those identified by other established clustering and network inference approaches using the TCGAterm BRCA and LUAD gene expression data (see Data Acquisition and Preprocessing in S1 Text for description of BRCA and LUAD data). Specifically, we considered two other types of coexpression networks including weighted co-expression networksterm·key (WGCN) and unweighted coexpression networks (FDRNterm·mth, based on the links at FDRterm·mth < 0.05)[ 14 , 16 ] (see S1 Text for details) and three established clustering techniques including infomap[ 38 ], walktrap[ 39 ], and leading eigenvector based spectral clustering[ 40 ]. Note that these clustering methods detect coherent clusters in complex network by optimizing for Newman’s modularityterm·key Qterm·mth. Two different similarity measures, MI and PCC, were used for constructing coexpression networks.

Weighted co-expression network analysisterm·mth (WGCNAterm·mth) uses its own clustering method which is not suitable for un-weighted networks like PFNsterm·mth and FDRNsterm·mth. Towards this end, we compared MEGENAterm·mth (as a combination of PFNterm·mth and MCAterm·mth), WGCNAterm·mth and the following 6 combinations of networks (PFNterm·mth and FDRNterm·mth) and clustering methods (infomap, walktrap, leading eigenvector), PFNterm·mth + infomap, PFNterm·mth + walktrap, PFNterm·mth + leading eigenvector, FDRNterm·mth + infomap, FDRNterm·mth + walktrap, and FDRNterm·mth + leading eigenvector. Furthermore, both MI and PCC were used to construct PFNsterm·mth, WGCNs and FDRNsterm·mth.

As there are a few oncogenic signatures available in LUAD, the evaluation of LUAD networks is less comprehensive than that of the BRCA networks. Therefore, we focus on the results from the BRCA data in the main text and report the results from the LUAD data in S1 Text .

Functional analysis of multiscale clusters We collected a large number of gene sets associated with known molecular functions and pathways from Molecular Signature DataBase (MSigDB) across GO-BP (Gene Ontology–Biological Processes), GO-CC (Gene Ontology–Cellular Components), GO-MF (Gene Ontology–Molecular Functions), KEGG (Kyoto Encyclopedia of Genes and Genomes) and REACTOME (Reactome database) categories. The significance of the overlap between a given cluster and an annotated gene set was calculated by the Fisher Exact Test (FET), and corrected for multiple testing via Bonferroni correction for the total number of comparisons (i.e. number of clusters Х number of gene sets tested). The detailed results from the enrichment analysis were included in S1 Data .

As shown in Fig 6A and 6B , the multiscale clustering analysisterm·mth (MCAterm·mth) has the best performance since the resulting clusters are enriched for the largest number of the annotated gene sets with respect to all significance levels in both MI- and PCC-based networks. More importantly, the MCA-derived clusters show the largest fold enrichment of the BRCA oncogenic signatures. Similar results are observed in PCC-based LUAD networks (see S3A and S3B Fig ). MCAterm·mth consistently outperforms the established co-expression network analysisterm·mth method, WGCNAterm·mth in both the BRCA and LUAD cases.

Comparison of MEGENAterm·mth (as a combination of the multiscale clustering analysisterm·mth and PFNterm·mth) and various combinations of the established clustering techniques (eigenvector, infomap, walktrap, WGCNAterm·mth) and the networks (PFNterm·mth, FDRNterm·mth, WGCN) using the TCGAterm BRCA gene expression data. Two different similarity measures (MI and PCC) were used to perform analyses to compare robustness with respect to difference in measures to evaluate interactions. A) The number of significantly enriched functional/pathway signatures (Bonferroni corrected FET p-values) from MSigDB at various p-value thresholds against. B) Number of significantly enriched functional/pathway signatures from MSigDB at the various odds ratio thresholds. C) Number of clusters predictive of patient survival (based on FDRterm·mth corrected Cox p-values) at various significance levels. D) Number of clusters predictive of patient survival (based on FDRterm·mth corrected Cox p-values) and associated to at least one significantly under-represented signatures with Bonferroni corrected FET p-value < 0.05.

Interestingly, the methods that directly optimize for Qterm·mth (infomap and eigenvector) consistently show better performance on the FDRNsterm·mth than the PFNsterm·mth. Qterm·mth assumes that the random networks follow the configurational modelterm·mth, where edges are shuffled while maintaining the degree sequence [ 29 ]. In other words, the underlying random model of Qterm·mth assumes that the chance of connecting two nodes by an edge is equal for all pairs to 2|E|/(|V|(|V|-1)) , regardless of the degrees of connected nodes. This implies that there is no correlation between the degrees of the nodes sharing an edge in the random network model. However, this assumption does not hold for geometrical networks such as PFNsterm·mth. We have previously shown that, these random geometrical networks do possess significant degree correlations [ 23 ], therefore optimization of Qterm·mth alone cannot properly address the optimal partition of PFNterm·mth. On the other hand, FDRNterm·mth is relatively free from these constraints, and therefore the resulting clusters by optimizing for Qterm·mth in FDRNterm·mth reflect the underlying biology better than the case of PFNterm·mth. On the other hand, the walktrap methodterm·mth resolves this problem by leveraging local random walks instead of Qterm·mth. Particularly, the walktrap methodterm·mth on the PFNsterm·mth outperforms the infomap and eigenvector methods on the FDRNsterm·mth. Altogether, the results imply that the modular structures in PFNsterm·mth can be better identified by clustering methods that capture local clustering structure, supporting the use of LPIterm·mth in MEGENAterm·mth, as described in Methods .

Prognostic analysis of multiscale clusters We then examined how each cluster identified from various approaches is associated with the overall survival based on principal component analysis (for details, see Section 4 Cluster-Trait Association Analysisterm·mth in Methods ). As shown in Fig 6C , MCAterm·mth (multiscale) identifies the largest number of significant clusters across a wide range of FDRterm·mth corrected Cox p-value thresholds, and this is commonly observed in both MI- and PCC-based networks. Although there are relatively more clusters from MCAterm·mth than other clustering methods due to the hierarchical divisive nature, several clusters from MCAterm·mth are most predictive of survival. This outstanding performance is also observed in the case of PCC-based LUAD data set (see S3C Fig ).

We further identified the clusters that were not enriched for any functions/pathway signatures based on Bonferroni corrected FET p-value < 0.05 and the odds ratio > 1. As these clusters have unknown functions, they are termed as unknown-functional clusters (UNCterm·dom). The UNCsterm·dom from various approaches are listed in S2 Data . We then checked how many of the UNCsterm·dom identified by each approach were predictive of survival at various thresholds for FDRterm·mth corrected Cox p-values. As shown in Fig 6D , MCAterm·mth identified the most number of UNCsterm·dom that were predictive of survival. Such an unambiguous trend was observed in both PCC- and MI-based networks.

In summary, MEGENAterm·mth as a combination of PFNterm·mth and MCAterm·mth identifies the largest number of functional clusters across wide range of cluster sizes in both BRCA and LUAD datasets (see S4 Fig ). In particular, MEGENAterm·mth outperforms WGCNAterm·mth in terms of more significant enrichment for known pathways and better predictive power of survival.

Highlight of an adipocytokine-enriched cluster in BRCA To highlight the findings by MEGENAterm·mth, we identified 9 functionally annotated clusters (FACs) detected only by MEGENAterm·mth. These clusters are significantly enriched for some GO-BP/KEGG/REACTOME gene sets that are not enriched in any clusters detected by any other aforementioned approaches based on a threshold of 0.05 for the multiple-test (Bonferroni) corrected p-values. The 9 FACs and their enriched gene sets are shown in S1 Table .

Among these MEGENA-specific FACs, one cluster (comp1_56) is enriched for the genes in adipocytokine signaling pathway (corrected FET p = 0.02, 2.7 fold). The hub genesterm·key of this cluster are all significantly associated with the overall survival of the BRCA Luminal-B patients. Luminal-B, a molecular subtype of hormone-receptor positive breast cancers, is associated with higher grade and increased proliferation rate, and has a poorer overall prognosis than its hormone-receptor positive counterpart, Luminal-A [ 41 ] . Fig 7A and 7B show the localization of genes with univariate Cox p-value < 0.05 for the Luminal B patients’ overall survival at the adipocytokine-enriched cluster. The hub genesterm·key of this cluster are all predictive of the Luminal B patients’ overall survival: AQP7 (Cox p-value < 1.5e-3), C14orf180 (Cox p-value < 4.4e-3), CIDEC (Cox p-value < 1.6e-3), CIDEA (Cox p-value < 2.1e-2) and MRAP (Cox p-value < 2.2e-2). We further examined the significance of the survival difference between expression median defined subgroups for each hub geneterm·key. Fig 8 shows the Kaplan-Meier plots for the two most predictive hub genesterm·key, AQP7 and CIDEC .

Identification of the adipocytokine-enriched cluster, comp1_56, which was specifically identified by MEGENAterm·mth. A) The Global BRCA PFNterm·mth. The nodes in red represent the genes that is predictive of overall survival of LumB patients (Cox p-value <0.05). The blue circle indicates the location of the cluster comp1_56. B) A magnified view of the cluster comp1_56. The nodes with labels are the hubs of the cluster.

Kaplan-Meier plots of subgroups separated by median expressions of two hub genesterm·key AQP7 (A) and CIDEC (B), showing significant logrank p-values. Blue curves showing lower risks correspond to lower expressions, and red curves showing higher risks correspond to higher expressions.

These genes were significantly repressed in the tumor samples in comparison to the matched normal samples for Luminal-B patients: CIDEC (fold change (FC) = 2.5E-2, T-test p-value (p) = 3.5E-4), CIDEA (FC = 4.9E-2, p = 1.4E-3), AQP7 (FC = 4.7E-2, p = 7.1E-4), MRAP (FC = 4.3E-2, p = 1.2E-3) and C14orf180 (FC = 4.7E-2, p = 5E-4). This is in agreement with the previous findings that AQP7 expression was down-regulated in breast tumor[ 42 ] and CIDEA and CIDEC are cell death-inducing DFFA-like effectors to activate apoptosis[ 43 ]. Particularly, CIDEA and CIDEC are involved in adipose tissue loss in cancer cachexia[ 44 ], an important, negative prognostic marker that has been linked to systemic inflammation and cell death[ 45 ].

In summary, these findings suggest that adipocytokine signaling pathways may play an important role in Luminal B subtype of breast cancer though the exact mechanisms need further experimental validation. Interestingly, the expression of the hub genesterm·key in the adipocytokine cluster/subnetwork were not significantly associated with the overall survival outcome of the Luminal A patients. Therefore, MEGENAterm·mth is capable of capturing finer-scale functional subnetworks to stratify a breast cancer subtype into the subgroups with prognostic significance.

Multiscale Organizations in PFNsterm·mth In this section, we will explore the multiscale clusteringterm·key structures in PFNsterm·mth constructed by MEGENAterm·mth. Here, we mainly focus on the PCC-based BRCA PFNterm·mth as it showed slightly better performance than the MI-based network ( Fig 3B–3E ).

Multiscale organization of functions/pathways We performed MHAterm·mth on the BRCA PFNterm·mth to identify the groups of scales that had similar interaction patterns and shared highly connected hubs across different scales. Six distinctive scale groups were identified: S1 (0.03 ≤ α ≤ 0.48), S2 (0.5 ≤ α ≤ 0.82), S3 (0.87 ≤ α ≤ 1), S4 (1.01 ≤ α ≤ 1.29), S5 (1.3 ≤ α ≤ 1.82) and S6 (1.83 ≤ α ≤ 6.8). Biological relevance of each scale group was evaluated by the number of significantly enriched MSigDB gene sets. We compared the performance of the clusters at each scale group and that of the clusters across all scale groups. Fig 9A shows that the combination of all the clusters across the different scale groups consistently outperforms the individual scale groups across almost the entire range of significance levels. Interestingly, the clusters at the scale S5 (1.30 ≤ α ≤ 1.82) show the best performance when compared against other scale groups. The clusters identified at the finest scale of S5 (α = 1.3) are shown in Fig 4 .

Hierarchical organization of functions and signaling pathways corresponding to the multiscale clusters identified by MEGENAterm·mth. A) Comparison of number of significantly enriched functions and pathway signatures across clusters identified at different scale groups. The scale groups identified from MHAterm·mth are colored according to the legend, and “all” denotes collection of clusters across the scale groups. B) Multiscale organization of clusters in PFNterm·mth. Each node is a cluster identified by multiscale clusteringterm·key in PFNterm·mth, where the node size is proportional to the cluster size, node color coincides with the cluster group color scheme in A, and node labels indicate most enriched function/signaling pathway for individual clusters. A directed link a→b indicates b is a sub-cluster of a.

To further understand the biology of the clusters at the different scales, we annotate each cluster with its most enriched function or pathway and then build up a cluster hierarchy based on the parent-child relationships (e.g., a child cluster is a subset of the parent cluster). As shown in Fig 9B , the cluster hierarchy clearly displays biological relevance of the multiscale/multisresolution organization patterns in the BRCA PFNterm·mth, by showing hierarchical aggregation of specific functions and pathways to more general terms. For instance, the cluster annotated as SKELETAL_DEVELOPMENT (FET p-value<2.7e-08, 2.48 fold enrichment) is the parent of the clusters with more specific functional categories associated with development/differentiation: HEMOSTASIS (P<2.49e-05, 5.35 fold), ADIPOCYTE_DIFFERENTIATION (P<3.95e-08, 19.42 fold), SYNAPTOGENESIS (P<4.32e-08, 7.18 fold), and MUSCLE_CONTRACTION (P<1.503e-08, 15.22 fold).

In the case of the PCC-based LUAD PFNterm·mth, there are twelve distinct scale groups (0.2 ≤ α ≤ 3.2) (see S5A and S5B Fig ). Again we observed the superior performance of the clusters across all scale groups over any individual scale (see S5C Fig ).

Identification of novel key drivers via MHAterm·mth Having validated that some known interactions are present in the BRCA PFNterm·mth via the oncogenic signatures, we hypothesize that hub genesterm·key of the BRCA PFNterm·mth are more likely to be key drivers of BRCA etiology, and further explore biological significance of novel key drivers of the network.

To verify biological significance of hub genesterm·key detected by MHAterm·mth, we compared expression fold changes of network hubs between different cancer stages, to those of the non-hub genes. We first identified the hub genesterm·key at each scale, and then intersected the hub geneterm·key sets at the different scales to identify a more stringent hub set, denoted as “multiscale hubsterm·dom”. We then evaluated the significance of the difference between the fold change distributions from each hub set and the corresponding non-hub genes using the Kolmogorov-Smirnov (KS) test. Fig 10 compares the distributions of expression fold changes in the two groups with respect to different stages of breast cancer. S6 Fig shows p-values from the KS test.

Comparison of expression fold changes (FC) of the hub genesterm·key and non-hub genes between different cancer stages in BRCA, against lists of genes identified by mutiscale hub analysis, where fc denotes expression fold change. The numeric labels on x-axis represent the ranges of α values defining the resolution levels of the hubs, “multiscale” represents intersection of hub genesterm·key across different scales, and “non.hub” represents the rest of genes.

In many cases, the hubs show significantly higher expression fold changes than the non-hub genes, and the average fold changes are greater for hubs at larger scales. Particularly, the multiscale hubsterm·dom show the highest average fold changes, suggesting that they may be transcriptomic drivers of breast cancer progression. There are 14 multiscale hubterm·dom genes including ROPN1 , TPX2 , TEKT1 , FOXA1 , ESR1 , CCNT1 , SDPR , THSD4 , MLPH , TBC1D9 , FOXC1 , SPDEF , SFRP1 , and AR . Many of these genes are known to play important roles in breast cancer etiology. For instance, AR and ESR1 are well-established endocrine receptors in breast cancer [ 46 , 47 ]. FOXA1 is a pioneer transcription cofactor for ESR1 , opening chromatin containing ESR1 target genes[ 48 ]. In ER- breast cancers, FOXA1 may promote androgen signaling and may contribute to the development of resistance to anti-androgen therapy [ 46 ]. SPDEF , a target gene of ESR1 , overexpressed in breast and other solid tumors, was shown to be associated with worse outcomes in patients with ER+ breast cancers and to be critical for the survival of ER+ breast cancer cells in vitro [ 49 ]. SPDEF was upregulated in ER + breast cancer cell models of estrogen-deprivation resistance and tamoxifen resistance[ 50 ]. TPX2 is an established regulator in spindle assembly and DNA damage response across many solid tumors[ 51 ], and indeed is a key regulator of the cell-cycle cluster identified by MEGENAterm·mth. FOXC1 , a prognostic biomarker of basal-like breast cancer patients[ 52 ] is a key regulator of NF - kappaB signaling pathways in basal-like breast cancer cells[ 53 ] and promotes breast cancer invasion[ 54 ]. SFRP1 is a known inhibitor of Wnt pathway and tumor suppressor gene, which is epigenetically silenced in a variety of tumors including breast cancer[ 55 ].

This multiscale hubterm·dom set also includes a number of novel genes as promising targets of breast cancer for further studies. One example is ROPN1 , which is an immediate neighborhood of FOXC1 and FOXA1 . ROPN1 , also known as ropporin , is a cancer-testis antigen[ 56 ] and a potential immune-therapeutic target for multiple myeloma as Chiriva-Internati et al . generated human leukocyte antigen class I-restricted cytotoxic lymphocytes to kill autologous multiple myeloma cells[ 57 ]. ROPN1 is significantly associated with the overall survival for all the BRCA patients in TCGAterm (Cox p-value < 2.6e-2, logrank p-value < 5.8 e-2), PR+ (Cox p-value < 9.7e-5, logrank p-value < 3.0e-4), and ER+ (Cox p-value < 3.3e-4, logrank p-value < 2.7e-4). Higher expression of ROPN1 was associated to better prognosis. Furthermore, ROPN1 is mostly down-regulated in tumor samples in comparison to normal samples (FDRterm·mth corrected test p-value < 8.4e-14, the ratio of the average expression in the tumor samples to that in the adjacent normal samples = 0.14). Fig 11 shows a significant survival difference between the ROPN1 -low and -high groups for all, ER+ and PR+ patients. Moreover, ROPN1 is predictive of survival of the Her2 subtype identified by PAM50 biomarkers (Cox p-value < 2.2e-3) [ 58 ]. These results suggest ROPN1 as a desirable immunotherapeutic target for further validation.

Kaplan-Meier plots of the subgroups defined by median expression of ROPN1 in A) all the patients, B) the ER+ patients and C) the PR+ patients. Blue and red curves correspond to the lower and higher expression levels of ROPN1 , respectively.

We also performed MHAterm·mth on the PFNterm·mth clusters based in the TCGAterm LUAD data set to verify the findings in the BRCA data set shown in S7 and S8 Figs . We again observed a similar pattern to the BRCA PFNterm·mth, i.e., the multiscale hubsterm·dom of the LUAD PFNterm·mth have significant expression fold changes. TPX2 , C16orf89 and GJB3 emerged as the multiscale hubsterm·dom present at all scale groups in the LUAD PFNterm·mth. Notably, GPR116 and HOPX directly connect with C16orf89 in the LUAD PFNterm·mth and are hubs at several scales. In particular, HOPX (HOP homeobox) is a lineage-specific transcriptional regulator of differentiation and has been shown to control the fate of LUAD progression where the cooperative expression of HOPX with GATA6 limits metastatic competence of LUAD cells[ 59 ]. GPR116 , an essential regulator of lung surfactant homeostasis, was shown to play a crucial role in preventing alveolar collapse through its ability to reduce surface tension[ 60 , 61 ]. We further showed that the DEG signatures from HOPX and GATA6 double-knockdown in the lung adenocarcinoma cell lines[ 59 ] and GPR116 knockdown in the murine type II alveolar epithelial cells from the Gpr116 knockout mice ( GSE41417 )[ 62 ] were significantly enriched in HOPX and Gpr116’s neighborhoods, respectively, (see S9 Fig ). Altogether, these results suggest C16orf89 as a novel therapeutic target in preventing lineage-specific metastasis in LUAD for murine type II alveolar epithelial tissue.

Discussion We developed a novel framework, Multiscale Embedded Gene Co-expression Network Analysisterm·mth (MEGENAterm·mth), to infer gene co-expression networksterm·key, by implementing a parallelized algorithm for embedding co-expression networksterm·key on topological sphere and a new clustering analysis algorithm to detect coherent clusters at various compactness scales. MEGENAterm·mth constructs a co-expression networkterm·key by enforcing an objective criterion of “embeddabilityterm·mth” of candidate connections on topological sphere, and thus reduces the inherent redundancy of pair-wise interactions. MEGENA-derived networks possess the hallmarks of complex networks [ 21 , 23 , 31 ]. Application of MEGENAterm·mth and the state-of-the-art network inference approaches to the simulated data with gold standard networks shows that the MEGENA-derived PFNsterm·mth have the best and most stable performance. The outstanding performance of MEGENAterm·mth was further demonstrated in the significant overlap between the siRNA knock-down signatures of a large number of functionally important TFs and their inferred network signatures.

We showed that the novel multiscale clusteringterm·key technique in MEGENAterm·mth can identify biologically more meaningful and relevant coexpressed gene clusters than established network clustering methods such as infomap, walktrap, leading eigenvector spectral clustering, and WGCNAterm·mth.

We highlighted the novel insights from the multiscale approach to decipher gene-gene interaction networks. Key drivers/hubs of MEGENA-derived networks were further identified by MHAterm·mth, a key procedure in MEGENAterm·mth. We identified ROPN1 and C16orf89 as the novel candidate drivers for BRCA and LUAD, respectively.

As an alternative approach to analyzing big Omics data, MEGENAterm·mth has demonstrated competitive performances and will have a great potential in unraveling novel pathways and key regulators in complex diseases. There are rooms to further improve MEGENAterm·mth. Currently, we are extending the network construction algorithm to higher genus to account for more complex interaction patterns as hyperbolic surfaces with higher genus are able to accommodate for more complexity [ 31 ].

Software Availability The MEGENAterm·mth R package for Windows can be downloaded from here: http://research.mssm.edu/multiscalenetwork/packages/MEGENAterm·mth_1.1.zip The MEGENAterm·mth R package for Linux can be downloaded from here: http://research.mssm.edu/multiscalenetwork/packages/MEGENAterm·mth_1.1.tar.gz Methods Fast Planar Filtered Network Constructionterm·mth A key component of MEGENAterm·mth is the construction of Planar Filtered Networksterm·mth (PFNsterm·mth). Here we developed a new procedure named FPFNCterm·mth to substantially improve the existing PMFGterm·mth in terms of efficiency and scalability. Specifically, we first introduced a parallelization process for testing planarity, and then implemented early termination options to construct ‘nearly maximal’ embedded networks which prevent inclusion of less informative but computationally expensive links. The procedure for constructing a PFNterm·mth is detailed below.

Compute similarities between gene expression profiles Given a gene expression dataset with N genes and M samples, we first compute the similarity between any two genes. A number of similarity measurements such as correlation, Euclidean distance and mutual information can be employed to compute the similarity between expression profiles. MEGENAterm·mth can take as input similarities from any similarity measurement. Comparison of various similarity measurements is beyond the scope of this paper. Gene-gene similarities are then filtered by False Discovery Rate (FDRterm·mth) to minimize the impact of false positives. FDRterm·mth is computed by permuting gene expression matrix across the samples (global FDRterm·mth), or by directly calculating pairwise nominal p-values by Fisher’s Z-transformation. In the current implementation of FPFNCterm·mth, we set FDRterm·mth < 0.05 as the default threshold for filtering similarities.

Construct PFNterm·mth The existing PMFGterm·mth algorithm embeds an input network onto a topological sphere by the following steps [ 21 ] ( Fig 12 ) : It begins from the empty network G o ( V o , E o ) where nodes are completely disconnected, i.e. E o = ∅.

Then, it rank-orders gene pairs by their co-expression similarities, and iteratively test planarity of each pair ij by Boyer-Myrvold algorithm.

If a pair ij passes the planarity test, this results to update G o by adding the pair as a link in the network, i.e. E f = E o ∪ { ij }. This is equivalent to embedding ij onto the sphere.

i) ~ iii) is repeated in a serial manner until the maximal number of edges |E max | = 3(|V o |-2) is reached, and results in G f (V f , E f ) where V o = V f and |E f | = |E max | .

Fast PFNterm·mth construction. A parallelized screening procedureterm·mth is developed to extract a subset of gene pairs which are highly likely to be embedded. A) FPFNCterm·mth begins with a rank-ordered list of association pairs. B) Then a subset of N c pairs undergo parallelized quality control by their embeddabilityterm·mth on a single platform of G o to identify the pairs which are more likely embedded in the subsequent network construction steps. C) These screened set of N ’ c pairs are then tested on the growing embedded network subsequently. D) A final updated network G ’ , which will be used as G o on the next cycle. The whole processes are repeated until the defined criterion for termination is met.

The resulting network G f is an embedded network on a topological sphere where every edge can be drawn without crossing another.

To speed up the planarity test, we designed a parallelized screening procedureterm·mth (PCPterm·mth) to extract a subset of gene pairs that are more likely to be embedded. The underlying premise of PCPterm·mth is explained as follows. If G ′( V ′, E ′) with E ′ = E o ∪ { ij } is embeddable, then G ″( V ′, E ″) with ( E ″ − { ij }) ⊆ E o is embeddable. Intuitively speaking, we are utilizing the fact that any subnetwork of a planar network is always planar. Therefore, if ij is one of edges included in the final network G f , the combination of any subnetwork in G f or G o with the edge ij must be embeddable since the combination is also a subnetwork of a planar network.

This allows to use G o as a platform to test the quality of N c individual pairs prior to the link embedding step (the step (iii) in PMFGterm·mth) in PCPterm·mth where each of N c pairs are combined with G o and tested for planarity in parallel. The qualified pairs after PCPterm·mth proceed to the serial embedding step in PMFGterm·mth, finally yielding to the updated network G f . This process is also shown in Fig 12 .

Currently, we implemented N c = 1000 x N core where N core is the number of cores used in parallel computation. Furthermore, we have set PCPterm·mth to start when the acceptance rateterm·var of each edge into G f falls below 10%, where acceptance rateterm·var is defined as the number of accepted pairs by the number of tested pairs. We recommend that users follow the current parameters as we show that PCPterm·mth with current parameters improves the overall computation time effectively as discussed in Section 1 in Results.

Terminate PFNterm·mth construction Early termination conditions are set up to further bypass unnecessary computations to embed less informative pairs that are not filtered out by Step 1.1 while keeping sufficient information for network construction. By setting up reasonable termination options, the resulting PFNterm·mth still harbors sufficient information. The construction process will be terminated if one of the following conditions is satisfied: (i) a network is maximally embedded (identical to PMFGterm·mth when Step 1.1 is bypassed)[ 21 ]; (ii) an embedded network reaches a certain saturation where the mathematically permitted maximum number of links is 3(|V|-2) ; (iii) the number of rejected pairs per one embedded link reaches a certain threshold. For the condition (iii), we set the threshold for the number of rejected links as 20|V| , corresponding to FDRterm·mth < 0.05. Empirically, we found out that the condition (ii) with 90~95% of the maximum number of links or the condition (iii) with the default number of rejected links can provide sufficient results with minimal information loss.

Multiscale Clustering Analysisterm·mth (MCAterm·mth) on PFNterm·mth As most biological networks exhibit highly modular and yet hierarchical organizations[ 3 , 7 ], we next set about to identify coherent modular structures in PFNsterm·mth. It has been well-known that characterization of the organization patterns in complex networks cannot be done by a single perspective, but requires a combination of multiple distinctive and diverse features[ 8 ].

Measures of compactness and local clusteredness We perform multiscale clustering analysisterm·mth of PFNsterm·mth by optimizing a number of key network features including within-cluster compactness, local clustering structures, overall modularityterm·key, and other network-theoretic measures.

Specifically, within-cluster compactness is measured via the shortest path distancesterm·mth (SPDterm·mth) on PFNterm·mth[ 28 ], local clustering structure via Local Path Indexterm·mth (LPIterm·mth)[ 63 ], and overall modularityterm·key via Qterm·mth[ 29 ]. SPDterm·mth is defined as distances of shortest geodesic paths among all nodes in a given network[ 64 ], thus can be naturally used to represent the degree of compactness for a given set of nodes by summarizing overall distances between the nodes on the network. However, SPDterm·mth tends to over-emphasize hierarchical organization of nodes, and thus ignores clustering structures[ 65 ].

In order to mitigate this problem, we chose LPIterm·mth to incorporate local clustering structures in PFNterm·mth. LPIterm·mth is defined as A 2 + ϵA 3 , where A is the binary adjacency matrix of the network, ϵ is a free parameter, and [A n ] ij is the number of paths of n steps linking nodes i and j [ 63 ]. LPIterm·mth is a quasi-local measure that evaluates structural similarity between two nodes on a network by accounting for all paths of lengths 2 and 3 between any two nodes. LPIterm·mth with ϵ = 0.01 in the context of link prediction has been shown to perform superior to some of already established network theoretic similarity measures in many real-world cases[ 63 , 66 ]. Particularly, the abundance of 3-cliques and 4-cliques in PFNterm·mth works synergistically with LPIterm·mth for the following reasons: a) a majority of the 2- and 3-step paths are confined within these cliques, and b) organization patterns of adjacent 3-cliques and 4-cliques provide rich clustering information of PFNterm·mth[ 22 , 23 , 27 ], and c) the adjacency information can be readily captured by LPIterm·mth by paths that span between different cliques. Furthermore, the computation time for LPIterm·mth is linear with the size of network[ 63 ], making it suitable for large-scale co-expression network analysisterm·mth. Lastly, we leveraged the capability of Qterm·mth to evaluate significance of within-cluster connectivity in tandem with between-cluster connectivity to guide proper splitting of networks into sub-clusters.

We introduced a network compactness measure ν(α) as a function of a resolution parameter α, to capture the cohesiveness of modular structures in PFNterm·mth at different resolutions (or scales) defined by α. Particularly, we leveraged geometrical characteristics of PFNterm·mth that retains similarity between two nodes by means of geodesic path distance such that coherent clusters exist as connected subgraphs with higher local clustering coefficients[ 21 ]. Furthermore, these geometrical networks exhibit a broad spectrum of compactness characteristics from loose-world to semi-ultra small-worldterm·key[ 23 ], and thus they allow different degrees of compactness to co-exist. Fig 1B presents an overview of the procedure in MCAterm·mth and the details of MCAterm·mth are described in the following subsections.

Network split: k-split At each iteration, a nested split is performed on each cluster to derive an optimal split with respect to the initial cluster via k-medoids clustering which detects k optimal clusters by minimizing the SPDterm·mth based intra-cluster distance[ 67 ]. However, SPDterm·mth based evaluation of the partitions at each k does not take into account clusteredness of local topology, leading to misclassification of the nodes at the boundaries between adjacent clusters. To remedy this problem, we designed a procedure to update the boundaries of the clusters derived from a partition at each k. Boundary nodes are those with immediate neighbors from at least two different clusters. Let V bnd be a set of boundary nodes. Then, for a node i ϵ V bnd , its cluster membership to V l is defined via the following equation, μ V l i = max V l ′ { μ V l ′ i | μ V l ′ i = ∑ j ∈ V l ′ L P I i j } (2) where membership of i is assigned to V l .

The Eq ( 1 ) assigns a given node i to a new cluster in which the node i has the maximal number of interactions as determined by LPIterm·mth. This boundary detection procedure (BDP) is iterated until there is no change in cluster membership, therefore leading to a stable partition. This process is illustrated in Fig 13 in the step “update boundary”, showing a toy example where misclassified red nodes in “Before” panel are correctly assigned in the “After” panel after BDP.

Flow chart of the clustering analysis procedure for each value of compactness resolution parameter, α. The upper panel illustrates the k-split procedure within each cluster to detect optimal sub-clusters. The lower panel describes the compactness evaluation procedure (CEP) after k-split. CEP compares the parent cluster prior to k-split with the sub-clusters after k-split by means of the compactness measure, ν l , and updates the partition accordingly. On the left, each step is illustrated by a graphical toy example. From the top, the pictures correspond to: the initial network subject to clustering, correct classification of boundary nodes by BDP (Before: before BDP, After: correction after BDP), identification of the optimal k via modularityterm·key Qterm·mth k , final clusters, and comparison between initial network and sub-clusters via compactness. These steps are iterated for all clusters from the newly updated partition until no further update can be made.

As network split requires determination of the optimal k value, we utilized an established measure of clusteredness on networks, Newman’s modularityterm·key, Qterm·mth[ 29 ]. Specifically, we repeat the k-medoids partition and BDP for a range of k values until no further optimal solution by Qterm·mth is available in the interval [k’—dk, k’], where dk is set as 10. In other words, we search for the optimal partition determined by Qterm·mth within 10 further partitions, and repeat the search until no better solutions are detected. Although we do not have the data for systematic evaluation on setting the appropriate dk value, we did not observe much differences in resulting clusters in the exemplified BRCA and LUAD PFNsterm·mth for dk ≥ 10, and dk = 10 often succeeds to explore up to k ≈ 60 in these cases. Given that the clustering is a hierarchical divisive method and 2 ≤ dk ≤ 6 are often the number of clusters explored for each split in hierarchical divisive methods [ 15 , 40 ], we recommend dk = 10 is sufficient to effectively identify these clusters.

Identification of significantly compact sub-clusters Upon splitting a network G o ( V o , E o ) into clusters, each cluster is evaluated by the compactness measure defined as the average shortest path distanceterm·mth within a cluster V l , normalized by the logarithm of the cluster size, υ l = ∑ i , j ∈ V l S P D i j log ( | V l | ) α ( | V l | ( | V l | − 1 ) / 2 ) = S P D ¯ log ( | V l | ) α , (3) where α is a scaling parameter to control the degree of compactness. Normalization by log(|V|) α is in accordance to the scaling effect from the diameter, S P D ^ , of geometrical networks such as PFNsterm·mth. In the thermodynamic limit of |V| → ∞, geometrical networks can be modeled by means of a stochastic network model where small clusters of “bubbles” aggregate to give characteristic features of real-world complex networks, i.e., S P D ^ ∼ log ⁡ ( | V | ) α , where α ≥ 1, a hallmark feature of small-worldterm·key property[ 23 ]. However, S P D ^ does not reflect network structural information since it relies on the minimum of the longest pair-wise shortest path lengths. To mitigate this, we utilize S P D ¯ to incorporate all information between all pairs of nodes, and take the scaling effect from S P D ^ as an upper bound of network compactness.

For a given network subject to split, a split is obtained by maximizing for Qterm·mth that is independent of α. Then, at a given α, the following criterion determines the acceptance or rejection of the split: { Accept : if ν o > υ l for any l , Reject : if ν o < υ l for all l , (4) where, ν o is the compactness of the parent cluster. The criterion implies that, for a given α , a split that results in at least one cluster l with improved ν l in comparison with G o is considered as valid.

In order to efficiently search for valid splits across all α, we leveraged the independence between the split detection via Qterm·mth, and the split acceptance criterion. For a given network G o , we define a characteristic α value as the one that can identify a more compact cluster l than G o : α l c ≔ max { α l | ν o > S P D ¯ log ( | V l | ) α l with α o c ≥ α l } (5) Calculation of α l c for all the clusters from a split, allows us to identify α values at which the split is accepted according to the criterion (4).

Individual sub-clusters from an accepted split are evaluated by calculating the statistical significance of ν l at a scale α = α l c . Specifically, a random PFNterm·mth is first generated by shuffling the link weights of the parent cluster. Then, |V l | nodes are sampled for 100 times from the random parent, and the compactness values for these random networks, ν l ′ , are calculated to estimate significance p-value of ν l . Random node insertion moves (namely T2 moves) are utilized to generate random PFNsterm·mth which exhibit scale-freeterm·key degree distributions and small-worldness in thermodynamic limits [ 23 ].

Termination of algorithm The search for sub-clusters is performed iteratively until no valid sub-clusters can be further identified under the following conditions: - no sub-clusters can be further identified as more compact than respective parent clusters at any α, or - no further sub-clusters show significant compactness (P>0.05).

Multiscale Hub Analysisterm·mth (MHAterm·mth) We developed a Multiscale Hub Analysisterm·mth (MHAterm·mth) procedure to identify highly connected nodes at each scale defined by α and across all the scales. MHAterm·mth identifies the nodes with significantly high connectivity within each significant clusters previously identified through the following steps: 1) Group the scales that show similar within-cluster connectivity patterns, 2) Identify hubs at each scale , and 3) identify multiscale hubsterm·dom by combining significance scores of individual nodes across all different scales. The procedure is detailed in the following subsection.

Grouping similar scales At each scale α, we define within-cluster connectivity of node v i as c w ( v i , α ) = ∑ v j ∈ V l α A ( v i , v j ) (6) where V l α is the set of nodes in the cluster l at the scale α , and A ( ν i , ν j ) = w is the adjacency matrix denoting the weight of link connecting v i and v j . Combining c w ( ν i , α ) across all nodes and α values, we obtain C w ( V ,A) where V is the set of all nodes in PFNterm·mth, and Φ is a set of α values that have at least one significant split.

Then, k-medoids clustering is performed for each α ∈ Φ. The optimal number of clusters, k, is then estimated by summarizing results from multiple established internal validity indices since a single validity index may emphasize only a specific criterion and may not provide good quality clusters [ 68 ]. Specifically, we utilized a number of internal validity indices such as average silhouette width, Normalized Gamma statistics, Dunn’s index and separation index [ 68 , 69 ] to evaluate quality of clustering results across k ∈ [2,|A| − 1].

We then summarized the results by calculating combined normalized ranks across these cluster quality indices by the following formula, score ( k = k ′ ) = ∑ m ′ ∈ M log ( rank ( k = k ′ , m = m ′ ) ( | Φ | − 1 ) ) (7) where M is the set of cluster quality indices, rank( k , m ) is the rank of k by the cluster quality index m , and score( k ) is the summarized score of k .

Fig 14 shows the grouping process in analyzing the BRCA PFNterm·mth where 6 distinct scales were identified.

Identification of hubs at various scale (defined by α) groups in the breast cancer PFNterm·mth. A) Plots of various internal validity indices used for selecting the optimal number of clusters to group α values. B) Barplot showing summarized scores from normalized ranks by internal validity indices from A). C) A heatmap of the pairwise Euclidean distances between any two vectors of the within-cluster connectivity (determined by C w (V , A)) of all the nodes at the corresponding scales. The color bar on the top of heatmap represents the distinct scale clusters identified by MHAterm·mth.

Identification of hubs at each scale For each subnetwork at a scale α, n s random networks are generated as described previously in MCAterm·mth description. The significance of within-cluster connectivity c w ( ν i , α ) is then evaluated by taking the average number of nodes that have higher connectivity than c w ( ν i , α ) from n s random networks with n s = 100 as default.

Identification of multiscale hubsterm·dom We adopt Fisher’s inverse Chi-square approach to compute the combined statistics for each node i across α values grouped together, termed a scale group. That is, S i = ∑ α ∈ A l − log 10 p i ( α ) (8) where, p i ( α ) is significance p-value of c w ( ν i , α ), and A l is the scale group including the set of α values grouped together by clustering for C w ( V ,A). In order to evaluate the significance of S i , we generate the null distribution of S i by randomly shuffling for α and i in the equation above. This shuffling is performed N s times to generate stable null distribution, where N s = 100 by default. Nominal p-values for each of S i is calculated from the null distribution, and is corrected for multiple testing by Bonferroni correction for the number of nodes in the PFNterm·mth. We have set Bonferroni corrected p-value < 0.05 as the default threshold to identify significant hubs for each scale group A l .

Cluster-Trait Association Analysisterm·mth (CTAterm·mth) To relate each cluster with clinical outcomes, principal component analysis (PCA) is first performed for each cluster and then the correlation between the first (or multiple) principal component(s) and each trait is computed as cluster relevance to the trait.

For patient survival data, the association is examined by multivariate Cox proportional hazards regression model that regresses patient survival onto the first (or multiple) principal component(s) of a given module, and Cox p-value is calculated to evaluate the significance. To further investigate the prognostic power of each cluster, logrank p-value is calculated to characterize the difference between the survival curves of two molecular subtypes defined by the median expression of the first PC of each cluster. The logrank p-values and Cox p-values are then corrected for multiple testing by Benjamini–Hochberg FDRterm·mth correction.

Computational Complexity Analysis Among the four major steps of MEGENAterm·mth including FPFNCterm·mth, MCAterm·mth, MHAterm·mth and CTAterm·mth, FPFNCterm·mth is most time consuming. The complexity of the existing serial PMFGterm·mth algorithm has a complexity of O(|V| γ ), 2 ≤ γ ≤ 3, where the worst case of O(|V| 3 ) is due to performing O(|V|) Myrvold-Boyer planarity test on O(|V| 2 ) correlation pairs. FPFNCterm·mth circumvents this problem by reducing the number of correlation pairs subject to the planarity test by means of testing significance of every correlation pair, taking O(|V| 2 ). Assuming that a certain threshold such as FDRterm·mth < 0.05 leaves a faction of nodes correlation to every node, we can approximate the number of remaining pairs to construct PFNterm·mth as ϵ|V|. Combining these two, the overall complexity is O(ϵ'|V| 2 ), which is a substantial improvement over the previous algorithm. Furthermore, the parallelization via PCPterm·mth with several cores allows to handle for the multiplicative factor ϵ', leading to O(ϵ''|V| 2 ) with ϵ' > ϵ'' ≥ 1. As a result, FPFNCterm·mth achieves a scalable computation of PFNterm·mth of |V| ~ 20,000 with moderate computational resources within a few days, while the exiting serial PMFGterm·mth algorithm takes over a week to handle a network with |V| ~ 5000. For instance, using FPFNCterm·mth on 16 cores (3.5 GHz Intel Ivy Bridge), the LUSC PFNterm·mth with |V| = 20523 ( Fig 2A and 2B ) was constructed in less than 36 hours and the THCA PFNterm·mth with |V| = 16639 ( Fig 2C and 2D ) took less than 18 hours. Construction of such large-scale PFNsterm·mth is not feasible for the existing serial PMFGterm·mth algorithm as we estimate that it will take over a month.

MCAterm·mth is governed by computation of shortest-path distance (SPDterm·mth) for all pairs of nodes, and iterative k-medoids clustering. We adopted the Bellman-Ford algorithm to compute SPDterm·mth[ 70 ] which has a computation complexity O(|V||E|). Given |E| ≤ 3(|V|-2) in embedded networks on surface with g = 0, the time complexity of computing SPDterm·mth is O(|V| 2 ). SPDterm·mth is calculated for multiple times in MCAterm·mth. It is calculated first from the global PFNterm·mth as the input dissimilarity matrix for k-medoids clustering, and then from multiple random planar networks to calculate statistics for cluster compactness to evaluate each split. Therefore, the initial computation of SPDterm·mth from PFNterm·mth dominates the running time since computing SPDsterm·mth for candidate clusters become relatively negligible due to dramatic decrease in cluster size. Therefore, the overall complexity involving all SPDterm·mth calculations becomes O(ϵ|V| 2 ), where ϵ corresponds to the number of clusters with sizes comparable to the PFNterm·mth. Additionally, the computational complexity of k-medoids clustering is O(|V| 2 /k)[ 67 ], where k is tested from k = 2,…,k max with k max reaching around 50 in practical cases with current implementation of k-split. Therefore, the overall time complexity of MCAterm·mth is O(ϵ'|V| 2 ). In the current implementation, the overall computation time of MCAterm·mth for a PFNterm·mth with |V| = 15402 took less than 2 hours on a single core (3.5 GHz Intel Ivy Bridge).

Lastly, the computational complexity of MHAterm·mth is dictated by calculation of significance of within-cluster connectivity for each node, across a range of α values. Given that we generate n s (= 100 by default) random planar networks for each unique cluster, and we calculate degree of each node for each random network, the time complexity for performing the statistical test for each cluster is O(|V l |), and for all clusters is O(Σ l |V l |). Since Σ l |V l | ~ ϵ|V| where ϵ is the mean number of instances that a single node appears in different clusters, the overall time complexity for MHAterm·mth is fairly linear with |V|. Indeed, MHAterm·mth for the BRCA and LUAD PFNsterm·mth in this manuscript took only few minutes.

Overall, the computation complexity of MEGENAterm·mth is O(β|V| 2 ), where β largely depends on the number of cores to perform parallelized computations. The space (memory) complexity of MEGENAterm·mth is O(|V| 2 ) due to a |V|x|V| similarity matrix. Based upon 16 cores of 3.5 GHz Intel Ivy Bridge, FPFNCterm·mth just needed less than 3 hours to construct the BRCA and LUAD PFNsterm·mth with 6999 and 7562 nodes, respectively while the existing PMFGterm·mth algorithm took over a week. The whole MEGENAterm·mth took less than 4 hours for both cases.

Supporting Information (TXT) Click here for additional data file.

(TXT) Click here for additional data file.

(DOCX) Click here for additional data file.

(DOCX) Click here for additional data file.

The error bars show standard deviation across 10 simulated data sets from a golden standard network from DREAM challenge.

(TIF) Click here for additional data file.

Node border colors represent different clusters identified at a scale α = 1. Node size and label size are proportional to node degree. Node fill colors are proportional to tumor expression fold changes in comparison to matched normal samples. Hub genesterm·key identified at any scale are labeled by gene symbols with font sizes proportional to the node degree.

(TIF) Click here for additional data file.

A) The number of significantly enriched functional/pathway signatures (Bonferroni corrected FET p-values) from MSigDB at various p-value thresholds. B) Number of significantly enriched functional/pathway signatures from MSigDB at the various odds ratio thresholds. C) Number of clusters predictive of patient survival (based on FDRterm·mth corrected logrank p-values) at various significance levels.

(TIF) Click here for additional data file.

Point shapes denote different clustering methods and point colors represent different co-expression networksterm·key, and point sizes represent significance of logrank p-values with FDRterm·mth corrected p-value < 0.05 threshold.

(TIF) Click here for additional data file.

A) Summarization scores from normalized ranks by various internal validity indices for clustering solutions across k. B) A heatmap of the pairwise Euclidean distances between α values. The distance was computed from the within-cluster connectivity matrix, Cw(V,A). The colorbar on top of the heatmap labels the scale clusters. C) The number of GO/KEGG/MSigDB gene sets enriched in the clusters at each scale group across a spectrum of Bonferroni corrected FET p-values.

(TIF) Click here for additional data file.

Each subplot compares two different stages of breast cancer. The y-axis represents–log10(Kolmogorov-Smirnov test p-value) and the x-axis represents different scales. The horizontal red line corresponds to KS p-value = 0.05.

(TIF) Click here for additional data file.

The x-axis shows the non-hub gene set and the hub geneterm·key sets at different scales. The category “multiscale” represents the hub geneterm·key set across all the different scales and the one “non.hub” represents the rest of genes. The y-axis shows the absolute values of log2(fc) between different cancer stages in LUAD.

(TIF) Click here for additional data file.

Each subplot compares two different stages of breast cancer. The y-axis represents–log10(Kolmogorov-Smirnov test p-value) and the x-axis represents different scales. The horizontal red line corresponds to KS p-value = 0.05.

(TIF) Click here for additional data file.

Red horizontal lines show Bonferroni corrected FET p-value = 0.05 for the number of nodes in LUAD PFNterm·mth.

(TIF) Click here for additional data file.