Vol. 19 Suppl. 1 2003, pages i273–i282
DOI: 10.1093/bioinformatics/btg1038

BIOINFORMATICS

Genome-wide discovery of transcriptional
modules from DNA sequence and gene
expression
E. Segal ∗, R. Yelensky and D. Koller
Computer Science Department of Stanford University, Stanford, CA 94305-9010,
USA
Received on January 6, 2003; accepted on February 20, 2003

ABSTRACTIn this paper, we describe an approach for understandingtranscriptional regulation from both gene expression andpromoter sequence data. We aim to identify transcriptional
modules—sets of genes that are co-regulated in a setof experiments, through a common motif profile. Using
the EM algorithm, our approach refines both the module
assignment and the motif profile so as to best explain
the expression data as a function of transcriptional motifs.
It also dynamically adds and deletes motifs, as required
to provide a genome-wide explanation of the expression
data. We evaluate the method on two Saccharomyces
cerevisiae gene expression data sets, showing that ourapproach is better than a standard one at recovering
known motifs and at generating biologically coherent
modules. We also combine our results with binding
localization data to obtain regulatory relationships withknown transcription factors, and show that many of theinferred relationships have support in the literature.
Contact: eran@cs.stanford.edu
Keywords: probabilistic models, gene expression, transcriptional regulation.

INTRODUCTION
Many cellular processes are regulated at the transcriptional
level, by one or more transcription factors that bind to
short DNA sequence motifs in the upstream regions of theprocess genes. These co-regulated genes then exhibit similar patterns of expression. Given the upstream regions of
all genes, and measurements of their expression under various conditions, we could hope to ‘reverse engineer’ theunderlying regulatory mechanisms and identify transcriptional modules—sets of genes that are co-regulated underthese conditions through a common motif or combinationof motifs.
In this paper, we take a genome-wide approach for discovering this modular organization, based on the premise∗ To whom correspondence should be addressed.

c Oxford University Press 2003; all rights reserved.
Bioinformatics 19(Suppl. 1) that transcriptional elements should ‘explain’ the observed
expression patterns as much as possible. We define a probabilistic graphical model (Pearl, 1988) that integrates both
the gene expression measurements and the DNA sequence
data into a unified model. The model assumes that genes
are partitioned into modules, which determine the gene’s
expression profile. Each module is characterized by a motif profile, which specifies the relevance of different sequence motifs to the module. A gene’s module assignmentis a function of the sequence motifs in its promoter region. However, our model does not assume that all motifs
are necessarily active. In fact, as motifs are usually short,there are many genes where a motif is randomly present
but does not play a role. Furthermore, our goal is to discover motifs that play a regulatory role in some particular
set of experiments; a motif that is active in some settings
may be completely irrelevant in others. Our model identifies motif targets—genes where the motif plays an activerole in affecting regulation in a particular expression data
set. These motif targets are genes that have the motif andthat are assigned to modules containing the motif in theirprofile.
Our algorithm is outlined in Figure 1. It begins byclustering the expression data, creating one module fromeach of the resulting clusters. As the first attempt towards
explaining these expression patterns, it searches for acommon motif in the upstream regions of genes assigned
to the same module. It then iteratively refines the model,
trying to optimize the extent to which the expression
profile can be predicted transcriptionally. For example,
we might want to move a gene g whose promoter region
does not match its current module’s motif profile, toanother module whose expression profile is still a good
match, and whose motif profile is much closer. Given theseassignments, we could then learn better motif models andmotif profiles for each module. This refinement processarises naturally within our algorithm, as a byproduct of the
expectation maximization (EM) algorithm for estimating
the model parameters.
i273

E.Segal et al.

data selection

Expression data

Upstream sequences data selection

Clustering
MotifsearchGene partition
E-stepM-stepMotif set

Transcriptional modules

X

Add/Delete motifs
Transcriptional module discovery procedure

Annotation
analysis

Gene annotations (GO)Visualization & analysis

Protein complexes

Post-processing

Fig. 1. Schematic flow diagram of our proposed method. The pre-processing step includes selecting the input gene expression and upstream
sequence data. The model is then trained using EM, and our algorithm for dynamically adding and deleting motifs. It is then evaluated on
additional data sets.

In general, the motifs learned will not suffice to characterize all of the modules. As our goal is to provide agenome-wide explanation of the expression behavior, ouralgorithm identifies poorly explained genes in modules
and searches for new motifs in their upstream regions. The
new motifs are then added to the model and subsequently
refined using EM. As part of this dynamic learning
procedure, some motifs may become obsolete and areremoved from the model. The algorithm iterates untilconvergence, adding and deleting motifs, and refining
motif models and module assignments.
Our algorithm has several important advantages overother attempts to relate upstream sequences and expression data. First, we use both expression and sequence
data together, requiring that modules display a coherent
profile for both. This approach allows us to refine both thecluster assignments and motifs within the same algorithm.
In contrast, many approaches (e.g. Brazma et al., 1998;
Liu et al., 2001; Roth et al., 1998; Sinha and Tompa,
2000; Tavazoie et al., 1999) use gene expression measurements to define clusters of genes that are potentiallyco-regulated, and then search for common motifs in
the upstream regions of the genes in each cluster. The
expression analysis and motif finding are thus decoupled,and neither the clusters nor the motifs are re-evaluatedonce they are learned. Other approaches (e.g. Bussemaker
et al., 2001; Pilpel et al., 2001) work in the opposite
direction, first identifying a set of candidate motifs, andthen trying to explain the expression using these motifs.
However, these approaches use a prespecified set of
motifs, which are never adapted during the algorithm.
Our approach is based on the framework of Segal etal. (2002) but extends it in several important directions.
i274

