BACKGROUND
over the last decade, the availability of fully sequenced genomes and the development of high-throughput technologies such as dna microarray-based transcript profiling have fuelled an exponential increase in the volume of functional genomics data.
this has led to a renewed interest in the study of molecular biology at the system level  <cit> .
the central paradigm in systems theory is that one can learn about a system by perturbing it and measuring the response.
this principle also applies to biological systems.
since mrna levels can nowadays easily be measured on a genome-wide scale, expression profiling has become a first method of choice for assessing the molecular response to experimental perturbation .
considerable efforts are put into creating compendia of expression profiles under genetic, chemical or environmental perturbations  <cit>  or in different tissues  <cit> .
such data compendia basically constitute a series of snapshots of expression states under a variety of conditions, and they contain a wealth of information concerning the underlying transcriptional network structure of an organism.
however, devising methods to efficiently and reliably extract that information is still a challenging task.
clustering of gene expression data allows the inference of functional correlations between genes through what was dubbed the 'guilt-by-association' principle  <cit> .
a classical clustering process consists of two steps  <cit> .
first, a matrix of distances between expression profiles is calculated using a distance or similarity measure, such as pearson's centered correlation coefficient .
based on this distance matrix, the actual clustering algorithm, for instance average linkage hierarchical clustering, groups similar profiles together.
traditional clustering methods are well suited for analyzing time-series expression data, but they fall short when applied to perturbational data, because the underlying similarity measures, such as pcc, primarily capture global correlation tendencies.
however, in compendia of perturbed expression profiles, genes do not necessarily show similar behavior under all experimental conditions: they may be coexpressed under some conditions and follow different expression regimes under other conditions.
one of the consequences is that genes may be coexpressed with multiple expression modules depending on the conditions, or in other words, expression modules may overlap.
these observations stimulated the development of alternative clustering strategies.
the process of detecting subsets of genes that exhibit similar expression behavior across a subset of conditions is known as biclustering.
several biclustering strategies exist today, each using its own heuristic approach to tackle this complex problem .
some biclustering methods use a greedy iterative search strategy to uncover biclusters, progressively subdividing, or adding and removing rows and columns from the biclusters obtained in a previous iteration in order to maximize a local score function  <cit> .
others use linear algebra  <cit>  or adopt a graph-theoretic approach to biclustering  <cit> .
yet other methods identify biclusters by proposing a statistical model and estimating the distribution parameters that minimize a certain model fit criterion  <cit> .
a feature that most biclustering methods share is that they do not explicitly define similarity measures on the global space of expression profiles, but rely on the emergent properties of groups of genes and conditions in order to identify significant subpatterns in the data.
evidently, a wide variety of biclustering algorithms exist, each of them having their own strengths and weaknesses.
for example, some of these methods are intrinsically less suited to find overlap between biclusters because they mask previously found biclusters with random noise  <cit> , or because they partition the data  <cit> .
others require extensive parameter tweaking, require the user to specify the desired number of biclusters in advance or generate very small or large  biclusters or highly redundant biclusters .
some have no publicly available implementation or are rather cumbersome to use, and most of them, notable exceptions being samba/expander  <cit> , genomica  <cit>  and cmonkey  <cit> , do not integrate or overlay other types of biological data, hampering their use as biological discovery tools.
also, to our knowledge, none of the existing biclustering methods uses the variational information in replicated expression experiments.
this information is routinely and successfully used to detect genes that are differentially expressed under a given perturbation  <cit> .
the main reason why biclustering methods do not use differential expression information is that they do not specifically focus on the analysis of perturbational data.
discretization-based biclustering methods such as samba  <cit>  and bimax  <cit>  could probably easily be modified to assess up- and downregulation of gene expression based on p-values for differential expression.
in their current implementation, however, these methods use rather arbitrary log-ratio or percentage cutoffs for this purpose.
in this study, we present a novel method, called enigma, that addresses some of these issues.
our goal was to build a method that:  leverages differential expression analysis results to extract co-differential expression networks and expression modules from perturbational gene expression data,  is able to detect significant partial coexpression relationships between genes and overlap between modules,  depends on parameters that can be automatically optimized or set on reasonably objective grounds.
produces a realistic amount of modules, and  visually integrates the expression modules with other data types such as gene ontology  information  <cit> , transcription factor  binding data, protein and genetic interactions, in order to facilitate the biological interpretation of the results.
below, we outline the enigma algorithm, test our methodology on artificial expression data and compare its performance to other methods.
we also apply enigma to a perturbational microarray compendium for budding yeast  <cit>  in order to assess its potential to generate testable hypotheses on real biological data.
RESULTS
algorithm
a global overview of the methodology is given in figure  <dig>  briefly, enigma takes as input a set of perturbational expression data, externally calculated p-values for differential expression  and other data types if available.
enigma uses a novel combinatorial statistic to assess which pairs of genes are significantly co-differentially expressed .
the resulting coexpression p-values are corrected for multiple testing and translated to edges in a coexpression network, which is clustered into expression modules  using a graph-based clustering algorithm inspired on the mcode algorithm  <cit> .
the clustering procedure depends on two parameters that control the density of individual modules and the overlap between modules.
the main reason why we chose a two-tier clustering approach  is that it allows simulated annealing-based optimization of the clustering parameters to obtain optimal coverage of the coexpression network, in terms of module overlap and redundancy.
the graph clustering method we use is very fast, which allows the parameters to be optimized in a reasonable amount of time.
in the postprocessing phase, enigma determines relevant condition sets for each module, visualizes their substructure and overlap with other modules, screens the modules for enriched go categories, suggests potential regulators for the modules based on regulator-module coexpression links and enrichment of tf binding sites, and overlays protein and genetic interaction data.
combinatorial statistic
consider the expression profiles of two genes a and b under n perturbations .
each gene is represented by a profile of n fields.
the gene expression values are discretized into three categories  based on their differential expression p-value.
if the gene is significantly upregulated in a given experiment , the corresponding field is labeled blue.
experiments in which the gene is significantly downregulated are similarly labeled yellow, and the remaining fields are labeled black.
let us now assume that the profiles of a and b contain ax and bx blue fields respectively, as well as ay and by yellow fields, and that they have x blue and y yellow fields in common.
we want to assess whether this overlap is statistically significant.
if the response of the genes a and b to the perturbations were uncorrelated , the blue and yellow fields would be independently distributed on both profiles.
under this hypothesis, the probability that the profiles overlap on exactly x blue and y yellow positions is given by the following recursive formula:
p=−∑x′=xmin⁡∑y′=ymin⁡≠p
the probability of observing an overlap of at least x blue and y yellow fields by chance is then expressed by the cumulative distribution:
pc=∑x′=xmin⁡∑y′=ymin⁡p
equation  <dig> can be understood by assuming that profile a is given, and that we randomly distribute bx blue and by yellow positions on profile b. the denominator of the first term then represents the total number of possible profiles b. the numerator represents the combinations in which x blue and y yellow matching positions are selected, and the residual positions are chosen at random.
however, in this manner, a number of combinations are selected while having more than exactly x blue and/or y yellow matching positions.
moreover, combinations with x' > x blue and/or y' > y yellow matching positions are counted c·c times, hence the second term .
although the probabilistic question formulated above can be cast in terms of contingency tables, the hypothesis tested by our statistic is different from that tested by standard contingency table analysis methods such as the χ <dig> test.
for example, situations in which a large amount of blue  fields in profile a are perfectly mapped onto the black fields  in profile b would yield a significant χ <dig> p-value, whereas they would not yield a significant p-value using equation  <dig>  our statistic only considers mappings of up- and down-regulation of the expression of a gene to up- or down-regulation of another gene to be meaningful for assessing coregulation, a premise which is motivated by the perturbational nature of the data we aim to analyze.
black fields are considered less informative from the perspective of coregulation.
multiple testing correction of coexpression p-values
in our probabilistic setup, each comparison of two profiles can be considered an individual test.
for n genes, n/ <dig> tests are performed to fish for co-differential expression relationships.
consequently, the obtained p-values have to be adjusted in order to control the type i error rate.
the raw p-values are corrected for multiple testing with the benjamini & hochberg procedure, which controls the false discovery rate   <cit> .
graph-based clustering
the set of significant coexpression relationships at a certain fdr threshold  is translated to a network, with nodes and edges representing genes and significant coexpression relationships, respectively.
enigma identifies coexpression modules from this network using a graph clustering technique inspired by the mcode algorithm  <cit> .
to identify potential module seeds, all nodes v are weighted based on the density of the highest k-core of the node neighborhood nv, denoted as the kmax-core of v .
analogous to bader and hogue  <cit> , the core-clustering coefficient ccore,v is defined as the density of the kmax-core of v, and the weight wv = ccore,v·kmax,v.
the kmax-core of the node with the highest weight is taken as the first module seed.
this module seed then grows by accreting nodes on which it exerts a pull above a certain threshold ν <dig>  the pull of a module with seed s on a node v outside the module is defined as |nv ∩ s|/|s|.
the next module is then initiated by taking the kmax-core of the node with the highest weight in the remaining graph.
an additional constraint is set by requiring that the overlap between the new seed s and any existing module m does not exceed ν1·min.
while the threshold ν <dig> controls the size and density of individual modules, ν <dig> controls the spacing or overlap between modules.
both parameters are optimized automatically.
clustering parameter optimization
in order to optimize the clustering parameters, the quality of the clustering for a given  is assessed by comparing the known input coexpression network  with the output coexpression network inferred by the modules.
the latter is constructed by translating the modules to fully connected components in the output network .
if we consider true/false positives  to be coexpression edges inferred by the clustering that are present/absent in the input coexpression network, and false negatives  as edges present in the input network that are not inferred by the clustering, we can define the precision p' = tp/ and the recall r' = tp/ of the clustering result.
enigma uses the f'-measure, i.e.
the harmonic mean of recall  and precision , f' = 2p'r'/, as a measure for the quality of the clustering.
we use the notation p', r', f' instead of the more commonly used p, r, f in order to distinguish between two different flavours of the f-measure used in this study for different purposes.
in contrast to the regular f-measure , the f'-measure penalizes overpredicted edges in order to avoid unnecessary overlap between the expression modules: an edge  that is inferred multiple times from the clustering, because the genes a and b belongs to the intersection of multiple  modules, is counted as  <dig> tp and x -  <dig> fp.
this is equivalent to drawing x edges between the genes a and b in the output coexpression network.
since there is only one edge in the input network, the x -  <dig> remaining edges can be considered false.
this penalization strategy has the intuitively pleasing property of not affecting the recall, but lowering the precision of the clustering result when the amount of edges 'explained' by multiple modules increases.
the parameters ν <dig> and ν <dig> are now optimized by monte-carlo simulated annealing   <cit>  using f' as the optimization criterion.
starting from an random initial guess for the parameters , random steps are taken in parameter space.
a step is accepted if
rand <eΔf'/t
with rand a random number drawn uniformly from the interval  <cit> , Δf' the change in f'-measure and t the simulated annealing parameter or 'temperature', which gradually decreases during the course of the optimization according to an exponential scheme ti = rcti- <dig>  with rc the cooling rate.
enigma uses a two-stage mcsa procedure.
in the first stage, a rough mcsa search of the clustering parameter space is performed in order to identify the most interesting parameter region .
in the second stage, a finer mcsa search is performed starting from the optimum obtained in the first stage .
at the end of each stage, an additional gradient descent is performed toward the nearest local optimum of f'.
by default, enigma performs  <dig> mcsa runs, starting from randomly chosen .
the convergence of the solutions of multiple runs can be used as a check on the adequacy of the mcsa parameter settings.
postprocessing of modules
for each gene module, enigma determines a condition set by selecting those conditions that show enrichment of up- or downregulated genes in the module .
thus, for a given module, the condition set contains the experimental conditions that elicit a significant and specific response in the module  and, by consequence, have been most influential in shaping the module.
the resulting 'bicluster' does not necessarily have a uniform expression pattern over all genes, but may show subpatterns for some genes under certain conditions, possibly indicating involvement in other expression modules.
these subpatterns are visualized by hierarchically clustering the module's expression data in both dimensions, using the cosine correlation coefficient  as a similarity measure.
the clustering tree can optionally be separated into leafs to make the subdivision more clear .
although conditions that show differential patterning within one module might appear to be irrelevant for the module as a whole, they are important for at least part of the module and may provide insight into inter-module connections or further substructure within the module.
in an attempt to provide the user with clues on how the expression modules are regulated, enigma searches for 'regulators' that are significantly more connected to a module, through positive or negative coexpression edges, than expected at random .
potential regulators are selected from a user-defined list or a user-defined set of go classes.
when chromatin immunoprecipitation  or tf motif data are available, enigma also screens the modules for enriched tf binding sites .
the expression profiles of significantly coexpressed or binding regulators are visualized on top of the modules.
significantly enriched go terms for both the gene and condition sets of the modules are determined using the bingo  <cit>  software, which is incorporated in enigma .
finally, enigma visually maps the available protein interaction data and genetic interaction data on the modules.
testing on artificial data
generating artificial expression data
to assess the performance of our method and compare enigma to other methods, we performed tests on artificial gene expression data.
we generated two types of artificial expression data, namely expression data containing overlapping biclusters  and expression data containing partially coexpressed genes but no biclusters .
in both cases, we built  <dig> expression datasets of  <dig> genes by  <dig> experiments .
for each dataset, artificial background expression data were randomly sampled from a normal distribution with mean μ =  <dig> and variance σ <dig> =  <dig> .
for the modular datasets, we implanted  <dig> biclusters in this background, each encompassing between 1–5% of all genes and 10–50% of all conditions.
bicluster sizes, member genes and conditions are chosen at random, with the restriction that at most 30% of the genes and 50% of the conditions overlap between any  <dig> biclusters .
except for a noise component , all genes in a bicluster share the same expression profile over the bicluster conditions.
however, a bicluster can be partially overwritten by other biclusters.
the bicluster profiles are sampled from a bimodal distribution consisting of  <dig> normal modes with means μ <dig> = - <dig>  and μ <dig> =  <dig>  and variances σ12=σ22= <dig> .
the expression profiles of individual genes in a bicluster are noisified by adding normally distributed noise .
the variances, bicluster size and overlap parameters are chosen so that the overall distribution of the simulated log ratio expression values mimicks the distribution of log ratio expression values in the rosetta compendium  <cit>  up to a scale factor .
note that, apart from the distribution of expression ratios, the structure of these toy datasets does not necessarily bear any resemblance to real biological data.
for the non-modular datasets, we implanted  <dig> pairs of partially coexpressed genes  in the background.
the expression profiles are constructed as described above.
the resulting expression value distribution again mimicks the rosetta distribution .
unlike for real data , we used log <dig> ratio thresholds to discretize the expression values of the artificial datasets, because the generation of meaningful artificial differential expression p-values proved to merit further study in its own right.
therefore, the artificial data cannot be used to assess the advantage of including variational information in enigma's discretization step .
on the other hand, we can still compare the performance of enigma with other methods that do not use variational information.
we used a log <dig> ratio threshold of  <dig> for upregulation and - <dig> for downregulation, corresponding to the means of the distributions used to generate the bicluster profiles.
in other words, half of the datapoints in the biclusters are presumed not to be significantly over- or underexpressed.
performance of enigma on artificial data and comparison with other methods
the performance of enigma on these toy datasets was compared with that of two commonly used similarity measures, namely pcc and the χ2-statistic, and two established biclustering methods, samba  <cit>  and isa  <cit> .
pcc was chosen as a representative of the global similarity measures used in traditional clustering algorithms, while we included the χ2-statistic because of its relation to the combinatorial statistic used by enigma .
the selection of biclustering methods was based on the following criteria:  the methods should be non-partitioning in nature,  they should have the capacity to generate overlapping biclusters,  a suitable implementation should be publicly available, and  they should produce a reasonable amount of biclusters  on the modular toy datasets.
we used the version of samba  <cit>  incorporated in the expander  <dig>  package  <cit> , and the implementation of isa  <cit>  available as part of the biclustering tool bicat  <cit> , both with default parameter settings.
the isa trajectories from randomly chosen starting points  converge to a limited number of 'fixed point' biclusters.
to prune nearly identical modules, we merged isa biclusters that overlap for more than 80%.
the clustering performance of all methods is only assessed in the gene dimension.
standard internal criteria for partitional clustering performance, such as the silhouette width or dunn's index  <cit> , cannot be used to assess the performance of algorithms that generate overlapping clusters.
instead, we use the f-measure and introduce a derivative, the f'-measure , to compare the performance of different clustering methods on artificial datasets.
in both cases, the coexpression network generated by a method  is compared to the artificial input coexpression network in terms of true and false positive edges and false negative edges, from which the different flavors of the f-measure are calculated .
the difference between the f-measure and the f'-measure is that the f-measure does not take into account the multiplicity of the inferred edges.
in other words, the f'-measure penalizes overpredicted  edges, whereas the f-measure does not.
this entails that the f'-measure is more useful to compare methods that generate overlapping clusters, whereas the f-measure can be used more generally to compare methods that generate both overlapping or non-overlapping clusters or pair-wise coexpression networks.
the performance of enigma is tested on two levels by assessing the overlap between the artificial input correlation network and  the network of significant correlations obtained in the first step of the enigma algorithm ;  the modules inferred by enigma .
the output networks for enigma-m and the biclustering methods samba and isa are obtained by converting the obtained modules/biclusters to fully connected network components.
the χ <dig> network is constructed by translating significant χ <dig> correlation p-values between the discretized expression profiles to edges in the output network.
we used the same discretization threshold  and fdr level  for the χ <dig> and enigma methods.
the performance of pcc was measured for different thresholds .
using the f-measure, enigma outperforms all other methods on the modular artificial data .
the performance of enigma-m was consistently higher than the χ <dig> performance  and the optimal pcc performance .
the global similarity measure pcc appears to perform surprisingly well.
however, the performance of pcc critically depends on the choice of the pcc threshold, and determining the optimal pcc threshold on real data is problematic.
in contrast, enigma has the advantage of having an easily tunable significance threshold: the false discovery rate  level.
to illustrate this, we plotted the performance curve of enigma for different non-corrected p-value thresholds , on figure 2a and 2b.
for all artificial datasets, the performance of enigma-n at fdr =  <dig>   is close to the optimum of this curve, indicating that fdr control at a reasonable level gives near-optimal performance.
among the biclustering methods, the rather poor performance of the isa algorithm  may seem somewhat surprising.
prelić et al  <cit> , using the same implementation of isa but other methods to generate artificial data and to assess biclustering performance, previously established that the performance of isa decreases with increasing overlap between biclusters.
our results seem to confirm that isa is not the optimal method in case there is substantial overlap between modules.
the performance of isa did not change significantly when using  <dig> starting points instead of the default  <dig> .
the performance gain of enigma-m over samba is substantially smaller , and on two out of  <dig> datasets, the performance of samba was slightly higher than that of enigma-m .
a more tangible advantage of enigma over samba  is that enigma nearly always recovered the correct number of modules , whereas samba consistently predicted more modules than there were in the input data .
isa predicted only one extra module on average, but with a higher variance than enigma .
in other words, samba and to a lesser extent isa produce more fragmented and/or more redundant modules.
redundancy makes the module output much harder to interpret, but it is not taken into account by the standard f-measure.
to quantify the effect of redundancy on the clustering quality, we compared samba, isa and enigma-m using the f'-measure.
as in the calculation of the f'-measure used in the clustering optimization procedure , edges that are inferred by multiple modules are counted multiple times, but in the present case, multiply defined edges may also occur in the input network if they overlap between multiple artificial input modules .
specifically, edges that are inferred by x output modules and y input modules are now counted as y tp and x - y fp in case x ≥ y, or x tp and y - x fn in case x <y.
using the f' criterion, the performance of enigma-m  is substantially higher than that of samba  and isa .
on non-modular artificial data, the performance of enigma-m and the biclustering methods samba and isa is very low .
this is not surprising since there are no modules to be found in these datasets.
in this respect, a particularly attractive feature of enigma is that it finds very few modules in the non-modular data , in contrast to isa and samba, which recover  <dig> ±  <dig> modules  and  <dig> ±  <dig> modules , respectively.
among the pair-wise methods, enigma-n invariably featured the highest performance, indicating that our combinatorial statistic detects partial coexpression relationships more efficiently than pcc and χ <dig>  the fact that enigma efficiently uncovers coexpression relationships in non-modular data opens perspectives for the exploration of the less modular parts of expression datasets.
real datasets typically contain a limited number of perturbation experiments that target a few specific processes.
these processes can be expected to be rather well resolved in terms of their coexpression relationships, whereas other processes will probably give rise to more fragmented  regions in the network.
moreover, despite the success of the modularity concept in the analysis of expression data and systems biology in general, it is not inconceivable that transcriptional networks might also contain genuinely non-modular regions.
testing on real data: the rosetta gene expression compendium
p-value versus log-ratio based discretization
although useful for testing and comparing methods, artificial datasets do not capture the complexity of real biological systems.
consequently, good performance on artificial data does not guarantee good performance on real biological data.
in order to assess the use of enigma for analyzing real data, we applied our methodology to the rosetta compendium of expression profiles, representing data on  <dig> different experimental perturbations of s. cerevisiae  <cit> .
experiments on  <dig> strains that were marked as aneuploid in the original dataset were left out, because they can give rise to artificial expression correlations between genes on the aneuploid chromosomes.
the log-ratio expression data and differential expression p-values were downloaded in prenormalized and preprocessed form.
genome-wide chip data for  <dig> tfs were obtained from harbison et al  <cit> .
all genes that are bound with p <  <dig>  by a certain tf were considered reliable targets.
protein and genetic interactions for s. cerevisiae were obtained from the biogrid database  <cit> .
using a differential expression p-value threshold of  <dig>  in the discretization step and an fdr threshold of  <dig>  for defining coexpression edges, enigma identified a network of  <dig>  significant positive coexpression links and  <dig>  negative coexpression links involving  <dig>  genes.
the clustering parameters  =  were optimized by mcsa as described in the algorithm section.
to assess the efficiency of the mcsa procedure, we performed an exhaustive screen of the parameter space to locate the global maximum of the f'-measure criterion .
the mcsa procedure found back the global optimum with 100% efficiency.
enigma discovered  <dig> modules in the rosetta dataset encompassing  <dig> genes and  <dig> conditions .
these numbers seem reasonable given that  <dig> of the  <dig> conditions included in the compendium contain less than five differentially expressed genes, which entails that they have a small chance of contributing to a module.
given the low amount of informative conditions, it is not surprising that only a third of the s. cerevisiae genes can be included in modules.
according to the go enrichment results,  <dig> out of  <dig> modules have a significant degree of functional coherence.
fifty-four modules are enriched in targets of one or more tfs, and  <dig> modules show enrichment of both go biological process categories and tf binding sites.
together, 60% of the modules show enrichment of go categories and/or tf binding sites, indicating that our method is capable of identifying biologically relevant expression modules.
in order to qualitatively assess the effect of using a differential expression p-value cutoff in the discretization step instead of a fold-change cutoff, we repeated the analysis using a |log <dig> ratio| discretization threshold of  <dig> .
the resulting coexpression network contains  <dig>  positive and  <dig>  negative links between  <dig>  genes.
the clustered network contains  <dig> modules encompassing  <dig>  genes.
ninety-three modules exhibit go enrichment,  <dig> exhibit tf binding enrichment and  <dig> exhibit both.
despite the significantly lower amount of connections in the log-ratio network, the number of functionally coherent modules and the number of clustered genes is roughly similar, and the optimized clustering parameters  =  are identical, indicating that the general structure of the network and its strongest modules are fairly well preserved.
indeed, many highly functionally coherent modules  are present in both networks.
not incidentally, many of these modules are related to the processes that were targeted by hughes et al  <cit> , which can be expected to show a pronounced expression response.
however, modules that show less pronounced expression variations, for example the modules related to ribosome biosynthesis, are not recovered in the log-ratio network.
this illustrates the main disadvantage of using a fixed log-ratio threshold: different processes show different amplitudes of expression change upon perturbation, which cannot be captured by a single threshold.
one could argue that this can easily be remedied by standardizing the expression profiles to zero mean, unit variance before applying the threshold, as is done by some methods, e.g.
samba  <cit> .
however, in the case of perturbational data, this manipulation runs the risk of effiectively breaking the connection to the reference condition, thereby distorting the meaning of up- and downregulation and introducing serious artifacts .
topological characteristics of the enigma coexpression network
since many cellular functions are carried out in a highly modular manner  <cit> , most cellular networks, including protein interaction networks, metabolic networks and gene expression networks, are modular in nature  <cit> .
on the other hand, many cellular networks, including coexpression networks, have been claimed to exhibit a node degree  distribution of the power-law type, p ~ k-γ, indicative of scale-free properties  <cit> .
the coexistence of modularity and a scale-free degree distribution can be explained by assuming a hierarchical modular network organization  <cit> .
according to this view, the network consists of a hierarchy of nested topological modules of increasing size and decreasing coherence.
in other words, small coherent modules combine to form larger and less coherent modules in a hierarchical fashion.
at reasonable levels of module resolution, the modules consist mainly of sparsely connected but highly clustered nodes .
the modules are linked together through a small number of highly connected nodes with a low clustering coefficient , often referred to as hubs.
in the case of coexpression networks, these hubs represent genes that are linked to different expression modules depending on the experimental conditions.
a few papers  <cit>  have cast doubt on the ubiquity of power-law degree distributions in biological networks, claiming that some of the supposed power-laws actually turn out to be closer to exponentials when rigorously analyzed.
indeed, the degree distribution of the enigma co-differential expression network appears to be exponentially distributed , at least for lower k. for higher k, the picture is different.
relative to the distribution obtained for lower degrees, the most highly connected nodes  seem to be underconnected.
this observation is exactly the opposite of what would be expected for a power-law  degree distribution , indicating that the coexpression hubs are not nearly as central in the network as would be expected in a scale-free network.
however, from the plots of the clustering coefficient c versus the degree k , it is apparent that the highly connected nodes still possess hub-like characteristics: they generally have a lower clustering coefficient and are assigned to multiple modules.
thus, highly connected nodes act more as local hubs that hold together a few modules.
these hubs, by virtue of their polytomous expression behavior, may represent genes that function at the interface of several processes.
an example of genes that probably interface between the cell cycle, mating pheromone response and cell wall biosynthesis is given below.
overall,  <dig> genes are linked to  <dig> or more modules and  <dig> are linked to  <dig> or more modules, indicative of extensive crosstalk at the transcriptional level.
comparison between enigma, samba and isa
rigorously comparing the performance of clustering algorithms on real data is extremely difficult, because of the lack of adequate gold standards and the subjectivity of the available external performance criteria  <cit> .
therefore, we limit ourselves to a more qualitative comparison of the enigma, samba and isa modules obtained on the rosetta dataset.
samba was run with default parameter settings, for isa we used the bicat implementation  <cit>  with parameters tg =  <dig> , tc =  <dig>  and  <dig>  starting points .
the isa biclusters were pruned by merging biclusters that showed more than 80% overlap.
the isa and samba biclusters were put through the enigma postprocessing pipeline to functionally annotate them and to screen them for tf binding enrichment.
samba identified  <dig> modules containing  <dig>  genes and  <dig> conditions.
<dig> modules were enriched in one or more go biological process categories,  <dig> modules were enriched in binding sites for one or more tfs, and  <dig> modules showed both go and tf binding enrichment.
isa identified  <dig> modules containing  <dig>  genes and  <dig> conditions.
eighty-one modules were enriched in one or more go biological process categories,  <dig> modules were enriched in binding sites for one or more tfs, and  <dig> modules showed both go and tf binding enrichment.
these numbers are not directly comparable between methods, because of the differing degrees of overlap  between modules in the three formalisms.
samba generates a lot of biclusters with largely overlapping gene content , whereas the gene overlap between the enigma modules and especially the pruned isa modules is more limited.
for instance, samba identified  <dig> modules enriched in conjugation-related genes, containing a total of  <dig> genes annotated to 'conjugation' in go .
in contrast, enigma and isa identified fewer conjugation modules , but containing similar amounts of known conjugation genes .
instead of comparing general properties such as the overall coverage of genes and conditions by biclusters, the proportion of go-enriched modules or the average specificity  of the enriched modules, we focused our comparison on the biological processes that were mainly targeted by hughes et al  <cit>  , namely mating , ergosterol biosynthesis, cell wall biogenesis, oxidative phosphorylation and iron ion homeostasis.
all three formalisms uncover modules that are highly enriched for these processes.
we used two criteria to assess the module representation of a given go class a, namely the overall recall, or proportion of genes annotated to a found across all modules enriched for a, and the top module precision, or the proportion of genes in the most significantly enriched module that belong to a. samba generally detects slightly more true positive genes than enigma , but at the expense of a lower top module precision and a higher amount of modules .
isa generally features a lower recall than samba and enigma, but frequently exhibits better top modules in terms of precision.
in short, the main distinction between the formalisms seems to be a difference in balance between precision and recall.
moreover, the interpretation of the criteria defined above is not always straightforward.
for instance, a lower top module precision is not always caused by a lack of functional coherence, but may be caused by the presence of genes involved in closely related processes.
if we look at the overlap between the gene sets identified by the three methods , it becomes clear that all three formalisms add extra information to the global picture.
for all  <dig> processes in table  <dig>  a sizeable core of genes is identified by all three methods, but the different methods also have substantial idiosyncrasies.
for instance, only  <dig> out of a total of  <dig> identified conjugation-related genes are found by all three formalisms.
eleven genes are found by enigma and samba but not isa, two genes are found by enigma and isa but not samba, and four are shared between samba and isa but not enigma.
five genes are enigma-specific,  <dig> are samba-specific and strikingly,  <dig> are isa-specific, although isa identifies the least number of conjugation genes in total and has the 'worst' top module.
this illustrates that different methods have different focuses and biases, and that integration of the results of different analysis methods often leads to a better global picture.
comparison of enigma, samba and isa results for selected biological processes targeted by hughes et al.
<cit> .
the three middle columns give the number of modules enriched for the go class in the first column, the total number of genes annotated to that go class in these modules  and the corresponding recall .
the three last columns contain the enrichment p-value of the top module, the number of true positives  and the proportion of genes in the top module annotated to the respective go class .
pheromone response modules
in order to further assess the capacity of enigma to discover biologically relevant connections between genes and processes, we took a closer look at the mating-related enigma modules.
the rosetta compendium contains expression data on at least  <dig> mating-related perturbations, and consequently the mating pheromone response system is well resolved in the enigma network.
several mating-related modules were uncovered .
module  <dig> is most strongly related to mating .
twenty-three of the  <dig> genes in this module are annotated to the go category 'conjugation' .
four tfs exhibit binding enrichment in module 28: ste <dig>  dig <dig>  mcm <dig> and tec <dig>  all of these function in the regulation of the mating pheromone response .
two regulators show significant coexpression links with the module: ste <dig>  an important regulator of the mating response  and tec <dig>  a transcription factor involved in the regulation of haploid invasive and diploid pseudohyphal growth.
the mating and invasive/pseudohyphal growth signaling pathways share many of the same components, and tec <dig> is believed to mediate an invasive growth response upon low levels of pheromone signaling  <cit> .
whereas ste <dig> appears to be the main regulator for module  <dig>  tec <dig> is mainly coexpressed with genes that are shared between module  <dig> and modules  <dig>   <dig> or  <dig>  modules  <dig> and  <dig> are smaller pheromone response-related modules .
both modules contain tec <dig> as a member gene, suggesting that these modules might be more related to pseudohyphal growth than to the conjugation process.
module  <dig> exhibits a more complicated substructure, with five major patterns  in the condition dimension and five in the gene dimension .
most of the known mating-related genes in module  <dig> reside in the gene leafs e and f. several genes in these leafs overlap with the mating module  <dig>  in contrast, most of the genes in the leafs b and c overlap with module  <dig>  a module enriched in cell wall biogenesis genes .
nevertheless, most of these genes contain binding sites for ste <dig> and dig <dig> and some for tec <dig>  justifying their presence in a pheromone response-related module.
while the genes in leaf c appear to be genuinely related to cell wall biogenesis, none of the genes in leaf b is annotated as such.
compared to the genes in leaf c, the genes in leaf b show a distinctive subpattern in condition leaf  <dig>  which mainly contains perturbations that affect the cell cycle, dna maintenance and dna repair.
interestingly, the genes in leaf b also distinguish themselves from the ones in leaf c  by the presence of tf binding sites for swi <dig> and swi <dig>  which together form the sbf complex that regulates transcription at the g1/s transition  <cit> .
additionally, the genes in leaf b show strong coexpression links with the cyclins cln <dig> and cln <dig>  both swi <dig> and swi <dig> are potential substrates of the protein kinase cdc <dig>  which is activated by cln <dig> and cln <dig>  <cit> .
together, these data suggest that the genes in leaf b function at the interface of cell wall biogenesis, the g1/s transition and mating/filamentous growth.
such a link makes sense since upon activation of the pheromone signaling pathway, the yeast cell cycle is arrested in g <dig> and extensive cell wall rearrangements take place  <cit> .
together with the genes in leafs b and f, most genes in leaf e are strongly repressed under the pheromone response-related perturbations in condition leaf  <dig>  unlike leafs b and f, only a few genes in leaf e  feature bona fide ste <dig> or tec <dig> binding sites.
however, the expression of the other genes in leaf e  is specifically and strongly downregulated upon haploid tec <dig> deletion , suggesting that these genes are somehow transcriptionally regulated by tec <dig>  further investigation made apparent that several of the genes in leaf e are flanked by or overlapping with an antisense ty <dig> retrotransposon long terminal repeat  on the 3' side  or the 5' side .
the presence of these ty elements is highly relevant, since tec <dig> was originally described as a gene required together with ste <dig> for full ty <dig> expression  <cit> .
three of these genes  were found to be directly or indirectly associated with tec <dig> in a previous study in which the rosetta compendium was analyzed using a bayesian network framework  <cit> .
a peculiar member gene of leaf e is hap <dig>  a transcription factor involved in the regulation of respiratory metabolism in response to levels of heme and oxygen.
interestingly, hap <dig> also contains a 3' ty <dig> insertion in the yeast strain used by hughes et al   <cit> , which helps explain its puzzling presence in a pheromone response module and strengthens our belief that the ty <dig> elements are responsible for the link between leaf e genes and mating genes.
the coexpression of ndj <dig> with tec <dig> can be directly explained by the presence of a 5' ty <dig> ltr in antisense direction .
for gas <dig>  ylr334c, hap <dig> and yol106w, the situation is different given the 3' location of the flanking ty <dig> ltrs.
tec <dig> and ste <dig> activation of these ty <dig> elements could in theory cause the production of antisense transcripts of these loci.
since the probes spotted on the microarray used by hughes et al  <cit>  contained both strands of the gene sequences, such antisense transcripts might be responsible for the observed coexpression pattern.
we did not test the antisense hypothesis; the analysis we present here is merely intended as a use case to show that enigma can generate hypotheses that can be tested in the lab.
we did however briefly investigate whether the ty1-associated genes  could be functionally related to the mating process.
only two genes in leaf e  are known to be involved in mating.
neither of them is flanked by a ty <dig> ltr.
one gene overlapping with an antisense ty <dig> ltr, yol106w, was previously reported to elicit a mating-related phenotype upon deletion  <cit> .
we performed mating experiments, halo assays and growth assays  for two other 3' ty1-associated genes in leaf e, namely ylr334c  and gas <dig> , in addition to a wild type  strain and sst2Δ, a mutant that is supersensitive to mating factor-induced g1-arrest.
the ylr334cΔ deletion strain did not yield an interesting phenotype in any of the assays.
the gas2Δ deletion strain exhibited an interesting phenotype in the halo assay, characterized by extensive colony formation inside the halo , which indicates that deletion of gas <dig> somehow facilitates the recovery from α-factor induced growth arrest.
in the mating and growth assays, we did not observe any effect of gas <dig> deletion on the mating ability .
gas <dig> is homologous to gas <dig>  which encodes a  <dig> -β-glucanosyltransferase required for cell wall assembly.
in a recent study, gas <dig> was found to be involved in spore wall assembly  <cit> .
ectopic expression of gas <dig> under control of the gas <dig> promoter was found to complement the gas1Δ phenotype only partially, and only at ph =  <dig>   <cit> .
it is therefore unlikely that gas <dig> directly functions in regular cell wall assembly or maintenance.
in one hypothetical scenario, antisense transcripts of gas <dig>  produced under control of tec <dig>  might interfere with the expression of its homolog gas <dig> and hence indirectly with the formation and maintenance of the cell wall.
an altered cell wall morphology might influence the efficiency with which α-factor is inactivated, which could explain the observed gas2Δ phenotype.
obviously, this is only a hypothesis and much more detailed experimentation is needed to unravel if and how gas <dig> is linked to the pheromone response pathway.
this is however outside the scope of the present study.
implementation
enigma is implemented as a command-line java application that is open-source and freely available  from  <cit> .
enigma can be used for any organism for which there is sufficient gene expression data available.
the only organism-specific part of the enigma algorithm is the functional annotation module, which is based on the bingo tool  <cit> .
enigma can be used out-of-the-box for  <dig> organisms, including yeasts, invertebrates, plants and mammals .
furthermore, enigma allows the use of custom go annotation files and go consortium files to accommodate other organisms.
CONCLUSIONS
we have developed a novel method, called enigma, to analyze perturbational microarray data.
one of the innovations of our methodology is the use of a combinatorial statistic that is capable of detecting significant partial coexpression relationships between genes.
in this respect, our method can be considered similar in purpose to biclustering methods, although enigma assesses coexpression links between individual genes rather than expression coherence in a group of genes under a group of conditions.
our method produces both a detailed network of significant pair-wise coexpression links and a high-level representation of the modularity in the expression network.
tests on artificial data have shown that enigma outperforms other methods, although enigma wins from samba on points rather than by knockout.
similar near-draws with samba were reported earlier for cmonkey  <cit>  and bimax  <cit> .
this indicates that the clustering field has matured to a point at which it becomes increasingly difficult to easily improve on the performance of existing methods.
however, enigma does have some specific advantages.
first, in contrast to other discretization-based methods such as samba, enigma discretizes the expression data based on differential expression p-values.
second, enigma efficiently retrieved the correct number of modules from artificial datasets and actively avoids generating redundant modules, which greatly improves the interpretability of the results.
third, enigma's clustering parameters are automatically optimized or can be set on relatively objective grounds.
a fourth advantage that is more obvious on real data is the use of enigma's expression module concept for biological discovery.
in contrast to the coherent biclusters generated by most methods, an enigma expression module may contain distinctive subpatterns.
from our analysis of the rosetta data, it became apparent that these subpatterns frequently represent more tightly coregulated gene clusters involved in biological processes related to a common functional theme.
in our view, the grouping of such different but statistically and functionally connected patterns in one module aids greatly in the biological interpretation of the data and in the assessment of crosstalk between biological processes.
the interpretation of a module's substructure is further facilitated by the integration of other data types.
this is illustrated in our analysis of module  <dig>  a pheromone response module which shows links to the cell cycle, cell wall biosynthesis and ty <dig> ltr-associated genes.
although numerous approaches have already been used to mine the rosetta compendium  <cit> , enigma offers yet another perspective on the data.
this mainly illustrates that the ideal clustering method does not exist  <cit> , and that no single approach can extract all the information hidden in large compendium datasets.
the elucidation of the regulatory networks governing the many different aspects of cellular function will therefore not only require the integration of different types of data, but also the integrated use of several complementary methods to analyze these data.
we believe that enigma constitutes a valuable addition to the existing repertoire of analysis methods.

 over the last decade, the availability of fully sequenced genomes and the development of high-throughput technologies such as dna microarray-based transcript profiling have fuelled an exponential increase in the volume of functional genomics data. this has led to a renewed interest in the study of molecular biology at the system level. the central paradigm in systems theory is that one can learn about a system by perturbing it and measuring the response. this principle also applies to biological systems. since mrna levels can nowadays easily be measured on a genome-wide scale, expression profiling has become a first method of choice for assessing the molecular response to experimental perturbation. considerable efforts are put into creating compendia of expression profiles under genetic, chemical or environmental perturbations or in different tissues. such data compendia basically constitute a series of snapshots of expression states under a variety of conditions, and they contain a wealth of information concerning the underlying transcriptional network structure of an organism. however, devising methods to efficiently and reliably extract that information is still a challenging task. clustering of gene expression data allows the inference of functional correlations between genes through what was dubbed the'guilt-by-association' principle. a classical clustering process consists of two steps. first, a matrix of distances between expression profiles is calculated using a distance or similarity measure, such as pearson's centered correlation coefficient. based on this distance matrix, the actual clustering algorithm, for instance average linkage hierarchical clustering, groups similar profiles together. traditional clustering methods are well suited for analyzing time-series expression data, but they fall short when applied to perturbational data, because the underlying similarity measures, such as pcc, primarily capture global correlation tendencies. however, in compendia of perturbed expression profiles, genes do not necessarily show similar behavior under all experimental conditions: they may be coexpressed under some conditions and follow different expression regimes under other conditions. one of the consequences is that genes may be coexpressed with multiple expression modules depending on the conditions, or in other words, expression modules may overlap. these observations stimulated the development of alternative clustering strategies. the process of detecting subsets of genes that exhibit similar expression behavior across a subset of conditions is known as biclustering. several biclustering strategies exist today, each using its own heuristic approach to tackle this complex problem. some biclustering methods use a greedy iterative search strategy to uncover biclusters, progressively subdividing, or adding and removing rows and columns from the biclusters obtained in a previous iteration in order to maximize a local score function. others use linear algebra or adopt a graph-theoretic approach to biclustering. yet other methods identify biclusters by proposing a statistical model and estimating the distribution parameters that minimize a certain model fit criterion. a feature that most biclustering methods share is that they do not explicitly define similarity measures on the global space of expression profiles, but rely on the emergent properties of groups of genes and conditions in order to identify significant subpatterns in the data. evidently, a wide variety of biclustering algorithms exist, each of them having their own strengths and weaknesses. for example, some of these methods are intrinsically less suited to find overlap between biclusters because they mask previously found biclusters with random noise, or because they partition the data. others require extensive parameter tweaking, require the user to specify the desired number of biclusters in advance or generate very small or large biclusters or highly redundant biclusters. some have no publicly available implementation or are rather cumbersome to use, and most of them, notable exceptions being samba/expander, genomica and cmonkey, do not integrate or overlay other types of biological data, hampering their use as biological discovery tools. also, to our knowledge, none of the existing biclustering methods uses the variational information in replicated expression experiments. this information is routinely and successfully used to detect genes that are differentially expressed under a given perturbation. the main reason why biclustering methods do not use differential expression information is that they do not specifically focus on the analysis of perturbational data. discretization-based biclustering methods such as samba and bimax could probably easily be modified to assess up- and downregulation of gene expression based on p-values for differential expression. in their current implementation, however, these methods use rather arbitrary log-ratio or percentage cutoffs for this purpose. in this study, we present a novel method, called enigma, that addresses some of these issues. our goal was to build a method that: leverages differential expression analysis results to extract co-differential expression networks and expression modules from perturbational gene expression data, is able to detect significant partial coexpression relationships between genes and overlap between modules, depends on parameters that can be automatically optimized or set on reasonably objective grounds. produces a realistic amount of modules, and visually integrates the expression modules with other data types such as gene ontology information, transcription factor binding data, protein and genetic interactions, in order to facilitate the biological interpretation of the results. below, we outline the enigma algorithm, test our methodology on artificial expression data and compare its performance to other methods. we also apply enigma to a perturbational microarray compendium for budding yeast in order to assess its potential to generate testable hypotheses on real biological data.  algorithm a global overview of the methodology is given in figure briefly, enigma takes as input a set of perturbational expression data, externally calculated p-values for differential expression and other data types if available. enigma uses a novel combinatorial statistic to assess which pairs of genes are significantly co-differentially expressed. the resulting coexpression p-values are corrected for multiple testing and translated to edges in a coexpression network, which is clustered into expression modules using a graph-based clustering algorithm inspired on the mcode algorithm. the clustering procedure depends on two parameters that control the density of individual modules and the overlap between modules. the main reason why we chose a two-tier clustering approach is that it allows simulated annealing-based optimization of the clustering parameters to obtain optimal coverage of the coexpression network, in terms of module overlap and redundancy. the graph clustering method we use is very fast, which allows the parameters to be optimized in a reasonable amount of time. in the postprocessing phase, enigma determines relevant condition sets for each module, visualizes their substructure and overlap with other modules, screens the modules for enriched go categories, suggests potential regulators for the modules based on regulator-module coexpression links and enrichment of tf binding sites, and overlays protein and genetic interaction data. combinatorial statistic consider the expression profiles of two genes a and b under n perturbations. each gene is represented by a profile of n fields. the gene expression values are discretized into three categories based on their differential expression p-value. if the gene is significantly upregulated in a given experiment, the corresponding field is labeled blue. experiments in which the gene is significantly downregulated are similarly labeled yellow, and the remaining fields are labeled black. let us now assume that the profiles of a and b contain ax and bx blue fields respectively, as well as ay and by yellow fields, and that they have x blue and y yellow fields in common. we want to assess whether this overlap is statistically significant. if the response of the genes a and b to the perturbations were uncorrelated, the blue and yellow fields would be independently distributed on both profiles. under this hypothesis, the probability that the profiles overlap on exactly x blue and y yellow positions is given by the following recursive formula: p=−∑x′=xmin⁡∑y′=ymin⁡≠p the probability of observing an overlap of at least x blue and y yellow fields by chance is then expressed by the cumulative distribution: pc=∑x′=xmin⁡∑y′=ymin⁡p equation can be understood by assuming that profile a is given, and that we randomly distribute bx blue and by yellow positions on profile b. the denominator of the first term then represents the total number of possible profiles b. the numerator represents the combinations in which x blue and y yellow matching positions are selected, and the residual positions are chosen at random. however, in this manner, a number of combinations are selected while having more than exactly x blue and/or y yellow matching positions. moreover, combinations with x' > x blue and/or y' > y yellow matching positions are counted c·c times, hence the second term. although the probabilistic question formulated above can be cast in terms of contingency tables, the hypothesis tested by our statistic is different from that tested by standard contingency table analysis methods such as the χ test. for example, situations in which a large amount of blue fields in profile a are perfectly mapped onto the black fields in profile b would yield a significant χ p-value, whereas they would not yield a significant p-value using equation our statistic only considers mappings of up- and down-regulation of the expression of a gene to up- or down-regulation of another gene to be meaningful for assessing coregulation, a premise which is motivated by the perturbational nature of the data we aim to analyze. black fields are considered less informative from the perspective of coregulation. multiple testing correction of coexpression p-values in our probabilistic setup, each comparison of two profiles can be considered an individual test. for n genes, n/ tests are performed to fish for co-differential expression relationships. consequently, the obtained p-values have to be adjusted in order to control the type i error rate. the raw p-values are corrected for multiple testing with the benjamini & hochberg procedure, which controls the false discovery rate. graph-based clustering the set of significant coexpression relationships at a certain fdr threshold is translated to a network, with nodes and edges representing genes and significant coexpression relationships, respectively. enigma identifies coexpression modules from this network using a graph clustering technique inspired by the mcode algorithm. to identify potential module seeds, all nodes v are weighted based on the density of the highest k-core of the node neighborhood nv, denoted as the kmax-core of v. analogous to bader and hogue, the core-clustering coefficient ccore,v is defined as the density of the kmax-core of v, and the weight wv = ccore,v·kmax,v. the kmax-core of the node with the highest weight is taken as the first module seed. this module seed then grows by accreting nodes on which it exerts a pull above a certain threshold ν the pull of a module with seed s on a node v outside the module is defined as |nv ∩ s|/|s|. the next module is then initiated by taking the kmax-core of the node with the highest weight in the remaining graph. an additional constraint is set by requiring that the overlap between the new seed s and any existing module m does not exceed ν1·min. while the threshold ν controls the size and density of individual modules, ν controls the spacing or overlap between modules. both parameters are optimized automatically. clustering parameter optimization in order to optimize the clustering parameters, the quality of the clustering for a given is assessed by comparing the known input coexpression network with the output coexpression network inferred by the modules. the latter is constructed by translating the modules to fully connected components in the output network. if we consider true/false positives to be coexpression edges inferred by the clustering that are present/absent in the input coexpression network, and false negatives as edges present in the input network that are not inferred by the clustering, we can define the precision p' = tp/ and the recall r' = tp/ of the clustering result. enigma uses the f'-measure, i.e. the harmonic mean of recall and precision, f' = 2p'r'/, as a measure for the quality of the clustering. we use the notation p', r', f' instead of the more commonly used p, r, f in order to distinguish between two different flavours of the f-measure used in this study for different purposes. in contrast to the regular f-measure, the f'-measure penalizes overpredicted edges in order to avoid unnecessary overlap between the expression modules: an edge that is inferred multiple times from the clustering, because the genes a and b belongs to the intersection of multiple modules, is counted as tp and x - fp. this is equivalent to drawing x edges between the genes a and b in the output coexpression network. since there is only one edge in the input network, the x - remaining edges can be considered false. this penalization strategy has the intuitively pleasing property of not affecting the recall, but lowering the precision of the clustering result when the amount of edges'explained' by multiple modules increases. the parameters ν and ν are now optimized by monte-carlo simulated annealing using f' as the optimization criterion. starting from an random initial guess for the parameters, random steps are taken in parameter space. a step is accepted if rand <eΔf'/t with rand a random number drawn uniformly from the interval, Δf' the change in f'-measure and t the simulated annealing parameter or'temperature', which gradually decreases during the course of the optimization according to an exponential scheme ti = rcti- with rc the cooling rate. enigma uses a two-stage mcsa procedure. in the first stage, a rough mcsa search of the clustering parameter space is performed in order to identify the most interesting parameter region. in the second stage, a finer mcsa search is performed starting from the optimum obtained in the first stage. at the end of each stage, an additional gradient descent is performed toward the nearest local optimum of f'. by default, enigma performs mcsa runs, starting from randomly chosen. the convergence of the solutions of multiple runs can be used as a check on the adequacy of the mcsa parameter settings. postprocessing of modules for each gene module, enigma determines a condition set by selecting those conditions that show enrichment of up- or downregulated genes in the module. thus, for a given module, the condition set contains the experimental conditions that elicit a significant and specific response in the module and, by consequence, have been most influential in shaping the module. the resulting'bicluster' does not necessarily have a uniform expression pattern over all genes, but may show subpatterns for some genes under certain conditions, possibly indicating involvement in other expression modules. these subpatterns are visualized by hierarchically clustering the module's expression data in both dimensions, using the cosine correlation coefficient as a similarity measure. the clustering tree can optionally be separated into leafs to make the subdivision more clear. although conditions that show differential patterning within one module might appear to be irrelevant for the module as a whole, they are important for at least part of the module and may provide insight into inter-module connections or further substructure within the module. in an attempt to provide the user with clues on how the expression modules are regulated, enigma searches for'regulators' that are significantly more connected to a module, through positive or negative coexpression edges, than expected at random. potential regulators are selected from a user-defined list or a user-defined set of go classes. when chromatin immunoprecipitation or tf motif data are available, enigma also screens the modules for enriched tf binding sites. the expression profiles of significantly coexpressed or binding regulators are visualized on top of the modules. significantly enriched go terms for both the gene and condition sets of the modules are determined using the bingo software, which is incorporated in enigma. finally, enigma visually maps the available protein interaction data and genetic interaction data on the modules. testing on artificial data generating artificial expression data to assess the performance of our method and compare enigma to other methods, we performed tests on artificial gene expression data. we generated two types of artificial expression data, namely expression data containing overlapping biclusters and expression data containing partially coexpressed genes but no biclusters. in both cases, we built expression datasets of genes by experiments. for each dataset, artificial background expression data were randomly sampled from a normal distribution with mean μ = and variance σ =. for the modular datasets, we implanted biclusters in this background, each encompassing between 1–5% of all genes and 10–50% of all conditions. bicluster sizes, member genes and conditions are chosen at random, with the restriction that at most 30% of the genes and 50% of the conditions overlap between any biclusters. except for a noise component, all genes in a bicluster share the same expression profile over the bicluster conditions. however, a bicluster can be partially overwritten by other biclusters. the bicluster profiles are sampled from a bimodal distribution consisting of normal modes with means μ = - and μ = and variances σ12=σ22=. the expression profiles of individual genes in a bicluster are noisified by adding normally distributed noise. the variances, bicluster size and overlap parameters are chosen so that the overall distribution of the simulated log ratio expression values mimicks the distribution of log ratio expression values in the rosetta compendium up to a scale factor. note that, apart from the distribution of expression ratios, the structure of these toy datasets does not necessarily bear any resemblance to real biological data. for the non-modular datasets, we implanted pairs of partially coexpressed genes in the background. the expression profiles are constructed as described above. the resulting expression value distribution again mimicks the rosetta distribution. unlike for real data, we used log ratio thresholds to discretize the expression values of the artificial datasets, because the generation of meaningful artificial differential expression p-values proved to merit further study in its own right. therefore, the artificial data cannot be used to assess the advantage of including variational information in enigma's discretization step. on the other hand, we can still compare the performance of enigma with other methods that do not use variational information. we used a log ratio threshold of for upregulation and - for downregulation, corresponding to the means of the distributions used to generate the bicluster profiles. in other words, half of the datapoints in the biclusters are presumed not to be significantly over- or underexpressed. performance of enigma on artificial data and comparison with other methods the performance of enigma on these toy datasets was compared with that of two commonly used similarity measures, namely pcc and the χ2-statistic, and two established biclustering methods, samba and isa. pcc was chosen as a representative of the global similarity measures used in traditional clustering algorithms, while we included the χ2-statistic because of its relation to the combinatorial statistic used by enigma. the selection of biclustering methods was based on the following criteria: the methods should be non-partitioning in nature, they should have the capacity to generate overlapping biclusters, a suitable implementation should be publicly available, and they should produce a reasonable amount of biclusters on the modular toy datasets. we used the version of samba incorporated in the expander package, and the implementation of isa available as part of the biclustering tool bicat, both with default parameter settings. the isa trajectories from randomly chosen starting points converge to a limited number of'fixed point' biclusters. to prune nearly identical modules, we merged isa biclusters that overlap for more than 80%. the clustering performance of all methods is only assessed in the gene dimension. standard internal criteria for partitional clustering performance, such as the silhouette width or dunn's index, cannot be used to assess the performance of algorithms that generate overlapping clusters. instead, we use the f-measure and introduce a derivative, the f'-measure, to compare the performance of different clustering methods on artificial datasets. in both cases, the coexpression network generated by a method is compared to the artificial input coexpression network in terms of true and false positive edges and false negative edges, from which the different flavors of the f-measure are calculated. the difference between the f-measure and the f'-measure is that the f-measure does not take into account the multiplicity of the inferred edges. in other words, the f'-measure penalizes overpredicted edges, whereas the f-measure does not. this entails that the f'-measure is more useful to compare methods that generate overlapping clusters, whereas the f-measure can be used more generally to compare methods that generate both overlapping or non-overlapping clusters or pair-wise coexpression networks. the performance of enigma is tested on two levels by assessing the overlap between the artificial input correlation network and the network of significant correlations obtained in the first step of the enigma algorithm; the modules inferred by enigma. the output networks for enigma-m and the biclustering methods samba and isa are obtained by converting the obtained modules/biclusters to fully connected network components. the χ network is constructed by translating significant χ correlation p-values between the discretized expression profiles to edges in the output network. we used the same discretization threshold and fdr level for the χ and enigma methods. the performance of pcc was measured for different thresholds. using the f-measure, enigma outperforms all other methods on the modular artificial data. the performance of enigma-m was consistently higher than the χ performance and the optimal pcc performance. the global similarity measure pcc appears to perform surprisingly well. however, the performance of pcc critically depends on the choice of the pcc threshold, and determining the optimal pcc threshold on real data is problematic. in contrast, enigma has the advantage of having an easily tunable significance threshold: the false discovery rate level. to illustrate this, we plotted the performance curve of enigma for different non-corrected p-value thresholds, on figure 2a and 2b. for all artificial datasets, the performance of enigma-n at fdr = is close to the optimum of this curve, indicating that fdr control at a reasonable level gives near-optimal performance. among the biclustering methods, the rather poor performance of the isa algorithm may seem somewhat surprising. prelić et al, using the same implementation of isa but other methods to generate artificial data and to assess biclustering performance, previously established that the performance of isa decreases with increasing overlap between biclusters. our results seem to confirm that isa is not the optimal method in case there is substantial overlap between modules. the performance of isa did not change significantly when using starting points instead of the default. the performance gain of enigma-m over samba is substantially smaller, and on two out of datasets, the performance of samba was slightly higher than that of enigma-m. a more tangible advantage of enigma over samba is that enigma nearly always recovered the correct number of modules, whereas samba consistently predicted more modules than there were in the input data. isa predicted only one extra module on average, but with a higher variance than enigma. in other words, samba and to a lesser extent isa produce more fragmented and/or more redundant modules. redundancy makes the module output much harder to interpret, but it is not taken into account by the standard f-measure. to quantify the effect of redundancy on the clustering quality, we compared samba, isa and enigma-m using the f'-measure. as in the calculation of the f'-measure used in the clustering optimization procedure, edges that are inferred by multiple modules are counted multiple times, but in the present case, multiply defined edges may also occur in the input network if they overlap between multiple artificial input modules. specifically, edges that are inferred by x output modules and y input modules are now counted as y tp and x - y fp in case x ≥ y, or x tp and y - x fn in case x <y. using the f' criterion, the performance of enigma-m is substantially higher than that of samba and isa. on non-modular artificial data, the performance of enigma-m and the biclustering methods samba and isa is very low. this is not surprising since there are no modules to be found in these datasets. in this respect, a particularly attractive feature of enigma is that it finds very few modules in the non-modular data, in contrast to isa and samba, which recover ± modules and ± modules, respectively. among the pair-wise methods, enigma-n invariably featured the highest performance, indicating that our combinatorial statistic detects partial coexpression relationships more efficiently than pcc and χ the fact that enigma efficiently uncovers coexpression relationships in non-modular data opens perspectives for the exploration of the less modular parts of expression datasets. real datasets typically contain a limited number of perturbation experiments that target a few specific processes. these processes can be expected to be rather well resolved in terms of their coexpression relationships, whereas other processes will probably give rise to more fragmented regions in the network. moreover, despite the success of the modularity concept in the analysis of expression data and systems biology in general, it is not inconceivable that transcriptional networks might also contain genuinely non-modular regions. testing on real data: the rosetta gene expression compendium p-value versus log-ratio based discretization although useful for testing and comparing methods, artificial datasets do not capture the complexity of real biological systems. consequently, good performance on artificial data does not guarantee good performance on real biological data. in order to assess the use of enigma for analyzing real data, we applied our methodology to the rosetta compendium of expression profiles, representing data on different experimental perturbations of s. cerevisiae. experiments on strains that were marked as aneuploid in the original dataset were left out, because they can give rise to artificial expression correlations between genes on the aneuploid chromosomes. the log-ratio expression data and differential expression p-values were downloaded in prenormalized and preprocessed form. genome-wide chip data for tfs were obtained from harbison et al. all genes that are bound with p < by a certain tf were considered reliable targets. protein and genetic interactions for s. cerevisiae were obtained from the biogrid database. using a differential expression p-value threshold of in the discretization step and an fdr threshold of for defining coexpression edges, enigma identified a network of significant positive coexpression links and negative coexpression links involving genes. the clustering parameters = were optimized by mcsa as described in the algorithm section. to assess the efficiency of the mcsa procedure, we performed an exhaustive screen of the parameter space to locate the global maximum of the f'-measure criterion. the mcsa procedure found back the global optimum with 100% efficiency. enigma discovered modules in the rosetta dataset encompassing genes and conditions. these numbers seem reasonable given that of the conditions included in the compendium contain less than five differentially expressed genes, which entails that they have a small chance of contributing to a module. given the low amount of informative conditions, it is not surprising that only a third of the s. cerevisiae genes can be included in modules. according to the go enrichment results, out of modules have a significant degree of functional coherence. fifty-four modules are enriched in targets of one or more tfs, and modules show enrichment of both go biological process categories and tf binding sites. together, 60% of the modules show enrichment of go categories and/or tf binding sites, indicating that our method is capable of identifying biologically relevant expression modules. in order to qualitatively assess the effect of using a differential expression p-value cutoff in the discretization step instead of a fold-change cutoff, we repeated the analysis using a |log ratio| discretization threshold of. the resulting coexpression network contains positive and negative links between genes. the clustered network contains modules encompassing genes. ninety-three modules exhibit go enrichment, exhibit tf binding enrichment and exhibit both. despite the significantly lower amount of connections in the log-ratio network, the number of functionally coherent modules and the number of clustered genes is roughly similar, and the optimized clustering parameters = are identical, indicating that the general structure of the network and its strongest modules are fairly well preserved. indeed, many highly functionally coherent modules are present in both networks. not incidentally, many of these modules are related to the processes that were targeted by hughes et al, which can be expected to show a pronounced expression response. however, modules that show less pronounced expression variations, for example the modules related to ribosome biosynthesis, are not recovered in the log-ratio network. this illustrates the main disadvantage of using a fixed log-ratio threshold: different processes show different amplitudes of expression change upon perturbation, which cannot be captured by a single threshold. one could argue that this can easily be remedied by standardizing the expression profiles to zero mean, unit variance before applying the threshold, as is done by some methods, e.g. samba. however, in the case of perturbational data, this manipulation runs the risk of effiectively breaking the connection to the reference condition, thereby distorting the meaning of up- and downregulation and introducing serious artifacts. topological characteristics of the enigma coexpression network since many cellular functions are carried out in a highly modular manner, most cellular networks, including protein interaction networks, metabolic networks and gene expression networks, are modular in nature. on the other hand, many cellular networks, including coexpression networks, have been claimed to exhibit a node degree distribution of the power-law type, p ~ k-γ, indicative of scale-free properties. the coexistence of modularity and a scale-free degree distribution can be explained by assuming a hierarchical modular network organization. according to this view, the network consists of a hierarchy of nested topological modules of increasing size and decreasing coherence. in other words, small coherent modules combine to form larger and less coherent modules in a hierarchical fashion. at reasonable levels of module resolution, the modules consist mainly of sparsely connected but highly clustered nodes. the modules are linked together through a small number of highly connected nodes with a low clustering coefficient, often referred to as hubs. in the case of coexpression networks, these hubs represent genes that are linked to different expression modules depending on the experimental conditions. a few papers have cast doubt on the ubiquity of power-law degree distributions in biological networks, claiming that some of the supposed power-laws actually turn out to be closer to exponentials when rigorously analyzed. indeed, the degree distribution of the enigma co-differential expression network appears to be exponentially distributed, at least for lower k. for higher k, the picture is different. relative to the distribution obtained for lower degrees, the most highly connected nodes seem to be underconnected. this observation is exactly the opposite of what would be expected for a power-law degree distribution, indicating that the coexpression hubs are not nearly as central in the network as would be expected in a scale-free network. however, from the plots of the clustering coefficient c versus the degree k, it is apparent that the highly connected nodes still possess hub-like characteristics: they generally have a lower clustering coefficient and are assigned to multiple modules. thus, highly connected nodes act more as local hubs that hold together a few modules. these hubs, by virtue of their polytomous expression behavior, may represent genes that function at the interface of several processes. an example of genes that probably interface between the cell cycle, mating pheromone response and cell wall biosynthesis is given below. overall, genes are linked to or more modules and are linked to or more modules, indicative of extensive crosstalk at the transcriptional level. comparison between enigma, samba and isa rigorously comparing the performance of clustering algorithms on real data is extremely difficult, because of the lack of adequate gold standards and the subjectivity of the available external performance criteria. therefore, we limit ourselves to a more qualitative comparison of the enigma, samba and isa modules obtained on the rosetta dataset. samba was run with default parameter settings, for isa we used the bicat implementation with parameters tg =, tc = and starting points. the isa biclusters were pruned by merging biclusters that showed more than 80% overlap. the isa and samba biclusters were put through the enigma postprocessing pipeline to functionally annotate them and to screen them for tf binding enrichment. samba identified modules containing genes and conditions. modules were enriched in one or more go biological process categories, modules were enriched in binding sites for one or more tfs, and modules showed both go and tf binding enrichment. isa identified modules containing genes and conditions. eighty-one modules were enriched in one or more go biological process categories, modules were enriched in binding sites for one or more tfs, and modules showed both go and tf binding enrichment. these numbers are not directly comparable between methods, because of the differing degrees of overlap between modules in the three formalisms. samba generates a lot of biclusters with largely overlapping gene content, whereas the gene overlap between the enigma modules and especially the pruned isa modules is more limited. for instance, samba identified modules enriched in conjugation-related genes, containing a total of genes annotated to'conjugation' in go. in contrast, enigma and isa identified fewer conjugation modules, but containing similar amounts of known conjugation genes. instead of comparing general properties such as the overall coverage of genes and conditions by biclusters, the proportion of go-enriched modules or the average specificity of the enriched modules, we focused our comparison on the biological processes that were mainly targeted by hughes et al, namely mating, ergosterol biosynthesis, cell wall biogenesis, oxidative phosphorylation and iron ion homeostasis. all three formalisms uncover modules that are highly enriched for these processes. we used two criteria to assess the module representation of a given go class a, namely the overall recall, or proportion of genes annotated to a found across all modules enriched for a, and the top module precision, or the proportion of genes in the most significantly enriched module that belong to a. samba generally detects slightly more true positive genes than enigma, but at the expense of a lower top module precision and a higher amount of modules. isa generally features a lower recall than samba and enigma, but frequently exhibits better top modules in terms of precision. in short, the main distinction between the formalisms seems to be a difference in balance between precision and recall. moreover, the interpretation of the criteria defined above is not always straightforward. for instance, a lower top module precision is not always caused by a lack of functional coherence, but may be caused by the presence of genes involved in closely related processes. if we look at the overlap between the gene sets identified by the three methods, it becomes clear that all three formalisms add extra information to the global picture. for all processes in table a sizeable core of genes is identified by all three methods, but the different methods also have substantial idiosyncrasies. for instance, only out of a total of identified conjugation-related genes are found by all three formalisms. eleven genes are found by enigma and samba but not isa, two genes are found by enigma and isa but not samba, and four are shared between samba and isa but not enigma. five genes are enigma-specific, are samba-specific and strikingly, are isa-specific, although isa identifies the least number of conjugation genes in total and has the'worst' top module. this illustrates that different methods have different focuses and biases, and that integration of the results of different analysis methods often leads to a better global picture. comparison of enigma, samba and isa results for selected biological processes targeted by hughes et al.. the three middle columns give the number of modules enriched for the go class in the first column, the total number of genes annotated to that go class in these modules and the corresponding recall. the three last columns contain the enrichment p-value of the top module, the number of true positives and the proportion of genes in the top module annotated to the respective go class. pheromone response modules in order to further assess the capacity of enigma to discover biologically relevant connections between genes and processes, we took a closer look at the mating-related enigma modules. the rosetta compendium contains expression data on at least mating-related perturbations, and consequently the mating pheromone response system is well resolved in the enigma network. several mating-related modules were uncovered. module is most strongly related to mating. twenty-three of the genes in this module are annotated to the go category'conjugation'. four tfs exhibit binding enrichment in module 28: ste dig mcm and tec all of these function in the regulation of the mating pheromone response. two regulators show significant coexpression links with the module: ste an important regulator of the mating response and tec a transcription factor involved in the regulation of haploid invasive and diploid pseudohyphal growth. the mating and invasive/pseudohyphal growth signaling pathways share many of the same components, and tec is believed to mediate an invasive growth response upon low levels of pheromone signaling. whereas ste appears to be the main regulator for module tec is mainly coexpressed with genes that are shared between module and modules or modules and are smaller pheromone response-related modules. both modules contain tec as a member gene, suggesting that these modules might be more related to pseudohyphal growth than to the conjugation process. module exhibits a more complicated substructure, with five major patterns in the condition dimension and five in the gene dimension. most of the known mating-related genes in module reside in the gene leafs e and f. several genes in these leafs overlap with the mating module in contrast, most of the genes in the leafs b and c overlap with module a module enriched in cell wall biogenesis genes. nevertheless, most of these genes contain binding sites for ste and dig and some for tec justifying their presence in a pheromone response-related module. while the genes in leaf c appear to be genuinely related to cell wall biogenesis, none of the genes in leaf b is annotated as such. compared to the genes in leaf c, the genes in leaf b show a distinctive subpattern in condition leaf which mainly contains perturbations that affect the cell cycle, dna maintenance and dna repair. interestingly, the genes in leaf b also distinguish themselves from the ones in leaf c by the presence of tf binding sites for swi and swi which together form the sbf complex that regulates transcription at the g1/s transition. additionally, the genes in leaf b show strong coexpression links with the cyclins cln and cln both swi and swi are potential substrates of the protein kinase cdc which is activated by cln and cln. together, these data suggest that the genes in leaf b function at the interface of cell wall biogenesis, the g1/s transition and mating/filamentous growth. such a link makes sense since upon activation of the pheromone signaling pathway, the yeast cell cycle is arrested in g and extensive cell wall rearrangements take place. together with the genes in leafs b and f, most genes in leaf e are strongly repressed under the pheromone response-related perturbations in condition leaf unlike leafs b and f, only a few genes in leaf e feature bona fide ste or tec binding sites. however, the expression of the other genes in leaf e is specifically and strongly downregulated upon haploid tec deletion, suggesting that these genes are somehow transcriptionally regulated by tec further investigation made apparent that several of the genes in leaf e are flanked by or overlapping with an antisense ty retrotransposon long terminal repeat on the 3' side or the 5' side. the presence of these ty elements is highly relevant, since tec was originally described as a gene required together with ste for full ty expression. three of these genes were found to be directly or indirectly associated with tec in a previous study in which the rosetta compendium was analyzed using a bayesian network framework. a peculiar member gene of leaf e is hap a transcription factor involved in the regulation of respiratory metabolism in response to levels of heme and oxygen. interestingly, hap also contains a 3' ty insertion in the yeast strain used by hughes et al, which helps explain its puzzling presence in a pheromone response module and strengthens our belief that the ty elements are responsible for the link between leaf e genes and mating genes. the coexpression of ndj with tec can be directly explained by the presence of a 5' ty ltr in antisense direction. for gas ylr334c, hap and yol106w, the situation is different given the 3' location of the flanking ty ltrs. tec and ste activation of these ty elements could in theory cause the production of antisense transcripts of these loci. since the probes spotted on the microarray used by hughes et al contained both strands of the gene sequences, such antisense transcripts might be responsible for the observed coexpression pattern. we did not test the antisense hypothesis; the analysis we present here is merely intended as a use case to show that enigma can generate hypotheses that can be tested in the lab. we did however briefly investigate whether the ty1-associated genes could be functionally related to the mating process. only two genes in leaf e are known to be involved in mating. neither of them is flanked by a ty ltr. one gene overlapping with an antisense ty ltr, yol106w, was previously reported to elicit a mating-related phenotype upon deletion. we performed mating experiments, halo assays and growth assays for two other 3' ty1-associated genes in leaf e, namely ylr334c and gas, in addition to a wild type strain and sst2Δ, a mutant that is supersensitive to mating factor-induced g1-arrest. the ylr334cΔ deletion strain did not yield an interesting phenotype in any of the assays. the gas2Δ deletion strain exhibited an interesting phenotype in the halo assay, characterized by extensive colony formation inside the halo, which indicates that deletion of gas somehow facilitates the recovery from α-factor induced growth arrest. in the mating and growth assays, we did not observe any effect of gas deletion on the mating ability. gas is homologous to gas which encodes a -β-glucanosyltransferase required for cell wall assembly. in a recent study, gas was found to be involved in spore wall assembly. ectopic expression of gas under control of the gas promoter was found to complement the gas1Δ phenotype only partially, and only at ph =. it is therefore unlikely that gas directly functions in regular cell wall assembly or maintenance. in one hypothetical scenario, antisense transcripts of gas produced under control of tec might interfere with the expression of its homolog gas and hence indirectly with the formation and maintenance of the cell wall. an altered cell wall morphology might influence the efficiency with which α-factor is inactivated, which could explain the observed gas2Δ phenotype. obviously, this is only a hypothesis and much more detailed experimentation is needed to unravel if and how gas is linked to the pheromone response pathway. this is however outside the scope of the present study. implementation enigma is implemented as a command-line java application that is open-source and freely available from. enigma can be used for any organism for which there is sufficient gene expression data available. the only organism-specific part of the enigma algorithm is the functional annotation module, which is based on the bingo tool. enigma can be used out-of-the-box for organisms, including yeasts, invertebrates, plants and mammals. furthermore, enigma allows the use of custom go annotation files and go consortium files to accommodate other organisms.  we have developed a novel method, called enigma, to analyze perturbational microarray data. one of the innovations of our methodology is the use of a combinatorial statistic that is capable of detecting significant partial coexpression relationships between genes. in this respect, our method can be considered similar in purpose to biclustering methods, although enigma assesses coexpression links between individual genes rather than expression coherence in a group of genes under a group of conditions. our method produces both a detailed network of significant pair-wise coexpression links and a high-level representation of the modularity in the expression network. tests on artificial data have shown that enigma outperforms other methods, although enigma wins from samba on points rather than by knockout. similar near-draws with samba were reported earlier for cmonkey and bimax. this indicates that the clustering field has matured to a point at which it becomes increasingly difficult to easily improve on the performance of existing methods. however, enigma does have some specific advantages. first, in contrast to other discretization-based methods such as samba, enigma discretizes the expression data based on differential expression p-values. second, enigma efficiently retrieved the correct number of modules from artificial datasets and actively avoids generating redundant modules, which greatly improves the interpretability of the results. third, enigma's clustering parameters are automatically optimized or can be set on relatively objective grounds. a fourth advantage that is more obvious on real data is the use of enigma's expression module concept for biological discovery. in contrast to the coherent biclusters generated by most methods, an enigma expression module may contain distinctive subpatterns. from our analysis of the rosetta data, it became apparent that these subpatterns frequently represent more tightly coregulated gene clusters involved in biological processes related to a common functional theme. in our view, the grouping of such different but statistically and functionally connected patterns in one module aids greatly in the biological interpretation of the data and in the assessment of crosstalk between biological processes. the interpretation of a module's substructure is further facilitated by the integration of other data types. this is illustrated in our analysis of module a pheromone response module which shows links to the cell cycle, cell wall biosynthesis and ty ltr-associated genes. although numerous approaches have already been used to mine the rosetta compendium, enigma offers yet another perspective on the data. this mainly illustrates that the ideal clustering method does not exist, and that no single approach can extract all the information hidden in large compendium datasets. the elucidation of the regulatory networks governing the many different aspects of cellular function will therefore not only require the integration of different types of data, but also the integrated use of several complementary methods to analyze these data. we believe that enigma constitutes a valuable addition to the existing repertoire of analysis methods. 
