BACKGROUND
current biology research generates an ever-increasing deluge of omics derived data.
each type of omics data pertains to a single level of the biological system under investigation .
while detailed knowledge of the individual genes, transcripts, proteins, metabolites and other cellular components remains important, understanding a biological system requires considering the networks into which these components are embedded and of their functioning as a  whole.
a major challenge in systems biology lies thus on developing effective and efficient methods to optimally extract the information contained in the aggregate of these datasets.
the increasing availability of genome-scale expression data has boosted the development of methods to infer the underlying regulatory networks.
a broad range of alternative methods are available, see  for reviews, and each of them uses different mathematical tools and/or biological assumptions.
a class of methods use differential equations to express transcript changes as a function of the transcript levels of the corresponding transcription factors  <cit> .
a second class of methods rely on bayesian networks to analyse the joint probability distributions obtained from the experimental data  <cit> .
other methods use the similarity between gene expression profiles to detect associations and to reconstruct a genome scale transcriptional regulatory network .
another class of methods use a combination of machine learning techniques to produce prioritized lists of transcription factors regulating each target gene  <cit> .
each method has different strengths and weaknesses, even methods using similar conceptual tools.
for example c3net  <cit>  uses mutual information  to reconstruct the core regulatory interactions in the network, that is, to recover the strongest interactions.
within this core, c3net was shown to outperform other methods also using mi such as clr  <cit>  and arance  <cit> .
each method provides different results, therefore much effort has been devoted to generate consensus networks from the multiple solutions.
it has been shown that in many instances, an integrative approach combining the outcome of each algorithm produces the best result  <cit> , however, a detailed analysis shows that for some interactions individual methods perform better than the consensus network  <cit> .
co-expression networks and module identification algorithms
co-expression networks contain genes as nodes and the edges of the network represent significant co-expression levels across the studied data set.
an open problem is still how these connexions are to be defined and how an adequate threshold is to be imposed  <cit> .
in differential network analysis two networks obtained using the same algorithm but alternative datasets are compared to identify interactions appearing only under a subset of conditions .
for example, such an approach was used to analyse prostate cancer datasets by comparing the networks inferred using c3net from normal or tumour samples datasets and was able to successfully identify cancer specific interactions  <cit> .
here, we construct and analyse co-expression networks extracted from the same dataset but using alternative algorithms to identify relevant interactions.
biclustering or module learning algorithms  <cit>  aim at the identification of functionally related genes showing co-expression patterns.
in some cases, such as the cmonkey algorithm  <cit> , additional biological information  is also considered.
the results of the genome-scale reconstruction methods can be displayed and visualized as one network whereas an identified module or bicluster contains a set of functionally related genes that might be differently regulated in different conditions.
therefore, the genes within a module might not form a cluster when considering all the conditions and would remain undetectable for network inference methods.
network inference methods and module prediction algorithms are highly complementary and new information can be obtained by combining the outcome of both approaches.
analysis of chip-seq data
there are experimental methods to directly reconstruct regulatory networks.
the knowledge of transcriptional regulatory events and specifically the transcription factors  binding sites can be greatly improved by chromatin immunoprecipitation and sequencing, chip-seq.
powerful techniques have been developed to process the sequencing data and to isolate the binding sites from the background noise generated by non-specific sequences, such as statistical tests, use of controls, techniques for signal de-convolution or lag-analysis among others  <cit> .
in the classical model of transcription regulation in prokaryotes a binding site in the promoter region of a gene is linked to a regulatory interaction between the corresponding tf and the gene.
however, the collection of chip-seq data for mycobacterium tuberculosis  hosted in the tb database   <cit>  shows that not all identified tf binding sites can be linked to regulatory interactions.
for example, in some instances the tf binding is much weaker than expected, in other cases this association is complicated by the occurrence of divergons, which are pairs of divergently transcribed operons or genes, or by the presence of binding sites within coding regions.
this lack of a one-to-one relationship between tf binding sites and regulatory events can be due to multiple reasons such as non-specific binding, cumulative effects of sites with weaker binding to regulate overall promoter affinity, specific binding generated in non-biological or non natural conditions, false positives as a result of the experimental procedure, or, simply errors in the genome annotation.
we propose that it is through the integration with additional data sources, specifically through the integration with expression data, that this challenge can be overcome.
integration of heterogeneous molecular networks
each inference algorithm has its weaknesses and strengths and each network holds its own intrinsic value and can be used to gain more specific insights on various aspects of the biological system  <cit> .
in a sense, using multiple algorithms to extract the networks and presenting them to the user is a similar approach to that followed by annotation pipelines, such as microscope  <cit> , diya  <cit>  or basys  <cit>  among others.
these pipelines present the user a list of gene centred information, with different annotation sub-fields, to enable supervised annotation.
integration efforts require a common layout for the data, therefore, we have chosen to represent the available layers of information, such as operon structure, known interactions between genes or proteins, enzymatic activity  or functional similarity, among others through network representations.
the networks themselves are represented in a common format, xgmml , that allows their simultaneous exploration.
in addition to the simultaneous visualization of different networks, additional analysis tools such as motif search and identification, go-enrichment analysis, and tools to overlay expression data or analyse expression profiles of multiple genes across different conditions can be used.
we have developed a pipeline for the generation of co-expression networks that is easily tunable to produce alternative networks .
the pipeline is based on using the similarity between gene expression profiles to detect associations and contains a higher order extension of the data processing inequality method  <cit>  to reduce the number of redundant or possible spurious links.
this pipeline is presented in figure  <dig> figure  <dig>
schematic of the pipeline to obtain co-expression networks.
from top to bottom the following steps are applied:  calculation of the similarity matrix,  z-transformation,  combine i,  threshold setting  inequality simplifier and  combine ii.
note that when applying the inequality simplifier to the zrc network the result will be a tree.
through the exploration of the multiple networks, the user gains new knowledge of the biological system, but it can also lead to the discovery of new strategies to integrate the information stored in the networks.
these newly gained strategies can then be translated into new pipelines.
as result of our exploration, we have developed a pipeline to uncover multiple layers of regulation, presented in figure  <dig> and another one to analyse chip-seq data and assign regulatory interactions to the detected binding sites, presented in figure  <dig> figure  <dig>
pipeline to uncover additional regulatory layers.
step 1: identify conditions linked to the main regulatory event for the initial gene set.
this can be done using biclustering techniques or by direct comparison with the expression levels of the regulator .
step 2: build co-expression networks in the remaining conditions.
step 3: identify the closest neighbours of the selected genes in the new networks.
step 4: iterative round of motif identification/matching to identify the secondary motif and the set of genes with this motif in their upstream regions.
pipeline to analyze chip-seq data.
after the locations of the chip-seq binding sites have been retrieved, their genomic context is analysed.
a core set is defined by selecting targets with i) literature evidence or ii) a hit in the upstream region of not divergently transcribed genes.
the expression levels of these genes are analysed and they are categorized through clustering.
finally the rest of the putative targets are assigned to these groups  based on the similarities of their expression patterns.
show case: deployment of the framework to unravel regulatory mechanisms in mycobacterium tuberculosis
we have analysed regulatory events in mtb.
due to its implications to human health this highly successful pathogen has been extensively studied and there is already a substantial body of information on mtb and its underlying regulatory networks, but still much remains to be learned.
we have analysed not only networks extracted from literature  but also networks extracted from publicly available databases such as string  <cit> , metacyc  <cit> , kegg  <cit> , tbdb  <cit>  and tuberculist  <cit> .
expression data from publicly available repositories and corresponding to  <dig> perturbations, have been used to generate the co-expression networks and to further analyse chip-seq data  <cit> .
we have finally combined the inferred networks with the outcome of biclustering algorithms, to further explain the functionalities of the modules.
we have analysed the regulation of dna repair systems within mtb, particularly we have focussed on two alternative regulatory mechanisms: the ‘reca-lexa dependent dna repair system’ and ‘reca non dependent regulation’.
we will show that the integration of expression data and its use to guide the exploration of upstream sequences allows to analyse alternative regulatory mechanisms for the same set of genes.
in addition, we have explored some of the regulatory mechanisms that allow mtb to survive within the host, such as devr  regulon, which is a key element for understanding the dormant state, and the regulation of the response to changes in zinc availability .
methods
gene expression data
<dig> two-colour microarrays for mtb  were retrieved from the gene expression omnibus database  <cit> .
<dig> of them aimed to capture the effect of  <dig> drugs targeting metabolic pathways  <cit>  whereas  <dig> captured stress induced dormancy in the wild-type and in devr activation genes knockout mutants  <cit> .
we followed boshoff et al.
<cit>  to classify the conditions in the compendium into  <dig> categories according to the experimental perturbation:  aromatic amides intracellularly hydrolyzed, low ph;  mutants of devr activation pathways;  translation inhibition;  acidified medium;  cell wall synthesis inhibition;  respiration inhibition ;  nutrient starvation;  dna damage;  transcription inhibition;  iron scavengers;  multiple stress sources applied simultaneously ;  minimal medium ;  not classified; and  conditions associated with devr upregulation.
a list of the used datasets is presented in additional file  <dig>  we applied a common normalization method, loess, to the arrays from each experiment and each of the six independently designed platforms.
linear models were constructed to consider biological and/or technical replicates ; within-array replicate probes were averaged; and an additional between-array quantile-quantile normalization was performed to ensure the comparability between experiments.
these manipulations were performed using the r limma package  <cit> .
a common locus tag format was introduced and missing values were filled up using knn-imputation from the impute r package  <cit>  with k =  <dig> and eliminating genes and experiments with more than 30% and 50% missing values respectively.
the resulting compendium contained information on  <dig> open reading frames across  <dig> different conditions .
biclustering
we adapted the original cmonkey r code to mtb data  <cit>  and we increased the number of biological information sources that the algorithm can consider.
for the biclustering process, we selected a subset of the interactions present at the string database  <cit> .
we selected specifically interactions obtained based on co-occurrence of linked proteins across different species, on curated databases, and on their association in the abstracts of scientific literature.
from prolinks  <cit>  we collected interactions obtained by the phylogenetic profile method.
in addition, we included the network based on the similarity among the annotated go terms , which was constructed using the sleipnir library  <cit> .
upstream sequences for mtb genes were retrieved using rsat  <cit>  and used for the motif detection and motif matching steps.
for the automatic biclustering algorithm, we used the 1000 bp region upstream of the translation start site, avoiding overlapping with upstream neighbour genes when present.
in the initial rounds we used a reduced matrix where only the leaders of the operons, as defined by roback and co-workers  <cit> , were kept.
multiple runs of the algorithm were performed, which used the default parameters except for: initial size of seed clusters ; number of iterations ; maximal number of clusters .
we obtained a set of  <dig> biclusters.
the jaccard similarity coefficient between each pair of biclusters was computed .
to reduce the redundancy in the set of biclusters, we merged the pairs that showed a jaccard similarity higher that  <dig> .
the biclusters in the merged set were enlarged by adding genes based on the predicted operon structure combined with the expression level measurements.
this new set was used to seed the biclustering algorithm in two subsequent optimization rounds one biased towards the detection of concurrent motifs in the upstream regions and a second one biassed towards the identification of sets linked to highly related biological processes .
we performed an additional manual merging step, that also considered the similarity between the detected motifs using the matrix comparison tools from rsat  <cit>  and we obtained a final set of  <dig> biclusters.
expression plots for these biclusters are in additional file  <dig>
pipeline for the generation of co-expression networks
the pipeline to generate co-expression networks is presented in figure  <dig>  departing from an nxm matrix containing the expression profile of n genes in m conditions, we construct a symmetric nxn similarity matrix s. the pipeline allows to choose between correlation  and mi as similarity measurements.
mi is computed using an estimator based on the entropy of the empirical probability distribution with initial data discretization into n  equal sized bins from the bioconductor package minet  <cit> .
in the second step, a z-score transformation is performed on the distribution of the scores from the similarity matrix.
the transformation is done by rows, to obtain the zr matrix; by columns, to obtain the zc matrix; or by both rows and columns simultaneously, to obtain the zrc matrix.
the z-score transformation allows for each possible interaction to be weighted regarding the background of interactions in which each member of the interacting pair is involved  <cit> .
therefore, in the zc case each element sij becomes zcij, which is the z-score of sij regarding the distribution si <dig>  si <dig>  …, sin.
in the zr case each element sij becomes zrij, which is the z-score of sij regarding the distribution s1j, s2j, …, snj.
in the zrc case each element sij becomes  <dig> zcij2+zrij2/ <dig>
in the following step the matrices are combined into two new matrices: cλ and rλ: cijλ=maxzcij,λ·zrcijrijλ=maxzrij,λ·zrcij with λ  a positive real number and i and j denoting the rows and columns of each matrix.
λ allows to fine-tune the results, since higher values of λ will lead to cλ and rλ matrices that are similar to each other and to zrc.
considering λ = √ <dig> means that, for those cases where one of the elements in eq.
<dig> is zero, both λ zrcij and zcij will be identical.
however, using the default λ =  <dig> will ensure that when the values of zcij and zrij differ then the highest one will be selected either through zcij  or zrcij.
in the following step a threshold can be applied to remove interactions with low weight .
any value below the threshold is set to  <dig>  this step produces a more sparse network, which might be needed to obtain a neat visualization.
the default threshold value is chosen so that the number of non zero edges is equal to a predefined value .
afterwards the inequality simplifier can be used to remove, possibly spurious, links from the network.
the inequality simplifier is an extended version of the data processing inequality .
in simple terms, the dpi states that given two interdependent random variables a <dig> and a <dig> and a third one a <dig> that only depends on one of them, for example a <dig>  then a <dig> cannot contain more information about a <dig> than a <dig> does.
this statement is mathematically represented by the following inequality among mi values mia <dig> a3≤minmia <dig> a <dig> mia <dig> a <dig>
the dpi must hold whenever a <dig> does not depend on a <dig>  so it can be used to remove spurious interactions from the network  <cit> , possibly caused by feed forward loops or mutually dependent regulators.
the spurious links are removed according to: ifmia <dig> a3≤mia <dig> a2&mia <dig> a3≤mia <dig> a <dig> thensetmia <dig> a3= <dig>
the dpi can be extended to higher-order interactions, therefore allowing a recursive implementation  <cit> .
the dpi is derived from the triangle inequality satisfied by any metric or distance, such as the mi or kendall rank distance.
our extension to remove spurious edges by considering higher order interactions contains two parts: the first part is the identification of alternative pathways connecting the same nodes and the second part is to apply the inequality simplifier to decide whether a link is spurious and should be removed.
the first part is based on dijkstra's shortest route algorithm  <cit> .
here, instead of searching for the shortest route between two nodes , we use it to find the alternative path between those two nodes with the best throughput.
the throughput of a path is defined as the value of the edge with the lowest throughput  value.
given two connected nodes, the link between them is considered spurious and removed if an alternative path is found with a higher throughput value.
for a given pair of nodes, ai and aj, , so that the best alternative path of length  passes through a <dig>  a <dig>  …, an- <dig> and an, the following rule  is applied: ifmiai,aj≤minmiai,a <dig> mia <dig> a <dig> …,mian− <dig> an,mian,ajthensetmiai,aj=0wheren∈ <dig> …,n− <dig>
this rule is applied for each possible pair of nodes and for each possible value of n. the alternative pathway joins ai and aj and does not contain any cross-linking, so that therms such as a2a <dig> or a1a <dig>  do not need to be considered.
this is a consequence of using the path with the best throughput so that no triangle inequality of the metric is required to derive this property.
figure  <dig> shows how the inequality simplifier acts to remove, possibly spurious, links from the network when applied to higher order terms.
while applying this rule, an additional nxn matrix, t, is built to keep track of how many links a particular edge has caused to be removed.
in this case, the links aia <dig>  a1a <dig>  a2a <dig>  …, an-1an and anak have caused the removal of link aiaj, therefore, in the t matrix, the elements tai,a <dig>  ta <dig> a <dig> ta <dig> a <dig> …, tan- <dig> an and tan,aj are increased by 1∕.figure  <dig>
the inequality simplifier.
the similarity values among the nodes  connected by the different edges have been indicated.
dotted lines represent the spurious links removed by the inequality simplifier.
left application of order two, which is equivalent to a direct application of the dpi.
right higher order application, both dotted lines are removed.
the dpi would only remove the blue dotted line.
finally the two networks cij and rij are combined into one final network rcij.
rcijλ=maxzcij,zr
the presented pipeline can be used to obtain different outputs, such as the rc or the zrc networks and their variants through the application of thresholds or the inequality simplifier.
the final output is a square matrix that represents a network through a weighted adjacency matrix.
this adjacency matrix is finally exported into a tabular format that can be imported into cytoscape  <cit>  to obtain a graphical representation.
the weights in the adjacency matrix or the t matrix produced by the inequality simplifier can be imported as edge attributes and used to set the thickness of the edges.
finally, a force directed algorithm can be used to generate the appropriate layout and the graphical representation of the network, which then can be exported to an xgmml file.
applying the inequality simplifier to a symmetric matrix, such as zrc, removes any possible loop in the network, therefore it results in a forest graph as depicted in figure  <dig>  applying the inequality simplifier to an initial matrix containing around  <dig> nodes takes at most five minutes in a standard desktop computer .
the networks used for our exploration of gene regulation networks in mtb were the result of applying the inequality simplifier to zr, zc and zrc .
belcastro et al.
<cit>  have shown that from the total number of possible interactions a selection of 5%  would be a sensible choice according to the assumption that biological network are sparse.
here, one of our goals is to obtain a clear visualization of the network, therefore, we have chosen an even lower number  so that the number of nodes included in the network is as high as possible without having a too crowded visual representation of the network.
an additional threshold was imposed so the number on non zero edges of zr, zc and zrc networks is equal to  <dig>  the zrc derived network contained  <dig> nodes and  <dig> edges whereas the mixed network contained  <dig> nodes and  <dig> edges.
these manipulations were implemented in r and java, the corresponding scripts can be found, together with an example data set in additional file  <dig>
additional networks and analysis tools
metabolic networks
the nodes in the metabolic networks are of two different types: reactions and metabolites, and additional information, such as the stoichiometry, directionality and the catalysing proteins of the reactions, also has to be stored.
we kept this information by using a combination of two file types: a svg  file that contains the graphical representation of the networks and an rdf  file containing the additional information .
other networks
a network formed by multiple disjoint interlinked clusters was generated using the predicted mtb operon structure  <cit> .
in addition, the information stored in the string database  <cit>  was used to generate networks based on: known protein-protein interactions; distance of the genes in the genome and that of orthologous genes in other species; co-occurrence in other species; orthologous genes being fused together in other species; association in the abstracts of scientific literature; and other databases, such as metacyc  <cit>  and kegg  <cit> .
to identify regions duplicated in the genome or homologous genes, a network was built by considering the sequence similarity between each pair of genes.
the similarity was computed using megablast  <cit>  and an edge was created for matches with e-values lower than of 10− <dig>
sequence analysis tools
we have used meme  <cit>  for motif elicitation within a set of genes and fimo  <cit>  to find occurrences of the chosen motif.
a post processing step ranks, using the fimo q-value, the genes  for each identified motif .
an iterative cycle of motif identification and motif matching can be established to further refine each motif.
the computational time required by each iteration greatly depends on the length of the upstream sequences, which has to be selected for each particular analysis.
for the selected show cases the motifs accumulated in the 200 bp region upstream of the translation start site.
gene ontology  enrichment analysis
on selected set of genes was performed using a hypergeometric function to model the probability density, as implemented in the gohypergall r function  <cit> .
the go annotation for mtb was obtained from the uniprot-goa database  <cit> , the mtb-goa database , tbdb  <cit>  and the recent re-annotation of the mtb genome  <cit> .
operons
we developed a putative operon extension tool which is only based on the gene orientation, so that genes that are downstream of the selected gene  and with the same orientation are added to the putative operon until a gene with the reverse orientation is found.
the user is then free to decided whether the expression data sustains the extension of the operon.
this approach allows to easily combine expression information with genomic context information.
venn diagrams have been created using the utilities from the venndiagram r package  <cit> .
expression plots
generated with a solid line representing the mean expression of the selected genes in the conditions included in the compendium and dots marking the conditions in the bicluster where the genes show high correlation .
the classification of the experimental conditions, previously mentioned, was included in the expression plots through an horizontal colour line that allows to quickly associate the behaviour of the genes with the type of perturbation.
to select the conditions  in which the selected genes show co-expression there is an initial step where all conditions are included in the bicluster, then an iterative loop starts were each condition is removed at a time and the new correlation values are computed, finally the condition leading upon removal to the highest correlation values in the remaining set, is removed from the set.
this process is iterated until the correlation in the remaining set of conditions is higher or equal to  <dig> .
discovery of additional regulatory layers
whenever two alternative regulators regulate the expression of genes in a given set, it might be that one overshadows the detection of the other.
we have developed the work flow to analyse these events .
given a set of genes under control of two regulatory events, the first step is to identify the conditions were the first event is active or taking place.
how these conditions are identified depends on the particular example.
one way would be to identify conditions with up/down regulation of the corresponding transcription factor.
to identify conditions were reca expression is not regulated by the reca/lexa mediated mechanism we built a model linking the time dynamics of reca to lexa expression levels using the inferrelator algorithm  <cit> .
the output of the algorithm is a prediction or fit of reca levels based on lexa levels .
those conditions with a poor agreement between the measured and fitted levels are the conditions  without lexa mediated induction of reca.
once the new set of conditions has been selected, the following step is to reconstruct the co-expression network, using only the expression data corresponding to this subset of conditions and analyse the location of the original set of genes in the newly built network.
furthermore an iterative round of motif identification and motif matching can be run to identify a possible motif in the upstream region of the cluster.
this iterative round can be performed either using cmonkey  <cit>   or manually using alternatively meme  <cit>  and fimo  <cit> .
analysis of chip-seq data
publicly available chip-seq data were obtained from tbdb  <cit> .
we have developed the workflow, presented in figure  <dig>  to analyse these data.
the chip-seq data had already been analysed with a peak calling algorithm.
for each of the considered tfs a list of predicted targets was obtained from tbdb, together with information about the peak location.
a reduced list, or core set, was obtained by selecting those targets that show a hit in their adjacent upstream intergenic regions.
an additional filtering was done to select only those hits where the peak is not flanked by divergently transcribed genes/operons.
additionally, the genes for which literature evidence supported the regulatory interaction with the tf were included in the core set.
in the following step, the expression levels of the genes in the core set were analysed and the matrix of correlations across the conditions in the expression compendium was computed.
the appearance of negatively correlated genes in this set is a signal of a dual repressing/activation function of the regulator and therefore two  subsets can be defined by hierarchical clustering of the set.
this process can be further enlarged to encompass additional groups that would be linked to alternative regulatory mechanisms or to the effect of additional tfs.
once this/these group were defined, the rest of the putative targets predicted by chip-seq were assigned to either one of these groups based on their average correlation with the members of the group .
topological overlap
the topological overlap between two genes i, j in an unweighed network is defined as  <cit> : ti,j=ni∩njminni,nj+1−ai,ji≠j1i=j
where n is the set of direct neighbours of gene x ; |x| represents the number of elements in set x and ai,j is the adjacency matrix of the network .
the topological overlap is bounded between  <dig> and  <dig>  two genes will have high topological overlap if there exists a connection between them and if they are connected to the same group of genes.
for weighted networks, ai,j represents the weight of the interaction between genes i and j and takes continuous values between  <dig> and  <dig>  in these cases, the previous formula can be generalized to  <cit> : ti,j=Σkai,kak,j+ai,jminΣkai,k,Σkaj,k+1−ai,j
the topological overlap of a group of genes was defined as the average of their mutual topological overlap, which was computed using the r wgcna package  <cit> .
only groups with more than  <dig> genes were considered.
an empirical p-value for these scores was calculated by randomly sampling  gene groups of the same size in the respective network.
network visualization
we have developed a visualization tool that allows the simultaneous visualization of networks in xgmml format and the sharing of identifiers between them.
technical characteristics of this tool are available in additional file  <dig>
RESULTS
co-expression networks
we have developed a pipeline to generate co-expression networks that allows for a myriad of possible networks.
choosing a subset of them highly depends on the available data and the process the user wants to explore, since different clusters appear in each of them.
this reflects the inherent modularity of biological networks and the dependency of regulation on the chosen conditions.
the pipeline presented in figure  <dig> includes some already well established and tested methods for genome-scale network inference.
for example, the similarity matrix computed using either correlation or mi is commonly used to analyse gene expression data  <cit> , but also other methods are contained in the pipeline.
when working with mutual inference as a similarity measure the matrix denoted zrc corresponds to the output of the clr algorithm  <cit> .
additionally, applying the inequality simplifier only up to order two on a mi matrix amounts to using the dpi to prune spurious interactions from the networks, which is a key element of the aracne algorithm  <cit> .
a comparative analysis of methods to reverse-engineer transcriptional regulatory networks was done using the results of the dream <dig> challenge.
in this challenge the teams had to reconstruct genome-scale transcriptional regulatory networks from expression data.
the networks proposed by the different teams, were evaluated through their comparison with a gold standard, a set of experimentally verified regulatory interactions in the target organisms.
in addition, marbach et al.
<cit>  constructed a consensus network integrating the predictions from the multiple methods.
whole network performance estimators such as area under receiver operating characteristic curve  or area under precision recall curve  show that over the entire network the consensus network outperforms individual methods due to their inherent complementarity.
however, for some interactions individual methods perform better than the consensus network .
the loss of information when building the consensus network does not affect equally the different transcription factors .
for example, for purr or lexa, a significant fraction of the total number of known interaction is better recovered by the individual methods than by the consensus network.
to evaluate the performance of our algorithm we have generated co-expression networks corresponding to the synthetic, escherichia coli and saccharomyces cerevisiae datasets used in the dream <dig> challeng.
for each dataset, two networks  were built .
in the challenge the goal was to identify links between regulators and target genes.
in co-expression networks no special emphasis is done on the regulators.
to evaluate the networks we have calculated the topological overlap  <cit>  of the known targets of the transcriptions factors .
the topological overlap measures how similar the neighbourhoods of two genes are  <cit> .
for 95%  of e.coli tfs considered  the topological overlap of the target genes in the zrc network is significantly higher than for the overall network , which means that they form cohesive modules in the network.
the relative number of cohesive clusters identified for the other datasets  are lower  overall the zrc networks maintain a higher degree of cohesiveness than the mixed networks .figure  <dig>
topological overlap of tf targets in the
e. coli
zrc co-expression network.
the zrc co-expression network was reconstructed using our pipeline using e. coli expression data from the dream <dig> challenge  <cit> .
only the  <dig> tf with more than  <dig> experimentally verified targets  were considered.
dashed line represents the average topological overlap in this network .
systematic analysis of regulation in mtb.
integration and simultaneous exploration of heterogeneous networks
we considered co-expression networks obtained with different approaches.
additional information is represented through networks, such as information on operons, go annotation similarity, data base stored knowledge and sequence homology among others.
this integrative approach allows to systematically explore functional modules in the network and it is highly complementary to existing bi-clustering and other module identification methods.
even when additional biological information is included, biclustering methods have to face the challenge of interpreting the function of each groups of genes.
therefore, a general approach to detect and understand functional modules in a given organism is to simultaneously explore the location of the genes in a given bicluster in the co-expression networks.
we have done so a to analyze the output of the cmonkey algorithm  <cit>  on mtb data.
we have run multiple instances of the algorithm to bias the search towards the detection of both putative co-regulated and functionally related sets of genes.
the obtained bi-clusters are presented in additional file  <dig>  to further investigate these biclusters, we have projected them in the multiple networks and genes proceed to their in-depth analysis.
an example of how this is done is shown in the following section when we analyse the regulation of dna damage repair systems in mtb by lexa.
in a given co-expression networks it might very well happen that the effect of one regulator on a given set of genes overshadows the identification of other regulators.
however, these mechanisms can be uncovered when comparing different networks generated under different experimental conditions.
we have developed the pipeline presented in figure  <dig> to uncover additional layers of regulation based on the assumption that regulation proceeds differently upon different perturbations, therefore some regulatory events will only be detectable in subset of conditions.
this pipeline is only useful when the user suspects or has prior knowledge of an additional regulatory interaction affecting the same set of target genes.
therefore, we have used it to further explore the regulation of dna repair systems in mtb, where additional regulators are known to exist  <cit> .
in addition, our approach allows to systematically explore and interpret additional data such as chip-seq data, through the pipeline presented in figure  <dig>  among the available data we have chosen to focus on one of the key subsystems key for mtb survival in the host, the regulation of the response to hypoxia and the induction of dormancy program via devr.
furthermore, the simultaneous network visualization, allows the exploration of common trends in the networks.
instead of focussing on clusters of genes appearing in the co-expression networks and how these clusters are reflected in, for example, the network liked to the go biological process annotation.
the alternative approach can be taken and for each set of genes linked to a particular go or cog category we can trace back their location in the network.
this analysis can point to interesting effects.
for example, we explored the genes coding for ribosomal proteins and in the different co-expression networks these genes appear to form a highly interlinked cluster.
however, in all versions of the networks, the rpmb2-rpmg1-rpsn2-rpsr <dig> operon  formed by genes coding for ribosomal proteins  appear as a separated set .
the combination of this information with the network of interactions extracted from literature  <cit>  led us to the analysis of the zinc uptake regulator, zur, and its targets.
dna repair systems in mtb
reca-lexa dependent dna repair system: the sos box
lexa is a repressor known to be involved in the control of mechanisms for dna repair.
under neutral conditions lexa, dimerizes and binds to the sos box, repressing its own expression and of other genes related to dna repair.
the consensus sequence of the sos box for mtb had been identified as the palindromic motif tcgaac4gttcga  <cit> .
upon dna damage, reca binds single stranded dna .
the complex reca-ssdna, stimulates autocatalytic cleavage of lexa, so genes repressed by lexa are induced.
acidic conditions trigger a similar response, since lexa can no longer dimerize, effectively preventing it from binding its target dna sequence  <cit> .
the iteratively mapping of the genes in each of the biclusters presented in additional file  <dig> to the multiple co-expression networks, showed that bicluster  <dig> is a  <dig> gene module that shares many of the characteristics of the lexa regulon.
to further investigate this bicluster, we proceed to analyse the location of these genes in the multiple networks and to identify genes that might be linked to them through the systematic analysis of their expression patterns .
for each of the candidates to be included in the regulon we combined information from the different information layers: literature, databases, functional annotation, and upstream sequence among others.
we were able to identify a total of  <dig> genes putatively in the regulon, listed in additional file 9: table s <dig>  this includes the  <dig> genes reported by davis et al.
<cit>  plus  <dig> additional ones, resulting in an effective enlargement of the regulon size by 75%, .
based on this extended regulon, we identified for the sos box in mtb a more specific motif: mkwmtcgaamrywtgttcga  plot of the average expression level of the members of the lexa regulon across the different conditions.
red dots mark conditions with high  correlation between the genes in lexa regulon.
the horizontal bar and its different regions indicated by numbers refer to the classification of the conditions as described in materials and methods.
high expression levels are observed in conditions corresponding to low ph or uv light.
b) clusters of genes involved in dna repair mechanisms in the co-expression network .
genes regulated by lexa are marked red.
c) refined lexa identified binding motif, positions  <dig> and  <dig> were previously non specific.
d) number of genes identified to be regulated by lexa.
previous indicates genes previously reported in the literature as lexa regulated  <cit> , whereas automatic refers to the genes initially identified by the automatic biclustering algorithm.
to verify our predictions about the  <dig> genes not previously assigned to this regulon, we compared our predictions with the lexa binding regions identified by smollett et al.
<cit> .
in that work chip-seq analyses of lexa binding sites was complemented with experimental measurements of gene expression upon dna damage induced by mytomycin c. the  <dig> genes that we assign to the lexa regulon contain a site in their upstream regions where lexa binding was detected.
in addition, these genes also show dis-regulation in dna damaging conditions.
it is important to stress that the results from smollett et al.
were used in our analysis only to verify our predictions.
therefore, we can conclude that we have successfully reconstructed the regulatory interactions of lexa.
davis et al.
identified a single sos box in the upstream region of the divergently transcribed genes whib <dig> and fbia, however none of them were listed as likely to be regulated by lexa.
the analysis of the expression profiles of these two genes, fbia and whib <dig> led us to conclude that fbia, a probable 2-phospho-l-lactate transferase involved in coenzyme f <dig> biosynthesis, does not belong to the regulon, since no significant expression changes are observed upon conditions related with induction of the lexa regulon .
the expression pattern of whib <dig> gene shows anti-correlation with the rest of the genes in the lexa regulon in all the conditions in our compendium that show upregulation of lexa.
therefore, we concluded that lexa regulation of whib <dig> expression proceeds through a different mechanism than the previously described so that the dimerized form of lexa acts as an inductor of whib <dig> expression.
this has also been independently confirmed in experimental conditions were mytomycin c was added to the medium  <cit> .
whib <dig> has been hypothesised to be involved in the regulation of cell division  <cit> .
functional analysis of the genes co-expressed with whib <dig> further supports the idea that whib <dig> is involved in the regulation of cell division.
we also found that not only rv2719c, as reported by davis and co-workers, but the complete rv2719c-nrdr-rv2717c operon is under the control of lexa.
information from the string database  <cit>   allows us to functionally link rv2717c with dna repair and cell wall synthesis; nrdr is involved in the control of the synthesis of dntp needed for dna replication and/or dna repair  <cit> ; and the cell wall hydrolase rv2719c is involved in suppressing the cell cycle by altering localization of ftsz rings  <cit> .
therefore, our analysis has unveiled that repression of whib <dig> and induction of the rv2719c-nrdr-rv2717c operon are the lexa regulated mechanisms to temporary arrest cell division upon dna damage.
furthermore, through the re-annotation of the newly discovered lexa targets, we can identify a putative functional mismatch repair  system.
no mr system has been previously identified in mtb.
a typical bacterial mr system contains: the mismatch-recognition protein muts that contains an hnh endonuclease domain; muth, a nicking endonuclease; and mutl, which acts as a scaffold between these.
the mr system additionally contains a dna helicase; a dna exonuclease to remove the mismatching nucleotides; and a dna polymerase together with a dna ligase to repair and ligate the created gap  <cit> .
the hypothetical mr system that we have identified is formed by: i) some of the  <dig> hnh endonuclease domain containing genes that belong to the 13e <dig> family, ii) the holiday junction dna helicases ruva and ruvb together with ruvc, a crossover junction endo deoxyribonuclease, and iii) the dna polymerase dnae <dig>  together with imua/b that are essential for its function  <cit> .
the existence of this repair system, under the control of lexa, would solve the apparent inconsistency between the low mutation rates in mtb and the absence of an mr system  <cit> .
additional regulation: reca independent dna repair system
reca is key to the correct regulation of the lexa regulon and it is also lexa regulated.
however, additional dna repair mechanims have been described in mtb, particularly the reca non-dependent  dna repair system  <cit> , that also regulates reca expression in a lexa independent manner.
to analyse this additional regulatory layer, we have used the pipeline in figure  <dig>
the conditions not linked to the main regulatory event are those  <dig> where no relationship was found between reca and lexa levels .
using these conditions we built a new co-expression network and afterwards we compared the original and the newly built networks.
the members of the lexa regulon appear in the original network as a tight cluster, however, in the second network only a subset of them appear clustered.
we selected the genes in this small cluster and proceed to an iterative round of motif identification and motif matching, to finally identify the genes regulated by the reca_nd mechanism.
our approach does not allow us to identify the regulator of the set, but previous studies point to clpr   <cit> .
we have verified our predictions for the identified motif and the list of targets genes,  through comparison with literature data, since they match those previously described by gamulin et al.
<cit>  for the reca_nd dna repair system.
however, there are two striking differences between our results and those previously reported by gamulin et al.
we don’t find that sigg responds to dna repair, which matches the result reported by smollett et al.
<cit> , on the other hand, we do find another sigma factor, sigc, that seems to be regulated in response to dna damage.
the list of genes in the regulon, show that in this case, regulation of cell cycle arrest upon dna damage is not linked to whib <dig>  but only to the rv2719c-nrdr-rv2717c operon.
hypoxia and induction of dormancy program: devr regulon
one of the main characteristics of mtb is its ability to switch to a non replicating or ‘dormant’ state that allows it to survive for a long time within the host and renders it less susceptible to antibiotics.
the environment inside the host is hypoxic and it might have high concentration of co or no released by the host macrophages.
under either low oxygen concentrations or high concentrations of no or co, the heme iron from the kinases devs  and dost becomes ferrous, the kinases become co or no bound and they get activated.
in the active form, devs and dost autophosphorylate and induce phosphorylation of devr, which in turn, can bind its dna recognition sequence and induce expression of the devr regulon resulting in the activation of the dormancy program  <cit> .
all the options to build the co-expression networks that we have explored, result in a tight cluster for the known members of the devr regulon.
this is to be expected since dormancy and devr regulon induction have been an active research topic in the mtb field.
from the  <dig> distinct conditions present in the expression compendium, almost 10% of them  are conditions associated with expression of devr regulon.
however, the fact that these genes always appear forming a tight cluster, points to the absence of additional regulatory elements that might cause a differential expression of some members of this regulon.
we have selected this regulon to validate our methodology for chip-seq data analysis.
chip-seq data corresponding to over-expression of devr were obtained from tbdb  <cit> .
the processed results for devr contain  <dig> detected peaks, that correspond to  <dig> genes that could be possibly regulated by devr, although in our compendium expression data for only  <dig> of them were available.
we have analysed this dataset following the methodology shown in figure  <dig>  and defined a core-set containing  <dig> genes.
the analysis of the correlations in the expression profile among this set as compared to the overall distribution indicates a common regulatory influence over the selected genes, figure 7a.
further analysis of the genes in this core set and their behaviour across the conditions in the compendium lead us to identify five distinct groups of genes within the identified targets .
in four of these groups there is a high correlation among the genes whereas no clear pattern can be identified in the expression of the genes in the fifth group.
one of these groups, that from now on we will refer to as the devr regulon, contains  <dig> genes that were identified by chip-seq, has an expression pattern consistent with the previously described devr regulon, and  <dig> of these genes have been previously identified as devr regulated genes .
additionally, we have found among the list of targets from tbdb,  <dig> genes that have been previously reported as devr regulated  however their expression patterns suggest that either they are not regulated by devr or there is a secondary regulatory event altering their expression levels, therefore it is arguable whether they can be included in the regulon.
a detailed list of members of this regulon is provided in additional file 13: table s <dig>  the genes in this table have been assigned to functional categories: cell wall, transport elements, anaerobic respiration, translational machinery, regulatory elements and elements related to stress response.
these six elements are linked to some of the main changes observed during growth arrest and dormancy, such as the changes in the cell wall, the arrest in protein synthesis and the adaptation to a hypoxic environment with reactive nitrogen species  <cit> .figure  <dig>
devr regulon.
a)
left: histogram in blue represents the correlation among all the genes present in our compendium, whereas the histogram in green is based upon the correlations of the identified targets for which expression data is available within our compendium .
both show the same overall distribution.
right: histogram in blue represents the correlation among all the genes present in our compendium, whereas the histogram in green is based upon the correlations of the  <dig> genes selected in the core set.
note there is a shift towards positive correlation values, pointing to a common regulatory influence over the selected genes.
b) number of genes identified in the devr regulon compared to the number of targets identified through chip-seq experiments or the ones cited in literature .
c) group assignment of the  <dig> targets identified by chip-seq for devr.
g <dig> contains genes for which non discernible expression pattern has been found.
g <dig> correspond to the usually named devr regulon, whereas g2- <dig> contains genes that show correlated expression patterns, although these patterns are not consistent with the previously described behaviour of devr regulon.
a detailed list of these genes is available in additional file  <dig> and the output of the go-enrichment analysis is shown in additional file  <dig>
interestingly, we found a faint link between the type vii secretion system  genes rv0282-rv <dig> and devr.
although the correlation analysis shows that they are not members of devr regulon, we indeed find that, in a reduced set of  <dig> conditions they show a high  correlation with devr regulon .
these genes are also known to be regulated by zur   <cit>  and ider   <cit>  and are required for mycobactin-mediated iron acquisition  <cit>  in mtb.
rv1734c belongs to the set of genes that had been previously reported as devr regulated but for which our analysis shows that their expression pattern is not compatible with this assertion.
our analysis is based on transcript levels, so it could be that the regulation of rv1734c happens at the post-transcriptional level.
however we have found additional evidence, albeit indirect, of our prediction.
chauhan et al.
<cit>  analysed the effect of mutations in the devr binding motif and found that positions  <dig> ,  <dig>  and  <dig>  are essential and a substitution in any of them dramatically reduces the binding affinity .
in addition, our analysis shows that position  <dig> always contains an a. motif analysis complements the results from the analysis of expression data and further supports our prediction that rv1734c is not regulated by devr, since its putative binding site contains a mutation  at position  <dig>
among the targets identified by chip-seq four other groups emerge  however none of these groups show significant up or down regulation in the conditions associated with devr regulon induction.
no clear common behaviour can be detected among the genes in group g <dig>  these could be considered as false positives as a result of the devr over expression performed prior to the chip-seq procedure  <cit> .
however the genes within each of the other groups  show a consistent co-expression across the conditions in the compendium, therefore we believe that these hits should not be discarded as being caused by non-biological reasons, instead alternative explanations such as weaker binding to regulate overall or regulation through a transcription factor homologous to devr should be further explored.
in addition the genes in each group are functionally related , specifically, genes in g <dig> are mostly linked to translation and might be linked to the translation arrest observed during dormancy.
genes in g <dig> and g <dig> are functionally linked to ‘metabolism’, ‘stress’ and ‘cell wall formation’, which are significantly different in the non replicating state.
zinc uptake regulator zur
as previously stated, the systematic analysis of genes linked to the different cog categories, showed that in the different co-expression networks the rpmb2-rpmg1-rpsn2-rpsr <dig> operon  appeared forming a distinct cluster, separated from the rest of ribosomal protein coding genes .
in addition, this operon also appears linked to the ppe3-rv <dig> operon.
therefore, we concluded that these two operons should be regulated by a specific mechanism and respond to a specific type of perturbations.
the location of these genes in the network of known regulatory interactions  <cit> , pointed to the zinc uptake regulator zur  as the most likely regulator of both operons.
therefore we set forth to study the possible targets of this tf.
initially we selected from the network extracted from the literature, a list of regulatory targets of zur that had been verified either by more than one source or by an appreciable upregulation in a zur knock out mutant  <cit> .
once this core set was selected, we proceed to identify the subset of conditions where these genes are co-expressed.
the method to construct the expression plots  was used to select  <dig> conditions were the members of the core set were correlated ; similarly we selected a core subset set of  <dig> conditions  and another subset of  <dig> conditions .
the correlation among the different putative members of the zurr regulon  was used as a signature to identify the other members of the regulon.
the list of genes identified as belonging to this regulon is provided in additional file  <dig> and the zur consensus binding motif is depicted in figure 8b.
the correlation analysis shows the importance of the biclustering approach to select only those conditions where no additional regulatory influences hinder the discovery of the members of the regulon.
for example, if we analyse the correlation of rpmb <dig> with the members of the core set across all conditions present in our compendium, we obtain no correlation .
however when we compute the correlation between rpmb <dig> and the genes assigned to the core set but only considering the previously selected subsets of conditions the correlation values raise to  <dig> ,  <dig>  and  <dig>  respectively, showing that, as previously reported, expression of rpmb <dig> is indeed regulated by zur.figure  <dig>
zur regulon.
a) bicluster formed by members of the zur as reported in literature  <cit> .
the grey line represents the average expression levels of the members of zur regulon in the conditions in our compendium.
the numbers identify the  <dig> conditions that have been included in the bicluster.
the horizontal bar and the different regions indicated by numbers refer to the classification of the conditions as described in materials and methods.
notice the clear up-regulation of this set in conditions of type 9: transcription inhibition, in particular these values correspond to experiments were rifapentine was added to the medium.
for clarity, expression values have been scaled, so that the mean value for each gene when all conditions are considered is zero.
b) identified zur binding motif.
e-value  <dig> *10− <dig>  c) number of genes in the zur regulon  compared to the ones previously identified  <cit> , the set anticorr.
contains rv <dig> and lpqr, that show anticorrelation with the rest of the genes in the regulon.
our predictions closely match the list proposed by maciag et al.
<cit> , however some differences in the target assignment appear .
from the  <dig> previously assigned genes to the zur regulon, we can only confirm that 79%  of them show correlated expression patterns.
in addition, among our predicted targets there are  <dig> genes for which no experimental evidence can be found in the literature.
these predicted additional members of the regulon are: rv0223c, lpqr , rv <dig>  pe <dig> , ppe <dig> , rv2617c, rv <dig>  rv2619c, rv3018c and rv3018b.
rv <dig> and lpqr appear anti-correlated with the genes in the previously defined core set.
in the case of lpqr this anti-correlation clearly increases when restricting the set of conditions.
this together with the presence of zur binding motif in its upstream region, leads us to conclude that regulation of lpqr  expression proceeds through a different mechanism than that of the rest of the genes in the regulon.
additionally, it is striking to notice the clear up-regulation of this set of genes in conditions where transcription was inhibited by addition of rifapentine to the medium .
this shows that, at least regarding these genes, the inhibition of transcription cause by rifapentine has a similar effect as zinc limitation.
however, a detailed explanation of this phenomenon is most likely unreachable from the analysis of the present dataset.
CONCLUSIONS
we have shown that by integrating data from different sources and through the combined analysis of data, we are able to obtain insights into the biological system under investigation that go beyond the specific research questions of each experimental design.
our systems level approach has allowed us to analyse in depth publicly available data on mtb and has enabled us to extract valuable new information from the existing datasets.
the processes that we have studied  are paramount in allowing mtb to thrive within the hostile host environment.
we generated comprehensive lists of genes involved in the response to such environmental conditions.
these compendia summarize and substantially extend and modify the current knowledge.
we believe that any further research on these adaptation mechanisms will make extensive use of this new knowledge and the hypotheses generated herewith.
we have developed a framework that allows the user to integrate information from gene or protein expression experiments, genome annotations and existing databases together with analysis tools.
most of the already existing databases and tools provide human user interfaces that only allow for querying one gene at a time and only provide a subset of the total available information.
within our framework this integration is done at once on a genome-scale, by using the gene co-expression networks.
instead of obtaining a network with the majority vote from the algorithms we keep the individual output from each algorithm and rely on the expertise of the user to select which of them contains, for each particular problem, the highest amount of relevant information.
this way of presenting the information greatly speeds up the process of analysing new data, since after the study of the different networks and the modules that each of them show, expert users will be able to quickly understand the changes occurring during a new experiment just by overlaying the new expression data over the networks and looking for the networks were the differentially expressed genes or proteins appear clustered together.
our approach is not, in itself a new analytical approach but a general framework that facilitates the consideration of multiple pieces of evidence to pose a hypothesis on a biological system.
comparison with experimentally verified networks  and with in silico generated datasets shows that the co-expression networks generated using our pipeline preserve the desired modularity of transcriptional networks and the regulatory targets of a tf tend to appear in the networks as forming a closely interconnected module.
comparison with e. coli data, shows that in 97% of the cases the target genes of a tf form clusters in the network.
modularity was on average less preserved for s. cerevisiae networks , due most likely to the increased complexity of regulatory interactions in eukaryotes.
this poorer performance of the inference methods for the yeast dataset was also observed in the dream <dig> challenge were at most a  <dig>  recall was obtained by any of the tested methods.
on average the zrc networks perform better that the mixed ones, but our analysis also showed that both types of networks contain complementary information, for example the targets of a yeast's ypr199c are only identified as forming a module in the mixed network .
we have established a protocol to assign regulatory interactions to binding sites identified through chip-seq experiments via the integration of expression data.
we have used this approach to correctly identify the target genes of a given tf  as a response to an specific type of perturbation  among the hundred of candidates from the experimental dataset.
within our framework biclusters can quickly be analysed, interactions between the biclusters can be identified in the overall networks and biclustering methods can easily be transformed into a tool for automatic detection of functionally related modules or underlying layers of regulation.
moreover, we have developed a method that allows to identify additional regulatory layers.
the comparison of the networks generated by considering different subsets of conditions allowed us to distinguish various regulatory mechanisms for dna damage response.
the analysis of zur regulon shows the potential of an integrative approach.
in our compendium there were no data corresponding to experiments designed to analyse the effects of zinc limitation on mtb.
however, the analysis of the expression patterns of the genes in the different conditions and the analysis of their correlations, allowed us to select a set of conditions that would complement the bioinformatic analysis of the upstream sequences of the genes and would allow us to decided which of the regions similar to the motif actually represents a zur binding site and characterize the members of the regulon.
we were able to compare our predictions with dedicated experiments performed with knock-out mutants and we found a good agreement between computational predictions and experimental data.
our work extends the already existing knowledge and produces a comprehensive list of the members of the devr regulon.
among the targets identified by chip-seq, we have uncovered three additional sets of genes that show a consistent expression pattern across the conditions in the compendium and are functionally related.
additional studies are required to further understand the regulation of these genes and their possible link to the non-replicative state.
the discovery of new regulatory mechanisms involved in dormancy has the potential to deliver a new set of drug targets.
the automatic biclustering was the basis of our analysis of dna repair systems and specifically lexa regulon.
we were able to confirm our computational predictions through the comparison we literature data and recently performed chip-seq experiments.
in addition, our work has produced a more specific binding motif for lexa through the identification of new members of its regulon.
additionally, the re-annotation of the identified new targets allowed us to identify a putative mr system in mtb.
the analysis of mutations and mechanisms to avoid them is of the uttermost importance for the further understanding of the evolution of antibiotic resistances and pathogenicity.
previously existing data were used to verify our predictions on the regulatory mechanisms of dna damage response.
we have clarified some points previously in dispute, such as the lack of involvement of sigg in the response to dna damage and the regulation of the alternate sigma factor sigc in these conditions.
correct identification of the sigma factor up-regulated upon dna damage is key to understanding the systemic response of mtb to this damage type.
in addition, our work has provided further evidence on the mechanisms leading to cell cycle arrest upon dna damage in mtb.
in addition, we identified a new regulatory mechanism for zur, since the analysis of the upstream regions of its target genes and their expression patterns show that rv <dig> and lpqr belong to its regulon, although their regulation must proceed through a different mechanism.
here, we have presented the results obtained by applying our integrative approach to mtb and whenever additional data was available, we have found good agreement between predictions and experiments.
the basic underlying principle of this approach relies on the comparison of the networks obtained using different sources of information or methodologies.
this approach can be readily extended to different organisms and to the comparison between different species, by using global identifiers together with a database of orthologous genes between species.
this would allow to select a gene or group of genes in one organism and see how they are arranged in the networks corresponding to a different organisms.
in addition, other types of data, such as synteny or evolutionary information, and protein structure and families could improve the evolutionary comparison of functional modules.
additional files
additional file 1:
microarray data sets considered in this meta-analysis.
identifiers from the gene expression omnibus database have been used for each sample, experiment and platform.
additional file 2:
identified biclusters.
the following information is presented for each bicluster.
left.
top
: expression plot, the grey line represents the average expression levels of the genes in all the conditions in this meta-study, the numbers identify the conditions in the biclusters.
the horizontal coloured line shows the classification of the conditions .
middle interactions among the members of the biclusters either from protein-protein interactions ; operons or similarity of annotated cog terms.
the nodes in the network represent the genes in the biclusters.
the colours are linked to the cog annotation, whereas the shape is linked to the functional classification: square: information storage and processing; rhomboid: cellular processes and signalling; triangle: metabolism; and circle: poorly characterized.
bottom: identified motifs.
right: genes in the bicluster, their upstream sequences have been plotted together with a representation of the location of the detected motifs: red, green and blue for the first, second and third motif .
additional file 3:
pipeline for the generation of co-expression networks.
the r code that implements the pipeline as depicted in figure  <dig>  the inequality simplifier is implemented in the joined jar file.
additional file 4:
technical characteristics of the visualization tool.
technical description of the network visualization and comparison tool used.
additional file 5:
relationship between the expression levels of
reca
and
lexa.
the fit corresponds to a model, based on ordinary differential equations, of lexa mediated induction of reca.
the numbers signal the conditions where this model actually holds.
for clarity, expression values have been scaled, so that the mean value for each gene, when all conditions are considered, is zero.
the horizontal bar includes the overall categorization of the conditions .
notice that the conditions were lexa is regulating reca expression are mainly linked to 1: aromatic amides intra cellularly hydrolysed, low ph; 4: acidified medium; 5: cell wall synthesis inhibition; and 8: dna damage.
additional file 6:
comparison of predictions for individual methods to consensus network.
a) overall predictions.
the plot shows the number of interactions from the gold standard that were predicted worse, equally or better by each individual method than by the consensus network  <cit> .
the predictions were assessed comparing the rank of the predictions with a  <dig>  cutoff.
only the e. coli network has been considered.
b) predictions for each tf.
fraction of the interactions that are better predicted by each type of method than by the consensus network.
only tf with more than  <dig> interactions have been considered.
these plots have been built using the methods described in  <cit> .
additional file 7:
topological overlap of tf targets in the co-expression networks obtain for the dream <dig> challenge datasets.
for each dataset  the zrc and mixed networks were built.
only tf with more than  <dig> experimentally verified targets  were considered.
dashed line represents the average topological overlap in each network.
additional file 8:
relative locations of the cluster of ribosomal proteins and the operon of ribosomal proteins
rpmb2-rpmg1-rpsn2-rpsr2
in different co-expression networks.
additional file 9: table s <dig>  lexa regulon.
additional file 10:
genes in the reca_nd regulated dna damage repair set.
operons and divergons are marked by alternating  row colours.
the fifth column indicates whether the targets had been previously identified  <cit> .
additional file 11:
identified motif in the upstream regions of the genes in the reca_nd dna repair system.
logo generated using meme  <cit> .
e-value  <dig>  e- <dig>  the location of the hits is included in additional file  <dig>
additional file 12:
classification of devr targets identified by chip-seq.
the column tbdb indicates whether the target was listed in the original data set as a possible target for devr regulation, some genes have been added to the table based on operon/divergon analysis.
hit location indicates the genomic location of the hit.
the last column  refers to the initial classification of the genes within the core set.
the column references refers to whether they have been previously identified as members of devr regulon: a identified by park et al.
<cit> ; b identified by balázsi et al.
<cit> ; c identified by bartek et al.
<cit> ; d identified by chao and rubin  <cit> .
<dig> groups are present in the table.
"the group  <dig> contains those genes whose expression patterns showed no significant correlation with any others, and we believe are not regulated by devr, however they have been previously identified as belonging to devr regulon in literature, therefore we have included them in the table.
the first group corresponds to the members of devr regulon.
additional file 13: table s <dig>  devr regulon members.
additional file 14:
bicluster formed by members of devr regulon and the type vii secretion system  genes
rv282-rv <dig>  the grey line represents the average expression levels of these genes in the conditions in our compendium.
the numbers identify the  <dig> conditions that have been included in the bicluster .
additional file 15:
devr motifs top: previously identified primary and secondary motifs for devr
[
84
] bottom: newly proposed motifs.
main differences in primary motif can be seen at positions  <dig>   <dig>   <dig> and  <dig>  main differences in secondary motif can be seen at position  <dig> and  <dig> to  <dig>
additional file 16:
results of the go enrichment analysis performed in the classified devr targets identified by chip-seq.
node size reefers to the number of genes annotated to that particular go term in the whole genome annotation file, whereas sample match indicates the number of genes annotated to the go term in the group .
additional file 17:
members of zurr regulon.
operons and divergons are marked by alternating  row colours.
rows signalled in yellow refer to genes that have been previously identified as belonging to this regulon, but for which our analysis does not support this claim.
the fifth column indicates whether the targets had been previously identified by maciag et al.
and the sixth column reflects the reported fold change in a zurr knock out mutant  <cit> .
column  <dig> contains the definition of the core set of genes, whereas columns  <dig> to  <dig> contains the average correlation between the expression level of the each gene and the genes in the core-set across all conditions in the compendium or a subset of them.
abbreviations
tbdbtb database
tftranscription factor
mimutual information
mtbmycobacterium tuberculosis
mrmismatch repair
gogene ontology
dpidata processing inequality
ssdnasingle stranded dna
reca_ndreca non-dependent
aurocarea under receiver operating characteristic curve
auprarea under precision recall curve
competing interests
the authors declared that they have no competing interests.
authors’ contributions
jd and msd: study concept and design, methodology and tools development, manuscript preparation and data analysis.
ps and vmds: study concept and design and critical revision of the manuscript.
all authors read and approved the final manuscript.

 current biology research generates an ever-increasing deluge of omics derived data. each type of omics data pertains to a single level of the biological system under investigation. while detailed knowledge of the individual genes, transcripts, proteins, metabolites and other cellular components remains important, understanding a biological system requires considering the networks into which these components are embedded and of their functioning as a whole. a major challenge in systems biology lies thus on developing effective and efficient methods to optimally extract the information contained in the aggregate of these datasets. the increasing availability of genome-scale expression data has boosted the development of methods to infer the underlying regulatory networks. a broad range of alternative methods are available, see for reviews, and each of them uses different mathematical tools and/or biological assumptions. a class of methods use differential equations to express transcript changes as a function of the transcript levels of the corresponding transcription factors. a second class of methods rely on bayesian networks to analyse the joint probability distributions obtained from the experimental data. other methods use the similarity between gene expression profiles to detect associations and to reconstruct a genome scale transcriptional regulatory network. another class of methods use a combination of machine learning techniques to produce prioritized lists of transcription factors regulating each target gene. each method has different strengths and weaknesses, even methods using similar conceptual tools. for example c3net uses mutual information to reconstruct the core regulatory interactions in the network, that is, to recover the strongest interactions. within this core, c3net was shown to outperform other methods also using mi such as clr and arance. each method provides different results, therefore much effort has been devoted to generate consensus networks from the multiple solutions. it has been shown that in many instances, an integrative approach combining the outcome of each algorithm produces the best result, however, a detailed analysis shows that for some interactions individual methods perform better than the consensus network. co-expression networks and module identification algorithms co-expression networks contain genes as nodes and the edges of the network represent significant co-expression levels across the studied data set. an open problem is still how these connexions are to be defined and how an adequate threshold is to be imposed. in differential network analysis two networks obtained using the same algorithm but alternative datasets are compared to identify interactions appearing only under a subset of conditions. for example, such an approach was used to analyse prostate cancer datasets by comparing the networks inferred using c3net from normal or tumour samples datasets and was able to successfully identify cancer specific interactions. here, we construct and analyse co-expression networks extracted from the same dataset but using alternative algorithms to identify relevant interactions. biclustering or module learning algorithms aim at the identification of functionally related genes showing co-expression patterns. in some cases, such as the cmonkey algorithm, additional biological information is also considered. the results of the genome-scale reconstruction methods can be displayed and visualized as one network whereas an identified module or bicluster contains a set of functionally related genes that might be differently regulated in different conditions. therefore, the genes within a module might not form a cluster when considering all the conditions and would remain undetectable for network inference methods. network inference methods and module prediction algorithms are highly complementary and new information can be obtained by combining the outcome of both approaches. analysis of chip-seq data there are experimental methods to directly reconstruct regulatory networks. the knowledge of transcriptional regulatory events and specifically the transcription factors binding sites can be greatly improved by chromatin immunoprecipitation and sequencing, chip-seq. powerful techniques have been developed to process the sequencing data and to isolate the binding sites from the background noise generated by non-specific sequences, such as statistical tests, use of controls, techniques for signal de-convolution or lag-analysis among others. in the classical model of transcription regulation in prokaryotes a binding site in the promoter region of a gene is linked to a regulatory interaction between the corresponding tf and the gene. however, the collection of chip-seq data for mycobacterium tuberculosis hosted in the tb database shows that not all identified tf binding sites can be linked to regulatory interactions. for example, in some instances the tf binding is much weaker than expected, in other cases this association is complicated by the occurrence of divergons, which are pairs of divergently transcribed operons or genes, or by the presence of binding sites within coding regions. this lack of a one-to-one relationship between tf binding sites and regulatory events can be due to multiple reasons such as non-specific binding, cumulative effects of sites with weaker binding to regulate overall promoter affinity, specific binding generated in non-biological or non natural conditions, false positives as a result of the experimental procedure, or, simply errors in the genome annotation. we propose that it is through the integration with additional data sources, specifically through the integration with expression data, that this challenge can be overcome. integration of heterogeneous molecular networks each inference algorithm has its weaknesses and strengths and each network holds its own intrinsic value and can be used to gain more specific insights on various aspects of the biological system. in a sense, using multiple algorithms to extract the networks and presenting them to the user is a similar approach to that followed by annotation pipelines, such as microscope, diya or basys among others. these pipelines present the user a list of gene centred information, with different annotation sub-fields, to enable supervised annotation. integration efforts require a common layout for the data, therefore, we have chosen to represent the available layers of information, such as operon structure, known interactions between genes or proteins, enzymatic activity or functional similarity, among others through network representations. the networks themselves are represented in a common format, xgmml, that allows their simultaneous exploration. in addition to the simultaneous visualization of different networks, additional analysis tools such as motif search and identification, go-enrichment analysis, and tools to overlay expression data or analyse expression profiles of multiple genes across different conditions can be used. we have developed a pipeline for the generation of co-expression networks that is easily tunable to produce alternative networks. the pipeline is based on using the similarity between gene expression profiles to detect associations and contains a higher order extension of the data processing inequality method to reduce the number of redundant or possible spurious links. this pipeline is presented in figure  figure schematic of the pipeline to obtain co-expression networks. from top to bottom the following steps are applied: calculation of the similarity matrix, z-transformation, combine i, threshold setting inequality simplifier and combine ii. note that when applying the inequality simplifier to the zrc network the result will be a tree. through the exploration of the multiple networks, the user gains new knowledge of the biological system, but it can also lead to the discovery of new strategies to integrate the information stored in the networks. these newly gained strategies can then be translated into new pipelines. as result of our exploration, we have developed a pipeline to uncover multiple layers of regulation, presented in figure  and another one to analyse chip-seq data and assign regulatory interactions to the detected binding sites, presented in figure  figure pipeline to uncover additional regulatory layers. step 1: identify conditions linked to the main regulatory event for the initial gene set. this can be done using biclustering techniques or by direct comparison with the expression levels of the regulator. step 2: build co-expression networks in the remaining conditions. step 3: identify the closest neighbours of the selected genes in the new networks. step 4: iterative round of motif identification/matching to identify the secondary motif and the set of genes with this motif in their upstream regions. pipeline to analyze chip-seq data. after the locations of the chip-seq binding sites have been retrieved, their genomic context is analysed. a core set is defined by selecting targets with i) literature evidence or ii) a hit in the upstream region of not divergently transcribed genes. the expression levels of these genes are analysed and they are categorized through clustering. finally the rest of the putative targets are assigned to these groups based on the similarities of their expression patterns. show case: deployment of the framework to unravel regulatory mechanisms in mycobacterium tuberculosis we have analysed regulatory events in mtb. due to its implications to human health this highly successful pathogen has been extensively studied and there is already a substantial body of information on mtb and its underlying regulatory networks, but still much remains to be learned. we have analysed not only networks extracted from literature but also networks extracted from publicly available databases such as string, metacyc, kegg, tbdb and tuberculist. expression data from publicly available repositories and corresponding to perturbations, have been used to generate the co-expression networks and to further analyse chip-seq data. we have finally combined the inferred networks with the outcome of biclustering algorithms, to further explain the functionalities of the modules. we have analysed the regulation of dna repair systems within mtb, particularly we have focussed on two alternative regulatory mechanisms: the ‘reca-lexa dependent dna repair system’ and ‘reca non dependent regulation’. we will show that the integration of expression data and its use to guide the exploration of upstream sequences allows to analyse alternative regulatory mechanisms for the same set of genes. in addition, we have explored some of the regulatory mechanisms that allow mtb to survive within the host, such as devr regulon, which is a key element for understanding the dormant state, and the regulation of the response to changes in zinc availability. methods gene expression data two-colour microarrays for mtb were retrieved from the gene expression omnibus database. of them aimed to capture the effect of drugs targeting metabolic pathways whereas captured stress induced dormancy in the wild-type and in devr activation genes knockout mutants. we followed boshoff et al. to classify the conditions in the compendium into categories according to the experimental perturbation: aromatic amides intracellularly hydrolyzed, low ph; mutants of devr activation pathways; translation inhibition; acidified medium; cell wall synthesis inhibition; respiration inhibition; nutrient starvation; dna damage; transcription inhibition; iron scavengers; multiple stress sources applied simultaneously; minimal medium; not classified; and conditions associated with devr upregulation. a list of the used datasets is presented in additional file we applied a common normalization method, loess, to the arrays from each experiment and each of the six independently designed platforms. linear models were constructed to consider biological and/or technical replicates; within-array replicate probes were averaged; and an additional between-array quantile-quantile normalization was performed to ensure the comparability between experiments. these manipulations were performed using the r limma package. a common locus tag format was introduced and missing values were filled up using knn-imputation from the impute r package with k =  and eliminating genes and experiments with more than 30% and 50% missing values respectively. the resulting compendium contained information on open reading frames across different conditions. biclustering we adapted the original cmonkey r code to mtb data and we increased the number of biological information sources that the algorithm can consider. for the biclustering process, we selected a subset of the interactions present at the string database. we selected specifically interactions obtained based on co-occurrence of linked proteins across different species, on curated databases, and on their association in the abstracts of scientific literature. from prolinks we collected interactions obtained by the phylogenetic profile method. in addition, we included the network based on the similarity among the annotated go terms, which was constructed using the sleipnir library. upstream sequences for mtb genes were retrieved using rsat and used for the motif detection and motif matching steps. for the automatic biclustering algorithm, we used the 1000 bp region upstream of the translation start site, avoiding overlapping with upstream neighbour genes when present. in the initial rounds we used a reduced matrix where only the leaders of the operons, as defined by roback and co-workers, were kept. multiple runs of the algorithm were performed, which used the default parameters except for: initial size of seed clusters; number of iterations; maximal number of clusters. we obtained a set of biclusters. the jaccard similarity coefficient between each pair of biclusters was computed. to reduce the redundancy in the set of biclusters, we merged the pairs that showed a jaccard similarity higher that. the biclusters in the merged set were enlarged by adding genes based on the predicted operon structure combined with the expression level measurements. this new set was used to seed the biclustering algorithm in two subsequent optimization rounds one biased towards the detection of concurrent motifs in the upstream regions and a second one biassed towards the identification of sets linked to highly related biological processes. we performed an additional manual merging step, that also considered the similarity between the detected motifs using the matrix comparison tools from rsat and we obtained a final set of biclusters. expression plots for these biclusters are in additional file pipeline for the generation of co-expression networks the pipeline to generate co-expression networks is presented in figure  departing from an nxm matrix containing the expression profile of n genes in m conditions, we construct a symmetric nxn similarity matrix s. the pipeline allows to choose between correlation and mi as similarity measurements. mi is computed using an estimator based on the entropy of the empirical probability distribution with initial data discretization into n equal sized bins from the bioconductor package minet. in the second step, a z-score transformation is performed on the distribution of the scores from the similarity matrix. the transformation is done by rows, to obtain the zr matrix; by columns, to obtain the zc matrix; or by both rows and columns simultaneously, to obtain the zrc matrix. the z-score transformation allows for each possible interaction to be weighted regarding the background of interactions in which each member of the interacting pair is involved. therefore, in the zc case each element sij becomes zcij, which is the z-score of sij regarding the distribution si si …, sin. in the zr case each element sij becomes zrij, which is the z-score of sij regarding the distribution s1j, s2j, …, snj. in the zrc case each element sij becomes zcij2+zrij2/ in the following step the matrices are combined into two new matrices: cλ and rλ: cijλ=maxzcij,λ·zrcijrijλ=maxzrij,λ·zrcij with λ a positive real number and i and j denoting the rows and columns of each matrix. λ allows to fine-tune the results, since higher values of λ will lead to cλ and rλ matrices that are similar to each other and to zrc. considering λ = √ means that, for those cases where one of the elements in eq. is zero, both λ zrcij and zcij will be identical. however, using the default λ =  will ensure that when the values of zcij and zrij differ then the highest one will be selected either through zcij or zrcij. in the following step a threshold can be applied to remove interactions with low weight. any value below the threshold is set to this step produces a more sparse network, which might be needed to obtain a neat visualization. the default threshold value is chosen so that the number of non zero edges is equal to a predefined value. afterwards the inequality simplifier can be used to remove, possibly spurious, links from the network. the inequality simplifier is an extended version of the data processing inequality. in simple terms, the dpi states that given two interdependent random variables a and a and a third one a that only depends on one of them, for example a then a cannot contain more information about a than a does. this statement is mathematically represented by the following inequality among mi values mia a3≤minmia a mia a the dpi must hold whenever a does not depend on a so it can be used to remove spurious interactions from the network, possibly caused by feed forward loops or mutually dependent regulators. the spurious links are removed according to: ifmia a3≤mia a2&mia a3≤mia a thensetmia a3= the dpi can be extended to higher-order interactions, therefore allowing a recursive implementation. the dpi is derived from the triangle inequality satisfied by any metric or distance, such as the mi or kendall rank distance. our extension to remove spurious edges by considering higher order interactions contains two parts: the first part is the identification of alternative pathways connecting the same nodes and the second part is to apply the inequality simplifier to decide whether a link is spurious and should be removed. the first part is based on dijkstra's shortest route algorithm. here, instead of searching for the shortest route between two nodes, we use it to find the alternative path between those two nodes with the best throughput. the throughput of a path is defined as the value of the edge with the lowest throughput value. given two connected nodes, the link between them is considered spurious and removed if an alternative path is found with a higher throughput value. for a given pair of nodes, ai and aj,, so that the best alternative path of length passes through a a …, an- and an, the following rule is applied: ifmiai,aj≤minmiai,a mia a …,mian− an,mian,ajthensetmiai,aj=0wheren∈ …,n− this rule is applied for each possible pair of nodes and for each possible value of n. the alternative pathway joins ai and aj and does not contain any cross-linking, so that therms such as a2a or a1a do not need to be considered. this is a consequence of using the path with the best throughput so that no triangle inequality of the metric is required to derive this property. figure  shows how the inequality simplifier acts to remove, possibly spurious, links from the network when applied to higher order terms. while applying this rule, an additional nxn matrix, t, is built to keep track of how many links a particular edge has caused to be removed. in this case, the links aia a1a a2a …, an-1an and anak have caused the removal of link aiaj, therefore, in the t matrix, the elements tai,a ta a ta a …, tan- an and tan,aj are increased by 1∕.figure the inequality simplifier. the similarity values among the nodes connected by the different edges have been indicated. dotted lines represent the spurious links removed by the inequality simplifier. left application of order two, which is equivalent to a direct application of the dpi. right higher order application, both dotted lines are removed. the dpi would only remove the blue dotted line. finally the two networks cij and rij are combined into one final network rcij. rcijλ=maxzcij,zr the presented pipeline can be used to obtain different outputs, such as the rc or the zrc networks and their variants through the application of thresholds or the inequality simplifier. the final output is a square matrix that represents a network through a weighted adjacency matrix. this adjacency matrix is finally exported into a tabular format that can be imported into cytoscape to obtain a graphical representation. the weights in the adjacency matrix or the t matrix produced by the inequality simplifier can be imported as edge attributes and used to set the thickness of the edges. finally, a force directed algorithm can be used to generate the appropriate layout and the graphical representation of the network, which then can be exported to an xgmml file. applying the inequality simplifier to a symmetric matrix, such as zrc, removes any possible loop in the network, therefore it results in a forest graph as depicted in figure  applying the inequality simplifier to an initial matrix containing around nodes takes at most five minutes in a standard desktop computer. the networks used for our exploration of gene regulation networks in mtb were the result of applying the inequality simplifier to zr, zc and zrc. belcastro et al. have shown that from the total number of possible interactions a selection of 5% would be a sensible choice according to the assumption that biological network are sparse. here, one of our goals is to obtain a clear visualization of the network, therefore, we have chosen an even lower number so that the number of nodes included in the network is as high as possible without having a too crowded visual representation of the network. an additional threshold was imposed so the number on non zero edges of zr, zc and zrc networks is equal to the zrc derived network contained nodes and edges whereas the mixed network contained nodes and edges. these manipulations were implemented in r and java, the corresponding scripts can be found, together with an example data set in additional file additional networks and analysis tools metabolic networks the nodes in the metabolic networks are of two different types: reactions and metabolites, and additional information, such as the stoichiometry, directionality and the catalysing proteins of the reactions, also has to be stored. we kept this information by using a combination of two file types: a svg file that contains the graphical representation of the networks and an rdf file containing the additional information. other networks a network formed by multiple disjoint interlinked clusters was generated using the predicted mtb operon structure. in addition, the information stored in the string database was used to generate networks based on: known protein-protein interactions; distance of the genes in the genome and that of orthologous genes in other species; co-occurrence in other species; orthologous genes being fused together in other species; association in the abstracts of scientific literature; and other databases, such as metacyc and kegg. to identify regions duplicated in the genome or homologous genes, a network was built by considering the sequence similarity between each pair of genes. the similarity was computed using megablast and an edge was created for matches with e-values lower than of 10− sequence analysis tools we have used meme for motif elicitation within a set of genes and fimo to find occurrences of the chosen motif. a post processing step ranks, using the fimo q-value, the genes for each identified motif. an iterative cycle of motif identification and motif matching can be established to further refine each motif. the computational time required by each iteration greatly depends on the length of the upstream sequences, which has to be selected for each particular analysis. for the selected show cases the motifs accumulated in the 200 bp region upstream of the translation start site. gene ontology enrichment analysis on selected set of genes was performed using a hypergeometric function to model the probability density, as implemented in the gohypergall r function. the go annotation for mtb was obtained from the uniprot-goa database, the mtb-goa database, tbdb and the recent re-annotation of the mtb genome. operons we developed a putative operon extension tool which is only based on the gene orientation, so that genes that are downstream of the selected gene and with the same orientation are added to the putative operon until a gene with the reverse orientation is found. the user is then free to decided whether the expression data sustains the extension of the operon. this approach allows to easily combine expression information with genomic context information. venn diagrams have been created using the utilities from the venndiagram r package. expression plots generated with a solid line representing the mean expression of the selected genes in the conditions included in the compendium and dots marking the conditions in the bicluster where the genes show high correlation. the classification of the experimental conditions, previously mentioned, was included in the expression plots through an horizontal colour line that allows to quickly associate the behaviour of the genes with the type of perturbation. to select the conditions in which the selected genes show co-expression there is an initial step where all conditions are included in the bicluster, then an iterative loop starts were each condition is removed at a time and the new correlation values are computed, finally the condition leading upon removal to the highest correlation values in the remaining set, is removed from the set. this process is iterated until the correlation in the remaining set of conditions is higher or equal to. discovery of additional regulatory layers whenever two alternative regulators regulate the expression of genes in a given set, it might be that one overshadows the detection of the other. we have developed the work flow to analyse these events. given a set of genes under control of two regulatory events, the first step is to identify the conditions were the first event is active or taking place. how these conditions are identified depends on the particular example. one way would be to identify conditions with up/down regulation of the corresponding transcription factor. to identify conditions were reca expression is not regulated by the reca/lexa mediated mechanism we built a model linking the time dynamics of reca to lexa expression levels using the inferrelator algorithm. the output of the algorithm is a prediction or fit of reca levels based on lexa levels. those conditions with a poor agreement between the measured and fitted levels are the conditions without lexa mediated induction of reca. once the new set of conditions has been selected, the following step is to reconstruct the co-expression network, using only the expression data corresponding to this subset of conditions and analyse the location of the original set of genes in the newly built network. furthermore an iterative round of motif identification and motif matching can be run to identify a possible motif in the upstream region of the cluster. this iterative round can be performed either using cmonkey or manually using alternatively meme and fimo. analysis of chip-seq data publicly available chip-seq data were obtained from tbdb. we have developed the workflow, presented in figure  to analyse these data. the chip-seq data had already been analysed with a peak calling algorithm. for each of the considered tfs a list of predicted targets was obtained from tbdb, together with information about the peak location. a reduced list, or core set, was obtained by selecting those targets that show a hit in their adjacent upstream intergenic regions. an additional filtering was done to select only those hits where the peak is not flanked by divergently transcribed genes/operons. additionally, the genes for which literature evidence supported the regulatory interaction with the tf were included in the core set. in the following step, the expression levels of the genes in the core set were analysed and the matrix of correlations across the conditions in the expression compendium was computed. the appearance of negatively correlated genes in this set is a signal of a dual repressing/activation function of the regulator and therefore two subsets can be defined by hierarchical clustering of the set. this process can be further enlarged to encompass additional groups that would be linked to alternative regulatory mechanisms or to the effect of additional tfs. once this/these group were defined, the rest of the putative targets predicted by chip-seq were assigned to either one of these groups based on their average correlation with the members of the group. topological overlap the topological overlap between two genes i, j in an unweighed network is defined as: ti,j=ni∩njminni,nj+1−ai,ji≠j1i=j where n is the set of direct neighbours of gene x; |x| represents the number of elements in set x and ai,j is the adjacency matrix of the network. the topological overlap is bounded between and two genes will have high topological overlap if there exists a connection between them and if they are connected to the same group of genes. for weighted networks, ai,j represents the weight of the interaction between genes i and j and takes continuous values between and in these cases, the previous formula can be generalized to: ti,j=Σkai,kak,j+ai,jminΣkai,k,Σkaj,k+1−ai,j the topological overlap of a group of genes was defined as the average of their mutual topological overlap, which was computed using the r wgcna package. only groups with more than genes were considered. an empirical p-value for these scores was calculated by randomly sampling gene groups of the same size in the respective network. network visualization we have developed a visualization tool that allows the simultaneous visualization of networks in xgmml format and the sharing of identifiers between them. technical characteristics of this tool are available in additional file  co-expression networks we have developed a pipeline to generate co-expression networks that allows for a myriad of possible networks. choosing a subset of them highly depends on the available data and the process the user wants to explore, since different clusters appear in each of them. this reflects the inherent modularity of biological networks and the dependency of regulation on the chosen conditions. the pipeline presented in figure  includes some already well established and tested methods for genome-scale network inference. for example, the similarity matrix computed using either correlation or mi is commonly used to analyse gene expression data, but also other methods are contained in the pipeline. when working with mutual inference as a similarity measure the matrix denoted zrc corresponds to the output of the clr algorithm. additionally, applying the inequality simplifier only up to order two on a mi matrix amounts to using the dpi to prune spurious interactions from the networks, which is a key element of the aracne algorithm. a comparative analysis of methods to reverse-engineer transcriptional regulatory networks was done using the results of the dream challenge. in this challenge the teams had to reconstruct genome-scale transcriptional regulatory networks from expression data. the networks proposed by the different teams, were evaluated through their comparison with a gold standard, a set of experimentally verified regulatory interactions in the target organisms. in addition, marbach et al. constructed a consensus network integrating the predictions from the multiple methods. whole network performance estimators such as area under receiver operating characteristic curve or area under precision recall curve show that over the entire network the consensus network outperforms individual methods due to their inherent complementarity. however, for some interactions individual methods perform better than the consensus network. the loss of information when building the consensus network does not affect equally the different transcription factors. for example, for purr or lexa, a significant fraction of the total number of known interaction is better recovered by the individual methods than by the consensus network. to evaluate the performance of our algorithm we have generated co-expression networks corresponding to the synthetic, escherichia coli and saccharomyces cerevisiae datasets used in the dream challeng. for each dataset, two networks were built. in the challenge the goal was to identify links between regulators and target genes. in co-expression networks no special emphasis is done on the regulators. to evaluate the networks we have calculated the topological overlap of the known targets of the transcriptions factors. the topological overlap measures how similar the neighbourhoods of two genes are. for 95% of e.coli tfs considered the topological overlap of the target genes in the zrc network is significantly higher than for the overall network, which means that they form cohesive modules in the network. the relative number of cohesive clusters identified for the other datasets are lower overall the zrc networks maintain a higher degree of cohesiveness than the mixed networks.figure topological overlap of tf targets in the e. coli zrc co-expression network. the zrc co-expression network was reconstructed using our pipeline using e. coli expression data from the dream challenge. only the tf with more than experimentally verified targets were considered. dashed line represents the average topological overlap in this network. systematic analysis of regulation in mtb. integration and simultaneous exploration of heterogeneous networks we considered co-expression networks obtained with different approaches. additional information is represented through networks, such as information on operons, go annotation similarity, data base stored knowledge and sequence homology among others. this integrative approach allows to systematically explore functional modules in the network and it is highly complementary to existing bi-clustering and other module identification methods. even when additional biological information is included, biclustering methods have to face the challenge of interpreting the function of each groups of genes. therefore, a general approach to detect and understand functional modules in a given organism is to simultaneously explore the location of the genes in a given bicluster in the co-expression networks. we have done so a to analyze the output of the cmonkey algorithm on mtb data. we have run multiple instances of the algorithm to bias the search towards the detection of both putative co-regulated and functionally related sets of genes. the obtained bi-clusters are presented in additional file to further investigate these biclusters, we have projected them in the multiple networks and genes proceed to their in-depth analysis. an example of how this is done is shown in the following section when we analyse the regulation of dna damage repair systems in mtb by lexa. in a given co-expression networks it might very well happen that the effect of one regulator on a given set of genes overshadows the identification of other regulators. however, these mechanisms can be uncovered when comparing different networks generated under different experimental conditions. we have developed the pipeline presented in figure  to uncover additional layers of regulation based on the assumption that regulation proceeds differently upon different perturbations, therefore some regulatory events will only be detectable in subset of conditions. this pipeline is only useful when the user suspects or has prior knowledge of an additional regulatory interaction affecting the same set of target genes. therefore, we have used it to further explore the regulation of dna repair systems in mtb, where additional regulators are known to exist. in addition, our approach allows to systematically explore and interpret additional data such as chip-seq data, through the pipeline presented in figure  among the available data we have chosen to focus on one of the key subsystems key for mtb survival in the host, the regulation of the response to hypoxia and the induction of dormancy program via devr. furthermore, the simultaneous network visualization, allows the exploration of common trends in the networks. instead of focussing on clusters of genes appearing in the co-expression networks and how these clusters are reflected in, for example, the network liked to the go biological process annotation. the alternative approach can be taken and for each set of genes linked to a particular go or cog category we can trace back their location in the network. this analysis can point to interesting effects. for example, we explored the genes coding for ribosomal proteins and in the different co-expression networks these genes appear to form a highly interlinked cluster. however, in all versions of the networks, the rpmb2-rpmg1-rpsn2-rpsr operon formed by genes coding for ribosomal proteins appear as a separated set. the combination of this information with the network of interactions extracted from literature led us to the analysis of the zinc uptake regulator, zur, and its targets. dna repair systems in mtb reca-lexa dependent dna repair system: the sos box lexa is a repressor known to be involved in the control of mechanisms for dna repair. under neutral conditions lexa, dimerizes and binds to the sos box, repressing its own expression and of other genes related to dna repair. the consensus sequence of the sos box for mtb had been identified as the palindromic motif tcgaac4gttcga. upon dna damage, reca binds single stranded dna. the complex reca-ssdna, stimulates autocatalytic cleavage of lexa, so genes repressed by lexa are induced. acidic conditions trigger a similar response, since lexa can no longer dimerize, effectively preventing it from binding its target dna sequence. the iteratively mapping of the genes in each of the biclusters presented in additional file to the multiple co-expression networks, showed that bicluster is a gene module that shares many of the characteristics of the lexa regulon. to further investigate this bicluster, we proceed to analyse the location of these genes in the multiple networks and to identify genes that might be linked to them through the systematic analysis of their expression patterns. for each of the candidates to be included in the regulon we combined information from the different information layers: literature, databases, functional annotation, and upstream sequence among others. we were able to identify a total of genes putatively in the regulon, listed in additional file 9: table s this includes the genes reported by davis et al. plus additional ones, resulting in an effective enlargement of the regulon size by 75%,. based on this extended regulon, we identified for the sos box in mtb a more specific motif: mkwmtcgaamrywtgttcga plot of the average expression level of the members of the lexa regulon across the different conditions. red dots mark conditions with high correlation between the genes in lexa regulon. the horizontal bar and its different regions indicated by numbers refer to the classification of the conditions as described in materials and methods. high expression levels are observed in conditions corresponding to low ph or uv light. b) clusters of genes involved in dna repair mechanisms in the co-expression network. genes regulated by lexa are marked red. c) refined lexa identified binding motif, positions and were previously non specific. d) number of genes identified to be regulated by lexa. previous indicates genes previously reported in the literature as lexa regulated, whereas automatic refers to the genes initially identified by the automatic biclustering algorithm. to verify our predictions about the genes not previously assigned to this regulon, we compared our predictions with the lexa binding regions identified by smollett et al.. in that work chip-seq analyses of lexa binding sites was complemented with experimental measurements of gene expression upon dna damage induced by mytomycin c. the genes that we assign to the lexa regulon contain a site in their upstream regions where lexa binding was detected. in addition, these genes also show dis-regulation in dna damaging conditions. it is important to stress that the results from smollett et al. were used in our analysis only to verify our predictions. therefore, we can conclude that we have successfully reconstructed the regulatory interactions of lexa. davis et al. identified a single sos box in the upstream region of the divergently transcribed genes whib and fbia, however none of them were listed as likely to be regulated by lexa. the analysis of the expression profiles of these two genes, fbia and whib led us to conclude that fbia, a probable 2-phospho-l-lactate transferase involved in coenzyme f biosynthesis, does not belong to the regulon, since no significant expression changes are observed upon conditions related with induction of the lexa regulon. the expression pattern of whib gene shows anti-correlation with the rest of the genes in the lexa regulon in all the conditions in our compendium that show upregulation of lexa. therefore, we concluded that lexa regulation of whib expression proceeds through a different mechanism than the previously described so that the dimerized form of lexa acts as an inductor of whib expression. this has also been independently confirmed in experimental conditions were mytomycin c was added to the medium. whib has been hypothesised to be involved in the regulation of cell division. functional analysis of the genes co-expressed with whib further supports the idea that whib is involved in the regulation of cell division. we also found that not only rv2719c, as reported by davis and co-workers, but the complete rv2719c-nrdr-rv2717c operon is under the control of lexa. information from the string database allows us to functionally link rv2717c with dna repair and cell wall synthesis; nrdr is involved in the control of the synthesis of dntp needed for dna replication and/or dna repair; and the cell wall hydrolase rv2719c is involved in suppressing the cell cycle by altering localization of ftsz rings. therefore, our analysis has unveiled that repression of whib and induction of the rv2719c-nrdr-rv2717c operon are the lexa regulated mechanisms to temporary arrest cell division upon dna damage. furthermore, through the re-annotation of the newly discovered lexa targets, we can identify a putative functional mismatch repair system. no mr system has been previously identified in mtb. a typical bacterial mr system contains: the mismatch-recognition protein muts that contains an hnh endonuclease domain; muth, a nicking endonuclease; and mutl, which acts as a scaffold between these. the mr system additionally contains a dna helicase; a dna exonuclease to remove the mismatching nucleotides; and a dna polymerase together with a dna ligase to repair and ligate the created gap. the hypothetical mr system that we have identified is formed by: i) some of the hnh endonuclease domain containing genes that belong to the 13e family, ii) the holiday junction dna helicases ruva and ruvb together with ruvc, a crossover junction endo deoxyribonuclease, and iii) the dna polymerase dnae together with imua/b that are essential for its function. the existence of this repair system, under the control of lexa, would solve the apparent inconsistency between the low mutation rates in mtb and the absence of an mr system. additional regulation: reca independent dna repair system reca is key to the correct regulation of the lexa regulon and it is also lexa regulated. however, additional dna repair mechanims have been described in mtb, particularly the reca non-dependent dna repair system, that also regulates reca expression in a lexa independent manner. to analyse this additional regulatory layer, we have used the pipeline in figure  the conditions not linked to the main regulatory event are those where no relationship was found between reca and lexa levels. using these conditions we built a new co-expression network and afterwards we compared the original and the newly built networks. the members of the lexa regulon appear in the original network as a tight cluster, however, in the second network only a subset of them appear clustered. we selected the genes in this small cluster and proceed to an iterative round of motif identification and motif matching, to finally identify the genes regulated by the reca_nd mechanism. our approach does not allow us to identify the regulator of the set, but previous studies point to clpr. we have verified our predictions for the identified motif and the list of targets genes, through comparison with literature data, since they match those previously described by gamulin et al. for the reca_nd dna repair system. however, there are two striking differences between our results and those previously reported by gamulin et al. we don’t find that sigg responds to dna repair, which matches the result reported by smollett et al., on the other hand, we do find another sigma factor, sigc, that seems to be regulated in response to dna damage. the list of genes in the regulon, show that in this case, regulation of cell cycle arrest upon dna damage is not linked to whib but only to the rv2719c-nrdr-rv2717c operon. hypoxia and induction of dormancy program: devr regulon one of the main characteristics of mtb is its ability to switch to a non replicating or ‘dormant’ state that allows it to survive for a long time within the host and renders it less susceptible to antibiotics. the environment inside the host is hypoxic and it might have high concentration of co or no released by the host macrophages. under either low oxygen concentrations or high concentrations of no or co, the heme iron from the kinases devs and dost becomes ferrous, the kinases become co or no bound and they get activated. in the active form, devs and dost autophosphorylate and induce phosphorylation of devr, which in turn, can bind its dna recognition sequence and induce expression of the devr regulon resulting in the activation of the dormancy program. all the options to build the co-expression networks that we have explored, result in a tight cluster for the known members of the devr regulon. this is to be expected since dormancy and devr regulon induction have been an active research topic in the mtb field. from the distinct conditions present in the expression compendium, almost 10% of them are conditions associated with expression of devr regulon. however, the fact that these genes always appear forming a tight cluster, points to the absence of additional regulatory elements that might cause a differential expression of some members of this regulon. we have selected this regulon to validate our methodology for chip-seq data analysis. chip-seq data corresponding to over-expression of devr were obtained from tbdb. the processed results for devr contain detected peaks, that correspond to genes that could be possibly regulated by devr, although in our compendium expression data for only of them were available. we have analysed this dataset following the methodology shown in figure  and defined a core-set containing genes. the analysis of the correlations in the expression profile among this set as compared to the overall distribution indicates a common regulatory influence over the selected genes, figure 7a. further analysis of the genes in this core set and their behaviour across the conditions in the compendium lead us to identify five distinct groups of genes within the identified targets. in four of these groups there is a high correlation among the genes whereas no clear pattern can be identified in the expression of the genes in the fifth group. one of these groups, that from now on we will refer to as the devr regulon, contains genes that were identified by chip-seq, has an expression pattern consistent with the previously described devr regulon, and of these genes have been previously identified as devr regulated genes. additionally, we have found among the list of targets from tbdb, genes that have been previously reported as devr regulated however their expression patterns suggest that either they are not regulated by devr or there is a secondary regulatory event altering their expression levels, therefore it is arguable whether they can be included in the regulon. a detailed list of members of this regulon is provided in additional file 13: table s the genes in this table have been assigned to functional categories: cell wall, transport elements, anaerobic respiration, translational machinery, regulatory elements and elements related to stress response. these six elements are linked to some of the main changes observed during growth arrest and dormancy, such as the changes in the cell wall, the arrest in protein synthesis and the adaptation to a hypoxic environment with reactive nitrogen species.figure devr regulon. a) left: histogram in blue represents the correlation among all the genes present in our compendium, whereas the histogram in green is based upon the correlations of the identified targets for which expression data is available within our compendium. both show the same overall distribution. right: histogram in blue represents the correlation among all the genes present in our compendium, whereas the histogram in green is based upon the correlations of the genes selected in the core set. note there is a shift towards positive correlation values, pointing to a common regulatory influence over the selected genes. b) number of genes identified in the devr regulon compared to the number of targets identified through chip-seq experiments or the ones cited in literature. c) group assignment of the targets identified by chip-seq for devr. g contains genes for which non discernible expression pattern has been found. g correspond to the usually named devr regulon, whereas g2- contains genes that show correlated expression patterns, although these patterns are not consistent with the previously described behaviour of devr regulon. a detailed list of these genes is available in additional file and the output of the go-enrichment analysis is shown in additional file interestingly, we found a faint link between the type vii secretion system genes rv0282-rv and devr. although the correlation analysis shows that they are not members of devr regulon, we indeed find that, in a reduced set of conditions they show a high correlation with devr regulon. these genes are also known to be regulated by zur and ider and are required for mycobactin-mediated iron acquisition in mtb. rv1734c belongs to the set of genes that had been previously reported as devr regulated but for which our analysis shows that their expression pattern is not compatible with this assertion. our analysis is based on transcript levels, so it could be that the regulation of rv1734c happens at the post-transcriptional level. however we have found additional evidence, albeit indirect, of our prediction. chauhan et al. analysed the effect of mutations in the devr binding motif and found that positions, and are essential and a substitution in any of them dramatically reduces the binding affinity. in addition, our analysis shows that position always contains an a. motif analysis complements the results from the analysis of expression data and further supports our prediction that rv1734c is not regulated by devr, since its putative binding site contains a mutation at position among the targets identified by chip-seq four other groups emerge however none of these groups show significant up or down regulation in the conditions associated with devr regulon induction. no clear common behaviour can be detected among the genes in group g these could be considered as false positives as a result of the devr over expression performed prior to the chip-seq procedure. however the genes within each of the other groups show a consistent co-expression across the conditions in the compendium, therefore we believe that these hits should not be discarded as being caused by non-biological reasons, instead alternative explanations such as weaker binding to regulate overall or regulation through a transcription factor homologous to devr should be further explored. in addition the genes in each group are functionally related, specifically, genes in g are mostly linked to translation and might be linked to the translation arrest observed during dormancy. genes in g and g are functionally linked to ‘metabolism’, ‘stress’ and ‘cell wall formation’, which are significantly different in the non replicating state. zinc uptake regulator zur as previously stated, the systematic analysis of genes linked to the different cog categories, showed that in the different co-expression networks the rpmb2-rpmg1-rpsn2-rpsr operon appeared forming a distinct cluster, separated from the rest of ribosomal protein coding genes. in addition, this operon also appears linked to the ppe3-rv operon. therefore, we concluded that these two operons should be regulated by a specific mechanism and respond to a specific type of perturbations. the location of these genes in the network of known regulatory interactions, pointed to the zinc uptake regulator zur as the most likely regulator of both operons. therefore we set forth to study the possible targets of this tf. initially we selected from the network extracted from the literature, a list of regulatory targets of zur that had been verified either by more than one source or by an appreciable upregulation in a zur knock out mutant. once this core set was selected, we proceed to identify the subset of conditions where these genes are co-expressed. the method to construct the expression plots was used to select conditions were the members of the core set were correlated; similarly we selected a core subset set of conditions and another subset of conditions. the correlation among the different putative members of the zurr regulon was used as a signature to identify the other members of the regulon. the list of genes identified as belonging to this regulon is provided in additional file and the zur consensus binding motif is depicted in figure 8b. the correlation analysis shows the importance of the biclustering approach to select only those conditions where no additional regulatory influences hinder the discovery of the members of the regulon. for example, if we analyse the correlation of rpmb with the members of the core set across all conditions present in our compendium, we obtain no correlation. however when we compute the correlation between rpmb and the genes assigned to the core set but only considering the previously selected subsets of conditions the correlation values raise to, and respectively, showing that, as previously reported, expression of rpmb is indeed regulated by zur.figure zur regulon. a) bicluster formed by members of the zur as reported in literature. the grey line represents the average expression levels of the members of zur regulon in the conditions in our compendium. the numbers identify the conditions that have been included in the bicluster. the horizontal bar and the different regions indicated by numbers refer to the classification of the conditions as described in materials and methods. notice the clear up-regulation of this set in conditions of type 9: transcription inhibition, in particular these values correspond to experiments were rifapentine was added to the medium. for clarity, expression values have been scaled, so that the mean value for each gene when all conditions are considered is zero. b) identified zur binding motif. e-value *10− c) number of genes in the zur regulon compared to the ones previously identified, the set anticorr. contains rv and lpqr, that show anticorrelation with the rest of the genes in the regulon. our predictions closely match the list proposed by maciag et al., however some differences in the target assignment appear. from the previously assigned genes to the zur regulon, we can only confirm that 79% of them show correlated expression patterns. in addition, among our predicted targets there are genes for which no experimental evidence can be found in the literature. these predicted additional members of the regulon are: rv0223c, lpqr, rv pe, ppe, rv2617c, rv rv2619c, rv3018c and rv3018b. rv and lpqr appear anti-correlated with the genes in the previously defined core set. in the case of lpqr this anti-correlation clearly increases when restricting the set of conditions. this together with the presence of zur binding motif in its upstream region, leads us to conclude that regulation of lpqr expression proceeds through a different mechanism than that of the rest of the genes in the regulon. additionally, it is striking to notice the clear up-regulation of this set of genes in conditions where transcription was inhibited by addition of rifapentine to the medium. this shows that, at least regarding these genes, the inhibition of transcription cause by rifapentine has a similar effect as zinc limitation. however, a detailed explanation of this phenomenon is most likely unreachable from the analysis of the present dataset.  we have shown that by integrating data from different sources and through the combined analysis of data, we are able to obtain insights into the biological system under investigation that go beyond the specific research questions of each experimental design. our systems level approach has allowed us to analyse in depth publicly available data on mtb and has enabled us to extract valuable new information from the existing datasets. the processes that we have studied are paramount in allowing mtb to thrive within the hostile host environment. we generated comprehensive lists of genes involved in the response to such environmental conditions. these compendia summarize and substantially extend and modify the current knowledge. we believe that any further research on these adaptation mechanisms will make extensive use of this new knowledge and the hypotheses generated herewith. we have developed a framework that allows the user to integrate information from gene or protein expression experiments, genome annotations and existing databases together with analysis tools. most of the already existing databases and tools provide human user interfaces that only allow for querying one gene at a time and only provide a subset of the total available information. within our framework this integration is done at once on a genome-scale, by using the gene co-expression networks. instead of obtaining a network with the majority vote from the algorithms we keep the individual output from each algorithm and rely on the expertise of the user to select which of them contains, for each particular problem, the highest amount of relevant information. this way of presenting the information greatly speeds up the process of analysing new data, since after the study of the different networks and the modules that each of them show, expert users will be able to quickly understand the changes occurring during a new experiment just by overlaying the new expression data over the networks and looking for the networks were the differentially expressed genes or proteins appear clustered together. our approach is not, in itself a new analytical approach but a general framework that facilitates the consideration of multiple pieces of evidence to pose a hypothesis on a biological system. comparison with experimentally verified networks and with in silico generated datasets shows that the co-expression networks generated using our pipeline preserve the desired modularity of transcriptional networks and the regulatory targets of a tf tend to appear in the networks as forming a closely interconnected module. comparison with e. coli data, shows that in 97% of the cases the target genes of a tf form clusters in the network. modularity was on average less preserved for s. cerevisiae networks, due most likely to the increased complexity of regulatory interactions in eukaryotes. this poorer performance of the inference methods for the yeast dataset was also observed in the dream challenge were at most a recall was obtained by any of the tested methods. on average the zrc networks perform better that the mixed ones, but our analysis also showed that both types of networks contain complementary information, for example the targets of a yeast's ypr199c are only identified as forming a module in the mixed network. we have established a protocol to assign regulatory interactions to binding sites identified through chip-seq experiments via the integration of expression data. we have used this approach to correctly identify the target genes of a given tf as a response to an specific type of perturbation among the hundred of candidates from the experimental dataset. within our framework biclusters can quickly be analysed, interactions between the biclusters can be identified in the overall networks and biclustering methods can easily be transformed into a tool for automatic detection of functionally related modules or underlying layers of regulation. moreover, we have developed a method that allows to identify additional regulatory layers. the comparison of the networks generated by considering different subsets of conditions allowed us to distinguish various regulatory mechanisms for dna damage response. the analysis of zur regulon shows the potential of an integrative approach. in our compendium there were no data corresponding to experiments designed to analyse the effects of zinc limitation on mtb. however, the analysis of the expression patterns of the genes in the different conditions and the analysis of their correlations, allowed us to select a set of conditions that would complement the bioinformatic analysis of the upstream sequences of the genes and would allow us to decided which of the regions similar to the motif actually represents a zur binding site and characterize the members of the regulon. we were able to compare our predictions with dedicated experiments performed with knock-out mutants and we found a good agreement between computational predictions and experimental data. our work extends the already existing knowledge and produces a comprehensive list of the members of the devr regulon. among the targets identified by chip-seq, we have uncovered three additional sets of genes that show a consistent expression pattern across the conditions in the compendium and are functionally related. additional studies are required to further understand the regulation of these genes and their possible link to the non-replicative state. the discovery of new regulatory mechanisms involved in dormancy has the potential to deliver a new set of drug targets. the automatic biclustering was the basis of our analysis of dna repair systems and specifically lexa regulon. we were able to confirm our computational predictions through the comparison we literature data and recently performed chip-seq experiments. in addition, our work has produced a more specific binding motif for lexa through the identification of new members of its regulon. additionally, the re-annotation of the identified new targets allowed us to identify a putative mr system in mtb. the analysis of mutations and mechanisms to avoid them is of the uttermost importance for the further understanding of the evolution of antibiotic resistances and pathogenicity. previously existing data were used to verify our predictions on the regulatory mechanisms of dna damage response. we have clarified some points previously in dispute, such as the lack of involvement of sigg in the response to dna damage and the regulation of the alternate sigma factor sigc in these conditions. correct identification of the sigma factor up-regulated upon dna damage is key to understanding the systemic response of mtb to this damage type. in addition, our work has provided further evidence on the mechanisms leading to cell cycle arrest upon dna damage in mtb. in addition, we identified a new regulatory mechanism for zur, since the analysis of the upstream regions of its target genes and their expression patterns show that rv and lpqr belong to its regulon, although their regulation must proceed through a different mechanism. here, we have presented the results obtained by applying our integrative approach to mtb and whenever additional data was available, we have found good agreement between predictions and experiments. the basic underlying principle of this approach relies on the comparison of the networks obtained using different sources of information or methodologies. this approach can be readily extended to different organisms and to the comparison between different species, by using global identifiers together with a database of orthologous genes between species. this would allow to select a gene or group of genes in one organism and see how they are arranged in the networks corresponding to a different organisms. in addition, other types of data, such as synteny or evolutionary information, and protein structure and families could improve the evolutionary comparison of functional modules. additional files additional file 1: microarray data sets considered in this meta-analysis. identifiers from the gene expression omnibus database have been used for each sample, experiment and platform. additional file 2: identified biclusters. the following information is presented for each bicluster. left. top: expression plot, the grey line represents the average expression levels of the genes in all the conditions in this meta-study, the numbers identify the conditions in the biclusters. the horizontal coloured line shows the classification of the conditions. middle interactions among the members of the biclusters either from protein-protein interactions; operons or similarity of annotated cog terms. the nodes in the network represent the genes in the biclusters. the colours are linked to the cog annotation, whereas the shape is linked to the functional classification: square: information storage and processing; rhomboid: cellular processes and signalling; triangle: metabolism; and circle: poorly characterized. bottom: identified motifs. right: genes in the bicluster, their upstream sequences have been plotted together with a representation of the location of the detected motifs: red, green and blue for the first, second and third motif. additional file 3: pipeline for the generation of co-expression networks. the r code that implements the pipeline as depicted in figure  the inequality simplifier is implemented in the joined jar file. additional file 4: technical characteristics of the visualization tool. technical description of the network visualization and comparison tool used. additional file 5: relationship between the expression levels of reca and lexa. the fit corresponds to a model, based on ordinary differential equations, of lexa mediated induction of reca. the numbers signal the conditions where this model actually holds. for clarity, expression values have been scaled, so that the mean value for each gene, when all conditions are considered, is zero. the horizontal bar includes the overall categorization of the conditions. notice that the conditions were lexa is regulating reca expression are mainly linked to 1: aromatic amides intra cellularly hydrolysed, low ph; 4: acidified medium; 5: cell wall synthesis inhibition; and 8: dna damage. additional file 6: comparison of predictions for individual methods to consensus network. a) overall predictions. the plot shows the number of interactions from the gold standard that were predicted worse, equally or better by each individual method than by the consensus network. the predictions were assessed comparing the rank of the predictions with a cutoff. only the e. coli network has been considered. b) predictions for each tf. fraction of the interactions that are better predicted by each type of method than by the consensus network. only tf with more than interactions have been considered. these plots have been built using the methods described in. additional file 7: topological overlap of tf targets in the co-expression networks obtain for the dream challenge datasets. for each dataset the zrc and mixed networks were built. only tf with more than experimentally verified targets were considered. dashed line represents the average topological overlap in each network. additional file 8: relative locations of the cluster of ribosomal proteins and the operon of ribosomal proteins rpmb2-rpmg1-rpsn2-rpsr2 in different co-expression networks. additional file 9: table s lexa regulon. additional file 10: genes in the reca_nd regulated dna damage repair set. operons and divergons are marked by alternating row colours. the fifth column indicates whether the targets had been previously identified. additional file 11: identified motif in the upstream regions of the genes in the reca_nd dna repair system. logo generated using meme. e-value e- the location of the hits is included in additional file additional file 12: classification of devr targets identified by chip-seq. the column tbdb indicates whether the target was listed in the original data set as a possible target for devr regulation, some genes have been added to the table based on operon/divergon analysis. hit location indicates the genomic location of the hit. the last column refers to the initial classification of the genes within the core set. the column references refers to whether they have been previously identified as members of devr regulon: a identified by park et al.; b identified by balázsi et al.; c identified by bartek et al.; d identified by chao and rubin. groups are present in the table. "the group contains those genes whose expression patterns showed no significant correlation with any others, and we believe are not regulated by devr, however they have been previously identified as belonging to devr regulon in literature, therefore we have included them in the table. the first group corresponds to the members of devr regulon. additional file 13: table s devr regulon members. additional file 14: bicluster formed by members of devr regulon and the type vii secretion system genes rv282-rv the grey line represents the average expression levels of these genes in the conditions in our compendium. the numbers identify the conditions that have been included in the bicluster. additional file 15: devr motifs top: previously identified primary and secondary motifs for devr [ 84 ] bottom: newly proposed motifs. main differences in primary motif can be seen at positions and main differences in secondary motif can be seen at position and to additional file 16: results of the go enrichment analysis performed in the classified devr targets identified by chip-seq. node size reefers to the number of genes annotated to that particular go term in the whole genome annotation file, whereas sample match indicates the number of genes annotated to the go term in the group. additional file 17: members of zurr regulon. operons and divergons are marked by alternating row colours. rows signalled in yellow refer to genes that have been previously identified as belonging to this regulon, but for which our analysis does not support this claim. the fifth column indicates whether the targets had been previously identified by maciag et al. and the sixth column reflects the reported fold change in a zurr knock out mutant. column contains the definition of the core set of genes, whereas columns to contains the average correlation between the expression level of the each gene and the genes in the core-set across all conditions in the compendium or a subset of them. abbreviations tbdbtb database tftranscription factor mimutual information mtbmycobacterium tuberculosis mrmismatch repair gogene ontology dpidata processing inequality ssdnasingle stranded dna reca_ndreca non-dependent aurocarea under receiver operating characteristic curve auprarea under precision recall curve competing interests the authors declared that they have no competing interests. authors’ contributions jd and msd: study concept and design, methodology and tools development, manuscript preparation and data analysis. ps and vmds: study concept and design and critical revision of the manuscript. all authors read and approved the final manuscript. 