First, their approach made use of DNA localization
data, which are not widely available for all organisms
and under multiple growth conditions. In contrast, weconstruct models that are based solely on sequence andexpression data, which are much easier to obtain. Second,
their approach used a predetermined number of motifs toconstruct the model. To allow a genome-wide analysis, ouralgorithm dynamically removes and adds motifs as needed
to explain the expression data as a whole. Finally, while
the models of Segal et al. (2002) allowed for detection ofcontext-specific regulation, the resulting structure is hardto interpret. Our model assigns each gene to one module,
facilitating interpretability.
We tested our method on two distinct Saccharomyces
cerevisiae expression datasets. We show that our learned
models find motifs that account for a much larger fraction of the observed expression patterns in comparison to
standard approaches that first cluster the expression profiles and then search for motifs in the upstream regionsof the genes in each cluster. Our approach also recovers amuch larger number of known motifs. We evaluated thefunctional coherence of our transcriptional modules using a gene functional annotation database and two protein complex databases that were not given to the model
as input. We found enrichment for many more groups inour models compared to standard approaches, suggestingthat our transcriptional modules are biologically more accurate. Finally, we used the recent binding assays of Lee
et al. (2002) to relate the actual transcription factors to themodules they regulate, resulting in a regulatory network;
we show that many of the regulatory relationships discovered have support in the literature.

Genome-wide discovery of transcriptional modules

PROBABILISTIC MODEL
The basic entities in our model are the genes in some set
G. We assume that the genes are partitioned into a setof K mutually exclusive and exhaustive transcriptional
modules. Thus, each gene is associated with an attributeM ∈ {1, . . . , K } whose value represents the module
to which the gene belongs. We now describe how thesemodules are related to expression profiles and to motif
profiles.
Gene expression model For each gene g in G, wehave expression measurements g.E 1 , . . . , g.E J , whereg.E j represents the log ratio mRNA expression levelmeasured for gene g in experiment j. We assume that
all of the genes in a single module exhibit the samegene expression behavior, and use the simple yet powerful
Naive Bayes model (Cheeseman and Stutz, 1995) torepresent this behavior. In this model, as applied in our
setting, we assume that the expression measurements areconditionally independent given the module assignment:P(E 1 , . . . , E J | M) =J
P(E j | M).

j=1

As the expression measurements are real-valued, we
model each conditional probability distribution P(E j |M = m) using a Gaussian distribution N (µ jm ; σ jm ).
Motif model The second key component in our model
is a set of variables that represent the regulation of thegene by motifs. For each gene g, we have a set ofbinary-valued Regulates variables R = {R1 , . . . , R L },where g.Ri takes the value true if motif i appears in the
promoter region of gene g, allowing the motif to play aregulatory role in the gene’s expression. We model the
motif using a standard position specific scoring matrix(PSSM Bailey and Elkan, 1994; Roth et al., 1998), whichassigns a weight to each position in the motif and eachnucleotide  ∈ {A, C, G, T }; this weight represents theextent to which the nucleotide’s presence in this positionis associated with the motif.
We use the discriminative PSSM approach of Segal et al.
(2002), which trains the PSSM weights to discriminate asmuch as possible between the presence and the absenceof the motif. This approach provides better predictions,and entirely avoids the problems arising from highfrequency but meaningless motifs that are common in
many upstream sequences. This model is specified using
a standard binary logistic model. We have p positionspecific weights w j [], one for each position j and eachletter  ∈ {A, C, G, T }, and a threshold w0 . For apromoter sequence of length N , we assume that binding
occurs once, and with equal probability at each of the N −p + 1 possible positions in the sequence. The probabilityP(g. R |g.S)
g. M P(g.E | g.M)

P(g.R |g .S )

g.S 1

g.S 2

g.S 3

1

g. M R

R

1
0

g.R1

2

g.R2

2

0

g.M

3
CPD 1

0
CPD 2

g.E1

g.E2

3
CPD 3

g.E3

Fig. 2. Illustration of our unified probabilistic model, for a simple
example with upstream regions of length three, two sequencemotifs, three possible module assignments and three expression
measurements for each geneof binding given the sequence is then specified simply as:
P(g.R = true | S1 , . . . , Sn ) =

 
n−
p+1
p