we show that enigma constitutes a valuable addition to the repertoire of methods to analyze such data. compendia of gene expression profiles under chemical and genetic perturbations constitute an invaluable resource from a systems biology perspective. finally, we apply enigma to the rosetta compendium of expression profiles for saccharomyces cerevisiae and we analyze one pheromone response-related module in more detail, demonstrating the potential of enigma to generate detailed predictions. in contrast to the biclusters produced by most other methods, enigma modules may show internal substructure, i.e. enigma leverages differential expression analysis results to extract expression modules from perturbational gene expression data. the grouping of these related patterns in one module greatly aids in the biological interpretation of the data. we introduce a novel method, called enigma, that addresses some of these issues. it is increasingly recognized that perturbational expression compendia are essential to identify the gene networks underlying cellular function, and efforts to build these for different organisms are currently underway. we show that enigma outperforms other methods on artificial datasets, using a quality criterion that, unlike other criteria, can be used for algorithms that generate overlapping clusters and that can be modified to take redundancy between clusters into account. however, the perturbational nature of such data imposes specific challenges on the computational methods used to analyze them. for instance, some biclustering methods are less suited to identify overlapping biclusters, while others generate highly redundant biclusters.
 compendia of gene expression profiles under chemical and genetic perturbations constitute an invaluable resource from a systems biology perspective. however, the perturbational nature of such data imposes specific challenges on the computational methods used to analyze them. in particular, traditional clustering algorithms have difficulties in handling one of the prominent features of perturbational compendia, namely partial coexpression relationships between genes. biclustering methods on the other hand are specifically designed to capture such partial coexpression patterns, but they show a variety of other drawbacks. for instance, some biclustering methods are less suited to identify overlapping biclusters, while others generate highly redundant biclusters. also, none of the existing biclustering tools takes advantage of the staple of perturbational expression data analysis: the identification of differentially expressed genes.  we introduce a novel method, called enigma, that addresses some of these issues. enigma leverages differential expression analysis results to extract expression modules from perturbational gene expression data. the core parameters of the enigma clustering procedure are automatically optimized to reduce the redundancy between modules. in contrast to the biclusters produced by most other methods, enigma modules may show internal substructure, i.e. subsets of genes with distinct but significantly related expression patterns. the grouping of these related patterns in one module greatly aids in the biological interpretation of the data. we show that enigma outperforms other methods on artificial datasets, using a quality criterion that, unlike other criteria, can be used for algorithms that generate overlapping clusters and that can be modified to take redundancy between clusters into account. finally, we apply enigma to the rosetta compendium of expression profiles for saccharomyces cerevisiae and we analyze one pheromone response-related module in more detail, demonstrating the potential of enigma to generate detailed predictions.  it is increasingly recognized that perturbational expression compendia are essential to identify the gene networks underlying cellular function, and efforts to build these for different organisms are currently underway. we show that enigma constitutes a valuable addition to the repertoire of methods to analyze such data. 