as a show case, we have analysed the gene co-expression networks of mycobacterium tuberculosis generated using over expression experiments. our workflow is based on representing the multiple types of information as network representations and presenting these networks in a synchronous framework that allows their simultaneous visualization while keeping specific associations from the different networks. however each separate set has an intrinsic value that is diluted and partly lost when building a consensus network. here we present a methodology to generate co-expression networks and, instead of a consensus network, we propose an integration framework where the different networks are kept and analysed with additional tools to efficiently combine the information extracted from each network. as show cases we implemented a pipeline to integrate chip-seq data and another method to uncover multiple regulatory layers. regarding dna damage repair, we identified sigc as a key control element, new targets for lexa, an updated lexa binding motif, and a potential mismatch repair system. meta-studydata integrationregulatory networksreverse engineeringmycobacterium tuberculosisdormancydevrzinc uptake regulationdna repairissue-copyright-statement© the author 2014 we discovered new genes linked to zinc uptake and a new regulatory mechanism for zur. by simultaneously exploring these networks and metadata, we gained insights into regulatory mechanisms in m. tuberculosis that could not be obtained through the separate analysis of each data type. our methodology relies on providing the user the means to simultaneously visualise and analyse the coexisting networks generated by different algorithms, heterogeneous datasets, and a suite of analysis tools. each algorithm produces different networks and efforts have been devoted to automatically integrate them into consensus sets. we expanded the devr regulon with genes while identifying targets wrongly assigned to this regulon. different methods have been developed to infer regulatory networks from heterogeneous omics datasets and to construct co-expression networks. the use of co-expression networks to perform system level analysis allows the development of custom made methodologies.
 different methods have been developed to infer regulatory networks from heterogeneous omics datasets and to construct co-expression networks. each algorithm produces different networks and efforts have been devoted to automatically integrate them into consensus sets. however each separate set has an intrinsic value that is diluted and partly lost when building a consensus network. here we present a methodology to generate co-expression networks and, instead of a consensus network, we propose an integration framework where the different networks are kept and analysed with additional tools to efficiently combine the information extracted from each network.  we developed a workflow to efficiently analyse information generated by different inference and prediction methods. our methodology relies on providing the user the means to simultaneously visualise and analyse the coexisting networks generated by different algorithms, heterogeneous datasets, and a suite of analysis tools. as a show case, we have analysed the gene co-expression networks of mycobacterium tuberculosis generated using over expression experiments. regarding dna damage repair, we identified sigc as a key control element, new targets for lexa, an updated lexa binding motif, and a potential mismatch repair system. we expanded the devr regulon with genes while identifying targets wrongly assigned to this regulon. we discovered new genes linked to zinc uptake and a new regulatory mechanism for zur. the use of co-expression networks to perform system level analysis allows the development of custom made methodologies. as show cases we implemented a pipeline to integrate chip-seq data and another method to uncover multiple regulatory layers.  our workflow is based on representing the multiple types of information as network representations and presenting these networks in a synchronous framework that allows their simultaneous visualization while keeping specific associations from the different networks. by simultaneously exploring these networks and metadata, we gained insights into regulatory mechanisms in m. tuberculosis that could not be obtained through the separate analysis of each data type. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords meta-studydata integrationregulatory networksreverse engineeringmycobacterium tuberculosisdormancydevrzinc uptake regulationdna repairissue-copyright-statement© the author 2014 