w0
exp{
wi [Si+ j−1 ]}.
logit logn − p + 1 j=1
i=1
Regulation model We define the motif profile of atranscriptional module to be a set of weights u mi , one foreach motif, such that u mi specifies the extent to which
motif i plays a regulatory role in module m. Roughly
speaking, the strength of the association of a gene g with aLmodule m is i=1
g.Ri u mi . The stronger the associationof a gene with a module, the more likely it is to be
assigned to it. We model this using a softmax conditional
distribution, a standard extension of the binary logistic
conditional distribution to the multi-class case:P(g.M = m̄ | R1 = r1 , . . . , R L = r L ) =L
exp{ i=1
u m̄i ri }.
LKm  =1 exp{ i=1 u m  i ri }As we expect a motif to be active in regulating onlya small set of modules in a given setting, we limit the
number of weights u 1i , . . . , u K i that are non-zero to someh  K . This restriction results in a sparse weightmatrix for P(M | R), and ensures that each regulator
affects at most h modules. In addition, for interpretability
considerations, we require all weights to be non-negative.
Intuitively, this means that a gene’s assignment to specifictranscriptional modules can only depend on features that
correpond to the presence of certain motifs and not on theabscence of motifs. For a module m, the set of motifs u mithat are non-zero are called the motif profile of m.
i275

E.Segal et al.

Unified model These three components, are put together
as a probabilistic graphical model, as shown in Figure 2;
the model defines the following joint distribution:
P(g.R, g.M, g.E | g.S) =LJP(g.Ri | g.S) · P(g.M | g.R) ·P(g.Ej | g.M),
i=1

j=1

where each of the above conditional probability distributions is parameterized as described in the previous sections.

LEARNING THE MODELSIn the previous section, we presented our probabilistic
model. We now turn to the task of learning this model from
data. Our data set D consists of a set of genes G, wherefor each gene g ∈ G we have a set of gene expression
measurements g.e j for j = 1, . . . , J and a DNA sequenceg.S in the upstream region of the transcription start site
for g. For this section, we restrict attention to a fixed
number of motifs, and address the problem of estimating
the model parameters to fit the data. The model parameters
to be estimated are: the means and variances of thenormal distributions of the expression model, the softmax
weights and structure of the module assignments (i.e.
which sequence motifs each module depends on), and the
PSSM weights for each sequence motif.
We follow the standard approach of maximum likelihood
estimation: we find the parameters θ that maximize
P(D | θ ). Our learning task is made considerably moredifficult by the fact that both the module assignment g.M
and the Regulates variables g.R are unobserved in the
training data. In this case, the likelihood function has
multiple local maxima, and no general method exists forfinding the global maximum. We thus use the Expectation
Maximization (EM) algorithm (Dempster et al., 1977),
which provides an approach for finding a local maximumof the likelihood function.
Starting from an initial guess θ (0) for the parameters,
EM iterates the following two steps. The E-step computes
the distribution over the unobserved variables given theobserved data and the current estimate of the parameters.
We use the hard assignment version of the EM algorithm,
where this distribution is used to select a likely completionof the hidden variables. The M-step then re-estimates theparameters by maximizing the likelihood with respect tothe completion computed in the E-step. This estimation
task differs for the different parts of the model.
E-step: inferring modules and regulation Our task in the
E-step is to compute the distribution over the unobserved
data, which in our setting means computing P(g.M, g.R |g.E, g.S). As genes are assumed to be independent,
i276

this computation can be done separately for each gene.
However, although the softmax distribution for P(g.M |g.R) has a compact parameterization, inference using this
distribution is still exponential in the number of Regulates
variables. Even if only a small number of these variables
is associated with any single module, for the purpose ofmodule assignment we need to consider all of the variables
associated with any module; this number can be quite
large, rendering exact inference intractable.
We devise a simple approximate algorithm for doing
this computation, which is particularly well-suited forour setting. It exploits our expectation that, while a large
number of sequence motifs determine module assignment,
only a small number of motifs regulate a particulartranscriptional module. Consequently, given the module
assignment for a gene, we expect a small number ofRegulates variables for that gene to take the value true.
Our approximate algorithm therefore searches greedily
for a small number of Regulates variables to activate foreach module assignment. For each gene g, it considers
every possible module assignment m, and finds a good
assignment to the Regulates variables given that g.M =m. This assignment is constructed in a greedy way,by setting g.R variables to true one at a time, aslong as P(g.M, g.R, g.E | g.S) improves. The joint
setting for g.M and g.R which gives the overall best
likelihood is then selected as the (approximate) most likely
assignment. For the remainder of this section, let g.m̄ andg.r̄1 , . . . , g.r̄ L represent the values selected for g.M andg.R1 , . . . , g.R L respectively by the E-step. Full details of
the algorithm are given in Figure 3a.
M-step: expression model Given the assignments ofgenes to modules as computed in the E-step, the maximum
likelihood setting for the parameters of the expression
model Gaussian distributions has a closed form solution.
Letting Nm be the number of genes assigned to module m,
we have that the mean and variance of the Gaussian forexperiment j given module assignment m areµm j =


1g.e j
Nm g∈G : g.m̄=mand
σm2 j =

1g.e2j − µ2jm .
Nm g∈G : g.m̄=mM-step: motif model We want the motif model tobe a good predictor of the assignment r̄ to the Regulates variables computed in the E-step. Thus, foreach Ri , we aim to find the values of the parameters
w j [] that maximize the conditional log probabilw0 , ityg∈G log P(g.r̄i | g.S1 , . . . , g.Sn ). Unfortunately,
this optimization problem has no closed form solution,

Genome-wide discovery of transcriptional modules

