Integrative Multi-omics Module Network Inference with Lemon-Tree

Module network inference is an established statistical method to reconstruct co-expression modules and their upstream regulatory programs from integrated multi-omics datasets measuring the activity levels of various cellular components across different individuals, experimental conditions or time points of a dynamic process. We have developed Lemon-Tree, an open-source, platform-independent, modular, extensible software package implementing state-of-the-art ensemble methods for module network inference. We benchmarked Lemon-Tree using large-scale tumor datasets and showed that Lemon-Tree algorithms compare favorably with state-of-the-art module network inference software. We also analyzed a large dataset of somatic copy-number alterations and gene expression levels measured in glioblastoma samples from The Cancer Genome Atlas and found that Lemon-Tree correctly identifies known glioblastoma oncogenes and tumor suppressors as master regulators in the inferred module network. Novel candidate driver genes predicted by Lemon-Tree were validated using tumor pathway and survival analyses. Lemon-Tree is available from http://lemon-tree.googlecode.com under the GNU General Public License version 2.0.


Introduction
Recent years have witnessed a dramatic increase in new technologies for interrogating the activity levels of various cellular components on a genome-wide scale, including genomic, epigenomic, transcriptomic, and proteomic information [1]. It is generally acknowledged that integrating these heterogeneous datasets will provide more biological insights than performing separate analyses. For instance, in 2005, Garraway and colleagues combined SNP-based genetic maps and expression data to identify a novel transcription factor involved in melanoma progression [2]. More recently, international consortia such as The Cancer Genome Atlas (TCGA) or the International Cancer Genome Consortium (ICGC) have launched large-scale initiatives to characterize multiple types of cancer at different levels (genomic, transcriptomic, epigenomic, etc.) on several hundreds of samples. These integrative studies have already led to the identification of novel cancer genes [3,4]. This is a PLOS Computational Biology Software Article

