DREAM5 - Gene Network Inference

Supplementary Notes Wisdom of crowds for robust gene network inference Daniel Marbach♣,1,2, James Costello♣,3, Robert Küffner♣,4, Nicci M Vega3, Robert Prill6, Diogo M Camacho3, Kyle R Allison3, the DREAM5 Consortium, Manolis Kellis1,2, James J Collins3,5, and Gustavo Stolovitzky6 The DREAM5 Consortium: Andrej Aderhold25,27, Richard Bonneau30,31,33, Yukun Chen9, Francesca Cordero12,13, Martin Crane34, Frank Dondelinger25,26, Mathias Drton8, Roberto Esposito12, Rina Foygel8, Alberto de la Fuente14, Jan Gertheiss7, Pierre Geurts21,22, Alex Greenfield31, Marco Grzegorczyk29, Anne-Claire Haury17,18,19, Benjamin Holmes1,2, Torsten Hothorn7, Dirk Husmeier25, Vân Anh Huynh-Thu21,22, Alexandre Irrthum21,22, Guy Karlebach35, Sophie Lèbre28, Vincenzo De Leo14,15, Aviv Madar30, Subramani Mani9, Fantine Mordelet17,18,19,20, Harry Ostrer32, Zhengyu Ouyang16, Ravi Pandya11, Tobias Petri4, Andrea Pinna14, Christopher S. Poultney30, Serena Rezny8, Heather J. Ruskin34, Yvan Saeys23,24, Ron Shamir35, Alina Sîrbu34, Mingzhou Song16, Nicola Soranzo14, Alexander Statnikov10, Paola Vera-Licona17,18,19, Jean-Philippe Vert17,18,19, Alessia Visconti12, Haizhou Wang16, Louis Wehenkel21,22, Lukas Windhager4, Yang Zhang16, and Ralf Zimmer4 1Computer Science and Artificial Intelligence Laboratory, Massachusetts Institute of Technology, Cambridge, MA, USA 2Broad Institute of MIT and Harvard, Cambridge, MA, USA 3Howard Hughes Medical Institute, Center for BioDynamics, and Dept. of Bioengineering, Boston University, Boston, MA, USA 4Ludwig-Maximilians University, Department of Informatics, Amalienstr. 17, 80333 Munich, Germany 5Wyss Institute for Biologically Inspired Engineering, Harvard University, Boston, MA, USA 6IBM T. J. Watson Research Center, Yorktown Heights, New York, NY, USA 7Ludwig-Maximilians University, Department of Statistics, Ludwigstr. 33, 80539 Munich, Germany 8University of Chicago, Department of Statistics, Chicago, IL, USA 9Department of Biomedical Informatics, Vanderbilt University, Nashville, TN, USA 10 Center for Health Informatics and Bioinformatics, New York University, New York, NY, USA 11 Microsoft Research, 1 Microsoft Way, Redmond, WA, USA 12 Department of Computer Science, Corso Svizzera 185, 10149, Torino, Italy 13 Department of Clinical and Biological Sciences, Regione Gonzole 10, Orbassano, Italy 14 CRS4 Bioinformatica, Parco Tecnologico della Sardegna, Edificio 3, Loc. Piscina Manna, 09010 Pula (CA), Italy 15 Linkalab, Complex Systems Computational Laboratory, 09100 Cagliari (CA), Italy 16 Department of Computer Science, New Mexico State University, Las Cruces, NM, USA 17 Mines ParisTech, CBIO, 35 rue Saint-Honoré, Fontainebleau, F-77300, France 18 Institut Curie, Paris, F-75248, France 19 INSERM, U900, Paris, F-75248, France 20 CREST, INSEE, Malakoff, F-92240, France 21 Department of Electrical Engineering and Computer Science, Systems and Modeling, University of Liège, Belgium 22 GIGA-Research, Bioinformatics and Modeling, University of Liège, Belgium 23 Department of Plant Systems Biology, VIB, Gent, Belgium 24 Department of Plant Biotechnology and Bioinformatics, Ghent University, Gent, Belgium 25 Biomathematics and Statistics Scotland, Edinburgh & Aberdeen, UK 26 School of Informatics, University of Edinburgh, UK 27 School of Biology, University of St Andrews, UK 28 LSIIT, UMR UdS-CNRS 7005, Illkirch - Université de Strasbourg, France 29 Department of Statistics, TU Dortmund University, Germany 30 Department of Biology, Center for Genomics & Systems Biology, New York University, New York, NY, USA 31 Computational Biology Program, New York University Sackler School of Medicine, New York, NY, USA 32 Human Genetics Program, Department of Pediatrics, New York University Langone Medical Center, New York, NY, USA 33 Computer Science Department, Courant institute of Mathematical Sciences, New York University, New York, NY, USA 34 Centre for Scientific Computing and Complex Systems Modelling, School of Computing, Dublin City University, Ireland 35 The Blavatnik School of Computer Science, Tel Aviv University, Israel ♣These authors contributed equally 1 Nature Methods: doi:10.1038/nmeth.2016 Supplementary Notes 1 DREAM5 network inference challenge description 5 2 Gene expression compendia 5 3 Gold standards 6 4 Assessment of network inference methods 9 4.1 Performance assessment ............................................ 9 4.2 Clustering of methods by PCA ........................................ 9 4.3 Network motif analysis ............................................. 12 5 Data information content 13 5.1 Estimation of inference difficulty ....................................... 13 5.2 Information content of different experiment types .............................. 16 6 Integration of predictions 17 6.1 Why community integration can outperform the best individuals ..................... 17 6.2 DREAM5 community networks ........................................ 20 7 E. coli and S. aureus community networks 21 7.1 Network construction ............................................. 21 7.2 S. aureus network evaluation using RegPrecise ............................... 22 7.3 Analysis of network modules ......................................... 23 8 Experimental validation 24 8.1 Transcription factor selection ......................................... 24 8.2 RNA Extraction and qPCR .......................................... 27 9 Methodological insights 27 10 Network inference methods 29 10.1 Regression 1 – Inferring gene regulatory networks with the stabilized Lasso ............... 30 10.2 Regression 2 – Network inference with regularized linear regression methods ............... 31 10.3 Regression 3 – Sparse piecewise linear regression based on changepoint processes ............ 32 10.4 Regression 7 – Simple L1-regularization ................................... 33 10.5 Regression 8 – Linear regression* ....................................... 34 10.6 Mutual Information 1 – Context Likelihood of Relatedness (CLR)* .................... 34 10.7 Mutual Information 2 – Mutual information* ................................ 35 10.8 Mutual Information 3 – Algorithm for the reconstruction of accurate cellular networks (Aracne)* . . . 35 10.9 Mutual Information 4 & 5 –The DREAM5 network inference challenge with a combination of fast tools 35 10.10 Correlation 2 – Pearson’s correlation* .................................... 36 10.11 Correlation 3 – Spearman’s correlation* ................................... 36 10.12 Bayesian 6 – Regulatory network inference with Bayesian networks .................... 36 10.13 Other 1 – Inferring regulatory networks using tree-based methods .................... 38 10.14 Other 2 – Inferring gene regulatory networks by ANOVA ......................... 39 10.15 Other 3 – Network inference through Boolean networks .......................... 40 2 Nature Methods: doi:10.1038/nmeth.2016 10.16 Other 6 – Network inference using quantitative modeling and evolutionary algorithms ......... 41 10.17 Other 7 – Detecting interactions by generalized logic ............................ 43 10.18 Other 8 – Finding gene-gene interactions by concurrence of change between conditions ......... 44 10.19 Meta 1 — Inferring gene regulatory networks using knockout data and resampling ........... 44 10.20 Meta 3 — Analyzing subsets of heterogeneous gene expression compendia to elucidate transcriptional regulatory networks .............................................. 47 10.21 Meta 5 – A naïve Bayes based approach to network inference ....................... 48 Figures S1 Evaluation of datasets and inference approaches ................................ 6 S2 AUPR and AUROC values using different gold standards for S. cerevisiae ................. 8 S3 PR and ROC curves for individual methods and community predictions .................. 10 S4 Comparison of the method rankings across compendia ............................ 11 S5 PCA for individual compendia ......................................... 12 S6 First principal component ............................................ 12 S7 Definition of motif types used to analyze prediction biases .......................... 13 S8 Mutual dependency of mRNA levels between transcription factors and target genes in different compendia 14 S9 Transcription factor-target gene dependency of mRNA levels for alternative yeast compendia and gold standards ..................................................... 15 S10 Information content of different experiment types for network inference ................... 17 S11 Significance of individual microarrays for network inference ......................... 18 S12 Why community integration can outperform the best individual inference methods ............ 19 S13 Community integration using unweighted and weighted voting ........................ 21 S14 Rank improvement of individual edges through community integration ................... 22 S15 Performance evaluation for S. aureus using the RegPrecise database .................... 23 S16 Functional enrichment of network modules in E. coli ............................. 25 S17 Functional enrichment of network modules in S. aureus ............................ 26 S18 Experimental support for newly predicted interactions ............................ 28 S19 Inference method design for Bayesian 6. .................................... 37 S20 Network inference pipelines tested for Meta 1. ................................. 46 S21 Inference method design for Meta 5. ...................................... 48 Tables S1 χ2contingency table for Other 8 ........................................ 45 Websites DREAM website (wiki.c2b2.columbia.edu/dream/index.php/D5c4) Current information on the DREAM5 network inference challenge, including anonymized data, challenge description, and evaluation scripts. M3D website (m3d.bu.edu/dream) The deanonymized expression compendia are made publicly available from the Many Microbe Microarrays Database (M3D). 3 Nature Methods: doi:10.1038/nmeth.2016 GP-DREAM web platform (dream.broadinstitute.org) The GP-DREAM platform offers a toolkit of network inference and consensus methods. Overview of the supplementary notes Supplementary Note 1,Challenge description. This chapter briefly introduces the general outline of the challenge including the provided datasets, participation requirements and evaluation criteria. These aspects are described in more detail in the subsequent three sections. Supplementary Note 2,Gene expression compendia. The four gene expression compendia provided to the participants are described here. This section also describes how the set of valid transcription factors was obtained. The set of transcription factors was also available to the participants. Supplementary Note 3,Gold standards. This section explains the compilation of the gene regulatrory networks, i.e. the gold standards used for the evaluation of the challenge participants. In case of S. cerevisiae, we describe the compilation of two alternative gold standards used to evaluate the poor performance for S. cerevisiae. Note that this and the previous section expand on the Methods: Expression data and gold standards section of the main manuscript. Supplementary Note 4,Assessment of network inference methods. This section describes the metrics used for the assessment of methods, i.e. the AUROC, AUPR as well as the overall score. It further details our PCA and Network motif analyses including additional results and figures. It thus expands on three Methods sections of the main manuscript, namely Performance metrics,Clustering of inference approaches by principal component analysis (PCA) and Network motif analysis. Supplementary Note 5,Data information content. This section provides the results and figures for two aspects that are mentioned briefly in the Discussion section of the main manuscript, namely the performance differences between the three organisms as well as which kinds of microarray experiments provide the most valuable information for network inference. In particular, the low performance of network inference in S. cerevisiae is analyzed and possible reasons are discussed with the corresponding references. Supplementary Note 6,Integration of predictions. Here, we discuss the details on our ensemble approach integrating 35 network inference approaches to obtain consensus predictions. This section also expands on Box 1 from the main manuscript and describes the simulation approach that generated the theoretical distributions shown in the box. Supplementary Note 7,E. coli and S. aureus community networks. This section describes the details of how the community networks were constructed along with further analysis of the S. aureus network using the RegPrecise database. Module detection and GO term enrichment methods are also fully explained. Supplementary Note 8,Experimental validation. Experiments were designed to test transcription factor to target gene regulation. The selection criteria for transcription factors along with the full experimental details are covered. Supplementary Note 9,Methodological insights. This section extends the Discussion section of the main manuscript and describes in detail the method specific advantages and disadvantages that we identified. We analyze and discuss the reasons of the network motif specific performance of certain methods and method classes. We also discuss performance improvements due to method specific enhancements as well as performance degradations due to method specific omissions. Supplementary Note 10,Network inference methods. This section contains the detailed descriptions of 15 network inference methods that were supplied by participants of the challenge. Further, six off-the-shelf methods (marked by an asterisk in the title) are described that were used to compare the performance of participants against state-of-the-art approaches. 4 Nature Methods: doi:10.1038/nmeth.2016 1 DREAM5 network inference challenge description The DREAM5 network inference challenge solicited predictions of genome-scale transcriptional regulatory networks for expression compendia in E. coli,S. cerevisiae,S. aureus, and an in silico compendium. Each compendium is represented as an expression matrix of ggenes by cchip measurements (Figure SS1). In addition to the gene expression data (see below), a number of descriptive features were supplied for each microarray experiment (e.g., temporal information if the experiment is part of a time series, or the deleted gene if it is a gene deletion experiment). Participants were further given a list of candidate transcription factors for each compendium. To fully anonymize the datasets, for the E. coli,S. cerevisiae, and S. aureus compendia, all genes were assigned an arbitrary gene identifier (e.g. G1, G2, . . . , GN, where there are Ngenes in a given compendium). Additionally, a set of decoy genes, representing roughly 5% of the compendium, were introduced by randomly selecting gene expression values from the compendium, itself. Gene expression profiles for the in silico network were derived from GeneNetWeaver.56 In order to participate in the challenge, predictions for all four compendia were required. For each network, a list of directed, unsigned edges had to be submitted ordered according to the participant’s confidence scores. The set of submitted predictions were compiled into the prediction matrix (Figure SS1) that contains 29 entries from the participating teams, 6 entries from publicly available methods, and an entry from integrating the 29 methods into an additional set of integrated community predictions. In the prediction matrix, each method is represented as an ordered list of predicted gene regulatory interactions, ranked for prediction confidence. Organism specific gold standards containing the known transcription factor to target gene (transcription factortarget gene) interactions (= true positives) were compiled for assessing the participating approaches. For the evaluation, we considered all transcription factor-target gene pairs that are not part of the gold standards as negatives, although, as the gold standards are based on incomplete knowledge, they might contain yet unknown true interactions. Expression compendia and gold standards are briefly described in the following sections. The full description of the challenge is available on the challenge websitea. ahttp://wiki.c2b2.columbia.edu/dream/index.php/D5c4 2 Gene expression compendia Expression compendia, gene lists, and transcription factor lists for all considered datasets are supplied in Supplementary Data 1. The data is also freely available from the DREAM websiteband the Many Microbe Microarrays Database (M3D)c. Staphylococcus aureus. A compendium of microarray data was compiled for S. aureus, where all chips are the same Affymetrix platform, the S. aureus Genome Arrayd. Chips were downloaded from Gene Expression Omnibus (GEO) (Platform ID: GPL1339), a publicly available web repository hosted by the National Center for Biotechnology Information (NCBI)e. In total, 160 chips with available raw data Affymetrix files (.CEL files) were compiled. Microarray normalization was done using Robust Multichip Averaging (RMA)9through the software RMAExpressf. All 160 chips were uploaded into RMAExpress and normalization was done as one batch. All arrays were background adjusted, quantile normalized, and probesets were summarized using median polish. Normalized data was exported as log-transformed expression values. Mapping of Affymetrix probeset ids to gene ids was done using the library files made available from Affymetrix. Control probesets and probesets that did not map unambiguously to one gene were removed, specifically probeset ids ending in _x, _s, _i were removed. Lastly, if multiple probesets mapped to a single gene, then expression values were averaged within each chip. Completion of these steps resulted in a total of 2,677 genes over the 160 microarrays. In addition to the gene ×condition matrix, DREAM5 participants were also supplied with a list of known or putative transcription factors. Transcription factors were identified based on their Gene Ontology (GO) annotation. We selected genes annotated with a biological process related to transcription, namely GO:0009299;mRNA transcription or GO:0006351;transcription, DNA dependent. Additionally, we also required that a gene be annotated with a molecular function of GO:0003677;DNA binding or any child terms. A total of 90 genes were designated as potential transcription factors. Escherichia coli. A compendium of microarray data was compiled for E. coli, where all chips are the same Affymetrix platform, the E. coli Antisense Genome Arrayg. Chips were downloaded from GEO (Platform ID: GPL199). In total, 805 chips with available raw data Affymetrix files (.CEL files) were compiled. bhttp://wiki.c2b2.columbia.edu/dream chttp://m3d.bu.edu/dream/ dhttp://www.affymetrix.com ehttp://www.ncbi.nlm.nih.gov/geo fhttp://rmaexpress.bmbolstad.com ghttp://www.affymetrix.com 5 Nature Methods: doi:10.1038/nmeth.2016 mRNA levels Section S1 mRNA levels Section S1 Gold standard Section S3 Inference methods Section S10 Evaluation of inference methods Section S4 Evaluation of data information content Section S5 Community integration Section S6 mRNA levels Section S2 confidence ranks Section S1 Data matrix: c chips x g genes Prediction matrix: m methods x i interactions Figure S1: Evaluation of datasets and inference approaches DREAM5 participants were solicited to infer gene regulatory interactions (prediction matrix, green) from three different expression compendia (data matrix, blue). An additional set of predictions was derived by integrating all teams into community predictions. Analysis steps discussed in the present paper differ according to which of the two matrices are used and whether they depend on the gold standard of known gene regulatory interactions. Microarray normalization was done as described for S. aureus above. Completion of microarray normalization and filtering resulted in a total of 4,297 genes over the 805 microarrays. The list of potential transcription factors was obtained from two sources. First, transcription factors defined by RegulonDB were used.31 Second, transcription factors were identified using GO terms as described for S. aureus above. A total of 296 genes were designated as potential transcription factors. Saccharomyces cerevisiae. A compendium of microarray data was compiled for S. cerevisiae, where all chips are the same Affymetrix platform, the Affymetrix Yeast Genome S98 Arrayh. Chips were downloaded from GEO (Platform ID: GPL). In total, 536 chips with available raw data Affymetrix files (.CEL files) were compiled. Microarray normalization was done as described for S. aureus above. Completion of microarray normalization and filtering resulted in a total of 5,667 genes over the 536 microarrays. We used the list of potential transcription factors defined by Zhu et al.92 We further added transcription factors identified using GO terms as described above. A total of 183 genes were designated as potential transcription factors. We also evaluated the performance of some algorithms on an independent yeast dataset of 904 chip measurements from the M3D database.27 This dataset was used for internal comparison and was not provided to the participants of the DREAM5 challenge. The preparation and hhttp://www.affymetrix.com preprocessing of this dataset was performed in the same way as the datasets of the challenge. 3 Gold standards In this section, we describe how the gold standards for E. coli and S. cerevisiae were compiled. S. aureus was not used for benchmark evaluation and no gold standard was constructed. For the in silico benchmark, the true network structure is known and was used as gold standard. All gold standards are supplied in Supplementary Data 1. Note that in contrast to the in silico network, the gold standards for E. coli and S. cerevisiae are obviously not perfect. For S. cerevisiae, we tested a range of alternative gold standards (see below). Of these, we chose the most stringent gold standards for both organisms, which include only interactions with strong experimental support. Thus, most of the edges contained in the gold standards are likely to be true (i.e., the gold standards are expected to have relatively few false positives). However, they contain only a subset of the true interactions (i.e., they have many false negatives). Thus, predicted interactions that are not part of the gold standard should not be considered incorrect — they may also be newly discovered interactions that are currently missing in the gold standard, as our experimental validation of such novel interactions demonstrates (Figure 4c of the main text). Consequently, the reported precision and false positive rate (Supplementary Note 4.1) of network predictions should be considered with caution. 6 Nature Methods: doi:10.1038/nmeth.2016 Escherichia coli. The model organism E. coli has a well-studied transcriptional regulatory network, which makes it well-suited as a benchmark for network inference. Known transcriptional interactions are collected in the manually curated EcoCyc42 and RegulonDB31 databases (the two databases are synchronized). Each interaction is annotated with a set of evidence codes, which are classified as either strong or weak evidencei. The E. coli gold standard was constructed from RegulonDB Release 6.8. Only transcriptional interactions with at least one strong evidence were included (2,066 interactions). Saccharomyces cerevisiae. Several large-scale studies and databases have produced genome-wide regulatory networks for S. cerevisiae. We have confirmed that our results are consistent and reproducible across several alternative gold standards. In particular, we find that the performance of inference methods is low for S. cerevisiae (see Discussion of the main text and Fig. SS8 &Fig. SS9) independently of the gold standard. In total, we have tested 16 gold standards derived from three sources.1,39,50 In what follows, we describe the different gold standards and discuss the measured performance (AUPR and AUROC, see Supplementary Note 4.1) of inference methods on these gold standards (Fig. SS2). The first set of gold standards was obtained from the study of MacIsaac et al.,50 which is based on a re-analysis of ChIP-chip data for 203 transcription factors from Harbison et al.33 Regulatory interactions were identified based on measured binding (ChIP) and/or presence of evolutionary conserved motifs of the transcription factor in the intergenic region upstream of target genes. By varying the thresholds required for binding and evolutionary conservation of motifs, different versions of the network were obtained. We considered all nine versions available from the author’s websitej(genes and transcription factors that are not part of our expression compendium were excluded). The gold standard based on the most stringent thresholds, which includes only interactions with strong evidence of binding and a strongly conserved motif, also leads to the “strongest signal” (highest AUROCs and biggest fold-improvements for the AUPRsk), i.e., it “agrees best” ihttp://regulondb.ccg.unam.mx/evidenceClassification.jsp jhttp://fraenkel.mit.edu/improved_map kThe AUPR of network predictions depends on the connectivity of the gold standard, i.e., the ratio of positives (present edges) to negatives (absent edges). The more densely connected the gold standard, the easier it is to correctly “guess” true edges and obtain a high AUPR. For example, the same predictions have an AUPR of ∼32% on the densely connected (low confidence) gold standard of the 1st row, and ∼2% on the loosely connected (high confidence) gold standard of the 9th row in Fig. SS2. Thus, the absolute value of the AUPR should be considered with caution and always be compared to the expected AUPR of a random prediction, for instance. with the inferred networks. In contrast, the gold standard based on the loosest thresholds, which requires no ChIP binding and no evolutionary conservation (only the presence of a motif), leads to the “weakest signal” — likely because it contains many false positives (motif instances that are not bound and/or not functional in vivo). Incidentally, our observations confirm the conclusion of MacIsaac et al. and others43 that evolutionary conservation of motifs can be used as a signal to improve the quality of ChIP-based networks. The second set of gold standards was derived from the study of Hu, Killion & Iyer,39 which conducted a comprehensive mRNA profiling experiment of 269 yeast transcription factor deletion mutants. Genes were considered a target of a given transcription factor if they exhibited a fold-change above a certain threshold. Independently of the considered threshold, the overlap of the inferred networks and this gold standard is not better than expected by chance (AUPR and AUROC values are similar to those expected for random predictions, Fig. SS2). Finally we evaluated a gold standard from the curated YEASTRACT database,1which compiles direct and indirect evidence for yeast gene regulatory interactions from more than 1,200 publications. While direct evidence (28,336 interactions, as of version 1.1503, Jun 26, 2010) was predominantly derived from ChIP experiments, indirect evidence (21,847 interactions) was gathered from expression measurements of transcription factor deletion or overexpression mutants. In the present document we will refer to the YEASTRACT gold standard as the intersection of interactions with both direct and indirect evidence (2,528 interactions, after filtering out transcription factors and genes not part of our expression compendium). The YEASTRACT gold standard leads to slightly better AUPR and AUROC values than the high-confidence network of MacIsaac et al, however, the better agreement between the inferred networks and YEASTRACT may be because YEASTRACT includes expression data as evidence and is thus not completely independent from the inferred networks, which are also expression-based. In contrast, the MacIsaac gold standard is based on orthogonal datasets (ChIP and conserved motifs) and does not incorporate expression-based evidence. Based on the above observations, we used the network based on the most stringent thresholds from MacIsaac et al. as gold standard for all method related assessments and all other analyses, unless noted otherwise. 7 Nature Methods: doi:10.1038/nmeth.2016 119 32.3 32.7 32.6 32.5 32.7 32.9 32.7 32.5 32.6 32.9 32.9 32.8 33.0 33.1 32.8 33.1 33.0 32.3 32.4 32.7 32.6 32.6 33.6 32.9 32.2 32.9 32.6 32.1 32.9 32.2 31.7 33.1 32.8 32.9 32.3 32.4 32.9 119 19.6 20.1 19.8 19.8 20.0 20.2 19.9 19.5 19.9 20.2 20.2 20.2 20.3 20.5 20.3 20.5 20.5 19.6 19.6 20.1 20.0 20.0 20.3 20.2 19.5 20.3 20.0 19.4 20.0 19.7 19.2 20.4 20.3 20.2 19.6 19.9 20.3 119 11.8 12.2 12.0 12.0 12.2 12.4 12.1 11.7 12.1 12.4 12.3 12.3 12.4 12.6 12.4 12.5 12.5 11.8 11.8 12.2 12.2 12.2 12.2 12.3 11.8 12.3 12.1 11.6 12.1 11.9 11.5 12.4 12.4 12.4 11.9 12.2 12.4 116 2.8 3.1 2.8 3.0 3.1 3.2 3.0 2.8 3.0 3.1 2.9 2.9 2.8 3.0 2.9 3.0 3.0 2.9 2.9 3.1 3.0 3.0 3.0 3.1 3.2 3.0 3.2 3.0 2.9 2.7 2.7 3.6 3.0 3.7 3.3 3.6 3.3 115 2.4 2.6 2.4 2.6 2.7 2.7 2.6 2.4 2.5 2.7 2.5 2.4 2.4 2.6 2.5 2.6 2.6 2.4 2.4 2.7 2.6 2.6 2.5 2.7 2.7 2.6 2.9 2.6 2.5 2.3 2.3 3.2 2.5 3.3 2.8 3.3 2.8 115 2.1 2.4 2.1 2.3 2.4 2.4 2.3 2.1 2.3 2.4 2.2 2.2 2.2 2.3 2.2 2.3 2.3 2.1 2.1 2.4 2.3 2.3 2.2 2.4 2.4 2.3 2.7 2.3 2.2 2.0 2.0 3.0 2.3 3.1 2.6 3.2 2.6 116 2.1 2.4 2.1 2.3 2.4 2.5 2.4 2.2 2.3 2.4 2.2 2.2 2.2 2.3 2.2 2.3 2.3 2.2 2.2 2.4 2.3 2.3 2.3 2.4 2.5 2.3 2.7 2.3 2.2 2.1 2.0 3.2 2.3 3.1 2.7 3.3 2.6 114 1.9 2.2 1.9 2.1 2.2 2.2 2.2 2.0 2.1 2.2 2.0 2.0 2.0 2.2 2.0 2.1 2.1 2.0 2.0 2.2 2.1 2.1 2.1 2.2 2.3 2.1 2.6 2.1 2.0 1.9 1.8 3.1 2.1 2.9 2.4 3.2 2.4 114 1.7 2.0 1.8 2.0 2.1 2.1 2.0 1.8 1.9 2.1 1.9 1.8 1.8 2.0 1.9 2.0 1.9 1.8 1.8 2.0 2.0 2.0 1.9 2.1 2.2 1.9 2.5 2.0 1.9 1.7 1.7 3.1 1.9 2.8 2.3 3.2 2.2 234 4.1 4.0 3.8 4.0 4.0 4.1 3.9 3.9 4.0 4.0 4.0 4.0 4.0 3.9 4.0 4.0 3.9 4.0 4.0 4.0 4.0 4.0 4.1 4.0 4.1 4.0 4.1 3.9 4.2 4.0 4.2 4.0 4.0 4.0 3.9 4.0 4.0 233 1.9 1.9 1.7 1.9 1.9 1.9 1.8 1.8 1.9 1.9 1.9 1.8 1.9 1.8 1.9 1.9 1.8 1.9 1.9 1.9 1.9 1.9 1.9 1.9 2.0 1.9 1.9 1.8 2.1 1.8 2.0 1.8 1.9 1.9 1.8 1.9 1.9 202 1.5 1.6 1.4 1.5 1.6 1.6 1.4 1.5 1.5 1.5 1.5 1.5 1.6 1.5 1.6 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.7 1.6 1.6 1.5 1.8 1.5 1.7 1.4 1.5 1.5 1.5 1.6 1.6 234 4.8 4.9 4.5 4.8 4.8 4.8 4.6 4.7 4.8 4.8 4.7 4.8 4.8 4.8 4.8 4.9 4.8 4.8 4.8 4.9 4.9 4.9 4.8 4.8 4.9 4.8 4.9 4.7 5.0 4.8 4.8 4.6 4.8 4.8 4.7 4.8 4.8 233 1.2 1.2 1.1 1.2 1.2 1.2 1.1 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.3 1.2 1.3 1.1 1.2 1.2 1.2 1.2 1.2 197 0.7 0.7 0.7 0.8 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.8 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.8 0.7 0.7 0.7 0.8 0.7 0.8 0.7 0.7 0.7 0.7 0.7 0.7 119 50.0 50.2 50.1 50.1 50.2 50.4 50.2 50.2 50.1 50.3 50.3 50.3 50.4 50.4 50.3 50.4 50.4 50.0 50.0 50.2 50.2 50.2 50.8 50.4 50.0 50.4 50.1 49.9 50.3 49.9 49.5 50.4 50.3 50.4 50.0 49.9 50.4 119 50.0 50.4 50.1 50.2 50.3 50.5 50.3 50.0 50.3 50.5 50.5 50.5 50.5 50.8 50.6 50.8 50.7 50.0 50.0 50.3 50.4 50.4 50.7 50.5 50.0 50.6 50.2 49.9 50.3 50.1 49.6 50.5 50.6 50.5 50.0 50.2 50.6 119 50.0 50.6 50.2 50.3 50.5 50.7 50.4 50.0 50.3 50.7 50.7 50.6 50.7 51.0 50.8 51.0 51.0 50.0 50.1 50.5 50.5 50.5 50.7 50.6 50.0 50.7 50.3 49.9 50.3 50.2 49.7 50.6 50.8 50.7 50.1 50.4 50.8 116 50.0 51.1 49.9 51.0 51.4 51.6 50.8 50.3 50.7 51.0 50.3 50.2 50.2 50.9 50.5 50.9 50.8 50.1 50.2 51.2 51.2 51.2 50.8 51.0 51.4 50.7 50.9 50.7 50.5 49.7 49.3 50.9 50.6 52.7 51.3 50.8 51.9 115 50.0 51.2 50.1 51.2 51.5 51.7 50.9 50.4 50.8 51.2 50.5 50.3 50.3 51.2 50.7 51.2 51.0 50.1 50.2 51.4 51.4 51.4 50.8 51.3 51.6 51.0 51.0 50.8 50.6 49.7 49.3 51.0 50.7 53.1 51.5 50.9 52.1 115 50.0 51.4 50.0 51.3 51.8 51.8 51.0 50.5 50.9 51.3 50.5 50.4 50.4 51.4 50.9 51.4 51.2 50.1 50.2 51.4 51.5 51.5 51.0 51.5 51.7 50.9 51.2 50.8 50.8 49.7 49.2 51.1 50.9 53.4 51.6 51.1 52.3 116 50.0 51.3 50.0 51.3 51.6 51.8 51.0 50.8 50.9 51.4 50.4 50.3 50.3 51.3 50.7 51.2 51.2 50.1 50.2 51.4 51.2 51.2 51.0 51.4 51.8 50.8 51.2 50.9 50.7 49.8 49.3 51.3 50.8 53.3 51.9 51.1 52.3 114 50.0 51.5 50.2 51.4 51.8 51.9 51.3 50.8 51.0 51.5 50.5 50.3 50.4 51.5 50.8 51.5 51.3 50.1 50.2 51.5 51.4 51.4 51.1 51.6 51.8 50.9 51.4 51.0 50.8 49.6 49.4 51.4 50.9 53.7 52.1 51.1 52.5 114 50.0 51.7 50.3 51.5 52.0 52.2 51.4 51.0 51.1 51.6 50.5 50.5 50.4 51.6 51.0 51.7 51.4 50.1 50.2 51.6 51.6 51.6 51.2 51.8 51.9 50.8 51.5 51.1 51.1 49.8 49.3 51.4 51.0 53.9 52.2 51.4 52.8 234 50.0 50.0 48.6 49.8 49.9 50.2 49.1 49.3 49.9 49.7 49.8 49.6 50.0 49.4 49.8 49.7 49.4 49.9 49.9 49.7 49.6 49.6 50.0 49.9 50.3 49.8 50.1 49.7 50.5 49.8 50.5 49.5 49.8 49.9 49.5 50.0 49.7 233 50.0 50.2 48.3 49.9 50.1 50.5 48.7 49.2 49.9 49.8 49.7 49.5 50.2 49.5 50.0 49.9 49.6 49.8 49.8 49.9 49.8 49.8 50.1 49.9 50.8 50.1 50.2 49.5 51.2 49.5 50.7 49.0 49.9 50.0 49.5 50.3 50.0 202 50.0 50.7 48.6 49.7 50.3 50.6 48.7 49.4 50.0 50.1 50.1 49.7 51.0 49.5 50.4 50.2 49.8 49.8 49.9 50.0 50.0 50.0 50.0 50.1 51.4 50.7 50.3 49.6 51.9 49.6 51.6 48.7 50.2 50.1 49.5 50.5 50.5 234 50.0 50.1 48.7 50.0 50.0 50.0 49.2 49.5 50.0 49.7 49.5 49.8 49.9 50.0 50.0 50.1 50.0 49.9 49.9 50.1 50.1 50.1 50.0 49.9 50.2 50.0 50.1 49.7 50.4 49.8 50.0 49.1 49.9 50.0 49.7 49.7 49.9 233 50.0 50.2 48.2 50.1 50.1 50.2 48.9 49.5 50.0 49.7 49.4 49.8 50.0 49.8 50.1 50.1 49.6 49.9 49.9 50.0 50.1 50.1 50.0 49.9 50.6 50.0 50.2 49.5 50.9 49.8 50.6 48.7 50.0 49.8 49.5 49.8 50.0 212216 117007 65104 11723 7990 5809 6981 5106 3940 46675 12295 4625 62013 14096 4047 212216 117007 65104 11723 7990 5809 6981 5106 3940 46675 12295 4625 62013 14096 4047 197 50.0 50.3 48.6 50.2 50.3 50.2 49.0 49.9 50.0 49.6 49.4 49.9 50.0 50.2 50.2 50.6 49.8 49.9 49.8 50.3 50.4 50.4 50.2 49.8 50.8 50.1 50.1 49.5 50.8 49.6 50.7 48.5 50.1 50.0 49.6 49.8 50.1 – weak strong – weak strong – weak strong > 2 > 3 > 4 > 2 > 3 > 4 – – – weak weak weak strong strong strong Fold change Transcription Factor deletions (Hu, Killion & Iyer) z-score ChIP-chip & conserved motifs (MacIsaac et al.) – weak strong – weak strong – weak strong > 2 > 3 > 4 > 2 > 3 > 4 – – – weak weak weak strong strong strong Fold change Transcription Factor deletions (Hu, Killion & Iyer) z-score ChIP-chip & conserved motifs (MacIsaac et al.) a b # Edges # TFs Random Community ChIP-chip Motif cons. Gold standard AUPR of inference methods AUROC of inference methods 7 6 5 4 3 2 1 8 Regression 7 6 5 4 3 2 1 8 Other 5 4 3 2 1 Meta 6 5 4 3 2 1 Bayesian network 5 4 3 2 1 Mutual information 3 2 1 Correlation Fold improvement compared to random Gold standard used in this paper 0 2 -2 AUROC 50 55 45 YEASTRACT YEASTRACT 2528 107 1.5 2.4 1.6 2.2 2.5 2.4 2.3 1.8 2.1 2.2 1.7 1.9 1.8 2.1 1.9 2.2 2.1 1.6 1.7 2.0 1.9 1.9 2.1 2.4 2.9 1.9 3.0 2.4 1.9 1.5 1.7 4.1 2.0 3.3 3.1 5.3 2.9 2528 107 50.0 53.3 50.2 53.1 54.1 54.6 53.3 52.4 52.4 53.0 50.9 51.7 51.5 53.5 52.5 53.7 53.2 50.4 50.5 52.7 52.4 52.4 53.4 54.0 54.4 51.6 53.4 52.8 52.1 49.6 50.9 52.2 52.5 55.9 54.8 53.7 56.5 Figure S2: AUPR and AUROC values using different gold standards for S. cerevisiae 8 Nature Methods: doi:10.1038/nmeth.2016 4 Assessment of network inference methods 4.1 Performance assessment We used a standard approach established in previous editions of the challenge to evaluate network predictions.55,69,81 Briefly, we evaluate network predictions as a binary classification task (edges are predicted to be present or absent) and use standard performance metrics from machine learning, specifically precision vs. recall (PR) and receiver operating characteristic (ROC) curves.23 We further evaluate predictions statistically and compute an overall score summarizing the performance across networks. Scoring metrics are described below. Evaluation results are reported in: •Figure SS3: PR and ROC curves; •Figure SS4: Comparison of the ranking across compendia; •Supplementary Data 2: Table with the area under the curves and scores; •Supplementary Data 1: Evaluation scripts. PR and ROC curves. To score the ranked lists of interactions (the prediction format is described in Supplementary Note 1) against a binary gold standard, performance was assessed by the area under the ROC curve (AUROC, true positive rate vs. false positive rate) and the area under the precision vs. recall curve (AUPR). Expressions for true positive rate (TPR), false positive rate (FPR), precision and recall as a function of the cutoff (k) in the edge list are as follows:69,81 recall(k) = TP(k) P, where TP(k) is the number of true positives in the top kpredictions in the edge list, and Pis the number of positives in the gold standard. precision(k) = TP(k) TP(k) + FP(k)=TP(k) k, where FP(k) is the number of false positives at cutoff kin the edge list. The true positive rate is equivalent to recall and is defined as: TPR(k) = TP(k) P. The false positive rate is the fraction of negatives that are incorrectly predicted at cutoff k FPR(k) = FP(k) N, where Nis the number of negatives in the gold standard. Note that the length of the prediction lists was limited to at most 100,000 edges (some teams also submitted shorter lists). Edges not included in the list are thus effectively predicted to be absent. We extended the PR and ROC curves in an analytical way after the end of the list by assuming a random ordering of the remaining edges, as described.81 Note that AUROC values can vary by up to 7.5 percentage points when considering truncated vs. full lists of predictions (the AUPR values, on the other hand, don’t vary significantly — the reason is discussed in the legend of Figure SS3). Note that predictions for transcription factors and genes that are not part of the gold standard, i.e., for which no experimentally supported interactions exist, were ignored in this evaluation. Empirical p-values and overall scores. AUROC and AUPR values were separately transformed into p-values by simulating a null distribution for a large number (25,000) of random networks. We fit the histogram of the randomly obtained AUROC and AUPR values using stretched exponentials as previously described81 to extrapolate the distribution to values beyond the immediate range of the histogram. Note that p-values obtained this way depend on how “random” is defined. We chose to construct random edge lists by sampling edges from the submitted edge lists of the participants. By design, an “average” performance should have a p-value of approximately 0.5 for the two metrics. We call this procedure “grading on a curve” since the p-values are designed to identify the best and worst relative performance but have no absolute interpretation. Overall scores are derived from the AUPR or the AUROC for the three different networks by calculating the geometric mean of the network specific p-values. We report the negative log10 of this value as a score, or equivalently, ROC score =1 3 3 X i=1 −log10 pROCi PR score =1 3 3 X i=1 −log10 pPRi. Finally, an overall score is obtained as the mean of the AUROC and AUPR derived scores. A high score corresponds thus to low (significant) p-values. 4.2 Clustering of methods by PCA To depict similarities and differences between inference approaches via Principal Component Analysis (PCA), 9 Nature Methods: doi:10.1038/nmeth.2016 tation (ChIP) experiments than the E. coli gold standards. The physical binding of a transcription factor to the promotor region of a target gene is a required but not sufficient condition for effective gene regulation. transcription factor-target gene pairs derived from ChIP thus contain many false positives. However, in addition to strong evidence of binding, we also required the presence of strongly conserved transcription factor-binding motifs, which should filter out at least part of the false positives from the ChIP data50 (Supplementary Note 3). The lesser quality of the S. cerevisiae gold standard compared to E. coli thus cannot explain the complete absence of dependency between transcription factor and target mRNA levels. The second reason for this absence of dependency is the greater complexity of gene regulation in eukaryotes. In particular, additional layers of regulation at the posttranscriptional and chromatin levels lead to a reduced dependency of mRNA concentration between regulators and targets, and thus an increased difficulty for expressionbased inference of eukaryotic networks.39,60 Hu et al.39 found gene regulatory network inference in yeast very challenging. Despite the availability of microarray measurements for all 270 S. cerevisiae transcription factors, the performance of inference was hardly better than guessing. They discussed four potential reasons for the higher difficulty of inference in eucaryotes, namely the lack of operons, the high functional overlaps of transcription factors, transcription factors regulated on the protein rather than the transcript level as well as the fact that transcription factors regulate a high number of targets, albeit only mildly. Michoel et al.60 found that the overall performance of network inference in S. cerevisiae is quite low and that better results can only be achieved on specific selected subsystems. According to Wu et al.,88 independent measurements of the actual activities of the transcription factors would be required in order to resolve the subtle expression dependencies between transcription factors and their targets. Also, Herrgard et al.37 pointed out that a lack of correlation might be due to the fact that most transcription factors themselves are not significantly transcriptionally regulated and their expression remains at a low constitutive level. Instead, many transcription factors such as Mig1 glucose repressor in yeast are regulated by phosphorylation and localization as well as other posttranscriptional regulatory mechanisms. 5.2 Information content of different experiment types We sought to evaluate how informative different experiment types, such as time courses, drug perturbations, and genetic perturbations, are for network inference. To evaluate the information content of different experiment types across the three compendia, we first computed a weight for each individual microarray chip using a machine learning framework (described below). Subsequently, individual weights were averaged across chips of particular experiment types, such as time series, drug perturbations, gene and transcription factor knockout or overexpression experiments, as well as combinations thereof (Fig. SS10). The chip weights were computed as follows. We applied feature subset selection approaches based on machine learning classifiers to the in silico,E. coli, and S. cerevisiae expression compendia to estimate the significance of individual features (here: microarray chips) with respect to their ability to correctly identify transcriptional regulatory interactions. Thus, an individual local binary classifier61 was constructed for each transcription factor to distinguish known target genes from non-targets as defined by the given gold standards. As classifiers, both decision trees and support vector machines (SVMs) were employed. Decision trees were constructed72 by selecting microarray chips as decision nodes based on information gain, which we interpret as chip specific weights. In addition, linear SVMs17 were trained by optimizing SVM coefficients (αweights). In the case of linear SVM models, chip specific weights can be computed by an α-weighted linear combination of the support vectors. Default parameters were used for the training of both classifiers. As suggested by Mordelet and Vert,61 a three-fold crossvalidation was employed. This analysis shows that direct transcription factor manipulations, i.e., knockout/deletion or overexpression, are most informative for the inference of interactions involving this transcription factor (Fig. SS10). In comparison to these transcription factor specific perturbations, the remaining categories have on average much lower weights. The difference between weights of transcription factor specific knockouts and the other kinds of experiments is more pronounced for decision trees (that optimize chip weights) as compared to SVMs (that optimize example weights). Note that the S. cerevisiae compendium included only three transcription factor-deletion measurements for a single transcription factor (GCN4): results are therefore not as reliable as for the other compendia. Transcription factor overexpression experiments were also very rare in the examined compendia and were thus included in the same category as the knockouts in Figure SS10. We compared the information content of knockout versus overexpression experiments in a different E. coli compendium from the M3D database27 (Supplementary Note 2) and found that chip weights from transcription factor specific overexpression are comparable to the knockout weights. We conclude that direct transcription factor perturbations (knockout/overexpression) are the most informative experiments on average. However, other kinds of experimental conditions can be similarly informative for the inference of 16 Nature Methods: doi:10.1038/nmeth.2016 0 0.5 1 1.5 0 0.02 0.04 0.06 In silico E. coli S. cerevisiae Decision trees Average weight Average weight SVMs Transcription factor KO/OE Other experiments 654321 7 8 Transcription factor KO/OE Transcription factor KO/OE Other experiments 654321 7 8 Other experiments 654321 7 8 Experiment types 1 Time series 2 Drug perturbation 3 Drug perturbation, time series 4 Gene KO/OE 5 Gene KO/OE, time series 6 Drug + gene KO/OE 7 Drug + gene KO/OE, time series 8 Transcription factor KO/OE KO/OE = knockout/overexpression 18.5 Figure S10: Information content of different experiment types for network inference Feature selection was employed to estimate the utility of microarray chips from different experiment types to correctly identify gene regulatory interactions. Microarray measurements were categorized into time series, drug perturbations, gene knockout/overexpression (KO/OE) experiments, as well as combinations thereof. The bar at position 8 only refers to the weight of transcription factor knock out/over-expression experiments for inferring interactions involving that transcription factor. Knockout or overexpression of a transcription factor is by far the most informative experiment to detect regulatory interactions, independently of the classifier and compendium used. Note that in the S. cerevisiae compendium there were only three transcription factor deletion experiments available (all for the same transcription factor), results are thus not as reliable as for the other two compendia. the targets of certain transcription factors (Fig. SS11), but here it will be much more difficult to decide a priori which kind of conditions that may be. 6 Integration of predictions 6.1 Why community integration can outperform the best individuals Of the many ways to integrate results from an ensemble of predictions,47 we have chosen the Borda count election method. It was originally developed by 18th-century political scientist Jean-Charles de Borda as a method to select candidates in a democratic election.12 In this method, voters rank candidates in order of preference, and the winner of the election is the candidate with the best average rank. Similarly, DREAM5 participants provide a ranked list of transcription factor-target gene predictions. Figure S12a exemplifies three possible prediction lists where we denote transcription factor-target gene interactions by the letters A, B, C, . . . , X, Y, Z. For example, Team 2 has high confidence that interaction Ais a true interaction, but Team 1 and Team 3 are less certain about the truth of this interaction, having placed it at position 5 and 20 of the ranked lists, respectively. The same is true at the other end of the list. Team 2 considered Yto be a very unlikely interaction, but Team 3 located it at position 10 of its ranked list. In this section, we wish to build some intuition to understand why the integration of predictions can outperform individual teams. Let us assume that we have P true transcription factor-target gene interactions (the positives) and Nnon-interacting transcription factor-target gene pairs (the negatives). In total, there are T=N+P possible predictions. Different applications of the same network reconstruction algorithm that attempt to infer a given network can create lists of predictions in which the same interaction Ican be placed in different rank positions. This will depend on the set of experiments used by the algorithm, random aspects of the method itself, biological variability, etc. Therefore, each method can be characterized by the probability ppos(r, I)that it places a given true interaction Iat rank r. We make the simplifying assumption that this probability is the same for all interactions (i.e.,ppos(r, I)is independent of I), and refer to this probability as ppos(r). By randomly ranking interactions, we expect ppos(r) = 1/T. However, if a method has better than random predictive power, there 17 Nature Methods: doi:10.1038/nmeth.2016 Figure S11: Significance of individual microarrays for network inference The heatmap depicts the internal node weights of decision trees trained to predict the target genes of E. coli transcription factors. A single weight thus represents the importance of an individual microarray for the detection of known target genes for a given transcription factor. If a microarray of a transcription factor deletion mutant is available (transcription factors marked by *), it often receives particularly strong weights (cyan arrows). Even tough other types of experiments have lower weights on average (Fig. SS10), they can be similarly informative in some conditions (bright spots scattered through the matrix). Hierarchical clustering was employed to sort rows and columns of the heatmap. will be some tendency for the method to have an enhanced probability of detecting true interactions at the top of the ranked predictions, such as: ppos(r) = 1 T+b P(1 −2r−1 T−1), where we introduce a parameter that accounts for bias, b. Note that if b= 0, the probability reduces to that of the random prediction method. There is also a corresponding probability of detecting a non-interacting transcription factor-target gene pair (a negative) at rank position r,pneg(r). Note that, as each ranked position in the prediction list must contain either an interacting or noninteracting transcription factor-target gene pair (a positive or a negative), the condition P·ppos(r) + N·pneg(r)=1 holds, and therefore: pneg(r) = 1 T+b N(1 −2r−1 T−1). In order for these probabilities to be positive, bhas to be smaller than the smaller of N/T and P/T.Figure S12b shows the results of these calculations where P= 30,N= 70, and b= 0.2. It can be shown that the average area under the precision recall curve (AUPR) for a method following these probability laws is 0.41 (for a random prediction we expect AUPR=P/T =0.3). Let us now explore what happens to the average rank of a true interaction if we integrate the predictions of a community of Kinference methods. The average rank assigned to a possible transcription factor-target gene interaction I, over the predictions of the Kinference methods, is computed as rBorda(I) = 1 K K X j=1 rj(I). As an example, the true interaction Ain Figure S12a has average rank 8.66, whereas non-interaction Zhas average rank 19.67. If a method has better than a random probability to predict true interactions, then the average rank for a true interaction will be different from the average rank for a transcription factor-target gene pair that doesn’t interact, given that the former is computed using the average over ppos(r)and the latter is computed using pneg(r). We will assume that all the teams in this commu18 Nature Methods: doi:10.1038/nmeth.2016 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 Team 2 Rank A B C X Y Z Team 1 B C A X Y Z Team 3 A B C X Y Z Positives (interacting gene pairs) Negatives (non-interacting pairs) a 0 0.02 0.04 Probability 0 0.02 0.04 Probability 1 20 40 60 80 100 0 0.04 0.08 Probability Average rank 1 method Integrating 5 methods Integrating 30 methods AUPR = 0.41 AUPR = 0.66 AUPR = 0.97 Interacting gene pairs Non-interacting gene pairs b c d Figure S12: Why community integration can outperform the best individual inference methods (a) A hypothetical example of predictions submitted by 3 separate teams. The challenge is to integrate predictions made by each team into a single ranked list. (b-d) Two sufficient conditions for integration to outperform individual inference methods are: (1) each of the inference methods must have better than random predictive power (i.e., on average interacting pairs are assigned better ranks than non-interacting pairs); and (2) predictions of different inference methods must be statistically independent.24 For illustrative purpose, we consider a simple scenario comprising T=100 candidate transcription factor-gene pairs, out of which P=30 interact and N=70 do not interact. For instance, individual inference methods that assign ranks to interacting and non-interacting pairs with the probabilities shown in Panel b suffice condition (1). Although interacting pairs are assigned better ranks on average, the probability of incorrect predictions (non-interacting pairs in the top-part or interacting pairs in the bottom-part of the prediction list) is considerable in this example, resulting in an AUPR of only 0.41 (for a random prediction, we expect AUPR=P/T=0.3). (c,d) If the assigned ranks are independent across methods (condition 2), the central limit theorem establishes that the average rank distribution will approach a Gaussian distribution. Its variance shrinks as more methods are integrated, thereby increasingly segregating interacting from non-interacting transcription factor-gene pairs (Panels b→d). Consequently, the probability that interacting pairs are ranked better than non-interacting pairs increases, resulting in an AUPR that tends to 1 (perfect prediction) as the number of integrated inference methods increases. Even though all predictions are based on the same datasets in the DREAM challenge, and are thus not statistically independent, we have shown that diversity arises due to method-specific capabilities to extract different kinds of information from the data (Fig. 2). Methods from different classes show greater levels of independence and thus contribute more to community performance (Figs. 2b and 3c). Since predictions are partially, but not completely independent, the AUPR increases as more methods are integrated (Fig. 3a) but tends to a value lower than 1 in practice. 19 Nature Methods: doi:10.1038/nmeth.2016 nity have the same ppos(r)and pneg(r). The distribution of the average ranks for interacting and non-interacting pair pairs are shown in Figure S12b-d for an individual inference method and the integration of 5 and 30 inference methods, respectively. As the number of integrated methods increases, a true interacting pair is more likely to have a better rank than a non-interacting pair. As this is true for any interacting and non-interacting pair, then the interacting segregate from the non-interacting pairs. This feature of integration results in a larger AUPR curve as the number of included methods increases (e.g., AUPR=0.66 for 5 teams and 0.97 for 30 teams). It is clear then that the aggregation of even 5 teams outperforms (according to the AUPR metric) each of the participating members of the aggregate, whose typical AUPR is 0.41. In order for the integration to outperform individual predictions, the methods being integrated need to be statistically independent, that is, the rank where an interaction is placed by a method has no statistical dependency on the rank where the same interaction is placed by any other method. If this assumption holds, then the central limit theorem of probability theory establishes that, as more predictions are averaged, the average rank distribution will approach a Gaussian distribution whose variance shrinks as the number of integrated teams increases. If the integrated methods have some predictive ability, then the mean ranks of the interacting transcription factor-target gene pairs will be better than the mean ranks of of the noninteracting transcription factor-target gene pairs, and the variances will eventually shrink to make all the interactions have average ranks that cluster tightly around their respective means. In our example it can be shown that more than 95% of the true interactions will have an integrated rank in the interval [T 2−bT2 6P−T √3k,T 2−bT2 6P+T √3K], whereas more than 95% of the non-interactions will have an integrated rank in the interval [T 2+bT2 6N−T √3k,T 2+bT2 6N+T √3K]. Thus, when K(the number of members in the community) is large enough, the 95% intervals will cease to overlap, and 95% of the positives will be ranked above 95% of the negatives, producing excellent precision and recall. It is clear that as Kincreases, having a positive in the interval where the negatives concentrate will be extremely unlikely. That is, under the assumption of independent predictions, the area under the precision recall curve (AUPR) tends to 1 as the number of integrated methods increases. In the real world scenario of the present DREAM5 network inference challenge, predictions from different methods are indeed considerably different (as shown in Figs. 2b,c in main text) but can obviously not be fully independent as all predictors use the same input data. Thus, as predictions are partially, but not completely independent, the AUPR increases when more predictors are integrated (see Fig. 3a in main text) but tends to a value lower than 1. 6.2 DREAM5 community networks All 29 DREAM5 submissions were used to construct the community-based transcriptional regulatory networks. For each compendium, integration of the individual team predictions into a single community network was done using the Borda count method described in the previous section. The resulting community-based networks consists of the reordered lists of transcription factor-target gene pairs. Each team was requested to submit a total of 100,000 predictions. Interactions not listed in the top 100,000 predictions were assigned a rank of 100,001. Upon completion of applying Borda’s method, and for each dataset (in silico, E. coli,S. cerevisiae,S. aureus), the top 100,000 predictions were selected and called the community predictions (Supplementary Data 3). Weighted voting. Borda count voting amounts to an unweighted rank average over the individual predictions of an interaction. To explore how the community predictions are affected when combining only the best-performing methods and/or giving methods with a better performance a higher weight, we tested several weighted voting schemes. We stress that these weighted voting methods are only applied to build an intuition of how the performance of community predictions is affected by good/poor predictors. Weighted voting cannot be used in practice when inferring an unknown regulatory network, as in the case of S. aureus here, because the performance of the inference methods is not known. The weighted average rank assigned to a possible transcription factor-target gene interaction I, over the predictions of the Kinference methods, is computed as r(I) = 1 PK j=1 wj· K X j=1 wjrj(I), where wjis a measure of performance of method j(e.g., the AUPR). To gain a sense of the performance of unweighted (Borda count) and weighted community predictions, we systemat20 Nature Methods: doi:10.1038/nmeth.2016 1 5 10 15 20 25 29 0 20 40 60 80 100 Team Overall score All three compendia 1 5 10 15 20 25 29 0 10 20 30 40 Team AUPR (%) In silico 1 5 10 15 20 25 29 0 5 10 15 Team AUPR (%) E. coli 1 5 10 15 20 25 29 0 1 2 3 4 Team AUPR (%) S. cerevisiae Teams Rank-average Confidence-avg. Weighted rank-average Weighted confidence-avg. Community networks: {1,2}, {1,2,3}, {1,2,3,4}, ...Individual methods: 1, 2, 3, ... Figure S13: Community integration using unweighted and weighted voting Community predictions were obtained by combining the two best teams, the three best teams, the four best teams, etc, using either unweighted (Borda count) or weighted voting. In addition to voting based on the edge ranks (diamonds), we also tested voting based on the edge confidence values (third column in the prediction format, Supplementary Note 1) assigned by the inference methods (squares). In accordance with Figure 3d of the main text, even unweighted voting is robust to inclusion of poor predictors. Although weighted voting based on edge confidence values performed best overall, the difference with the other approaches is relatively small on the three individual compendia. We stress that only the unweighted voting that combines all methods at hand (rightmost points marked with arrows) is truly unsupervised, i.e., can be applied when inferring an unknown regulatory network, where the performance of the individual methods is not known a priori. ically formed communities composed of the top two methods, the top three methods, the top four methods, etc., until the last community, which contains all 29 methods applied by the participants of the challenge (Fig. SS13). This analysis confirms an observation made in the main text (Figure 3d): adding poor predictions hardly degrades the consensus of the more accurate predictions, even when using unweighted voting. Although weighted voting improves the performance overall, the difference is rather small on the individual compendia. Therefore, and since the performance of inference methods is difficult to estimate when inferring an unknown regulatory network, integrating all inference methods at hand using unweighted voting seems to be a good choice. 7 E. coli and S. aureus community networks 7.1 Network construction Community predictions for E. coli and S. aureus were obtained using the Borda count method as described in the previous section. Note that these community predictions are weighted networks, as they assign a measure of confidence (the average rank) to edges. To obtain an unweighted network that classifies edges simply as present or absent, a confidence threshold must be chosen. Edges above the threshold are considered present, and those below absent. Choosing a threshold amounts to a trade-off between sensitivity and specificity. For the analysis of the E. coli and S. aureus community networks in the main text, we chose a cutoff of 1,688 edges, which corresponds to an estimated precision of 50% for E. coli. Precision was estimated based on the RegulonDB gold standard of interactions described in Supplementary Note 3. Note that many genes and (potential) transcription factors still have no experimentally supported interactions in RegulonDB, i.e., they are not part of our gold standard. Edges involving genes or transcription factors that are not part of the gold standard were ignored when computing the precision (i.e., they were neither counted as true nor as false positives, instead they were simply excluded from the calculation). For example, at the cutoff of 1,688 edges, only 200 edges are part of the gold standard. Of these 200 edges, 100 were true positives and 100 were false positives, resulting in the estimated precision of 50%. Note that this is a conservative estimate, because some of the novel edges that are considered false positives may in fact be newly discovered regulatory interactions that are currently missing in RegulonDB, as our independent experimental validation of such novel interactions shows (Figure 4 of the main text). As discussed, there is no gold standard set of interactions that exist for S. aureus, which presents the problem of not being able to directly ascribe a precision with a given nework size. In this case, we make the assumption that the S. aureus network predictions perform similarly to the E. 21 Nature Methods: doi:10.1038/nmeth.2016 0 500 1000 1500 2000 2500 3000 3500 4000 Community worse Community similar Community better Regression Mutual information Correlation Bayesian networks Other Meta Inference methods Number of predicted interactions Figure S14: Rank improvement of individual edges through community integration The plot shows the number of true interactions that were predicted worse, similarly, or better by the integrated community network than by the individual inference methods (the average over all methods is shown in grey). Interactions were counted as predicted better or worse than the individual methods if the difference rank(integrated)−rank(methodi) exceeds 10,000 or -10,000 ranks, respectively. Otherwise, the predictions were considered similar. The majority of true interactions were ranked better (by a large margin of over 10,000 ranks) in the community network than even in the best individual predictions. Note that the ordering of methods is the same as in the main document (Table 1 and Figure 2). coli community network. Under this assumption, we use the network size of 1,688 edges since this is the network size derived from the E. coli network at an estimated 50% precision. 7.2 S. aureus network evaluation using RegPrecise RegPrecise is a database of transcriptional regulatory interactions in prokaryotes that have been inferred using manually reviewed, homology based methods.63 For the DREAM challenge evaluations, we required all gold standard interactions to be experimentally supported and did not consider electronically inferred interactions. As the interactions in RegPrecise for S. aureus are largely electronically inferred, we did not consider RegPrecise as a gold standard for the overall evaluation. Nonetheless, the RegPrecise interactions represent a rich set of information that we have used to test our assumption that the 50% precision threshold derived in E. coli can be used as a proxy for performance in S. aureus. RegPrecise contains 517 interactions comprised of 38 transcription factors and 446 target genes that match up with the genes represented in the microarray compendium supplied in the DREAM5 network inference challenge. All individual methods and the community network performance were evaluated using AUPR and AUROC (Supplementary Note 4.1). We performed the same analysis using RegPrecise interactions and report the performance in S15. Using the AUPR, we find that the S. aureus community network ranks 3rd. The 50% precision reported for E. coli correponds to 1,688 interactions. We identified 50% precision in the S. aureus community network from the RegPrecise analysis and this threshold corresponds to a network with 988 interactions. For the E. coli gold standard interactions reported in RegulonDB, there were 2,066 experimentally supported interactions comprised of 144 transcription factors and 999 target genes. Considering that the E. coli gold standard is 4 times the size of RegPrecise and only experimentally supported interactions from RegulonDB were used, we feel that the RegPrecise interactions understimate the number of true positive relationships. Therefore, we report the S. aureus community network using the 50% threshold in E. coli and perform further analyses on this network of 1,688 edges. For completeness, we provide the S. aureus community network using the RegPrecise 50% precision threshold 22 Nature Methods: doi:10.1038/nmeth.2016 False positive rate 0.75 0.5 0.25 0 1 Precision 0.75 0.5 0.25 0 1 True positive rate Individual methods Best indiv. method Community network Random Regression 3 0 0.25 0.5 0.75 1 Recall 0 0.05 0.1 0.15 0.2 0.25 Correlation 2 Regression MI Corr. Bayesian Other Meta Community Random 1 2 6543 7 8 1 2 543 1 2 3 1 2 6543 1 2 6543 7 8 1 2 543 0 5 10 15 AUPR (%) a b Figure S15: Performance evaluation for S. aureus using the RegPrecise database (a)The Area Under the Precision Recall (AUPR) for all network inference algorithms (including the community integration) evaluated for S. aureus using the RegPrecise database as the “gold standard.” (b) The precision/recall (PR) and receiver operating characteristic (ROC) are show with the community network highlighted in red and the top performing algorithm in black. The remaining DREAM participant sumbissions and the off-the-shelf predictions are shown in grey. (988 edges) in Supplementary Data 5. 7.3 Analysis of network modules Module detection. We identified network modules, i.e., groups of transcription factors and genes that are more densely connected among themselves than expected in a randomized network with the same degree distribution, using Newman’s spectral method (including the greedy optimization step after the spectral decomposition).62 We found that both the E. coli and S. aureus community networks are highly modular, as shown in Figures 4a and bof the main text for the two community networks at the 50% precision cutoff (1,688 edges). Gene Ontology (GO) term enrichment. GO term enrichment analysis was performed on each of the identified network modules for both the E. coli and S. aureus networks. GO term gene annotations were downloaded from the Gene Ontology websitemfor E. coli. For S. aureus, gene annotations were taken from the Affymetrix annotation files.nGO terms under the biological process branch of the ontology were used. Genes were mapped directly to the ontology and propagated to the root node. This process ensures that all parent GO terms recursively inherit the annotations of their child terms. GO terms annotated with less than 3 genes and GO terms annotated with greater than 500 genes were removed. The remaining GO terms were used as input to the GO term enrichment calculation. Each identified network module can be represented by a set of genes. For each module, all associated GO term annotations were tested by counting the number of occurrences of the term in comparison to the number of occurrences of the term in the entire network. Statistical significance for each GO term was assessed using the hymwww.geneontology.org/GO.downloads.annotations.shtml nwww.affymetrix.com/support/support_result.affx 23 Nature Methods: doi:10.1038/nmeth.2016 pergeometric distribution. Estimated p-values were then multiple hypothesis corrected using the q-value calculation.82 We found that the network modules—which were identified solely based on network connectivity—are also strongly enriched for specific biological processes, i.e., network modules coincide with functional modules (Fig. SS16 and Fig. SS17 show module enrichments for the networks at the 50% precision cutoff, similar results were obtained at different cutoffs). The genes assigned to each module, as well as the p-values and q-values for the GO terms, are supplied in Supplementary Data 6. 8 Experimental validation All cultures were grown in 1 mL of indicated growth medium in 14 mL Falcon tubes. Incubation of all cultures was performed in darkened shakers (300 RPM) at 37 ◦C. Experimental conditions for growth and transcription factor induction were all designed to reproduce the conditions under which the induction of these transcription factors have previously been studied. The wild-type E. coli strain used was BW25113. Knockout strains were constructed via transduction from the KEIO knockout library.6 8.1 Transcription factor selection Transcription factors were selected from the E. coli community network with 50% predicted precision or greater. Selected transcription factors had at least 8 predicted, but untested target genes. The conditions under which these transcription factors are active are known and can be replicated during aerobic growth under laboratory conditions. Targets were chosen as follows. For each transcription factor, if confirmed target genes were available in the data set, 1-2 of these geness were chosen as positive controls. All unconfirmed target genes were used in qPCR unless suitable primers could not be obtained for the sequence (see Supplementary Note 8.2 for primer design specifications). Where multiple targets were encoded within a single operon, only the first gene in the operon was used. rhaR is the transcriptional activator of the rhamnose utilization operon. Native production of rhamnose activated genes is extremely low in the absence of a chemical inducer, and induction by rhamnose is slow, requiring 40-50 minutes to reach steady state.26 Overnight LB cultures of both wild-type and ∆rhaR E. coli were inoculated 1:500 in minimal salt media (M9) + 0.2% casamino acids + 50 µM thiamine + 0.4% glycerol. Cultures were grown to early exponential phase (A600 ≈0.2, 3.5 hours) before addition of 0.2% (w/v) L-rhamnose26 rhamnose and were incubated 45 minutes before stabilization for RNA extraction. RNA was also extracted from untreated samples of wild-type and ∆rhaR grown under the same conditions but without the addition of rhamnose. As rhaR is cotranscribed with its confirmed target rhaS, to avoid any problems with qPCR due to possible scarring at the C-terminal end of the rhaS transcript, the rhamnose transporter gene rhaT was used as a positive control. purR is active under conditions of purine nucleotide deficiency. Overnight LB cultures of wild-type and ∆purR E. coli were inoculated 1:500 in minimal salt media (M9) + 0.2% casamino acids + 6.6 µM thiamine + 0.4% glucose, with 100 µg/mL adenine added to activate purR and repress purine nucleotide biosynthesis.18,73 Cultures were grown to exponential phase (A600 ≈0.2, 3.5 hours) with or without adenine, as previously described in,18,36 before stabilization for RNA extraction. RNA was also extracted from untreated samples of wild-type and ∆purR grown under the same conditions but without the addition of adenine. gadE is a central transcriptional activator of the principal acid resistance system.38 Overnight LB cultures of wildtype and ∆gadE E. coli were inoculated 1:500 in minimal salt media (M9) + 0.2% casamino acids + 0.4% glucose + 30 µM thiamine (pH 7). Cultures were grown to exponential phase (A600 ≈0.2, 3.5 hours) before adjustment of pH to 5.4-5.7 by addition of 45 µL 1M HCl. Cultures were incubated for an additional 2 hours before stabilization for RNA extraction.16 RNA was also extracted from untreated samples of wild-type and ∆gadE grown under the same conditions but without the pH adjustment. mprA (emrR) is known to respond to toxic molecules such as salicyclic acid (5 mM), 2,4-dinitrophenol (DNP, 0.5 mM), carbonyl cyanide m-chlorophenylhydrazone (CCCP, 10 µM), and carbonyl cyanide p-(trifluoromethoxy)phenylhydrazone (FCCP).48 This transcription factor has been shown to interact directly with DNP, CCCP and FCCP, which reduces the capacity of mprA to bind to DNA.15,89 Overnight LB cultures of wild-type and ∆mprA E. coli inoculated 1:500 in LB and allowed to grow to exponential phase (A600 ≈0.3, 2.5 hours) before addition of the trancriptional inducer, 10 µM CCCP. Cultures were incubated 30 minutes before stabilization for RNA extraction. RNA was also extracted from untreated samples of wild-type and ∆mprA grown under the same conditions but without the addition of CCCP. cueR is a metal-binding transcription factor and the regulator of the primary copper homeostasis system in E. coli.83,90 Overnight LB cultures of wild-type and ∆cueR E. coli were inoculated 1:500 in minimal salt media (M9) + 0.2% casamino acids + 0.4% glucose and grown to ex24 Nature Methods: doi:10.1038/nmeth.2016 Modules 1 141312111098765432 15 2019181716 21 22 23 P-value <10-6 >10-2 10-4 GO terms Figure S16: Functional enrichment of network modules in E. coli Network modules are strongly enriched for very specific biological processes, with only few processes being enriched across more than one module. Only network modules that comprise at least five genes are shown. 25 Nature Methods: doi:10.1038/nmeth.2016 10.3 Regression 3 – Sparse piecewise linear regression based on changepoint processes The idea behind this approach is to apply a Bayesian piecewise linear regression model based on a multiple changepoint process to all genes, where the explanatory variables are taken from the set of candidate regulatory genes provided by the DREAM5 network inference challenge. The multiple changepoint process acts on a linear ordering of the chips, resulting from a preprocessing step that takes exogenous information about experimental conditions and perturbations into consideration. Inference is based on reversible jump Markov chain Monte Carlo (RJMCMC). Owing to the high computational costs, a preprocessing step based on L1-regularized linear regression (Lasso) is included. The individual steps are discussed in more detail below. Data preprocessing. All gene profiles were standardized to zero mean and unit variance. Chip ordering. The multiple changepoint processes (described below) are based on a chip ordering. This ordering is based on experimental conditions and perturbations, and has been obtained as follows: (1) Chips belonging to the same experimental conditions and perturbations were grouped together, and an average expression profile for that group was computed. (2) The first principal component in the space of average expression profiles was computed and used to initialize a one-dimensional selforganizing map (SOM).p(3) The SOM learning algorithm was applied to obtain an ordering of groups. (4) Within each group, chips were ordered by the same process: Initialization of a one-dimensional SOM with the first principal component, followed by the application of the SOM learning algorithm. (5) The final chip ordering was given by the ordering thus obtained, subject to a manual correction to enforce the natural ordering of time series and perturbation dilution series as a rigid constraint. Active interventions. Active interventions, like gene knockouts, were dealt with in the regression model by presenting inferred values only for the explanatory variables, but removing them for the target variable. Gene filtering. Due to the high computational costs of the RJMCMC simulations, we applied a filtering step based on the modified Lasso approach proposed by Ahmed and Xing,2which is implemented in the software package TESLA. This method can be regarded as a piecewise linear sparse regression approach that is based on L1-regularization with group-specific regression parameter vectors. In addition to the standard L1-norm penalty pWe used the R package som available on CRAN (http:cran.r-project.org) term an additional L1-norm regularization term was introduced, which penalizes deviations between vectors of regression parameters associated with different groups. The groups correspond to different experimental conditions and perturbations, as described above. The regularization constants were set so as to optimize the BIC score. We modified the original TESLA code to switch from logistic to linear regression. However, for lack of time we did not reprogram the function for computing the likelihood, and hence computed the BIC score with the original code based on logistic regression. For each gene, potential regulators were ranked according to the modulus of the regression parameter. The 20 highest-ranked regulators (corresponding to about 10% of of all regulators) were kept for the follow-up analysis. Bayesian piecewise linear regression model. We adapted the model presented in45 to create a piecewise linear regression model for the regulation of each gene, given its transcription factors. This model can learn the structure of the underlying regulatory network, as well as changepoints in the linear ordering of chips which separate different experimental conditions or perturbations. We incorporated sparse Poisson priors on the number of changepoints and the number of potential regulators for each gene. In contrast to,45 this model assumes that the structure of the transcription factor-target gene interactions do not change with different experimental conditions, only the parameters associated with these interactions do. We learned the structure and parameters using RJMCMC. RJMCMC simulations. The RJMCMC simulations for inferring changepoints and parameters of the Bayesian regression and multiple changepoint model were run on two high-performance computer clusters. For data set 2, simulations were started from two different initializations, and scatter plots of the marginal posterior probabilities of the edges were obtained for monitoring convergence. For 2×106RJMCMC steps these scatter plots indicated sufficient convergence, but it turned out that these simulations could not be finished in time. We therefore had to reduce the simulation lengths to 105RJMCMC steps, despite the fact that the scatter plots indicated a certain lack of convergence. The method could only be applied to data sets 1 and 2 for lack of time. Submitted edge orderings. The submitted edge orderings were obtained as follows. For each target node, a set of 20 potential regulators passed the filter based on TESLA, corresponding to roughly 10% of the potential incoming edges. For these regulators, the original intention was to order the edges on the basis of the marginal posterior probabilities from the RJMCMC simulations. However, owing to hardware-related problems with the computer cluster (see below), the RJMCMC simulations did not reach the desired convergence level by the DREAM 32 Nature Methods: doi:10.1038/nmeth.2016 deadline. To guard against false negatives resulting from potential entrapment of MCMC trajectories in metastable states, we computed the ranking from the arithmetic mean of the marginal posterior probabilities and the modulus of the corresponding regression coefficients, where the latter were rescaled to the unit interval. After ranking all the node-regulator pairs for those regulators that had passed the filter, the remaining edges were ranked on the basis of the modulus of the regression parameter obtained with TESLA. For those data sets for which no RJMCMC simulations could be run (see below), the rankings were solely based on the modulus of the regression parameters from TESLA. Discussion This method was among the top quarter of performers by overall performance, even though we had to grapple with some difficulties, which were as follows: Data preprocessing. We spent substantial time on the preprocessing and the chip ordering, where we made various worrying observations. The correlation of the regression parameters obtained with Lasso from data subjected to different preprocessing schemes were rather modest, indicating that the preprocessing of the data matters a lot. Without explicit knowledge about the biological processes it is difficult to decide on the appropriate choice, and in the absence of this information we opted for a standardization to transform all gene profiles to zero mean and unit variance. Chip ordering. In this approach, information about experimental conditions and perturbations resides in the chip ordering. It turned out that a one-dimensional ordering cannot group both the same experimental conditions and the same perturbation types together. This suggests that a two-dimensional multiple changepoint model would be a more appropriate approach. However, in the available time this model could not be implemented and tested, and we therefore decided to stick to the one-dimensional approach. Preliminary analysis based on TESLA showed poor correlation between the regression parameters obtained from different chip orderings. This indicates that the chip ordering matters a lot, but without extra biological knowledge it is impossible to decide what matters most: the experimental conditions, or the perturbations. RJMCMC simulations. The net time available for the RJMCMC simulations - following the preprocessing steps and associated investigations - was about two weeks. Unfortunately, external circumstances beyond the control caused extended downtime of one of the computer clusters. As a consequence, the RJMCMC simulations could only be completed for two out of the four data sets, and even here the convergence was sub-optimal. Lack of consistency between the methods. While we expect that the Bayesian regression and multiple changepoint model should achieve an improvement one TESLA, we were surprised by how low the edge rank correlations between the two methods were. This might be indicative of a certain paucity of true patterns in the data, possibly as a consequence of the choice of the pre-processing and chip ordering schemes. Future Work. The proposed method critically depends on the pre-processing step, which is rather heuristic. Besides looking into a more principled alternative to SOMs for obtaining a chip ordering, especially such that explicit biological knowledge is included, we are working on extending the changepoint process to more than one dimension. As discussed above, this would allow groups of perturbations in addition to experimental conditions to be kept together, and we would expect that to be reflected in an improved network reconstruction accuracy. 10.4 Regression 7 – Simple L1regularization In a graphical model, an edge between two variables means that these two quantities are conditionally dependent, given the remaining variables. In this case, these variables are genes, and for each gene we would like to find the transcription factors that (directly) influence a respective gene. This means we are looking for those target genes and transcription factors which are conditionally dependent, given the other transcription factors. For that purpose, target genes are regressed on the transcription factors. Non-zero regression coefficients indicate (conditional) dependencies, whereas (conditional) independence is indicated by zero coefficients. To identify zero and non-zero coefficients, we use L1-type regularization techniques. Data preprocessing. In order to detect sets of genes sharing common expression patterns we employed a clustering approach consisting of two steps: (1) computation of the pairwise similarity between genes and (2) the detection of clusters. As a similarity score (cf. (1)) we employed the Spearman correlation on the given raw data (i.e. using a data point per chip and gene from the provided data matrices). The pairs of genes with a Spearman correlation of above 0.8 were subjected to Markov clustering (cf. (2)) using an inflation parameter of 3.0. The data provided by DREAM5 shows a noticeable clustering structure. The identified clusters show very little between-gene variability suggesting to treat the gene clusters as a single entity when used as the response variable in the regression model (see below). L1-Regularization using the Lasso. We used an adaption of the method proposed by58 who dealt with estimating 33 Nature Methods: doi:10.1038/nmeth.2016 high-dimensional (undirected) graphs under the assumption of a multivariate normal distribution. In this special case, the method can be seen as an approximate version of a method proposed by.30 The latter approach estimates the inverse covariance matrix (P−1) of a high-dimensional multivariate normal distribution, with the sum of absolute values of the elements of this matrix being penalized. On one hand, the use of a penalty makes the matrix estimable; on the other hand, the chosen L1-type penalty causes that some elements of the estimate of P−1were set to zero. Zero elements (P−1)ij in the inverse covariance matrix mean that the variables (or genes in this case) iand jare conditionally independent given the remaining variables (genes). In the situation given in the DREAM5 competition, however, the method30 cannot be used, since no independent and identically distributed samples from a multivariate normal distribution are given. In addition, the estimated graph should be directed. In,58 by contrast, the authors try to find conditional dependencies by regressing each gene separately on the remaining ones. A similar procedure can be applied to the DREAM5 data. A zero regression coefficient indicates that the target gene (the response) is not influenced by the respective explanatory gene (the regressor), given the other genes. To locate zero coefficients, a Lasso84 penalty was used, so regression coefficients may be set to zero. Non-zero coefficients indicate conditional dependencies, and edges can be drawn from the corresponding explanatory genes to the target. This procedure was repeated with each gene serving once as the response. By construction, the resulting graph is directed. Since only transcription factors may regulate other genes, only transcription factors are considered as (potential) regressor genes (but all genes may serve as response). If any of the (potential) effector genes are deleted or over-expressed, this information is taken into account by adding indicator variables to the set of explanatory genes. If such an indicator is selected by the Lasso, an edge between the corresponding effector gene and the considered target is drawn as well. Also permutations can be taken into account by using indicator variables. If data come from time series, the considered target gene is also regressed on lag one of the transcription factors. Obtaining confidence scores. We did not use stability selection59 to compute confidence scores, but a more simple procedure. The Lasso, which is applied to the single regression problems to select variables (that is, effector genes), depends on a tuning parameter t, which determines the strength of penalization. More precisely, the residual sum of squares is minimized as a function of regression parameters b1, . . . , bp, subject to |b1|+. . .+|bp|< t. The smaller t, the higher the penalty and the less variables are selected. To derive a measure of confidence, we used the value of the tuning parameter where a considered edge is selected the first time. After normalization with respect to the most/least reliable edges (out of the first 100,000), these scores have values between zero and one. Discussion. Though L1-type regularization seems promising for the selection of edges in a gene regulatory network, the performance of this approach was rather bad on the DREAM5 data. This may be for two reasons. First, selection patterns of single regression models are rather unstable and hence less reliable. Second, it is doubtful that tuning parameters from different regression models can be directly compared. That means it is difficult to rank potential edges selected in regression models with different target genes. Tough stability selection is computationally much more expensive than the simple approach directly using the Lasso tuning parameter, it is apparently more reliable when judging on the relevance of edges in the network. 10.5 Regression 8 – Linear regression* Linear Regression is one of the 6 commonly-used, off-theshelf algorithms run by the DREAM organizers. A full description of this method can be found in28 and is made available online.qThis method was implemented as an attempt to estimate a Bayesian network. The computational challenge underlying Bayesian networks for gene regulatory network inference is to exhaustively search the space of possible regulatory relationships. Given the size of the DREAM5 datasets and number of regulatory relationships to be inferred a true implementation of a Bayesian network is intractable. By constraining the number of regulatory relationships considered and simplifying the mathematical model to a linear regression, a regulatory network can be inferred from the DREAM5 data. Given the underlying linear regression model and algorithm clustering results shown in Figure 2, this method was placed in the Linear Regression group. 10.6 Mutual Information 1 – Context Likelihood of Relatedness (CLR)* CLR is one of the 6 commonly-used, off-the-shelf algorithms run by the DREAM organizers. A full description of this method can be found in28 and is made available online.rThe input into the CLR algorithm is a gene ×condition matrix of expression values. The CLR algorithm progresses through 2 main steps. First, a matrix of mutual information values is calculated for all input pairs of genes iand j. As the authors suggest, mutual information was calculated using B-spline smoothing and the gene expression values were discretized into 10 qhttp://gardnerlab.bu.edu/data/PLoS_2007/ rhttp://gardnerlab.bu.edu/data/PLoS_2007/ 34 Nature Methods: doi:10.1038/nmeth.2016 bins. Second, CLR estimates the significance of a gene pair by comparing the mutual information value between gene iand jto a empirically defined background distribution of mutual information values. The significance of a gene pair is defined by a modified z-score. Transcription factor-target gene predictions are ranked according to the modified z-score and the highest scoring gene pairs are selected. 10.7 Mutual Information 2 – Mutual information* Mutual information is one of the 6 commonly-used, off-the-shelf algorithms run by the DREAM organizers. A full description of this method can be found in28 and is made available online.sFor two random variables, Xand Y, mutual Information is defined as: I(X;Y) = X i,j P(xi, yi) log p(xi, yi) p(xi)p(xj) In this implementation, variables Xand Yrepresent a transcription factor and target gene, respectively. Mutual information calculations require discretized data and as gene expression data from microarrays are continuous, a B-spline smoothing and discretization method is used.22 The number of bins was set to 10 as suggested by the authors. Final transcription factor-target gene predictions were ranked and selected based on the highest to lowest mutual information scores. 10.8 Mutual Information 3 – Algorithm for the reconstruction of accurate cellular networks (Aracne)* Aracne is one of the 6 commonly-used, off-the-shelf algorithms run by the DREAM organizers. A full description of this method can be found in57 and is made available online.tAs input, Aracne accepts a matrix of gene ×condition expression values and a list of defined transcription factors. The Aracne algorithm can be separated into 2 main steps. First, a matrix of mutual information values is calculated for all input pairs of genes iand j. Aracne estimates mutual information using the Gaussian Kernel estimator.8Statistically significant relationships are determined and the nonsignificant edges are removed. Second, an additional pruning step is performed based on the the theoretical property known as the data processing inequality (DPI).21 Given shttp://gardnerlab.bu.edu/data/PLoS_2007/ twiki.c2b2.columbia.edu/califanolab a gene interaction network, the DPI aims to eliminate any edges that can be explained through the remaining interactions in the network. Consider the small network gi↔gj, gj↔gk, gi↔gk, an edge between giand gj will be removed if I(gi, gj)≤min[(gi, gk),(gj, gk)]. After pruning, the remaining transcription factor-target gene relationships are then ranked based on their mutual information values. The Aracne software package has several parameters that can be adjusted by the user. As stated, default parameters were used to best estimate a naïve application of Aracne, however, it should be noted that better results are likely to be achieved after tuning of two parameters, 1) the mutual information p-value cutoff, and 2) the DPI tolerance parameter. 10.9 Mutual Information 4 & 5 –The DREAM5 network inference challenge with a combination of fast tools In the DREAM5 network inference challenge, the task is to discover relationships between genes from four gene expression datasets. We use mutual information and Fisher’s score as the scoring function to compute a probabilistic dependency between pairs of variables, and combine it with the BLCD-HITON-PC search algorithm52 to find the Yarcs in the network. An example of a Ysubstructure is A→C,B→C, and C→D. The Ysubstructures have special properties that make causal discovery possible under plausible assumptions. In this example the arc C→Dis a Yarc. The Yarcs in a Bayesian network represent unconfounded causal influences under assumptions.54 Since the Yarcs represent unconfounded causal influences we expect that they could discover arcs with high precision and can be used to orient some of the edges that pair-wise mutual information and statistical analysis (Fisher’s score) cannot accomplish. Mutual information. We used the implementation for fast calculation of mutual Information for all pairs of genes in the challenge.71 For the largest network (Network 4) in the challenge with 5950 genes and 536 chips, it took about 19 minutes to generate the pair-wise mutual information for the entire set of genes in the system. Fisher’s score. We used correlation coefficients (CC) implemented in MATLAB to find the correlation between all pairs of genes in the networks. This gave us additional information to access relationship for the whole set of pairs in the network. It is also very fast, which took 27 seconds to complete the whole run for Network 4. BLCD-HITON-PC. We implemented the Bayesian local causal discovery algorithm (BLCD)52,53 to efficiently identify unconfounded direct causal relationships of gene 35 Nature Methods: doi:10.1038/nmeth.2016 variables in the network with high precision. It consists of 2 steps: Markov Blanket step and Yarcs step. For the Markov Blanket generation step in BLCD, we used the computational causal discovery method HITONPC3,4which outputs the parents and children set for each target variable. For the Yarcs generation step, it searches for the Yarcs locally in the parents and children set. The second step can be run in parallel for efficiency. For Network 4, we separated the Yarc discovery task into 15 groups and finished running in less than 2 and a half hours for each group. BLCD outputs the probability of a Yarc which we converted to a value of 1 or 0 using a threshold of 0.5. We used binary output with threshold of 0.5 since precision of Yarcs is high and we put higher weight than mutual information and Fisher’s score. Combination. We first combined the results of mutual Information and Fisher/Correlation coefficients score (CC) based on the ranking. Two ranking arrays are averaged and normalized into one array (called MICC array) containing numbers in the range from [0,1] so that the edge with number closer to 1 represents a stronger dependent relationship. Secondly, this array is combined with the BLCD binary output such that the arcs with BLCD output ’1’ and their reverse arc with BLCD output ’0’ are averaged with the MICC array. This process could fix part of the direction problem from mutual information and Fisher’s score. For example, I(a, b) = 0.5represents the confidence of which acan regulate bby mutual information or correlation coefficient or their combination; I(a, b) == I(b, a); If BLCD outputs blcd(a, b)=1and blcd(b, a) = 0, we would average the result and the final output for these 2 arcs will be s(a, b)=0.75 and s(b, a)=0.25. This could output higher precision and recall than transcription factor or correlation coefficient or their combination if the BLCD output is correct. Discussion. The reason we combine transcription factor and correlation coefficient is based on the preliminary experimental result for similar gene expression networks. It shows that if both mutual information and correlation coefficient can output good AUPR (Area under Precision and Recall curve), their combination can produce a better one; but it is possible that the CC performance is worse than mutual information, and then the combination result can be worse than the mutual information itself. So we would rather submit 2 results: Submission DSL is based on mutual information + correlation coefficient + BLCDHITON-PC, and Submission DSL2 is based on mutual information + BLCD-HITON-PC. There are some search strategies in BLCD based algorithm we could optimize to improve the performance. The original BLCD searches for Yarcs from a Markov Blanket set of each node that is derived by greedy search; however, it is not efficient to run the greedy search for such a high-dimensional dataset for this challenge. Thus, we run an efficient algorithm HITON-PC to output a smaller set (parents and children set), but it may reduce the recall. The threshold of output would also affect the final performance. If the threshold is too high (for example, 0.9), the number of Yarcs would be very small, so it would not change the result by just using mutual information or correlation coefficient; if the threshold is too low (for example, 0.1), the number of Yarcs would be very large but the precision may be reduced, so it may not positively affect the final result. Therefore, we use the default 0.5 as the threshold. 10.10 Correlation 2 – Pearson’s correlation* Pearson’s correlation is one of the 6 commonly-used, off-the-shelf algorithms run by the DREAM organizers. Pearson’s correlation coefficient rwas calculated between all transcription factors xand all target genes yas follows: rxy =nPxiyy−PxiPyi pnPx2 i−(Pxi)2pnPy2 i−(Pyi)2, where nis the number of measurements of xand y.xto y relationships were ranked via rxy, i.e. positively correlated gene pairs receive the highers confidence. 10.11 Correlation 3 – Spearman’s correlation* Spearman’s correlation is one of the 6 commonly-used, off-the-shelf algorithms run by the DREAM organizers. Spearman’s correlation ρwas calculated between all transcription factors xand all target genes yas follows: ρxy = 1 −6Pd2 i n(n2−1), where nis the number of conditions that xand yhave been sampled and dis the difference in rank order between gene xand gene yover the nconditions. xto yrelationships were ranked on ρxy and the most correlated gene pairs were selected. 10.12 Bayesian 6 – Regulatory network inference with Bayesian networks We built a simple Bayesian network to model the influence of a potential transcription factor on a gene. The network includes observed variables for inputs, such as perturbations, knockouts, or overexpression applied in each 36 Nature Methods: doi:10.1038/nmeth.2016 experiment. It then relates these to hidden variables for unobserved quantities, such as the magnitude of the influence or whether the perturbation affected this relationship. These hidden variables are then related to random variables for the observed quantities of transcription factor and target gene expression levels. An alternate network for the absence of regulatory influence modeled each expression level as an independent random variable with random observation noise. For each possible pair of transcription factor to target gene relationships (transcription factor-target gene), we computed the relative probability of each model (no influence, A→B,B→A) given the observation data. The highest probability for each pair was selected and sorted to give the final prediction. Method. The Bayesian network used for the simplest case of a steady-state observation experiment involving a single transcription factor-target gene pair is shown in Figure SS19. The nodes of the network are: •Perturbation Enabled. If a perturbation is present for the experiment, and which one (combinations of perturbations were treated as a unique perturbation, to keep the model simple). •Affected. Bernoulli random variable for whether the target is affected by the perturbation •Perturbation Factor. Gaussian random variable for the amount of perturbation •Knockout/Overexpression Knockout or overexpression factor of the target for a given experiment •Edge Influence. Gaussian random variable for the magnitude of influence of the transcription factor on the target gene •Transcription factor Level, Target Level. Expression levels of the transcription factor and target genes, normalized to mean 1.0, variance 0.01 of observations from steady-state experiments. •Transcription factor Level Observed, Target Level Observed. Observations with Gaussian noise The null hypothesis of no relationship was a simple model where each expression level is an independent Gaussian random variable observed with Gaussian noise. A more complex model was used for time series experiments, a one-step dynamic Bayesian network modeling the change in the target as a function of its prior level and the prior level of the transcription factor. The null hypothesis modeled each level independently as returning to mean with some time constant. The data were analyzed using a Bayesian network library, to calculate the relative probability of the respective models (no edge, A→B, or B→A) given the observations. These calculations were then performed pairwise Figure S19: Inference method design for Bayesian 6. The organization and categorization of information for the design of the Bayesian classifier. for all possible transcription factor-transcription factor or transcription factor-target gene interactions in each network on an internal computing cluster, since the method is quite computationally intensive. Transcription factortarget gene interactions were modeled in both directions and the highest probability taken, on the assumption that the magnitude of influence was significant even if directionality was suspect. The 100,000 most likely edges were then selected and sorted. Discussion. This method did not perform well on any of the DREAM5 data sets. There are several arbitrary constant parameters in the model that define the prior distributions of the hidden variables. These were trained on a small (50-node) yeast network generated from GeneNetWeaver,56 using a conjugate simplex method to find the parameter values that gave predictions best matching the actual network. The resulting model and parameters were then validated on the DREAM4 dataset, on which they performed reasonably. However, it was clear during the training phase that the prediction accuracy was quite sensitive to these parameters. It is likely that the DREAM5 datasets were sufficiently different from the training set that the parameters were no longer suitable. There are a couple of ways in which this issue can be addressed in the future. Most simply, the parameters can be trained on the actual dataset; in the absence of a known network, we could optimize for fit of predicted probabilities to an expected distribution based upon general features of the target network (e.g. inand out-degree distributions). This would, however, be quite computationally 37 Nature Methods: doi:10.1038/nmeth.2016 intensive. Alternately, the parameters themselves can also be incorporated into the model as hidden variables. This would require more complex Bayesian network modeling algorithms, and the library we used did not support model selection over such second-order networks (e.g. a random variable whose mean and variance were themselves defined by other random variables). It might also be possible to reconfigure the Bayesian network itself to reduce its dependence on such arbitrary parameters. 10.13 Other 1 – Inferring regulatory networks using tree-based methods This algorithm, called GENIE3 (for “Gene Network Inference with Ensemble of Trees”), decomposes the prediction of a regulatory network between pgenes into pdifferent regression problems. In each of the regression problems, the expression pattern of one of the target genes is predicted from the expression patterns of all known transcription factors, using a tree-based ensemble method called Random Forests.13 The importance of a transcription factor in the prediction of the target gene expression pattern is taken as an indication of a putative regulatory link. Putative regulatory links are then aggregated over all genes to provide a ranking of interactions from which the whole network is reconstructed. GENIE3 does not make any assumption about the nature of gene regulation, can deal with combinatorial and nonlinear interactions, produces directed gene regulatory networks, and is fast and scalable. This method is described in more detail in40 and available online.u Network inference procedure. The GENIE3 procedure works as follows: 1. For gene j= 1 to p: -Generate the learning sample of input-output pairs for gene j:LSj={(xT F k, xj k), k = 1, . . . , N}, where the input xT F kis the vector of expression values of all known transcription factor genes (except gene jif it is a transcription factor) in the kth experiment, the output xj kis the expression value of gene jin the kth experiment, and N is the number of experiments. -Use the Random Forests method on LSjto compute confidence level wi,j, for each transcription factor gene i6=j. 2. Aggregate the pindividual gene rankings to get a global ranking of all regulatory links. Tree-based ensemble methods. The basic idea of treebased methods in regression is to recursively split the uhttp://www.montefiore.ulg.ac.be/~huynh-thu/software. html learning sample with binary tests based each on one input variable (here, the expression of one potential transcription factor). These tests are optimized in such a way as to reduce as much as possible the variance of the output variable (here, the expression of the target gene) in the resulting subsets of samples. Candidate splits for numerical variables typically compare the input variable values with a threshold which is determined during the tree growing. Single trees are usually very much improved upon by ensemble methods, which average the predictions of several trees. In the network inference procedure, we used the Random Forests method,13 i.e. each tree is built on a bootstrap sample from the original learning sample, and at each test node, Kattributes are selected at random among all candidate attributes before determining the best split. We used K=√n, where nis the number of input variables, and grow ensembles of 1000 trees. We selected the Random Forests method (with this value of K) because it leads to the best performance among other tree-based ensemble methods on the prediction of the E. coli regulation network.40 Variable importance measure. To associate a confidence level wi,j to the regulation of target gene jby the transcription factor gene i, we directly exploit the variable importance measure of the expression of ias derived from the tree-based model learned for j. While several variable importance measures have been proposed in the literature for tree-based methods, we consider in this procedure a measure that computes at each test node Nthe total reduction of the variance of the output variable due to the split, defined by:14 w(N)=#S.V ar(S)−#St.V ar(St)−#Sf.V ar(Sf) where Sdenotes the set of samples that reach node N, St(resp. Sf) denotes its subset for which the test is true (resp. false), V ar(.)is the variance of the output variable in a subset, and # denotes the cardinality of a set of samples. For a single tree, the overall importance of one transcription factor gene is then computed by summing the wvalues of all tree nodes where the expression of this transcription factor was used to split. For an ensemble, the importance is then obtained by averaging importance scores over all trees in the ensemble. Global regulatory link ranking. Each tree-based model thus yields a separate ranking of the transcription factors as potential regulators of a target gene in the form of weights wi,j computed as sums of total variance reductions. The sum of the importances of all input variables for a tree is equal to the total variance of the output variable explained by the tree, which in the case of unpruned trees (as they are in the case of Random Forests ensembles) is usually very close to the initial total variance of the output. 38 Nature Methods: doi:10.1038/nmeth.2016 As a consequence, if we trivially order the regulatory links according to the weights wi,j , this is likely to introduce a positive bias for regulatory links towards the more highly variable genes.40 To avoid this bias, we first normalized the gene expression values so that they all have a unit variance in the training set, before applying the tree-based ensemble method. This normalization indeed implies that the different weights inferred from different models predicting the different gene expressions are comparable. Discussion. We developed GENIE3, a gene network inference algorithm based on feature selection with Random Forests. This method was the overall top performer of the challenge. It was also the top performer on the in silico generated data but did not perform as well as some other teams on the in vivo microarray data. Interestingly, setting Random Forests parameter Kto nwould have significantly improved the performance on the in silico data but with a decrease of the performance on the in vivo networks. One potential reason for the better performance on the in silico benchmarks is that the in silico dataset probably contains more statistically useful experiments than the in vivo datasets in which there may be some redundancy among the experiments or some bias in their selection. Such differences in terms of the quality of the data might affect more the non parametric approach than alternative approaches that make stronger assumptions. Other potential reasons may originate from the fact that the E. coli and S. cerevisiae gold standard networks are not complete and noisy, as well as from the discrepancy that certainly exists between the simulation model used to generate the in silico data and the in vivo regulation mechanisms of E. coli and S. cerevisiae. In the future, we would like to investigate further these differences. In principle, any feature selection algorithm74 could be substituted to Random Forests within the general procedure. Method Regression 1 (see Supplementary Note 10.1) could be seen as an instance of this procedure. We have also carried out experiments with feature ranking based on linear models trained with support vector regression but these methods were not as successful as treebased ensemble methods. In the future, we plan to consider other feature selection techniques. 10.14 Other 2 – Inferring gene regulatory networks by ANOVA This approach to network inference is based on the assumption that transcription factors (TFs) and their corresponding target genes (TGs) exhibit mutual expression dependencies in at least a subset of the measured experimental conditions (time points, perturbations etc). Such candidate interactions, i.e. pairs of a TF and a TG, are ranked by a score s. The score scan be any measure of dependency between the expression of the TF and its TG. Frequently used measures of dependency are based on Pearsons or Spearmans correlation coefficients, mutual information, or in case of Bayesian network inference on conditional probability tables. We evaluate candidate interactions by η2,19 a nonparametric, non-linear correlation coefficient obtained from a two-way analysis of variance (ANOVA). It is fast, easy to apply and does not require discretization of the input data. Refinements of this approach also enable incorporation of additional information, e.g. the specific over-expression or deletion of genes in given experiments. This method will be described in more detail in.44 Data preprocessing. Basal gene levels can be quite different between experiments. To account for these differences, we transformed the absolute expression values into expression fold changes. Fold changes were computed by mapping each measured condition to one or more control conditions from the same experiment. Control conditions were defined via their set of treatments (such as knockout, over-expression, drug treatments) that is required to be a subset of the treatments applied in the given measured condition. More than one fold change may be computed if more than one control is available. For example, the control conditions for an experiment where two genes are deleted (double knockout) may be a single knockout and/or the wild type. In the case of time series, we required the time points for measured conditions to match the time points of the corresponding controls. Note that each control usually has several replicates. Fold changes were computed by subtracting the (logtransformed) gene expression values of the controls, averaged across the replicates, from the given chips. As mentioned above, the mapping is not unique, i.e. a given chip can have 0, 1 or more control replicate sets assigned. For instance, from the 805 chips (487 different replicate sets) of the E. coli compendium, we could compute gene fold changes for 599 chips (379 replicate sets). Because of the multiplicity of controls, we obtained a total of 935 fold changes per gene (602 replicate sets). Comparison of TFs and putative TGs by η2. We employed a two-way ANOVA to test the differential expression of TFs and their putative TGs. The ANOVA compares the means of populations that are classified in two different ways, or the mean responses in an experiment with two factors. The factors analyzed here are the expression of two different genes (factor or dimension A) across a range of conditions (factor B) that is represented as a 2 rows by 602 columns matrix. In this analysis we tested (1) if at least two experimental conditions exhibit significant differences in their population means (i.e., differential expression) and (2) if these 39 Nature Methods: doi:10.1038/nmeth.2016 differences exceed the differences between the expression of the TF and a TG. Phrased in terms of the two-way ANOVA, the strength of an association is proportional to the fraction of the variance across conditions (factor A) in the total variance. Such a fraction (F-value) follows the F-statistic, which can be used to derive the statistical significance of the involved factors as p-values. Refinement of the basic ANOVA. Genetic perturbations are valuable as they help to establish directed causal relationships between regulators and their TGs. The information whether a TF was subjected to genetic perturbations (deletion or overexpression) was taken from the provided chip-feature descriptions. Conditions that indicate a perturbation of the currently tested TF are given a higher weight than other conditions. Informally, the weight is processed by inserting (w-1) additional copies of such a condition into the ANOVA matrix. Note that conditions where non-TFs or TFs other than the currently tested TF are perturbed receive the standard weight. We set w= 50 based on an analysis we performed on expression data obtained from the M3D database27 and a gene regulatory network obtained from RegulonDB.31 Discussion. For the detection of dependencies we proposed the measure η2that is derived from an analysis of variance (ANOVA). To our knowledge, η2has not been widely applied to network inference or to other problems in Bioinformatics, although it has a number of features that can facilitate the detection of gene dependencies. Like Pearson’s correlation, but in contrast to Bayes conditional probability tables or mutual information, η2does not require the discretization of the input data. This increases the robustness of this method as inappropriate discretization might lead to loss of signal. In contrast to Pearson’s linear correlation coefficient, η2is a nonparametric, non-linear correlation coefficient. Some of the known E. coli interactions identified by this approach were quite interesting biologically. For instance, an interaction between the multiple antibiotic resistance (mar) genes marA and marB was active after antibiotic treatment but not in growth phase experiments.44 The measure η2allowed us to detect such local correlations arising from condition-specific interactions. This increased sensitivity is due to the effective utilization of replicated measurements to model the measurement error and to estimate the statistical significance of TF-TG dependencies. In the future, we intend to further improve this ANOVAbased inference approach by including a dedicated treatment of time series and by using conditional correlations to distinguish direct from indirect interactions. 10.15 Other 3 – Network inference through Boolean networks The underlying model in this method is a Boolean network where the topology and logic is inferred from continuous expression data. It is based on the information theoretic conditional entropy criterion,5but assumes that the Boolean state of the network is not directly observed. Instead we are given a dataset of continuous measurements that reflect probabilistically the Boolean state. For every gene X, we want to find the set of regulators Ythat give the best conditional entropy score H(X|Y) among all sets of putative transcription factors up to some size limit. Since conditional entropy is computed for discrete random variables, and we have continuous measurements, we first transform the continuous data into a set of discrete observations. The simplest way would be to discretize the continuous data, e.g. set all values above some threshold to 1 and all values below it to 0. Instead, we interpret a vector of continuous values as a probability distribution over all possible Boolean vectors of the same dimension. To put it simply, instead of creating one Boolean vector with probability 1 for every continuous vector, for every continuous vector we create all possible Boolean vectors of the same dimension, and assign each such vector a probability. The probabilities are chosen as follows: we first normalized the continuous expression values of every gene to have mean 0 and standard deviation 1.5 (a value determined empirically). After normalization, we set the probability that a single (one-dimensional) continuous value ccorresponds to the Boolean value 1 to 1 1+e−c(the logistic function with input c). The probability that a continuous vector ¯ccorresponds to a specific Boolean vector ¯ bthen becomes: p(¯ b|¯c) = Y bi=1 1 1 + e−ciY bi=0 1−1 1 + e−ci, where ci(bi) is the value of the ith entry of ¯c(¯ b). Note that by setting the standard deviation we avoid using any parameters in the logistic function. Given a continuous dataset of Nsamples (experiments/microarrays) that are assumed to be independent and identically distributed, the probability of seeing the Boolean vector ¯ bin this dataset is: P(¯ b) = X ¯ci∈samples p(¯ b|¯ci) N In other words, for each Boolean vector we sum the probabilities that each continuous vector corresponds to it. With the probability distribution over all Boolean vectors in hand, we can use information theory to evaluate different topologies of the network. Similar approaches were 40 Nature Methods: doi:10.1038/nmeth.2016 previously used for network reconstruction.46,49 Denote by HC(X|Y)the conditional entropy for a gene Xand a set Yof regulators as computed using continuous data. We selected for every gene Xthe set Yof regulators that gave the best HC(X|Y)score among all sets and transcription factors of size ≤3. For time-series, Xwas taken from time tand Yfrom time t−1. Otherwise, we assumed samples were in steady state and Xand Ywere taken from the same experiment. Experiments in which Xis perturbed were discarded when computing HC(X|Y). For the DREAM5 network inference challenge, the transcription factor-target gene relationships identified by the above procedure were ranked according to the conditional entropy of the regulator set to which the transcription factor belonged to. This constituted 4,738, 13,120 and 17,519 interactions for networks 1,3 and 4 respectively. Since the challenge allowed submitting up to 100,000 regulatory interactions, we added to the list of predictions those pairs of genes that had the highest correlation, ranked by their correlation. In case of time series, the correlation was computed between the level of the regulator at time t−1 and the level of the regulatee at time t. In order to cope with the large number of regulator sets in networks 3,4 we used a computer cluster to distribute the tasks. Note that this method is inherently incompatible with the scoring scheme of DREAM5, because it assigns a score to the set of regulators of every gene and not to every transcription factor-target gene pair separately. Another difference is the limit of three regulators per gene that we imposed. This limit allowed us to examine all sets of 3 regulators, a very large set given the DREAM5 network’s size. In practice, however, some genes have more than 3 regulators. A speedup version in which regulators are selected incrementally can allow more regulators per gene.34 Finally, a set of 3 regulators will always score better in practice (sometimes insignificantly better) than a set of 1 or 2 regulators, and we did not set a criterion that prevents the addition of regulators that do not improve the score significantly. An advantage in this method’s focus on sets of regulators is that there is a natural way to derive the regulatory logic given the set of best scoring HC(Xi|Y)for every gene i. This can be achieved by performing steepest descent on the function , i.e. on the total entropy of the network. The partial derivative with respect to every continuous variable can be computed exactly. A tool implementing an extension of this method is available for download at: acgt.cs.tau.ac.il/ modent/ 10.16 Other 6 – Network inference using quantitative modeling and evolutionary algorithms The method presented here for inference of DREAM5 gene networks is based on quantitative modeling using an integrative evolutionary approach. The model used is a single-layered artificial neural network (ANN) and allows for extraction of qualitative information on connections, with the advantage of maintaining the ability to simulate quantitative behavior. Given the high dimensionality of the four datasets, the networks are less suitable for direct quantitative modeling, as simulation for model evaluation is computationally expensive, and the gene interaction space is very large. Although reverse engineering is difficult, quantitative models are valuable for networks of this size, as they allow for large scale in-silico simulation of the real system. We have included in the workflow several mechanisms to reduce the search space required for reverse engineering. These include grouping of tightly correlated genes into modules and filtering of putative transcription factors for each gene, prior to model inference by evolutionary optimization. Additionally, models have been obtained for each gene at a time for non-transcription factor genes, while whole network analysis has been performed for the transcription factor subnetwork only. Data preprocessing. Due to the requirements of modeling and the inferential approach, which uses ANNs, the expression values in the datasets had to be scaled to values on the interval [0,1]. Scaling was performed differently for the in silico and in vivo data, as values in the former were more homogeneous than the latter, due to lack of experimental differences. For the former, all values were scaled by dividing by the maximum value in the dataset. However, in the in vivo datasets, individual genes had very different ranges of expression. Additionally, we have observed that, for some knockout (KO) genes, expression values were very large. In consequence, each gene vector (column) has been first translated so that the minimum value over all experiments is close to 0, after which the scaling of all values (as above) was performed. This ensured a common scale for gene values and, at the same time, brought KO genes to a low level of expression, while not affecting the oscillations seen in the data. Module computation. One dimensionality reduction mechanism was grouping tightly correlated genes (values from all available experiments) into modules, and considering these as a single gene in the network. These may correspond to operons, which are common functionality groupings in GRNs. Module computation was based on pair-wise Pearson correlation and thresholds used were 0.9 for Networks 1, 3 and 4 and 0.95 for Network 2 (the higher threshold used for the last dataset was due to the limited number of experiments available, compared to the others). 41 Nature Methods: doi:10.1038/nmeth.2016 not perturbed alone and then we could not subtract its effect, accepting to make some mistakes as a trade off for also identifying real targets of TFi. In case of triple knockouts, we proceeded in a similar way. S1 ij(the confidence in TFi→Tj) was obtained by averaging the absolute RT Fi→Tjvalues from single, double and triple knockout and over-expression experiments. Partial correlation analysis. This was applied to steady state data (we only used the last time point of the time series), with the goal to exploit the expression correlation between a TFiand its targets. We applied full order partial correlation75 using the GeneNet R package.vS2 ij is the absolute value of partial correlations ωT Fi,Tj. Co-deviation analysis. This was applied to chips with drug perturbations, and chips featuring non-TF single gene perturbations. The goal was to check if target Tjconsistently makes large deviations in expression when TFi makes large deviations. We first converted the gene expression values into z-scores. Then, for each transcription factor TFi: 1) The data was split into two subsets: one group of observations DH iwith zi> d (T Fiis ’high’) and another other group of observations DL iwith zi<−d(TFiis ’low’); best results were obtained with d= 0.5, 2) To identify the potential targets Tjof TFi, for each Tja two-sided t-test was performed to check whether its mean in DH iis significantly different from its mean in DL i.S3 ij is the absolute value of the t-statistic tij (test performed for Tjwhen datasets are formed based on deviations of TFi). Combining. The edge confidences S1 ij,S2 ij and S3 ij were combined into a single confidence score by using a weighted mean of the ranks. The weights were obtained through an optimization process on simulated data. Discussion. This method was among the top performers on the in vivo data, but had only average performance on the simulated data. These results reveal that the combination of different inference techniques is indeed useful, especially with heterogeneous compendia, provided that these data are correctly subdivided. Moreover, even if gene expression simulations are more and more realistic, it clearly emerges that inference methods having good performances on synthetic datasets cannot always be expected to obtain the same results on in vivo gene expression data. This emphasizes the need for realistic simulation, to generate data with properties similar to those observed in in vivo data, an issue in which we have put much effort when creating SysGenSIM. vstrimmerlab.org/software/genenet Figure S21: Inference method design for Meta 5. The organization of the data and the the algorithms that were used as input to the naïve Bayesian classifier. 10.21 Meta 5 – A naïve Bayes based approach to network inference In this methodology we start by splitting the data according to the experiment type, then we analyze them according to specific statistical methods. Each method is used to evaluate the plausibility of the edges. Finally, to combine the statistics we follow a naïve Bayes approach. Method. To analyze the data, the four statistics summarized below have been used: •Pearson correlation coefficient (PCC): a standard method used to get a basic guess about the relationship between each pair of genes. •Limma: a linear models approach to assess differential genes expression. This is a commonly used algorithm for analyzing experiments when the number of samples is limited.80 •maSigPro: a method for the analysis of experiments that include time series. In working with this tool a sequential experiment design matrix has been used.20 •z-score: the standard zstatistic. We used it in experiments when other tools could not be applied. We dealt with each type of experiment design using one (or two) of the above statistics. The statistic used for each kind of experiment is shown in Figure S10.21. In deletion experiments we put an edge between the deleted gene and all the significantly expressed genes. In the other experiments we find the clusters of co-expressed genes and put an edge between each transcription factor and all other genes. The combination of the statistics was evaluated using a naïve Bayes approach. The goal is to compute the probability that an edge belongs to the network given the experimental evidence. Let us denote with Xthe variable that is equal to 1 when a given edge belongs to the network and to 0 otherwise. Let us also define with Y1. . . Ymthe values of statistics assessing the given edge. We would like to compute the probability P(X= 1|Y1. . . Ym). By Bayes’s theorem, this 48 Nature Methods: doi:10.1038/nmeth.2016 quantity can be written as: P(X= 1|Y1. . . Ym) =P(Y1···Ym|X= 1)P(X= 1) PP(Yi···Ym|X=x)P(X=x) =QP(Yi|X= 1)P(X= 1) PQP(Yi|X=x)P(X=x) where the equality holds by assuming independence between the statistic values. In order to compute the given formula, we specify the probabilities for P(X= 1) and P(X= 0). We evaluated them by exploiting the fact that in scale free networks: k−γis the probability that a randomly chosen node has exactly kedges (where γ∈(2,...,3)). Here, as suggested by Barabasi and Albert,7we set γ= 3. It follows that the number of edges e(N)in a scale free network of size Ncan be computed as: e(N) = PN k=1 N×P(k)×k=NPN k=1 1 k2. We approximate the quantity PN k=1 1 k2with π2 6(its limit for N→inf). The apriori probability P(X= 1) of picking up an edge belonging to a network of size Nis given by e(N) N2=π2 6N. Finally, each of the P(Yi|X=x)distributions has been manually set by taking into consideration the characteristics of each statistic.87 Discussion. Gene network reverse engineering is a major challenge in computational biology and, the presented method is one of the simplest approaches that can be developed for a problem of this complexity. Indeed, simplicity has been one of the goals we strived to attain in its design. In fact, past DREAM contests emphasized that simpler methods could perform as well as others. Also, when in vivo networks are to be analyzed, data scarcity and its quality demand for classifiers built using a small number of well understood parameters. This method showed average performances in the DREAM contest. An interesting facet of this methodology is that it performed remarkably better in the case of in vivo networks than with in silico ones: it is among the top performers (it ranks third) when the in silico dataset is not considered (it ranks 15th otherwise). A number of interesting questions could be raised by this observation: what is in synthetic datasets that set them apart from natural ones? Should one strive to optimize new algorithm more aggressively on natural dataset? Could the culprit be found in the quality of in vivo data, so that most of these methods will perform much better when this quality increases? We believe that the answers to these questions may help in better understanding current tools and in developing new ones. 49 Nature Methods: doi:10.1038/nmeth.2016 References [1] D. Abdulrehman, P. T. Monteiro, M. C. Teixeira, N. P. Mira, A. B. Lourenco, S. C. dos Santos, T. R. Cabrito, A. P. Francisco, S. C. Madeira, R. S. Aires, A. L. Oliveira, I. Sa-Correia, and A. T. Freitas. Yeastract: providing a programmatic access to curated transcriptional regulatory associations in saccharomyces cerevisiae through a web services interface. Nucleic Acids Res, 39:136–140, 2010. [2] A. Ahmed and E. P. Xing. Recovering time-varying networks of dependencies in social and biological studies. PNAS of USA, 106(29):11878–11883, 2009. [3] C. F. Aliferis, A. Statnikov, I. Tsamardinos, S. Mani, and X. D. Koutsoukos. Local causal and markov blanket induction for causal discovery and feature selection for classification. part i: Algorithm and empirical evaluation. Journal of Machine Learning Research, 11:171–234, 2010. [4] C. F. Aliferis, A. Statnikov, I. Tsamardinos, S. Mani, and X. D. Koutsoukos. Local causal and markov blanket induction for causal discovery and feature selection for classification. part ii: Analysis and extensions. Journal of Machine Learning Research, 11:235–284, 2010. [5] C. Arndt. Information Measures: Information and its description in Science and Engineering. Springer, Berlin, 2001. [6] T. Baba, T. Ara, M. Hasegawa, Y. Takai, Y. Okumura, M. Baba, K. Datsenko, M. Tomita, B. Wanner, and H. Mori. Construction of escherichia coli k-12 in-frame, single-gene knockout mutants: the keio collection. Molecular Systems Biology, 2(2006.0008), 2006. [7] A. L. Barabasi and R. Albert. Emergence of scaling in random networks. Science, 286:509–512, 1999. [8] J. Beirlant, E. Dudewicz, L. Gyorfi, and E. van der Meulen. Nonparametric entropy estimation: An overview. Int J Math Stat Sci, 6(1):17–39, 1997. [9] B. M. Bolstad, R. A. Irizarry, M. Astrand, and T. P. Speed. A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinformatics, 19(2):185–193, 2003. [10] R. Bonneau, M. T. Facciotti, D. J. Reiss, A. K. Schmid, M. Pan, and et al. A predictive model for transcriptional control of physiology in a free living cell. Cell, 131:1354– 1365, 2007. [11] R. Bonneau, D. J. Reiss, P. Shannon, M. T. Facciotti, L. Hood, and et al. The inferelator: an algorithm for learning parsimonious regulatory networks from systemsbiology data sets. Genome Biology, 7, 2006. [12] J. C. Borda. Memoire sur les elections au scrutin, 1781. [13] L. Breiman. Random forests. Machine Learning, 45(1):5– 32, 2001. [14] L. Breiman, J. H. Friedman, R. A. Olsen, , and C. J. Stone. Classification and Regression Trees. Chapman & Hall, 1984. [15] A. Brooun, J. J. Tomashek, and K. Lewis. Purification and ligand binding of emrr, a regulator of a multidrug transporter. J. Bacteriol., 181(16):5131–5133, 1999. [16] N. A. Burton, M. D. Johnson, P. Antczak, A. Robinson, and P. A. Lund. Novel aspects of the acid response network of e. coli k-12 are revealed by a study of transcriptional dynamics. J. Mol. Biol., 401(5):726–742, 2010. [17] C.-C. Chang and C.-J. Lin. LIBSVM: a library for support vector machines, 2001. Software available at http://www. csie.ntu.edu.tw/~cjlin/libsvm. [18] K. Y. Choi and H. Zalkin. Regulation of escherichia coli pyrc by the purine regulon repressor protein. J. Bacteriol., 172(6):3201–3207, 1990. [19] J. Cohen. Eta-squared and partial eta-squared in fixed factor anova designs. Educational and Psychological Measurement, 33(1):107, 1973. [20] A. Conesa, M. J. Nueda, A. Ferrer, and M. Talon. maSigPro: a method to identify significantly differential expression profiles in time-course microarray experiments. Bioinformatics, 22(9):1096–1102, 2006. [21] T. M. Cover and J. A. Thomas. Elements of Information Theory. John Wiley & Sons, New York, NY, 1991. [22] C. O. Daub, R. Steuer, J. Selbig, and K. S. Estimating mutual information using b-spline functions–an improved similarity measure for analysing gene expression data. BMC Bioinformatics, 5:118, 2004. [23] J. Davis and M. Goadrich. The relationship between Precision-Recall and ROC curves. In ICML ’06: Proceedings of the 23rd international conference on Machine learning, pages 233–240, Pittsburgh, Pennsylvania, 2006. ACM. [24] T. G. Dietterich. Ensemble methods in machine learning. In Multiple Classifier Systems, volume 1857, pages 1–15. Springer Berlin Heidelberg, Berlin, Heidelberg, 2000. [25] B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani. Least angle regression. Ann. Stat., 32(2):407–499, 2004. [26] S. M. Egan and R. F. Schleif. A regulatory cascade in the induction of rhabad. J. Mol. Biol., 234(1):87–98, 1993. [27] J. J. Faith, M. E. Driscoll, V. A. Fusaro, E. J. Cosgrove, B. Hayete, F. S. Juhn, S. J. Schneider, and T. S. Gardner. Many microbe microarrays database: uniformly normalized affymetrix compendia with structured experimental metadata. Nucleic Acids Res, 36:866–870, 2008. [28] J. J. Faith, B. Hayete, J. T. Thaden, I. Mogno, J. Wierzbowski, and et al. Large-scale mapping and validation of escherichia coli transcriptional regulation from a compendium of expression profiles. PLoS Biology, 5:54– 66, 2007. [29] R. Foygel and M. Drton. Exact block-wise optimization in group lasso and sparse group lasso for linear regression. Technical report, University of Chicago, Department of Statistics, 2010. available at: http://arxiv.org/abs/1010.3320. [30] J. H. Friedman, T. Hastie, and R. Tibshirani. Sparse inverse covariance estimation with the graphical lasso. Biostatistics, 9:432–441, 2008. [31] S. Gama-Castro, H. Salgado, M. Peralta-Gil, A. SantosZavaleta, L. Muniz-Rascado, H. Solano-Lira, V. JimenezJacinto, V. Weiss, J. S. Garcia-Sotelo, A. Lopez-Fuentes, 50 Nature Methods: doi:10.1038/nmeth.2016 L. Porron-Sotelo, S. Alquicira-Hernandez, A. MedinaRivera, I. Martinez-Flores, K. Alquicira-Hernandez, R. Martinez-Adame, C. Bonavides-Martinez, J. MirandaRios, A. M. Huerta, A. Mendoza-Vargas, L. ColladoTorres, B. Taboada, L. Vega-Alvarado, M. Olvera, L. Olvera, R. Grande, E. Morett, and J. Collado-Vides. RegulonDB version 7.0: transcriptional regulation of escherichia coli k-12 integrated within genetic sensory response units (Gensor units). Nucleic Acids Research, 39(Database issue):D98–105, 2011. [32] A. Greenfield, A. Madar, H. Ostrer, and R. Bonneau. Dream4: Combining genetic and dynamic information to identify biological networks and dynamical models. PLoS one, 5:e13397, 2010. [33] C. T. Harbison, D. B. Gordon, T. I. Lee, N. J. Rinaldi, K. D. Macisaac, T. W. Danford, N. M. Hannett, J. Tagne, D. B. Reynolds, J. Yoo, E. G. Jennings, J. Zeitlinger, D. K. Pokholok, M. Kellis, P. A. Rolfe, K. T. Takusagawa, E. S. Lander, D. K. Gifford, E. Fraenkel, and R. A. Young. Transcriptional regulatory code of a eukaryotic genome. Nature, 431(7004):99–104, 2004. [34] R. F. Hashimoto, E. Dougherty, M. Brun, Z. Zhou, M. L. Bittner, and et al. Efficient selection of feature sets possessing high coefficients of determination based on incremental determinations. Signal Process, 83(4):695–712, 2003. [35] A. Haury, F. Mordelet, P. Vera-Licona, and J. Vert. TIGRESS: trustful inference of gene regulation using stability selection. arXiv:1205.1181, May 2012. [36] B. He and H. Zalkin. Regulation of escherichia coli pura by purine repressor, one component of a dual control mechanism. J. Bacteriol., 176(4):1009–1013, 1994. [37] M. J. Herrgard, M. W. Covert, and B. O. Palsson. Reconciling gene expression data with known genome-scale regulatory network structures. Genome Res., 13(11):2423– 2334, 2003. [38] F. Hommais, E. Krin, J. Y. Copp’ee, C. Lacroix, E. Yeramian, A. Danchin, and P. Bertin. Gade (yhie): a novel activator involved in the response to acid environment in escherichia coli.Microbiology, 150(1):61–72, 2004. [39] Z. Hu, P. J. Killion, and V. R. Iyer. Genetic reconstruction of a functional transcriptional regulatory network. Nat Genet, 39:683–687, 2007. [40] V. A. Huynh-Thu, A. Irrthum, L. Wehenkel, and P. Geurts. Inferring regulatory networks from expression data using tree-based methods. PLoS one, 5(9):e12776, 2010. [41] E. Keedwell and A. Narayanan. Discovering gene networks with a neural-genetic hybrid. IEEE/ACM Trans. Comput. Biol. Bioinformatics, 2:231–242, 2005. [42] I. M. Keseler, J. Collado-Vides, A. Santos-Zavaleta, M. Peralta-Gil, S. Gama-Castro, L. MuÃśiz-Rascado, C. Bonavides-Martinez, S. Paley, M. Krummenacker, T. Altman, P. Kaipa, A. Spaulding, J. Pacheco, M. Latendresse, C. Fulcher, M. Sarker, A. G. Shearer, A. Mackie, I. Paulsen, R. P. Gunsalus, and P. D. Karp. EcoCyc: a comprehensive database of escherichia coli biology. Nucleic Acids Research, 39(Database issue):D583–590, 2011. [43] P. Kheradpour, A. Stark, S. Roy, and M. Kellis. Reliable prediction of regulator targets using 12 drosophila genomes. Genome Research, 17(12):1919–1931, 2007. [44] R. Küffner, T. Petri, L. Windhager, and R. Zimmer. Inferring gene regulatory networks by anova, In preparation. [45] S. Lebre, J. Becq, F. Devaux, M. P. H. Stumpf, and G. Lelandais. Statistical inference of the time-varying structure of gene regulation networks. BMC Systems Biology, 4(1):130, 2010. [46] S. Liang, S. Fuhrman, and R. Somogyi. Reveal, a general reverse engineering algorithm for inference of genetic network architectures. In Pacific Symposium on Biocomputing, volume 3, pages 18–29, 1998. [47] S. Lin. Rank aggregation methods, 2010. [48] O. Lomovskaya, K. Lewis, and A. Matin. Emrr is a negative regulator of the escherichia coli multidrug resistance pump emrab. J. Bacteriol., 177(9):2328–2334, 1995. [49] F. M. Lopes, D. C. Martins, and R. M. Cesar. Feature selection environment for genomic applications. BMC Bioinformatics, 9:451, 2008. [50] K. D. MacIsaac, T. Wang, D. B. Gordon, D. K. Gifford, G. D. Stormo, and E. Fraenkel. An improved map of conserved regulatory sites for saccharomyces cerevisiae. BMC Bioinformatics, 7:113, 2006. [51] A. Madar, A. Greenfield, E. Vanden-eijnden, and R. Bonneau. Dream3: Network inference using dynamic context likelihood of relatedness and the inferelator. PloS one, 5(3):e9803, 2010. [52] S. Mani, G. Cooper, and A. Statnikov. Causal discovery algorithms based on y structures, Under review. [53] S. Mani and G. F. Cooper. A bayesian local causal discovery algorithm. In M. Fieschi and et al., editors, Proceedings of the World Congress on Medical Informatics, pages 731–735, 2004. [54] S. Mani, P. Spirtes, and G. Cooper. A theoretical study of y structures for causal discovery. In Proceedings of the Conference on Uncertainty in Artificial Intelligence, pages 314–323, 2006. [55] D. Marbach, R. J. Prill, T. Schaffter, C. Mattiussi, D. Floreano, and G. Stolovitzky. Revealing strengths and weaknesses of methods for gene network inference. Proceedings of the National Academy of Sciences of the United States of America, 107(14):6286–6291, 2010. [56] D. Marbach, T. Schaffter, C. Mattiussi, and D. Floreano. Generating realistic in silico gene networks for performance assessment of reverse engineering methods. Journal of Computational Biology, 16(2):229–239, 2009. [57] A. A. Margolin, K. Wang, W. K. Lim, M. Kustagi, I. Nemenman, and A. Califano. Reverse engineering cellular networks. Nature Protocols, 1:662–671, 2006. [58] N. Meinshausen and P. Bühlmann. High-dimensional graphs and variable selection with the lasso. The Annals of Statistics, 34:1436–1462, 2006. 51 Nature Methods: doi:10.1038/nmeth.2016 [59] N. Meinshausen and P. Bühlmann. Stability selection. J. Royal. Statist. Soc. B., 72(4):417–473, 2010. [60] T. Michoel, R. De Smet, A. Joshi, Y. Van de Peer, and K. Marchal. Comparative analysis of module-based versus direct methods for reverse-engineering transcriptional regulatory networks. BMC Systems Biology, 3(1):49, 2009. [61] F. Mordelet and J. P. Vert. SIRENE: supervised inference of regulatory networks. Bioinformatics, 24:76–82, 2008. [62] M. E. J. Newman. Modularity and community structure in networks. PNAS, 103(23):8577–8582, 2006. [63] P. Novichkov, O. Laikova, E. Novichkova, M. Gelfand, A. Arkin, I. Dubchak, and D. Rodionov. Regprecise: a database of curated genomic inferences of transcriptional regulatory interactions in prokaryotes. Nucleic Acids Research, 38(Database issue):D111–118, 2010. [64] J. Pearl. Causality. Models, Reasoning, and Inference. Cambridge University Press, Cambridge, UK, 2009. [65] M. W. Pfaffl. A new mathematical model for relative quantification in real-time rt-pcr. Nucl. Acids Res., 29(9):e45, 2001. [66] A. Pinna, N. Soranzo, and A. de la Fuente. From knockouts to networks: Establishing direct cause-effect relationships through graph analysis. PLoS one, 5(10):e12912, 2010. [67] A. Pinna, N. Soranzo, V. De Leo, and A. de la Fuente. Elucidating transcriptional regulatory networks from heterogeneous gene-expression compendia, In preparation. [68] A. Pinna, N. Soranzo, I. Hoeschele, and A. de la Fuente. Simulating system genetics data with SysGenSIM, In press. [69] R. J. Prill, D. Marbach, J. Saez-Rodriguez, P. K. Sorger, L. G. Alexopoulos, and et al. Towards a rigorous assessment of systems biology models: the dream3 challenges. PloS ONE, page e9202, 2010. [70] A. T. Puig, A. Wiesel, and A. O. Hero. A multidimensional shrinkage-thresholding operator. In IEEE Workshop on Statistical Signal Processing, pages 113–116, 2009. [71] P. Qiu, A. J. Gentles, and S. K. Plevritis. Fast calculation of pairwise mutual information for gene regulatory network reconstruction. Computer Methods and Programs in Biomedicine, 94(2):177–180, 2009. [72] J. R. Quinlan. Improved use of continuous attributes in c4.5. Journal of Artificial Intelligence Research, 4:77–90, 1996. [73] R. J. Rolfes and H. Zalkin. Purification of the escherichia coli purine regulon repressor and identification of corepressors. J. Bacteriol., 172(10):5758–5766, 1990. [74] Y. Saeys, I. Inza, and P. Larranaga. A review of feature selection techniques in bioinformatics. Bioinformatics, 23:2507–2517, 2007. [75] J. Schäfer and K. Strimmer. An empirical bayes approach to inferring large-scale gene association networks. Bioinformatics, 21:754–764, 2005. [76] A. Scheinine, W. Mentzen, E. Pieroni, G. Fotia, F. Maggio, G. Mancosu, and A. de la Fuente. Inferring gene networks: Dream or nightmare? part 2: Challenges 4 and 5. Annals of the New York Academy of Sciences, 1158:287– 301, 2009. [77] C. Shannon. A mathematical theory of communication. The Bell System Technical Journal, 27:379–423, 1948. [78] A. Sîrbu, H. J. Ruskin, and M. Crane. Stages of Gene Regulatory Network Inference: the Evolutionary Algorithm Role, chapter 27, pages 521–545. 2011. [79] M. E. Smoot, K. Ono, J. Ruscheinski, P. Wang, and T. Ideker. Cytoscape 2.8: new features for data integration and network visualization. Bioinformatics (Oxford, England), 27(3):431–432, 2011. [80] G. Smyth. limma: Linear models for microarray data. In R. Gentleman, V. Carey, S. Dudoit, R. Irizarry, and W. Huber, editors, Bioinformatics and Computational Biology Solutions Using R and Bioconductor, Statistics for Biology and Health,, pages 397–420. Springer, 2005. [81] G. Stolovitzky, R. J. Prill, and A. Califano. Lessons from the DREAM2 challenges. Ann N Y Acad Sci, 1158:159– 195, 2009. [82] J. D. Storey. A direct approach to false discovery rates. J. Royal. Statist. Soc. B, 64:479–498, 2002. [83] J. V. Stoyanov, J. L. Hobman, and N. L. Brown. Cuer (ybbi) of escherichia coli is a merr family regulator controlling expression of the copper exporter copa. Molecular Microbiology, 39(2):502–512, 2001. [84] R. Tibshirani. Regression shrinkage and selection via the lasso. J. Royal. Statist. Soc. B, 58(1):267–288, 1996. [85] V. G. Tusher, R. Tibshirani, and G. Chu. Signicance analysis of microarrays applied to the ionizing radiation response. PNAS of USA, 98:5116–21, 2001. [86] A. Untergasser, H. Nijveen, X. Rao, T. Bisseling, R. Geurts, and J. A. M. Leunissen. Primer3plus, an enhanced web interface to primer3. Nucl. Acids Res., 35(suppl 2):W71–W74, 2007. [87] A. Visconti, R. Esposito, and F. Cordero. Tackling the dream challenge for gene regulatory networks reverse engineering. In AI*IA 2011: Artificial Intelligence Around Man and Beyond, XIIth Internationall Conference of the Italian Association for Artificial Intelligence, Palermo, Italy, 2011 (to appear). [88] M. Wu and C. Chan. Learning transcriptional regulation on a genome scale: a theoretical analysis based on gene expression data. Briefings in Bioinformatics, 2011. [89] A. Xiong, A. Gottman, C. Park, M. Baetens, S. Pandza, and A. Matin. The emrr protein represses the escherichia coli emrrab multidrug resistance operon by directly binding to its promoter region. Antimicrobial Agents and Chemotherapy, 44:2905–2907, 2000. [90] K. Yamamoto and A. Ishihama. Transcriptional response of escherichia coli to external copper. Molecular Microbiology, 56(1):215–227, 2005. [91] M. Yuan and Y. Lin. Model selection and estimation in regression with grouped variables. J. Royal. Statist. Soc. B, 68(1):49–67, 2006. 52 Nature Methods: doi:10.1038/nmeth.2016 [92] C. Zhu, K. J. Byers, R. P. McCord, Z. Shi, M. F. Berger, D. E. Newburger, K. Saulrieta, Z. Smith, M. V. Shah, M. Radhakrishnan, A. A. Philippakis, Y. Hu, F. De Masi, M. Pacek, A. Rolfs, T. Murthy, J. LaBaer, and M. L. Bulyk. High-resolution DNA-binding specificity analysis of yeast transcription factors. Genome Research, 19(4):556– 566, 2009. [93] H. Zou and T. Hastie. Regularization and variable selection via the elastic net. J R Statist Soc B, pages 301–320, 2005. 53 Nature Methods: doi:10.1038/nmeth.2016