For each gene g ∈ G
Set g.M = 1
Set g.Ri = false for 1 ≤ i ≤ L
Set p = P(g.M, g.R | g.S, g.E)
For m = 1 to K // for all modules
Repeat // Find g.Ri that increases p
Set pbest = p
For i = 1 to L // for all regulates variables
Set g.Ri = truep  = P(g.M = m, g.R | g.S, g.E)
if p  > p
Set g.M = m
Set p = p else
Set g.Ri = false
Until pbest = p

Set U = {}Set iteration = 0
Let V = {vmi }1≤m≤K ,1≤i≤L
Set MaxScore = maxV Score[V ]// MaxScore = score of unconstrained fit
Set T = Threshold for closeness to MaxScore
RepeatSet iteration = iteration + 1
Let U  = {u mi }1≤m≤K ,1≤i≤L − U
Set U  = argmaxU  ≥0 Score[U  , U ]// Optimize weights not in U ; weights in U fixed
For i = 1 to L // for all regulates variables
Let m = argmax{u  }
 m mi 1≤m≤K
Set U = U {u mi } // Add new non-zero weightSet U = argmaxU ≥0 Score[U, 0]// Reoptimize weights in U ; other weights = 0
Until iteration = max iteration or Score[U ] >= MaxScore − T

Delete:For i = 1 to L // for all regulates variables
Set U  = U
Set u mi = 0 for 1 ≤ m ≤ K
If Score[U ] − Score[U  ] ≤ threshold
Delete Ri
Set U = U 
Add:
For m = 1 to K // for all modules
Let G = {}For each g s.t. g.m̄ = m
Set g.r̄ = argmaxr̄ P(g.r̄ | g.S)
Set g.m  = argmaxm P(g.M = m | g.R = g.r̄ )If m  = m

Set G = G {g}Learn motif with positive set G
Add new Regulates variable with learned PSSM

(a)

(b)

(c)

Fig. 3. (a) Search procedure for E-step of EM. (b) Learning the softmax distribution for P(g.M | g.R) in the M-step. (c) Procedurefor dynamically deleting and adding Regulates variables. In (b) and (c), U denotes the non-zero weights of PU (g.M | g.R), andScore[U ] = g∈G log PU (g.m̄ | g.r̄).

and there are many local maxima. We therefore use a
conjugate gradient ascent to find a local optimum in theparameter space. Conjugate gradient starts from an initial
guess of the weights w
 (0) . As for all local hill climbing
methods, the quality of the starting point has a huge
impact on the quality of the local optimum found by
the algorithm. We therefore initilize the weights using
the method of Barash et al. (2001), which efficiently
generates motif seeds of length 6–15 and then scores themusing the hypergeometric significance test. Each seedproduced by this method is then expanded to produce a
PSSM of the desired length, whose weights serve as an
initialization point for the conjugate gradient procedure.

motif profile U , which contains motifs that have non-zeroweight. The optimal setting for the weights in U is thenfound by optimizing these weights, under the constraintthat each weight in U is non-negative and the weights
not in U must be zero. This problem is also convex, andcan be solved using gradient methods. The algorithm thencontinues to search for additional motifs to include in
the profile U . It finds the optimal setting to all weights
while holding the weights in U fixed; it then selects thehighest weight motifs not in U , adds them to U , andrepeats. Weights are added to U until the sparseness limit
is reached, or until the addition of motifs to U does not
improve the overall score.

M-step: regulation model Finally, we consider the taskof estimating the parameters for the distribution P(g.M |g.R). Our goal is to find a setting for the softmax weights{u mi }1≤m≤K
,1≤i≤L so as to maximize the conditional log
probabilityg∈G log P(g.M = g.m̄ | g.R = g.r̄).
Although this optimization does not have a closed form
solution, the function is convex in the weights of thesoftmax. Thus, a unique global maximum exists, which
we can find using gradient ascent.
However, as we discussed in the previous section, we
also constrain this weight matrix to be sparse and eachweight to be non-negative. These constraints lead to moredesirable models, but also turn our task into a hardcombinatorial optimization problem. We use a greedy
selection algorithm, that tries to include non-zero weights
for the most predictive motifs for each Regulates variable
Ri . The algorithm, shown in Figure 3b, first finds the
optimal setting to the full weight matrix; as we discussed,
the optimal setting can be found using gradient ascent.
For each variable Ri , it then selects the most predictive
motif—the one whose weight is largest—and adds it to the

DYNAMICALLY ADDING AND REMOVINGSEQUENCE MOTIFSIn the previous section, we showed how to optimize themodel parameters given a fixed set of motifs. We now
wish to devise a dynamic learning algorithm, capable ofboth removing and adding sequence motifs as part of thelearning process. As we learn the models, some motifs
may not turn out to be predictive, or redundant given the
newly discovered motifs. Conversely, some modules may
not be well explained by sequence motifs, so that new
motifs should be added.
We add and remove motifs after each completion of the
EM algorithm. (Note that EM itself iterates several times
between the E-step and the M-step.) To determine whether
Ri should be deleted, we compute the conditional log
probability g∈G log P(g.m | g.r̄) both with and without
g.Ri , leaving the values of other Regulates variables
fixed. This computation tells us the contribution that Ri
makes towards the overall fit of the model. Variables
that contribute below a certain threshold are subsequently
removed from the model.
i277