Design and Implementation
Lemon-Tree is a platform-independent command-line tool written in Java which implements previously validated algorithms for model-based clustering [15] and module network inference [16]. The principal design difference between Lemon-Tree and other module network softwares (e.g. Genomica [5] or CONEXIC [12]) consists of the separation of module learning and regulator assignment. We have previously shown that running a two-way clustering algorithm until convergence, and thereafter identifying the regulatory programs that give rise to the inferred condition clusterings for each gene module results in higher module network model likelihoods and reduced computational cost compared to the traditional approach of iteratively updating gene modules and regulator assignments [14,16]. Hence Lemon-Tree is run as a series of tasks, where each task represents a self-contained step in the module network learning and evaluation process and the output of one task forms the input of another (a work flow representation of the different steps is illustrated in Fig. 1): Flow chart for integrative module network inference with Lemon-Tree. This figure shows the general workflow for a typical integrative module network inference with Lemon-Tree. Blue boxes indicate the pre-processing steps that are done using third-party software such as R or user-defined scripts. Green boxes indicate the core module network inference steps done with the Lemon-Tree software package. Typical postprocessing tasks (orange boxes), such as GO enrichment calculations, can be performed with Lemon-Tree or other tools. The Lemon-Tree task names are indicated in red (see main text for more details). doi:10.1371/journal.pcbi.1003983.g001 Task "ganesh" Run one or more instances of a model-based Gibbs sampler [15] to simultaneously infer co-expression modules and condition clusters within each module from a gene expression data matrix.
Task "tight_clusters" Build consensus modules of genes that systematically cluster together in an ensemble of multiple "ganesh" runs. Consensus modules are reconstructed by a novel spectral edge clustering algorithm which identifies densely connected sets of nodes in a weighted graph [28], with edge weight defined here as the frequency with which pairs of genes belong to the same cluster in individual "ganesh" runs. Details about the tight clustering algorithm are provided in S1 Text.
Task "regulators" Infer an ensemble of regulatory programs for a set of modules and compute a consensus regulator-to-module score. Regulatory programs take the form of a decision tree with the (expression level of) regulators at its internal nodes. The regulator score takes into account the number of trees a regulator is assigned to, with what score (posterior probability), and at which level of the tree [16]. An empirical distribution of scores of randomly assigned regulators is provided to assess significance. Regulator data need not come from the same data that was used for module construction but can be any continuous or discrete data type measured on the same samples. When multiple regulator types are considered, the "regulators" task is run once for each of them.
Task "experiments" For a fixed set of gene modules, cluster conditions separately for each module using a model-based Gibbs sampler [15] and store the resulting hierarchical condition trees in a structured XML file.
Task "split_reg" Assign regulators to a given range of one or more modules. This task allows parallelization of the "regulators" task and needs the output of the "experiments" task as an input.
Task "figures" Draw publication-ready visualizations for a set of modules in postscript format, consisting of a heatmap of genes in each module, organized according to a consensus clustering of the samples, plus heatmaps of its top-scoring regulators, separated according to the regulator type (cf. S1 Fig.).
Task "go_annotation" Calculate gene ontology enrichment for each module using the BiNGO [29] library.
While a typical run of Lemon-Tree will apply tasks "ganesh", "tight_clusters" and "regulators" in successive order, the software is designed to be flexible. For instance, the "tight_clusters" task can be equally well applied to build consensus clusters from the output of multiple third-party clustering algorithms, regulators can be assigned to the output of any clustering algorithm, etc. To facilitate this interoperability with other tools, input/output is handled via plain text files with minimal specification, the only exception being the storage of the regulatory decision trees which uses an XML format. Tasks also permit customization by changing the value of various parameters. We have purposefully provided default values for all parameters, based on our experience accrued over many years of developing and applying the software to a great variety of datasets from multiple organisms, and avoided mentioning any parameter settings in the Tutorial such that first-time users are presented with a simple workflow. Detailed instructions on how to integrate or extend (parts of) Lemon-Tree and a complete overview of all parameters and their default values are provided on the project website (http://lemon-tree.googlecode.com/).

Benchmark between Lemon-Tree and CONEXIC
We compared the performance of Lemon-Tree with CONEXIC (COpy Number and Expression In Cancer), a state-of-the-art module network algorithm designed to integrate matched copy number (amplifications and deletions) and gene expression data from tumor samples [12]. The general framework is the same for the algorithms, with modules of co-expressed genes associated to a list of regulators assigned via a probabilistic score. However, the probabilistic techniques used to build the modules and to assign regulators are different. We ran the two programs on the same large-scale reference data set to evaluate these differences. We used Gene Ontology (GO) enrichment and a reference network of protein-protein interactions to compare the co-expressed modules and the regulatory programs.
We downloaded gene expression and copy number glioblastoma datasets from the TCGA data portal [3] and we built an expression data matrix of 250 samples and 9,367 genes. We limited the number of samples for this benchmark study in order to save computational time. For the candidate regulators, we selected the top 1,000 genes that were significantly amplified or deleted as input genes for both CONEXIC and Lemon-Tree. To run CONEXIC, we followed the instructions of the manual and more specifically used the recommended bootstrapping procedure to get robust results. For Lemon-Tree, we generated an ensemble of two-way clustering solutions that were assembled in one robust solution by node clustering. Then we assigned the regulators using the same input list as with CONEXIC. A global score was calculated for each regulator and for each module and we selected the top 1% regulators as the final list (see S1 Text). The total run-time for the two software programs on the benchmark dataset was quite similar, with a small advantage for Lemon-Tree (S5 Table).
To compare the Gene Ontology (GO) categories between Lemon-Tree and CONEXIC, we built a list of all common categories for a given p-value threshold and converted the corrected p-values to −log 10 (p-value) scores. We selected the highest score for each GO category and we counted the number of GO categories having a higher score for Lemon-Tree or CONEXIC, and calculated the sum of scores for each GO category and each software. The results shown in Fig. 2 indicate that Lemon-Tree clusters have a higher number of GO categories with lower p-values than CONEXIC ( Fig. 2A), and that globally the p-values are lower for Lemon-Tree clusters (Fig. 2B). To benchmark the regulators' assignment of each software, we used a scoring scheme developed by Jornsten et al. [30]. For a given interaction distance in a reference protein-protein interaction network, we calculated the relative enrichment of known interactions in the networks inferred by Lemon-Tree and CONEXIC with respect to known interactions in networks where edges have been randomly re-assigned (see S1 Text). Fig. 2C shows the relative enrichments for interaction distances ranging from 1 (direct interaction) to 4. The Lemon-Tree inferred network is enriched for short or direct paths, a desired characteristic for well-estimated networks [30].
These results are consistent with a previous study conducted on bacteria and yeast data, where we showed a better performance in terms of enrichment in functional categories and known regulatory interactions of the algorithms underlying the Lemon-Tree software over Genomica (a software tool on which CONEXIC is based) [17]. Taken together, these results show that Lemon-Tree compares favorably with state-of-the-art module network inference algorithms.

Integrative analysis of TCGA glioblastoma expression and copy-number data
Lemon-Tree can be used to integrate various types of 'omics' data and generate new biological and biomedical insights. Here, we exemplify how to integrate copy-number and expression data for a large dataset of glioblastoma tumor samples and show that the results are enriched in known key players of canonical tumor pathways as well as novel candidates. Malignant gliomas are the most common subtype of primary brain tumors and are very aggressive, highly invasive and neurologically destructive. Glioblastoma multiforme (GBM) is the most malignant form of gliomas, and despite intense investigation of this disease in the past decades, most patients with GBM die within approximately 15 months of diagnosis [31]. Somatic copy-number alterations (SCNA) are extremely common in cancer and affect a larger fraction of the genome than any other types of somatic genetic alterations. They have critical roles in activating oncogenes and inactivating tumor suppressor genes, and their study has suggested novel potential therapeutic strategies [32,33]. However, distinguishing the alterations that drive cancer development from the passenger SCNAs that are acquired over time during cancer progression is a critical challenge. Here we use the module network framework implemented in the Lemon-Tree software tool to build a module network relating genes located in regions that are significantly amplified or deleted to modules of co-expressed genes. In other words, the module network selects and prioritizes copy-number altered genes that might play a role (direct or indirect) for clusters of co-expressed genes, performing important biological functions in glioblastoma. The resulting module network is used to prioritize SCNA genes that are amplified or deleted, and to provide novel hypotheses regarding drivers of glioblastoma.
We downloaded data from the TCGA project portal [3] and we selected 484 glioblastoma tumor samples from different patients (representing 91% of the available samples). We selected 7,574 gene expression profiles and generated an ensemble of two-way clustering solutions that were assembled in one robust solution by node clustering, resulting in a set of 121 clusters composed of 5,423 genes (S1 Text and S1 Table). We assembled a list of genes amplified and deleted in glioblastoma tumors from the most recent GISTIC run of the Broad Institute TCGA Copy Number Portal on glioblastoma samples. GISTIC [34] is the standard software tool used for the detection of peak regions significantly amplified or deleted in a number of samples from copy-number profiles. We also included in the list a number of key genes amplified or deleted from previous studies [34][35][36]. The final list is composed of 353 amplified and 2,007 deleted genes (with all genes present on sex chromosomes excluded). We did not use extremely stringent statistical thresholds for the selection, to avoid the exclusion of potentially interesting candidates. From this list we built SCNA gene copy-number profiles using TCGA data and used those profiles as candidate regulators for the co-expressed gene clusters. We assigned regulators independently for amplified and deleted genes, and we selected the top 1% highest scoring regulators as the final list (a cutoff well above assignment of regulators expected by chance), with 92 amplified and 579 deleted selected genes (S1 Text; S2 and S3 Tables). The resulting glioblastoma module network is composed of 121 clusters of co-expressed genes, together with associated prioritized lists of high-scoring SCNA genes (associated to amplified and deleted regions).
More than 60% of the clusters have a significant Gene Ontology (GO) enrichment (corrected p-value < 0.05, Table 2 and S4 Table). Several of those enriched clusters can be related to the hallmarks of cancers, ten distinctive and complementary capabilities that have been defined as the fundamental biological capabilities acquired during tumor development [37,38]. For instance, we have 11 clusters enriched for GO categories related to cell cycle processes and regulation (p-value < 0.05), with three of them having very strong enrichment (corrected pvalues 4×10 −18 , 6×10 −24 and 9×10 −71 , Table 2). The cell cycle is deregulated in most cancers and is at the heart of the "sustaining proliferative signaling" hallmark. Eight clusters are enriched for categories related to immune response, with two of them displaying strong enrichment (corrected p-values 6×10 −33 and 6×10 −45 , Table 2). Most tumor lesions contain immune cells present at various degrees of density. Intense recent research has shown that this immune response is linked to two phenomena. First, it is obviously an attempt by the immune system to eradicate the tumor, but secondly, there is now a large body of evidence showing that immune cells also have strong tumor-promoting effects, and both aspects are categorized as part of the hallmarks of cancer [38]. For instance, microglia are a type of glial cells that act as macrophages of the brain and the spinal cord and thus act as the main form of immune response in the central nervous system. They constitute the dominant form of glioma tumor infiltrating immune cells, and they might promote tumor growth by facilitating immunosuppression of the tumor microenvironment [39]. The development of blood vessels (angiogenesis) is another crucial hallmark of cancer, providing sustenance in oxygen and nutrients and a way to evacuate metabolic wastes and carbon dioxide [38]. Glioblastoma multiforme is characterized by a striking and dramatic induction of angiogenesis [31]. There are seven clusters enriched for GO categories related to angiogenesis and blood vessel development, with two of them having strong enrichment (corrected p-values 4×10 −6 and 9×10 −16 , Table 2). A recent large-scale integrative study of hundreds of glioblastoma samples has shown that chromatin modifications could potentially have high biological relevance for this type of tumor [40]. Interestingly, we have a cluster highly enriched in chromatin assembly and organization (corrected p-value 5×10 −17 and 9×10 −24 , Table 2). Taken together, these results show that the clusters of co-expressed genes in the module network are representative of the molecular functions and biological processes involved in tumor in general and more specifically in glioblastoma.
In the glioblastoma module network, we inferred a list of amplified and deleted SCNA genes linked to one or more clusters of co-expressed genes. Some of those SCNA genes are highly connected, representing potential master copy-number regulators for module activity. To identify and analyze those SCNA hub genes, we calculated for each high-scoring regulator the sum of the scores obtained in each module, and ranked them by decreasing score for amplified (Table 3) and deleted (Table 4) genes. Among these genes, we find many well-known oncogenes and tumor supressors that are frequently amplified, deleted or mutated in glioblastoma. Those genes include EGFR, PDGFRA, FGFR3, PIK3CA, MDM4, CDKN2A/B, PTEN and are all members of the core alterated pathways in glioblastoma controlling key phenotypes such as proliferation, apoptosis and angiogenesis (Fig. 3, [3,35,40,41]). Those genes and pathways are also frequently impaired in many other types of tumors [42][43][44]. In addition, we find in those lists of hub genes a number of interesting new candidates, both in amplified and deleted genes, that have not been associated with glioblastoma before. To better visualize the importance and role of both the well-known and novel SCNAs prioritized by Lemon-Tree, we represent those that are part of the three core pathways altered in glioblastoma as a network with edges representing activation or inhibition relationships, together with their levels of gene gains and losses in glioblastoma samples (Fig. 3).
Within the list of amplified gene hubs (Table 3), we find a number of genes that have been rarely or never associated before with glioblastoma. INSR is a gene encoding for the insulin receptor, a transmembrane receptor activated by insuline and IGF factors, member of the tyrosine receptor kinase family, and playing a key role in glucose homeostasis. INSR is selected as a high-scoring regulator in 15 modules and ranked in third position in the list of amplified gene hubs. It is found to be amplified as low-level gain or higher in 39% of the samples (Table 3). Beyond its well-known role in glucose homeostasis, INSR stimulates cell proliferation (Fig. 3) and migration and is often aberrantly expressed in cancer cells [45]. Consequently, amplification of INSR in glioblastoma may enhance proliferation. MYCN encodes a transcription factor (Nmyc) highly expressed in fetal brain and critical for normal brain development. It is also a wellknown proto-oncogene, and amplification of N-myc is associated with poor outcome in neuroblastoma [46]. MYCN is amplified as low-level gain or higher in 8% of the glioblastoma samples and is connected to 21 modules (Table 3). MYCN is part of the RB signaling pathway, and is also strongly connected to the RTK / PI3K and p53 pathways (Fig. 3), with a direct influence on proliferation. For that reason, its amplification may also favor proliferation in glioblastoma. KRIT1 (also known as CCM1) is a gene crucial for maintaining the integrity of the vasculature and for normal angiogenesis. Loss of function of this gene is responsible for vascular malformations in the brain known as cerebral cavernous malformations [47,48]. It is amplified as lowlevel gain or higher in 83% of the glioblastoma samples and it is listed in the top 10 hubs in our list ( Table 3). The consequences of KRIT1 amplification are not completely clear, but we may hypothesize that it is required for proper angiogenesis development, which is a hallmark of glioblastoma [31], and that it may also help decrease apoptosis (Fig. 3).
In the list of putative deleted genes, PAOX (polyamine oxidase) is ranked first, with a connection to 54 modules and the highest sum of scores value. It is classified as single loss (GISTIC call value of -1 or lower) in 89% of the samples. This is a very high value, comparable to the value obtained for the classical tumor suppressor CDKN2A (75%, Table 4). Amine oxidases are involved in the metabolism of polyamines, regulating their intracellular concentrations and elimination. The products of this metabolism (e.g. hydrogen peroxyde) are cytotoxic and have been considered as a cause for apoptotic cell death. Amine oxidases are considered as biological regulators for cell growth and differentiation, and a primary involvement of amine oxidases in cancer growth inhibition and progression has been demonstrated [49]. Therefore, PAOX might have a tumor suppressor activity and its deletion in many glioblastoma samples could provide a selective advantage to glioblastoma tumor cells. Interestingly, amino acids metabolism is not part of the standard alterated pathways in glioblastoma (explaining why we did not represent PAOX on Fig. 3), but targeting this pathway could lead to novel therapeutic treatments [50]. KLLN encodes a nuclear transcription factor and shares a bidirectional promoter with PTEN. It is activated by p53 and is involved in S phase arrest and apoptosis [51]. Recent studies show that KLLN has a tumor supressor effect and is associated with worse prognosis in prostate and breast carcinomas [52,53]. Consequently, the loss of KLLN that is observed in 88% of the glioblastoma samples (Table 4) would help the development of tumor cells by decreasing apoptosis and favoring proliferation (Fig. 3).
To assess the biological relevance of the amplified and deleted gene hubs in the module network, we analyzed the prognosis value of the top gene hubs by survival analysis, using the clinical data available for TCGA samples (survival time and status of the patient). We constructed Kaplan-Meier estimates using GISTIC putative calls to define genes having single or deep copy loss (i.e. GISTIC call value -1) and genes having low-level gains or high-level amplifications (i.e. GISTIC call value 1. The differences between groups were formally tested and a total of 3 amplified genes and 18 deleted genes from the lists displayed in Tables 3 and 4 have significant p-values < 0.05 ( Fig. 4 and S6 Table). Interestingly, among those genes we find the well-known glioblastoma oncogene EGFR and tumor suppressors CDKN2A and PTEN, but also novel candidates such as KRIT1 and PAOX described in the previous paragraph. Glioblastoma patients having copy-number alterations for those genes have a worse survival prognostic. This indicates the biological relevance of those genes that may be used as biomarkers.

Availability and Future Directions
The Lemon-Tree software is hosted at Google Code (http://lemon-tree.googlecode.com/). The source code, executables and documentation can be downloaded with no restrictions and no registration, and are released under the terms of the GNU General Public License (GPL) version 2.0. Developers and users can join the project by contacting the authors and there is a mailing list for discussions and news about module networks and the project. A step-by-step tutorial to learn how to install and use the software is available on the wiki section, together with the corresponding data sets.
In the future, we intend to extend Lemon-Tree's support for explicitly modelling causal relations between regulator types and to incorporate complementary algorithms available in the literature for integrating gene-based methods, physical interactions and cross-species data. Firstly, the current version of Lemon-Tree is able to associate co-expression modules to multiple 'regulator' types (e.g. expression regulators, structural DNA variants, phenotypic states, etc.) by assigning each of those independently as regulators of a module. We will extend the software with Bayesian methods to account for possible causal relations between regulator types, e.g. when the association between a module and expression regulator can be partly explained by a structural DNA variant. Secondly, a key long-term objective of the Lemon-Tree project is to provide a general open-source repository for module network inference tools with a consistent user interface. As a first step, the current version of Lemon-Tree implements algorithms previously developed by our group [14][15][16][17]. In the future, we intend to extend it with complementary algorithms developed by other groups, including algorithms to combine the strengths of module network methods with gene-based methods [18], to incorporate physical protein-protein or protein-DNA interactions as a prior in the regulator assignment procedure [13] or to infer module networks from multiple species simultaneously [19]. A document detailing guidelines to implement new functions in the Lemon-Tree Java codebase is available on the project wiki.
Supporting Information S1  Kaplan-Meier survival curves for a selection of top hub glioblastoma genes predicted by the Lemon-Tree algorithm. The top three panels are genes having low-levels gains or high-level amplifications (magenta) compared to normal (blue), the bottom three panels are genes having single-copy loss or homozygous deletions (green) compared to normal (blue). All genes display significant differences between the groups (p < 0.05, see S6 Table for a