We try to add motifs when the current set of motifs does
not provide a satisfactory explanation of the expression
data: when there are genes for which the sequence
predictions do not match the expression profile. We define
the residual for a transcriptional module m to be the setof genes that are assigned to module m in the E-step, butwould not be assigned to m based on the sequence alone.
We determine the sequence-only assignment of each geneby computingg.r̄ = argmaxr P(g.r | g.S)andg.m  = argmaxm P(g.M = m | g.R = g.r̄ ).
We then attempt to provide a better prediction for the
residual genes by adding a sequence motif that is trained
to match these genes. Once a new Regulates variables is
added, it becomes part of the model and its assignment
and parameterization is adapted as part of the next
EM iteration, as described in the previous section. This
process tests whether a new motif contributes to theoverall model fit, and may assign it a non-zero weight.
Importantly, a motif that was trained for the residualsof one module often gets non-zero weights for othermodules as well, allowing the same motif to participate inmultiple modules. Full details of the algorithm are given
in Figure 3c.

RESULTS
Models learned We evaluated our method separately on
two different S.cerevisiae gene expression datasets, one
consisting of 173 microarrays, measuring the responses tovarious stress conditions (Gasch et al., 2000), and anotherconsisting of 77 microarrays, measuring expression during
cell cycle (Spellman et al., 1998). We also obtained the
500bp upstream region of each gene (sequences wereretrieved from SGD (Cherry et al., 1998)).
The EM algorithm requires an initial setting to all parameters. We use the standard procedure for learning motifs from expression data to initialize the model parameters: we first cluster the expression profiles, resulting in apartition of genes to clusters, and then learn a motif foreach of the resulting clusters. For clustering the expression, we use the probabilistic hierarchical clustering algorithm of Segal et al. (2001). For learning motifs, we use
the motif finder described above. To specify the initial parameterization of our model, we treat these clusters andmotifs as if they were the result of an E-step, assigning avalue to all of the variables g.M and g.R, and learn themodel parameters as described above.
For the stress data, we use 1010 genes which showeda significant change in expression, excluding members of
the generic stress response cluster (Gasch et al., 2000).
i278

Num Genes Correctly Predicted

E.Segal et al.

350
330
310
290
270
250
230
210

Cell cycle modelStress model

190
170
150
1

2

3

4

5

6

Num Iterations

Fig. 4. Number of genes whose module assignment can be correctlypredicted based on sequence alone, where a correct prediction is
one that matches the module assignemnt when the expression is
included. Predictions are shown for each iteration of the learning
procedure.

We initialized 20 modules using standard clustering, andlearned the associated 20 sequence motifs. From thisstarting point, the algorithm converged after 5 iterations,of an EM step and a motif addition/deletion step, resulting
in a total of 49 motifs. For the cell cycle data, welearned a model with 15 clusters over the 795 cell cycle
genes defined in (Spellman et al., 1998). The algorithm
converged after 6 iterations, ending with 27 motifs.
Predicting expression from sequence Our approach aims
to explain expression data as a function of sequencemotifs. Hence, one metric for evaluating a model is its
ability to associate genes with modules based on theirpromoter sequence alone. Specifically, we compare themodule assignment of each gene when we consider only
the sequence data to its module assignment consideringboth expression and sequence data. Figure 4 shows thetotal number of genes whose expression-based module
assignment is correctly predicted using only the sequence,as the algorithm progresses and sequence motifs areadded. As can be seen, the predictions improve acrossthe learning iterations, and significantly outperform thestandard approach (which is iteration 1). Ultimately, ourmodel converges to 340 and 296 genes correctly predictedin the cell cycle and stress models, respectively, compared
to 158 and 152 for the standard approach.
Gene expression coherence These results indicate that
our model assigns genes to modules such that genes
assigned to the same module are generally enriched forthe same motifs. However, we can achieve such an organization by simply assigning genes to modules based only
on their sequence, while entirely ignoring the expresssion
data. To verify the quality of our modules relative to gene
expression data, we define the expression coherence ofa module to be the average Pearson correlation between
each pair of genes assigned to it, where the Pearson

14

0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

7

-Log (pvalue) for our method

0.8-Log(pvalue) for our method

Expression coherence (our method)

Genome-wide discovery of transcriptional modules

12
10
8
6
4
2
0
0

0.8

2

4

6

10

12

(a)

0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0.3

0.4

2
1
0
0

1

2

3

4

0.5

0.6

0.7

0.8

0.9

5

6

7

-Log (pvalue) for standard clustering(c)9-Log (pvalue) for our method-Log(pvalue) for our method

0.8

0.2

3

14

35

0.1

4

(b)

0.9

0

5

-Log(pvalue) for standard clustering

Expression coherence (standard clustering)

Expression coherence (our method)

8

6

30
25
20
15
10
5
0
0

5

10

15

20

25

30

35

8
7
6
5
4
3
2
1
0
0

1

2

3

4

5

6

7

8

Expression coherence (standard clustering)-Log(pvalue) for standard clustering-Log (pvalue) for standard clustering

(d)

(e)

(f)

9

Fig. 5. Comparison of standard clustering and the proposed method. (a)–(c) are for the cell cycle dataset (Spellman et al., 1998) and (d)–(f)are for the stress expression dataset (Gasch et al., 2000). (a), (d) Comparison of the expression coherence for each inferred module (or clusterin the standard clustering model). (b), (e) Comparison of enrichment of the targets of each motif for functional annotations from the GO
database. For each annotation, the largest negative log p-value obtained from analyzing the targets of all motifs is shown. (c), (f) Comparisonof enrichment of the targets of each motif for protein complexes. For each protein complex, shown is the largest negative log p-value obtained
from any of the motifs.

correlation isPearson(gi .E, g j .E) =L
1
(gi .El − µi ) (g j .El − µ j )
,
L l=1
σi
σj

where µi , σi are the mean and standard deviation of
the entries in gi .E. Figure 5a,d compares the expression
coherence of our modules to those built from standard
clustering for the cell cycle and stress data, showing
identical coherence of expression profiles. For the cell
cycle data, there was even a slight increase in thecoherence of the expression profiles for our model. Thus,
our model results in clusters that are more enriched for
motifs, while achieving the same quality of expression
patterns as standard clustering which only tries to optimize
the expression score.
Coherence of motif targets As we discussed, the motif
profile characterizing a module allows us to define a notionof motif targets—genes that contain the motif, and where
the motif plays a role in its expression profile, i.e. those

assigned to a module whose motif profile contains the
motif. In the standard clustering model, we can define thetargets of a motif to be those genes that have the motif andbelong to the cluster from which the motif was learned.
We tested whether our motif targets correspond tofunctional groups, by measuring their enrichment forgenes in the same functional category according to thegene annotation database of GO (Ashburner et al., 2000).
We used only GO categories with 5 or more genes
associated with them, resulting in 537 categories. For eachannotation and each motif, we computed the fraction ofgenes in the targets of that motif associated with thatannotation and used the hypergeometric distribution to
calculate a p-value for this fraction, and took p-value <0.05 to be significant. We compared, for both expression
data sets, the enrichment of the motif targets for GO
annotations between our model and standard clustering.
We found many annotations that were enriched in bothmodels. However, there were 24 and 29 annotations thatwere significantly enriched in our cell cycle and stressi279

E.Segal et al.

models, respectively, that were not enriched at all in thestandard clustering model, compared to only 4 and 14
categories only enriched in the standard clustering modelfor these respective models. Among those categories
enriched only in our model were carbohydrate catabolism,
cell wall organization and galactose metabolism, all ofwhich are processes known to be active in response tovarious stress conditions that we can now characterize by
sequence motifs. A full comparison of the GO enrichment
for both datasets is shown in Figure 5b,e.
Since functional categories do not necessarily correspond to co-regulation groups, we also tested theenrichment of our motif targets for protein complexes,
as compiled experimentally in the assays of Gavin et al.
(2002) and Ho et al. (2002), consisting of 590 and 493
complexes, respectively. The member genes of protein
complexes are often co-regulated and we thus expect
to find enrichment for them in our motif targets. We
associated each gene with the complexes it is assigned
to in each protein complex dataset and computed thep-value of the enrichment of the targets of each motif
for each complex, as we did above for the GO annotations. The results for the cell cycle and stress datasets
are summarized in Figure 5c,f, showing much greater
enrichment of our motif targets than the targets of the
motifs identified using the standard approach, with 63 and10 complexes significantly enriched only in our model,
and no complexes only enriched in the standard approach,for the stress and cell cycle models, respectively.
Motifs and motif profiles We compared the motifs we
identified to motifs from the literature. Of the 49 motifs
learned for the stress model, 22 are known, compared
to only 10 known motifs learned using the standard
approach. For the cell cycle model, 15 of the 27 learned
motifs are known, compared to only 8 known motifs
learned using the standard approach. Many of the known
motifs identified, such as the stress element STRE, theheat shock motif HSF and the cell cycle motif MCM1, arealso known to be active in the respective datasets.
A powerful feature of our approach is its ability tocharacterize modules by motif profiles. This ability isparticularly important for higher eukaryotes, in which
regulation often occurs through multiple distinct motifs.
To illustrate the motif profiles found by our approach, we
found for each motif all modules enriched for the presenceof that motif. This was done by associating each gene withthe motifs in its upstream region, and then computing thep-value of the enrichment of the member genes of each
module. Figure 6a shows all the module-motif pairs inwhich the module was enriched for the motif with p-value<0.05. In addition, the figure indicates (by red circles) all
pairs in which the motif appears in the module’s motif
profile. As can be seen, many of profiles contain multiple
i280

motifs, and many motifs were used by more than one
module. Even though modules share motifs, each module
is characterized by a unique combination of motifs.
Inferring regulatory networks Identifying the active motifs is a significant step towards understanding the regulatory mechanisms governing gene expression. However, we
would also want to know the identity of the transcription
factor (TF) molecules that bind to these sequence motifs.
We used the DNA binding assays of Lee et al. (2002),
that directly detect to which promoter regions a particular
TF binds in vivo, and associated TFs with the motifs welearned. For each motif, we computed the fraction, among
the motif targets, of genes bound by each TF, as measured
in the data of Lee et al. We used the hypergeometric distribution to assign a p-value to each such fraction and took
p-value <0.05 to be significant. Inspection of the significant associations showed that, in most cases, there was aunique motif that was significant for the TF and that a highfraction (>0.5) of the TF’s binding targets were among the
motif target genes.
Based on this strong association between TFs andmotifs, for each such TF-motif pair, we predicted that
the TF regulates all the modules that are characterized
by the motif. By combining all associations, we arrived
at the regulatory network shown in Figure 6b. Of the
106 transcription factors measured in Lee et al., 28 wereenriched in the targets of at least one motif and werethus added to the resulting network. Of the 20 modules,
16 were associated with at least one TF. To validate
the quality of the network, we searched the biologicalliterature and compiled a list of experimentally verified
targets for each of the 28 TFs in our network. We thenmarked each association between a TF and a module assupported if the module contains at least one gene that the
TF is known to regulate from biological experiments. As
current knowledge is limited, there are very few known
targets for most TFs. Nevertheless, we found support
for 21 of the 64 associations. We also computed the pvalue for each supported association between a TF and a
module, using the binomial distribution with probabilityof success p = t/N , where K is the total number ofknown targets for the TF and N is the total number ofgenes (1010). The p-value is then P(X ≥  | X ∼B( p, n)), where  is the total number of known targets of
the regulator in the supported module and n is the numberof genes in the supported module. The resulting p-values
are shown in Figure 6b by edge thickness and color.
We assigned a name to each module based on a concise
summary of its gene content (compiled from both geneannotation and literature). The regulatory network thus
contains predictions for the processes regulated by each
TF, where for each association the prediction includes
the motif through which the regulation occurs. In many

5 8
Fzf1

Swi4

Nrg1

Gat1

Ste12

Hap4

Dal82

Cha4

Hir2

Hms1

Gcr1

17 18 11
9 12

19

1

14

13

4

0

10

16

7

Supported (p<.001)
Supported (p<.05)

2

Pdr1

Dig1

Arg81

Met4

Aro80

Cbf1

Fkh1

Yap3

Dot6

Crz1

Rfx1

Mig1

Gln3

Swi6

3
Gcr2

19 – Glycolysis (80)
18 – Redox regulation (10)17 – Mixed I (74)16 – Protein folding (45)
15 – Chromatin remodeling (8)14 – Energy and TCA cycle (66)
13 – Oxidative phosphorylation (31)
12 – Unknown I (75)
11 – Nitrogen metabolism (35)
10 – Cell cycle (44)9 – Protein folding (47)
8 – Ergosterol biosynthesis (65)
7 – Sucrose metabolism (55)
6 – Translation (80)
5 – Heat shock (61)
4 – Pheromone response (93)
3 – Transport (30)
2 – Aldehyde metabolism (34)
1 – Galactose metabolism (15)0 – Amino acid metabolism (52)

Cad1

15

Swi5

Motif 15
Motif 5
Motif 12
Motif 2
Motif 17
Motif 7
Motif 14
Motif 19
Motif 11
Motif 9
Motif 1
Motif 16
Motif 6
Motif 13
Motif 3
Motif 10
Motif 8Added Motif 1
Added Motif 2
Added Motif 3
Added Motif 6
Added Motif 7
Added Motif 9
Added Motif 11
Added Motif 12
Added Motif 13
Added Motif 14
Added Motif 15
Added Motif 16
Added Motif 17
Added Motif 21
Added Motif 27
Added Motif 29
Added Motif 30
Added Motif 31
Added Motif 33
Added Motif 34
Added Motif 36
Added Motif 39
Added Motif 42
Added Motif 48
Added Motif 51
Added Motif 54
Added Motif 55
Added Motif 60
Added Motif 61
Motif 20
Motif 4
Added Motif 10

Genome-wide discovery of transcriptional modules

6

Supported (p>.05)
Not supported

Fig. 6. (a) Matrix of motifs vs. modules for the stress data, where a module-motif entry is colored if the member genes of that module wereenriched for that motif with p-value <0.05. The intensity corresponds to the fraction of genes in the module that had the motif. Entries in themodule’s motif profile are circled in red. Modules were assigned names based on a summary of their gene content. (b) Regulatory network
inferred from our model using the DNA binding assays of Lee et al.. Ovals correspond to transcription factors and rectangles to modules (see(a) for module names).

cases, our approach recovered coherent biological processes along with their known regulators. Examples ofsuch associations include: Hap4, the known activator ofoxidative phosphorylation, with the oxidative phosphorylation module (13); Gcr2, a known positive regulatorof glycolysis, with the glycolysis module (19); Mig1, a
glucose repressor, with the galactose metabolism module(1); Ste12, involved in regulation of pheromone pathways,
with the pheromone response module (4); and Met4, a
positive regulator of sulfur amino acid metabolism, withthe amino acid metabolism module (0).

CONCLUSIONSWe presented a unified probabilistic model over both gene
expression and sequence data, whose goal is to identify
transcriptional modules and the regulatory motif bindingsites that control their regulation within a given set of experiments. Our results indicate that our method discovers
modules that are both highly coherent in their expression
profiles and significantly enriched for common motif binding sites in upstream regions of genes assigned to the same
module. A comparison to the common approach of constructing clusters based only on expression and then learning a motif for each cluster shows that our method recovers
modules that have a much higher correspondence to external biological knowledge of gene annotations and protein
complex data.
ACKNOWLEDGEMENTS
This work was supported by the National Science Foundation, grant ACI-0082554. Eran Segal was also supported
by a Stanford Graduate Fellowship (SGF).

REFERENCES
Ashburner,M., Ball,C.A., Blake,J.A., Botstein,D., Butler,H.,
Cherry,J.M., Davis,A.P., Dolinski,K., Dwight,S.S., Eppig,J.T.,
et al. (2000) Gene ontology: tool for the unification of biology.
The gene ontology consortium. Nat. Genet., 25, 25–29.
Bailey,T.L. and Elkan,C. (1994) Fitting a mixture model by
expectation maximization to discover motifs in biopolymers.
Proc. Int. Conf. Intell. Syst. Mol. Biol. Volume 2. pp. 28–36.
Barash,Y., Bejerano,G. and Friedman,N. (2001) A simple hypergeometric approach for discovering putative transcription factor
binding sites. Algorithms in Bioinformatics, Number 2149 in
LNCS, pp. 278–293.
Brazma,A., Jonassen,I., Vilo,J. and Ukkonen,E. (1998) Predicting
gene regulatory elements in silico on a genomic scale. GenomeRes., 8, 1202–1215.
Bussemaker,H., Li,H. and Siggia,E. (2001) Regulatory element
detection using correlation with expression. Nat. Genet., 27,
167–171.
Cheeseman,P. and Stutz,J. (1995) Bayesian classification (AutoClass): Theory and results. Advances in Knowledge Discovery
and Data Mining. AAAI Press, Menlo Park, CA, pp. 153–180.
Cherry,J.M., Adler,C., Ball,C., Chervitz,S.A., Dwight,S.S., Hester,E.T., Jia,Y., Juvik,G., Roe,T., Schroeder,M., Weng,S., Botstein,D. (1998) Sgd: Saccharomyces genome database. Nucleic
Acid Res., 26, 73–79.
Dempster,A.P., Laird,N.M. and Rubin,D.B. (1977) Maximum likelihood from incomplete data via the EM algorithm. J. Roy. Stat.
Soc. B, 39, 1–39.
Gasch,A.P.,
Spellman,P.T.,
Kao,C.M.,
Carmel-Harel,O.,
Eisen,M.B., Storz,G., Botstein,D. and Brown,P.O. (2000)
Genomic expression program in the response of yeast cells
to environmental changes. Mol. Biol. Cell, 11, 4241–4257.
Gavin,A.C., Bosche,M., Krause,R., Grandi,P., Marzioch,M.,
Bauer,A., Schultz,J., Rick,J.M., Michon,A.M. and Cruciat,C.M.,
et al. (2002) Functional organization of the yeast proteome
by systematic analysis of protein complexes. Nature, 415,
141–147.

i281

E.Segal et al.

Ho,Y., Gruhler,A., Heilbut,A., Bader,G.D., Moore,L., Adams,S.L.,
Millar,A., Taylor,P., Bennett,K., Boutilier,K., et al. (2002) Systematic identification of protein complexes in Saccharomyces
cerevisiae by mass spectometry. Nature, 415, 180–183.
Lee,T., Rinaldi,N.J., Robert,F., Odom,D.T., Bar-Joseph,Z., Gerber,G.K., Hannett,N.M., Harbison,C.T., Thompson,C.M., Simon,I., et al. (2002) Transcriptional regulatory networks in Saccharomyces cerevisiae. Science, 298, 824–827.
Liu,X., Brutlag,D. and Liu,J. (2001) Bioprospector: discoveringconserved DNA motifs in upstream regulatory regions of coexpressed genes. Pac. Symp. Biocomput. pp. 127–138.
Pearl,J. (1988) Probabilistic Reasoning in Intelligent Systems.
Morgan Kaufmann.
Pilpel,Y., Sudarsanam,P. and Church,G. (2001) Identifying regulatory networks by combinatorial analysis of promoter elements.
Nat. Genet., 29, 153–159.
Roth,F., Hughes,P., Estep,J.D. and Church,G. (1998) Finding DNA
regulatory motifs within unaligned noncoding sequences clustered by whole-genome mRNA quantitation. Nat. Biotechnol.,i282

16, 939–945.
Segal,E., Taskar,B., Gasch,A., Friedman,N. and Koller,D. (2001)
Rich probabilistic models for gene expression. Bioinformatics,
17(Suppl 1), S243–S252.
Segal,E., Barash,Y., Simon,I., Friedman,N. and Koller,D. (2002)
From sequence to expression: A probabilistic framework. In
Proc. RECOMB. pp. 263–272.
Sinha,S. and Tompa,M. (2000) A statistical method for findingtranscription factor binding sites. In Proc. Int. Conf. Intell. Syst.
Mol. Biol., Volume 8, pp. 344–354.
Spellman,P.T., Sherlock,G., Zhang,M.O., Iyer,V.R., Anders,K.,
Eisen,M.B., Brown,P.O., Botstein,D. and Futcher,B. (1998)
Comprehensive identification of cell cycle-regulated genes ofthe yeast Saccharomyces cerevisiae by microarray hybridization.
Mol. Biol. Cell, 9(12), 3273–3297.
Tavazoie,S., Hughes,J.D., Campbell,M.J., Cho,R.J. andChurch,GM. (1999) Systematic determination of genetic
network architecture. Nat. Genet., 22, 281–285 Comment in:Nat. Genet. 22 213-215.