Module networks revisited: computational assessment and prioritization of model predictions

1 MB·English

by  Anagha Joshi

Module networks revisited: computational assessment and prioritization of
model predictions
Anagha Joshi,1,2 Riet De Smet,3 Kathleen Marchal,3,4 Yves Van de Peer,1,2 and Tom Michoel1,2,∗
1Department of Plant Systems Biology, VIB, Technologiepark 927, B-9052 Gent, Belgium
2Department of Molecular Genetics, UGent, Technologiepark 927, B-9052 Gent, Belgium
3CMPG, Department of Microbial and Molecular Systems, KULeuven, Kasteelpark Arenberg 20, B-3001 Leuven,
Belgium
4ESAT-SCD, KULeuven, Kasteelpark Arenberg 10, B-3001 Leuven, Belgium
Motivation:Thesolutionofhigh-dimensionalinferenceandpredictionproblemsincomputationalbi-
ology is almost always a compromise between mathematical theory and practical constraints such
as limited computational resources. As time progresses, computational power increases but well-
establishedinferencemethodsoftenremainlockedintheirinitialsuboptimalsolution.
Results:WerevisittheapproachofSegaletal.(2003)toinferregulatorymodulesandtheircondition-
specificregulatorsfromgeneexpressiondata. Incontrasttotheirdirectoptimization-basedsolution
weuseamorerepresentativecentroid-likesolutionextractedfromanensembleofpossiblestatistical
modelstoexplainthedata. Theensemblemethodautomaticallyselectsasubsetofmostinformative
genesandbuildsaquantitativelybettermodelforthem.Geneswhichclustertogetherinthemajority
9 ofmodelsproducefunctionallymorecoherentmodules.Regulatorswhichareconsistentlyassignedto
0 amodulearemoreoftensupportedbyliterature,butasinglemodelalwayscontainsmanyregulator
0 assignmentsnotsupportedbytheensemble. Reliablydetectingcondition-specificorcombinatorial
2 regulationisparticularlyhardinasingleoptimumbutcanbeachievedusingensembleaveraging.
Availability: Allsoftwaredevelopedforthisstudyisavailablefromhttp://bioinformatics.psb.
n ugent.be/software/.
a Supplementary information: Supplementary data and figures are available from http://
J bionformatics.psb.ugent.be/supplementary data/anjos/module nets yeast/.

2
1

I. INTRODUCTION tionalmethods. Theseimprovementsconcerntheuseof
]
M Monte Carlo (Liu, 2004) and ensemble strategies (Car-
valhoandLawrence,2007;Webb-Robertsonetal.,2008).
Q Oneofthecentralgoalsofthetop-downapproachto
systemsbiologyistoinferpredictivemathematicalnet- Following (Hartwell et al., 1999) a ‘module’ is to be
.
o work models from high-throughput data. Much of the viewed as a discrete entity composed of many types of
i
b drivingforceforthedevelopmentofnetworkinference moleculesandwhosefunctionisseparablefromthatof
- methodshascomefromtheavailabilityofvarioustypes other modules. Understanding the general principles
q
of large-scale data sets for particular model organisms that determine the structure and function of modules
[
likeS.cerevisiaeandE.coli. Incontrast,datageneration and the parts they are composed of can be considered

1 for other organisms has been much slower and mainly oneofthemainproblemsofcontemporarysystemsbiol-
v focusedongeneexpressiondata.Thesegeneexpression ogy(Hartwelletal.,1999). Themodulenetworkmethod
4 data sets for typically more complex organisms pose of(Segaletal.,2003)addressesthisproblemusinggene
4
theirownchallenges,suchasahighernumberofgenes, expressiondataasitsinput. Ithasyieldednovelbiolog-
5
1 limited number of experimental conditions, and sup- icalinsightsinanumberofcomplexeukaryoticsystems
. posedlyamorecomplexunderlyingtranscriptionalnet- (Lee et al., 2006; Li et al., 2007; Novershtern et al., 2008;
1
work. Therefore,improvementandrefinementofmeth- Segaletal.,2003,2007;Zhuetal.,2007)andhasbeenthe
0
ods for network inference from gene expression data source of inspiration for numerous computational ap-
9
0 continues to be of great interest. Several reviews on a proachestonetworkinferenceasevidencedbyitshigh
: varietyofmethodshavebeenwritten(Bansaletal.,2007; numberofcitations. Amodulenetworkisaprobabilis-
v
Bussemaker et al., 2007; Friedman, 2004; Gardner and tic graphical model (Friedman, 2004) which consists of
i
X Faith,2005),anddevelopmentofnewmethodsremains modulesofcoregulatedgenesandtheirregulatorypro-
r anactiveareaofresearch(AlterandGolub,2005;Basso grams. Aregulatoryprogramusestheexpressionlevel
a etal.,2005;Bonneauetal.,2006;Faithetal.,2007). Here of a set of regulators to predict the condition depen-
we revisit the module network method of (Segal et al., dent mean expression of the genes in a module. (Segal
2003) to infer regulatory modules and their condition- etal.,2003)usedadeterministicoptimizationalgorithm
specificregulatorsfromgeneexpressiondataandshow thatsearchessimultaneouslyforapartitionofgenesinto
that better and more refined module networks can be modulesandaregulationprogramforeachmodule.We
obtained by using advanced statistical and computa- consider both as separate tasks. When searching for
modules, often many local optima exist with partially
overlappingmodulesdifferingfromeachotherinafew
genes. WeuseaGibbssamplingapproachfortwo-way
∗Correspondingauthor,E-mail:[email protected] clustering of genes and conditions to generate an en-
2
sembleofpartiallyoverlappingpartitionsofgenesinto hypotheses are statistically supported by the ensemble
modules and produce an ensemble averaged solution method. This example illustrates the usefulness of an
(Joshietal.,2008). Thiscentroidsolutionconsistsofso- approachwhichcangenerateinternalsignificancemea-
calledtightclusters,subsetsofgeneswhichconsistently sures,inparticularifnootherdatasourcesareavailable
clustertogetherinalmostalllocaloptima. Wealsousea tovalidatehypothesesgeneratedbyasinglelocalopti-
probabilistic method for learning regulatory programs. mum.
These regulatory programs take the form of fuzzy de- Together all these results convincingly show that the
cisiontreeswithregulatorexpressionlevelsatthedeci- ensemblemethodforlearningmodulenetworkssignif-
sion nodes and generalize the regression tree approach icantlyimprovesthedirectoptimizationmethodof(Se-
of (Segal et al., 2003). By summing the strength with galetal.,2003). Unlikeasingleoptimum,ensembleav-
whicharegulatorparticipatesineachmemberofanen- eraging allows the assessment and prioritization of the
sembleofregulatoryprogramsforacertainmodule,we statistically most reliable modules and their condition-
obtain a regulator score which gives a statistical confi- specific regulators. Such high-confidence modules can
dencemeasurefortheassignmentofthatregulator. To- be used directly for generating experimentally verifi-
gether,theGibbssamplingclusteralgorithmandproba- ablehypothesesorcanbeintegratedwithother,perhaps
bilisticregulatoryprogramlearningprovideacomputa- smaller-scale,datasourcestocreateamorecomprehen-
tionallyefficientmethodtogenerateensemblesofmod- siveviewoftheunderlyingnetworks.
ulenetworksfromwhichacentroid-likesummarization
canbeconstructed.
II. RESULTS AND DISCUSSION
We have applied this ensemble method to the very
same data set as (Segal et al., 2003) and performed sev-
A. Data and procedure
eral comparison tasks. First, we considered the proba-
bilistic models and evaluated them on training as well
We obtained all data from the supplemental website
astestdata. Weshowthatthemodelinferredby(Segal
of (Segal et al., 2003), including expression data, gene
et al., 2003) is equivalent to a single instance of the en-
modules and regulatory programs. Using the Gibbs
semble of models inferred by our algorithm. The tight
sampler we generated 12 different partitions of genes
clustersobtainedfromtheensemblesolutiongeneratea
intomoduleswhichwerecombinedintoonesetoftight
quantitatively better model than each of the single in-
clusters. The number of clusters is determined auto-
stances, includingthemodelof(Segaletal.,2003). Sec-
matically by the Gibbs sampler and ranges from 65 to
ond, we compared the clustering of genes. Tight clus-
78 in the different runs, compared to the predefined
ters are in general more functionally coherent and im-
value of 50 of (Segal et al., 2003). 1892 of the 2355
prove the original modules in two ways. They can re-
genes in the data set could be assigned with high con-
movespuriousprofilesandfetchonlythecoreoftightly
fidence to 69 tight clusters. To generate regulator as-
coexpressed genes from a single module, or they can
signment scores, we learned 10 probabilistic regulation
merge separate but related modules into one cluster.
programspermodulewith100regulatorandsplitvalue
Third, we used the regulator score to analyze the net-
pairs sampled per regulation program node. More de-
work of modules and their associated regulators from
tails about these procedures are given in the Methods.
(Segaletal.,2003). Weshowthatthisnetworkcontains
Thisresultedinfourdifferentmodulenetworkmodels:
both high- and low-scoring regulators and that several
high-scoring regulators are missed by the solution of 1. SCSR: Segal clusters with Segal regulation pro-
(Segal et al., 2003). In general, regulator assignments grams,correspondingexactlytotheresultsof(Se-
whichcanbevalidatedbyexternalsourcessuchasChIP galetal.,2003).
data or literature are highly ranked. In combination
withthetightclusters,theprobabilisticmethodassigns 2. SCPR: Segal clusters with probabilistic regulation
more regulators supported by literature and the clus- programs.
terstowhichtheyareassignedcontainahigherratioof
3. GCPR: Gibbs sampler clusters (single run) with
knowntargetscomparedtothemodulenetworkof(Se-
probabilisticregulationprograms.
galetal.,2003). Fourth,weshowthattheregulatorscor-
ingschemecanalsobeusedtoinfercontext-specificand
4. TCPR:Tightclusters(multipleGibbssamplerruns
combinatorialregulationbyidentifyingpairsofregula-
combined)withprobabilisticregulationprograms.
torswhichoccursignificantlyoftentogetherinthesame
regulationprogram.
Finally we have applied the ensemble method to B. Model evaluation
a bHLH module network that was recently inferred
for mouse brain (Li et al., 2007). (Li et al., 2007) Amodulenetworkinfersaprobabilisticmodelwhich
used their module network to make several hypothe- explains relations between expression levels of a set
sesaboutmodesofcombinatorialregulationamongdif- of genes. More precisely, there is a probability dis-
ferent brain tissues. We show that only few of these tribution p(x ,...,x ) which computes the probability
1 N
3
(a) (b) 1
2.5 1.4
0.9
1.2 0.8
2
0.7
1
0.6
1.5 0.8 0.5
0.4
1 0.6 0.3
0.4 0.2
0.5 0.1
0.2
0
1111....−01224684 −3.5 −3 −2.5 lo(gc−(2p))/N −1.5 −1 −0.5 0 12..−02554 −3.5 −3 −2.5 l(ogd−(2p))/N −1.5 −1 −0.5 0 purine nucleotide anabolismamino acid transportprotein folding electron transport, energy conservationrespirationelectron transporttransported compounds amino acid metabolism regulationgeneral transcription activitiesprotease inhibitoraerobic respirationmitochondrial transportenergy conversion detoxification by modificationhomeostasis of protonsstress responsetranscriptional controlperoxidase reactionmetabolism of the aspartate familymetabolism of methioninefatty acid metabolismglyoxylate cyclelipid, fatty acid metabolismnuclear membranebiosynthesis of methioninepurin nucleotide metabolismnuclear transportstructural protein bindingtranscription activationnutrient starvation responsemetabolism of asparaginemembrane lipid metabolismglycolipid metabolismdegradation of asparagineoxygen proteasomal degradation sulfatetransporthomeostasis of sulfatemodification with sugar residues oxidative stress responsehomeostasis of anionscomplex bindingmetabolism of aspartateanion transportmetabolism of tyrosinemetabolism of the cysteine metabolism of phenylalaninefermentationdetoxification
0.8 1
0.6 FIG.2 Histogramofthehighestfractionofgenesinonemod-
0.4 0.5
uleinaMIPSfunctionalcategoryforTC(red)andSC(blue),
0.2
sortedbyratiodifference.
−04 −3.5 −3 −2.5 log−(2p)/N −1.5 −1 −0.5 0 −04 −3.5 −3 −2.5 log−(2p)/N −1.5 −1 −0.5 0
FIG.1 Modelevaluationexperiments. (a)HistogramofŁ =
1 logp(x ,...,x )forSCPR(blue)andnon-parametricfitof hasahighermeanŁthanSCPRwithaone-tailedt-test
N 1 N
the histogram for GCPR (black curve). (b) Histogram of Ł (α = 0.01). This shows that tight clusters are selecting
forSCSR(red)overlayedonhistogramofŁforSCPR(blue), a subset of genes which are the most informative and
with non-parametric fits (black curves). (c) Histogram of Ł thereforegenerateabettermodel.
for SCPR (blue) overlayed on histogram of Ł for TCPR (ma-
Nextwetestedhowwellthesemodelsexplainunseen
genta),withnon-parametricfits(blackcurves).(d)Histogram
data by performing a cross-validation experiment. We
andnon-parametricfit(leftblackcurve)ofŁforGCPRlearned
removed10%oftheconditionsatrandomfromthecom-
ontrainingdataandevaluatedontestdata(green)andnon-
pletedata(thetestset)andrantheGibbssampleronce
parametricfitofthesamemodelsevaluatedontrainingdata
(rightblackcurve).Allhistogramsandcurvesarenormalized on the remaining 90% (the training set). The resulting
tohaveareaequalto1. model was then evaluated on the test set. This proce-
dure was repeated 10 times and all test set evaluation
values were collected in one histogram and compared
(density)toobserveaparticularcombinationofexpres- tothetrainingsetvalues(Figure1(d)). Thecurveofthe
sion levels x for a set of N genes. This probabilis- test set is slightly shifted to the left with respect to the
i
tic model predicts the response in expression of genes trainingsetcurve,asonewouldexpect,butbothcurves
in a module upon perturbations of its regulators, such havethesamemeanwithaone-tailed t-test(α = 0.01).
as knock-out or overexpression, and thus yields bio- This shows that the probabilistic models indeed gener-
logically verifiable hypotheses. For a module network, alizetounseendata.
the distribution p(x ,...,x ) is a product of N factors
1 N
(see Methods), so we consider the normalized quan-
tity Ł = 1 logp(x ,...,x ) which can be compared C. Gene clustering improvement
N 1 N
between models with potentially different numbers of
genes. Higher values of Ł mean better explanation of We have shown in the previous section that SCPR is
the data by the model, i.e. more accurate prediction of equivalenttoGCPRbutTCPRgivesabettermodelover
theoutcomeofnewexperiments. SCPR.Wealsoobservethattightclusters(TC)areover-
Firstweperformedevaluationsoneachofthecondi- allmorefunctionallycoherentthantheclustersobtained
tionsintheoriginaldataset. Figure1(a)showsthatthe in(Segaletal.,2003)(SC).Figure2showsthefractionof
histogram of Ł-values for SCPR fits well within a non- genes in a cluster belonging to a MIPS functional cate-
parametric curve fit of the histogram for GCPR. This gorywhichissignificantlyoverrepresented(p < 0.001)
impliesthattheclustersfoundby(Segaletal.,2003)are in SC and TC. Several examples illustrate the general
equivalenttoonelocaloptimumidentifiedbytheGibbs trend seen in this figure. In TC-40, 4/7 genes are in-
samplerprocedure. Figure1(b)showsthehistogramof volvedinaminoacidtransportcomparedtoSC-27with
Ł-valuesforSCSR(red)overlayedonthehistogramfor 8/53 genes. In TC-27, 7/9 genes belong to purine nu-
SCPR (blue), both with non-parametric curve fits. The cleotideanabolismcomparedtoSC-11with6/53genes.
meanŁ-valuesobtainedbySCPRarehigherthanSCSR Segal cluster 1 (SC-1) contains 55 genes, out of these
byaone-tailedt-test(α = 0.01)provingthatprobabilis- 32(58%)arevalidatedtargetsofHap4, aglobalregula-
ticregulationprogramsgiveabetterexplanationofthe tor of respiratory genes, according to the YEASTRACT
data. InFigure1(c)wecomparedSCPRtoTCPR.TCPR database (Teixeira et al., 2006). This cluster has maxi-
4
(a) (b)
SC-8
SC-11 SC-9
Hap4 targets in Respiratory
Yeastract database genes
TC-11
TC-27
FIG.4 (a)TC-27fetchesthecoreoftightlycoexpressedgenes
fromSC-11;67%genesinTC-27areknowntobeBas1targets.
(b)SC-8andSC-9whichhavesimilarexpressionaremerged
into TC-11. SC-8, SC-9 and TC-11 all are enriched for Gat1
targets.
FIG. 3 TC-7 with Hap4 assigned as a top regulator. Genes
knowntoberegulatedbyHap4inYEASTRACTaremarkedin
blueandthoseinvolvedinrespirationaremarkedinorange.
wereseparateinSC(Figure4(b)).
mum overlap with tight cluster 7 (TC-7) with 30 genes
D. Regulator assignment prioritization
out of which 25 (83%) are known Hap4 targets. The
five remaining genes are Qcr6, Cox5a and Fum1, all
The ensemble approach generates multiple equally
located in mitochondrion and involved in respiration,
plausible regulatory programs for a single module in
and two unknown genes Ygl188c and Ygr182c. With
a probabilistic fashion. The regulator assignment score
24/30 respiratory genes (80%), TC-7 even improves on
which takes into account how often a regulator is as-
COGRIM (Chen et al., 2007) which combines multiple
signedtoamodule,withwhatscore,andatwhichlevel
data sources. Using expression data alone (the same
intheregulationtree,canthereforebeusedtoprioritize
datasetas(Segaletal.,2003)),(Chenetal.,2007)obtain
regulators(highestregulatorscoregetstopmostrank).
aclusterwith32/51(62%)genesbelongingtoMIPSres-
First we consider only the difference between prob-
pirationcategory. UsingbothChIPandexpressiondata,
abilistic regulator assignment and the original method
theyobtainaclusterwith23/34(68%)respiratorygenes,
bycomparingSCSRwithSCPR,hencekeepingthegene
significantlylowerthanTC-7.Figure3showsTC-7with
modules the same for both methods. Figure 5 shows
known Hap4 targets and respiratory genes marked in
regulator-module links in SCSR (cfr. Figure 5 in (Segal
blueandorangerespectively.
et al., 2003)). The edges colored red are the ones sup-
TC-27containsninegeneswhichformasubsetofSC-
ported by literature (data from (Segal et al., 2003)). To
11 containing 53 genes (Figure 4 (a)). Six genes (67%)
each edge we add the rank with which it is assigned
inthisclusterareknownBas1targetscomparedtoonly
in SCPR. Regulator-module links supported by litera-
18%Bas1targetsinSC-11.TC-28andTC-37contain70%
ture have often a higher rank. SCSR assigns Hap4, a
and 100% known targets of Msn4. These clusters have
globalregulatorofrespiratorygenes,toSC-1. Thisclus-
alargeoverlapwithSC-3andSC-41respectively,which
ter contains 58% known Hap4 targets and Hap4 has
have 55% and 93% known targets of Msn4. TC-1 con-
second highest rank in SCPR. SCSR also assigns Hap4
sists of 51 genes, out of which 28 (55%) are known to
to SC-10 which contains genes involved in amino acid
be Swi4 targets. This module merges genes from SC-
metabolism. SC-10 has only 2/37 (5%) known Hap4
10, 29 and 30. They have 4/37 (11%), 19/41 (46%) and
targets according to YEASTRACT and this assignment
8/30 (27%) Swi4 targets respectively. TC-11 contains
is ranked very low (rank 73) in SCPR. Several high-
genes of SC-8 and SC-9 whose highest ranked regula-
rankingSCPRassignmentswhichweremissedbySCSR
torisGat1(seeSectionII.D)(Figure4(b)). YEASTRACT
could also be validated using (Harbison et al., 2004)
data confirms 17% of these targets, while for SC-8 and
data (p-value < 0.005). We assign Gal80, a transcrip-
9 overall 15% targets are confirmed by YEASTRACT.
tionalregulatorinvolvedintherepressionofGalgenes
TC-35isoverrepresentedforgenesinvolvedinRNAex-
in the absence of glucose, with second rank to SC-6.
port from nucleus (p-value 10−8). It overlaps with SC-
ThisisaclusteroffourGalgenes,Gal1,Gal2,Gal7and
19, 31 and 36 (p-values ∼ 10−3). TC-31 contains genes
Gal10. Met32, a zinc-finger DNA-binding protein in-
mainlyinvolvedinribosomalbiogenesis(p-value10−13) volved in transcriptional regulation of the methionine
andcombinesrelevantgenesfromSC-13, 14and15(p- biosyntheticgenesassignedwiththirdranktoSC-8,and
values∼10−4). Gis1, a histone demethylase assigned to SC-3 with 5th
Weconcludethattightclustersimproveclusteringre- rank,aresupportedbyYEASTRACT(respectively5/29
sults obtained by (Segal et al., 2003) in two ways. They and6/31knowntargets).
canfetchonlythecoreoftightlycoexpressedgenesfrom Next we compared TCPR with SCSR to analyse the
a SC (Figure 4 (a)), or they can merge clusters which combined improvement made by ensemble averaging
5
FIG. 5 Module network inferred by (Segal et al., 2003) with edge-ranks computed by the ensemble method described in the
currentpaper.Rededgesmeanthemoduleisoverrepresentedinknowntargetsoftheconnectedregulator.
at the level of gene clustering as well as at the level of 1
regulatorassignment. ForTCPR,weselectedthetopsix
0.9
regulators for each cluster. This rank cutoff was deter-
minedasfollows. Wecomputedthesignificanceforthe 0.8
overlap between each tight cluster and each transcrip- 0.7
tion factor target set using the YEASTRACT database.
0.6
A reference module network was formed by keeping
alltranscriptionfactor-tightclusteredgesbelowacer- 0.5
tain p-value cutoff. By comparison with this reference
0.4
network we found that a rank cutoff of six gives the
0.3
bestoverallF-measurescoreatdifferent p-valuecutoffs
(seeSupplementaryinformation). Asimilaranalysisfor 0.2
SCSRshowsthattheF-measureforTCPRisconsistently
0.1
higher (see Supplementary information). To compare
TCPR and SCSR in more detail, we identified for each 0
rtaerggueltastoinrtYhEeAclSuTsRteArCwTi.thLitkheewhiisgehwesetffirnadctitohneobfesktncolwusn- YAP1(6)DAL80(1)HAP4(1)MGA1(5)GAL80(2)GAT1(2)MSN2(4)RFX1(3)STE12(2)UME6(3)GCN4(6)MGA1(6)GZF3(6)TOS8(5)PHD1(2)MET4(6)HAC1(5)GCR2(2)MET32(3)RLM1(3)MET28(2)PDR1(3)YAP5(5)CIN5(5)XBP1(3)MCM1(6)THI2(5)RGM1(4)SWI4(3)OAF1(1)PDR3(1)DAL82(4)SKO1YAP6MSN4GAL4(3)
terforeachregulatorinSCSR.Figure6showsthatTCPR FIG.6 Histogramofthehighestfractionofknowntargetsof
assignsmoreregulatorssupportedbyYEASTRACTand transcriptionfactorsinamoduleusingTCPR(red)andSCSR
alsothattheclusterscontainahigherratioofknowntar- (blue)accordingtotheYEASTRACTdatabase. Therankwith
gets.Therearesixregulatorsassignedbybothmethods, whicharegulatorisassignedtoamoduleinTCPRisindicated
four of which HAP1, GAT1, TOS8 and XBP1 all are as- inbrackets.
signedtoclustersmoreenrichedintheirknowntargets
intheTCPRsolution.
thesetofdecisiontreesforacertainmodule.InSCPR,59
regulatorassignmentsdividedover39(outof50)mod-
E. Context-specific and combinatorial regulation uleshaveasignificantscorecontribution(value > 100)
from a non-root level (see the Methods for the decom-
(Segal et al., 2003) used a decision tree approach to positionofthescorefunctionoverdifferenttreelevels).
model regulatory programs because it can represent, Combinatorial regulation means two (or more) regula-
at least in principle, context-specific and combinato- torsareconsistentlyassignedtogetheratdifferentlevels
rial regulation. In the ensemble language, context- in all decision trees. Although this form of combina-
specificitymeansaregulatorgetsahighoverallscoreby torial regulation may correspond to genuine biological
beingassignedconsistentlytoalower,non-rootlevelin combinatorialregulation, wetakeastrictlydatadriven
6
definitionhere: combinatorialregulationinthedecision genesinthesameGOcategory.
tree sense means the expression levels of both regula- Only 11/28 SC have a high-scoring regulator with
tors are needed together to explain the expression level a significant score contribution from a non-root level,
ofthemodule(‘AND’regulation).Alternatively,two(or compared to 39/50 for yeast. (Li et al., 2007) use the
more) high-scoring regulators may achieve their high co-occurrence of Neurod6 and Hey2 in the SR regula-
rankfromthesamedecisiontreelevel(usuallytheroot tion programs of SC-10, 15 and 27 to predict a cross-
level). In this case both regulators explain the module repression between Neurod6 and Hey2 with different
equallywellalone(‘OR’regulation).InSCPR,therearea modes of coregulation in different brain tissues. In
totalof100regulatorassignmentswithsignificantscore the probabilistic regulation programs (PR), Hey2 is the
contributionfromtherootlevel(ORregulation),which highest ranked regulator for SC-10, consistently as-
can be combined with the 59 assignments at level 1 for signed to the root level. However, at level 1, there
potentialANDcombinatorialregulation. arethreeequallygoodregulatorsHes5(overallrank4),
OnlyfewofthesignificantANDcombinatorialregu- Neurod6(overallrank5)andNpas4(overallrank2).For
lation pairs are present in the single-optimum solution SC-15, Neurod6isthehighestrankedregulator, consis-
of SCSR (see edge ranks in Figure 5). SC-47 has Gcn20 tently assigned to the root level, but the assignment of
as the highest ranked regulator at level 0 and Cnb1 at Hey2 at level 1 has a very low score (overall rank 4).
level 1, and both assignments are supported by liter- For SC-27, we find consistent assignments of Hey2 at
ature. SC-36 has two validated regulators Gcn20 and root level with overall rank 1 and Neurod6 at level 1
Not3rankedfirstandthirdrespectivelyinSCPR,butthe with overall rank 2. Thus the cross-repression mecha-
scoreofNot3islowandnotdeemedsignificant. SC-4is nism predicted by (Li et al., 2007) is supported only in
anexampleofORregulationwronglyassignedinSCSR. the case of SC-27 and not SC-10 and 15. This example
InSCSR,Ypl230wisassignedatlevel0andGac1atlevel underscores the usefulness of an ensemble method to
1,butinSCPRbothareassignedatlevel0withfirstand assessconfidencelevelsofpredictedinteractions, espe-
thirdrankrespectivelyandnohigh-scoringregulatoris cially in cases with limited amount of expression data
found at level 1. Some of the AND combinatorial reg- andnoothervalidationsourcesavailable.
ulationpairsinSCPRthatweremissedinSCSRcanbe
validated by YEASTRACT. SC-40 has Tos8 assigned at
level0(overallrank1)andYap1atlevel1(overallrank III. CONCLUSIONS
2). Tos8 has 3/15 known targets in this module while
Yap1 has all known targets (15/15). SC-26 has Gac1 at Wehavereexaminedthemodulenetworkmethodof
root level (overall rank 1) and Mal13 at level 1 (overall (Segal et al., 2003) and compared an ensemble-based
rank2).Mal13hastwoknowntargets(outofsixknown) strategytothestandarddirectoptimization-basedstrat-
inSC-26. egy. Ensembleaveragingselectsasubsetofmostinfor-
Duetothehighnumberofpossibleregulatorcombi- mative genes and builds a quantitatively better model
nations,identifyingstatisticallysignificantregulationof forthem. Itfindsfunctionallymorecoherenttightgene
AND-typeisanevenmorecomplexproblemthansim- clusters and is able to determine the statistically most
ple regulator assignment. These examples show that significant regulator assignments. The difficult prob-
also for this problem, the ensemble approach is well lem of identifying multiple regulators which explain
suited. together, but not separately, the expression of a mod-
ule can be addressed in a reliable way. The ensem-
ble method is thus able to deliver the promise to infer
F. Module network in mouse brain context-specific and combinatorial regulation through
theprobabilisticmodulenetworkmodel.
Recently, (Li et al., 2007) reconstructed a bHLH tran-
scription factor regulatory network in mouse brain by
adirectapplicationofthemethodof(Segaletal.,2003). IV. METHODS
Theyselectedasmalldatasetof198genesand22con-
ditions, built a module network using 22 bHLH tran- A. Bayesian two-way clustering
scription factors as candidate regulators and assigned
15differentregulatorsto28modules(denotedagainby Weassociatetoeachgene i acontinuousvaluedran-
SC),outofwhich12(43%)haveatleasttwogenesinthe domvariable X measuringthegene’sexpressionlevel.
i
sameGOcategory. Basedontheco-occurenceofregula- ForadatamatrixD = (x )withexpressionvaluesfor
im
tors in the regulation programs of individual modules, N genes in M conditions, the module network model
(Lietal.,2007)makehypothesesaboutdifferentmodes of (Segal et al., 2003) gives rise to a probabilistic model
ofcoregulationamongbraintissueswhicharecurrently for two-way clusters, where a two-way cluster k is de-
notconfirmedbyotherdatasources. Weappliedtheen- fined as a subset of genes A ⊂ {1,...,N} with a par-
k
semblemethodonthisdatasetandgot17tightclusters tition E of the set {1,...,M} into condition clusters.
k
(denotedbyTC),outofwhich11(65%)haveatleasttwo TheBayesianposteriorprobabilityforasetofcoclusters
7
(A ,E ),denotedC,isgivenby where
k k
Ppost(C) ∝ ∏ ∏ (cid:90)(cid:90) dµdτ p(µ,τ) ∏ ∏ p(xi,m | µ,τ), αE(cid:0){xr,r ∈ Rk}(cid:1) = ∏p(yt | xrt,zt,βt)
t
k E∈Ek i∈Akm∈E
with the values y determined by the unique path
t
where p(x | µ,τ) is a normal distribution with mean µ
throughthedecisiontreethatendsatleafE.
andprecisionτandp(µ,τ)isanormal-gammadistribu- Foracocluster (A ,E ) inferredfromadataset D =
k k
tion(see(Segaletal.,2005)or(Joshietal.,2008)formore (x ) by the method summarized in the previous sec-
im
details). We use the Gibbs sampler strategy developed
tion, we can derive a posterior probability function for
in(Joshietal.,2008)tosamplemultiplehigh-scoringco-
eachregulatorateachnodetasfollows. Firstnotethat
clusteringsfromthisposteriordistribution. Fromthese eachconditionmbelongstoexactlyonesetEinA and
k
multiple solutions we extract tight gene clusters using
hence determines a unique path through the decision
the procedure outlined in (Joshi et al., 2008). It consists
tree,orinotherwordsasetofvaluesy ateachnodet.
t,m
ofagraphspectralmethodextractingdenselyconnected
Furthermore,eachnodethasanassociatedconditionset
regions from the graph on the set of genes with edge-
E consisting of the union of all condition sets E which
t
weightsp ,thefrequencythatgeneiandjbelongtothe
ij canbereachedfromnodet.Hencewecandefineateach
samecoclusterineachofthesampledsolutions.
nodeaposteriorprobabilityby
(cid:16) (cid:17)
P [(r,z)] ∝max ∏ p(y | x ,z,β) , (3)
post t,m r,m
B. Probabilistic regulatory programs β m∈Et
where for computational simplicity we maximize over
ForeachsetofconditionsEintheconditionpartition
E for a given module k we have an associated normal βinsteadofmarginalizingoverapriordistribution. By
k
distribution with parameters (µ ,τ ) which can be es- allowing only a discrete set of split values, eq. (3) be-
E E
comes a discrete distribution from which it is easy to
timated from the posterior distribution. Hence such a
sample. Typically, we consider as possible split values
conditionsetcanbeinterpretedasadiscreteexpression
z the expression values x for m ∈ E , but simpler
state for the module. A regulatory program ‘predicts’ r,m t
schemes such as only allowing one or two split values
theexpressionstateofanyconditionintermsoftheex-
pressionlevelsofasmallsetR ofregulators,i.e.,there can be used to reduce computation time for large data
k
sets.
isaconditionaldistribution
The posterior probability eq. (3) measures how well
(cid:0) (cid:1)
p xi | {xr,r ∈ Rk} = p(xi | µE,τE). the expression values of a regulator ‘predict’ the parti-
tion into two sets of E induced by the condition parti-
The selection of an expression state is done by con- t
structing a decision tree with the states E ∈ E at the tionEk. Wedefinetheaveragepredictionprobabilityof
k
(r,z)atnodetbythegeometricaverage
leaves. Toeachinternalnodet,weassociatearegulator
rt and split value zt. In (Segal et al., 2003), the decision (cid:16) ∏ (cid:17)1/|Et|
at the node is based on the test xrt ≥ zt or xrt < zt. pt(r,z) = p(yt,m | xr,m,z,βmax) , (4)
Hereweextendthismodeltoallowfuzzydecisiontrees. m∈Et
Moreprecisely, wesorttheexpressionstates E ∈ Ek by whereβmaxisthemaximizerineq. (3).
their mean µ , and link this ordered set hierarchically.
E
Then we can associate to each internal node a binary
variable y = ±1, where y = −1 means ‘decrease ex- C. Regulator assignment score
t t
pression state’ (go ‘left’ in decision tree) and y = +1
t
means ‘increase expression state’ (go ‘right’ in decision ToassessthesignificanceZt(r)forassigningaregula-
tree). Again we also associate a regulator r and split torrtoanodetinacertainregulationprogram,weuse
t
valuez tonodet,andaconditionalprobability theaveragepredictionprobabilities(eq. (4))anddefine:
t
∑
1 Z (r) = w p (r,z). (5)
p(y | x ,z ,β ) = . (1) t t t
t rt t t 1+e−βtyt(xrt−zt) z
Given expression values xr for all r ∈ Rk, we traverse Atypicalchoicefortheweightfactorwt iswt = |EMt|,ex-
thedecisiontreeinaprobabilisticfashion,takingthede- pressing that we have more confidence in assignments
cision yt = ±1 at each node t by tossing a biased coin to nodes supported on more conditions. The sum ∑z
withbiaseq. (1). Theoriginalmodelwithharddecision runs over the discrete set of split values for regulator r
treesisrecoveredifβ = ±∞foreachnode. at node t. The overall significance Z(r) for assigning a
t
The conditional distribution or regulatory program regulator r to a module is defined by summing eq. (5)
nowbecomesanormalmixturedistribution over all nodesof all regulation programsfor that mod-
(cid:0) (cid:1) ∑ (cid:0) (cid:1) ule:
p x | {x ,r ∈ R } = α {x ,r ∈ R } p(x | µ ,τ )
i r k E r k i E E ∑ ∑
E∈Ek Z(r) = Zt(r).
(2) T∈T t∈T
8
D. Model evaluation References
Foranexperimentwithexpressionlevels(x ,...,x ), Alter,O.,andG.H.Golub,2005,PNAS102,17559.
1 N
wecanevaluatetheprobabilitydistribution Bansal, M., V. Belcastro, A. Ambesi-Impiombato, and
D.diBernardo,2007,MolSystBiol3,78.
Basso, K., A. A. Margolin, G. Stolovitzky, U. Klein, R. Dalla-
N
p(x ,...,x ) = ∏p(cid:0)x | {x ,r ∈ R }(cid:1), Favera,andA.Califano,2005,NatGenet37,382.
1 N i r k(i) Bonneau, R., D. J. Reiss, P. Shannon, M. Facciotti, L. Hood,
i=1
N.S.Baliga,andV.Thorsson,2006,GenomeBiol7,R36.
Bussemaker,H.J.,B.C.Foat,andL.D.Ward,2007,AnnuRev
withk(i)themoduletowhichgeneibelongsandR the
k BiophysBiomolStruct36,329.
regulatorsetofmodulek,usingtheconditionaldistribu-
Carvalho,L.E.,andC.E.Lawrence,2007,PNAS105,3209.
tions (2). We only consider genes for which the model
Chen, G., S. T. Jensen, and C. J. Stoeckert Jr, 2007, Genome
makes actual predicitions, i.e., genes belongingto clus- Biology8,R4.
ters with a regulation tree. For the cross-validation ex- Faith,J.J.,B.Hayete,J.T.Thaden,I.Mogno,J.Wierzbowski,
periment,weremoved10%oftheconditionsrandomly G. Cottarel, S. Kasif, J. J. Collins, and T. S. Gardner, 2007,
fromthetotalof173conditions.Welearnedmodulenet- PLoSBiol5,e8.
worksontheremaining90%dataandrepeatedthispro- Friedman,N.,2004,Science308,799.
cedure10times. Gardner,T.S.,andJ.J.Faith,2005,PhysLifeRev2,65.
Harbison, C. T., D. B. Gordon, T. I. Lee, N. J. Rinaldi, K. D.
Macisaac, T.W.Danford, N.M.Hannett, J.B.Tagne, D.B.
Reynolds, J. Yoo, E. G. Jennings, J. Zeitlinger, et al., 2004,
E. Data sets Nature431,99.
Hartwell, L. H., J. J. Hopfield, S. Leibler, and A. W. Murray,
Yeast expression data for 2355 differentially 1999,Nature402,C47.
expressed genes in 173 stress conditions, gene Joshi,A.,Y.VandePeer,andT.Michoel,2008,Bioinformatics
24(2),176.
clusters, their regulators, split values and re-
Lee, S., D. Pe’er, A. Dudley, G. Church, and D. Koller, 2006,
gression trees were downloaded from the sup-
Proc.Natl.Acad.Sci.U.S.A.103,14062.
plemental website of (Segal et al., 2003) at
Li, J., Z. Liu, Y. Pan, Q. Liu, X. Fu, N. Cooper, Y. Li, M. Qiu,
http://robotics.stanford.edu/∼erans/module nets/.
andT.Shi,2007,GenomeBiol8,R244.
MIPS functional catagories were downloaded from Liu, J. S., 2004, Monte Carlo strategies in scientific computing
ftp://ftpmips.gsf.de/catalogue/annotation data. For (Springer).
TC and SC we calculated the p-value whether the Novershtern, N., Z. Itzhaki, O. Manor, N. Friedman, and
overlapbetweenagivenclusterandagivenfunctional N.Kaminski,2008,AmJRespirCellMolBiol38,324.
catagory is statistically significant. We used data on Segal,E.,D.Pe’er,A.Regev,D.Koller,andN.Friedman,2005,
genome-wide binding and phylogenetically conserved JournalofMachineLearningResearch6,557.
Segal,E.,M.Shapira,A.Regev,D.Pe’er,D.Botstein,D.Koller,
motifs for 102 transcription factors from (Harbison
andN.Friedman,2003,NatGenet34,166.
etal.,2004). Foragiventranscriptionfactor,onlygenes
Segal, E., C.B.Sirlin, C.Ooi, A.S.Adler, J.Gollub, X.Chen,
that were bound with high confidence (significance
B.K.Chan,G.R.Matcuk,C.T.Barry,H.Y.Chang,andM.D.
level α = 0.005) and showed motif conservation in
Kuo,2007,NatBiotech25,675.
at least one other Saccharomyces species (besides S.
Su, A., T. Wiltshire, S. Batalov, H. Lapp, K. Ching, D. Block,
cerevisiae) were considered true targets. We also down- J.Zhang, R.Soden, M.Hayakawa, G.Kreiman, M.Cooke,
loadedallknownregulatortargetinteractionsfromthe J.Walker,etal.,2004,Proc.Natl.Acad.Sci.U.S.A.101,6062.
YEASTRACT database http://www.yeastract.com. We Teixeira, M., P. Monteiro, P. Jain, S. Tenreiro, A. Fernandes,
calculated the p-value whether the overlap between a N. Mira, M. Alenquer, A. Freitas, A. Oliveira, and I. Sa`-
given cluster and a given transcription factor target set Correia,2006,NucleicAcidsRes.34,D446.
isstatisticallysignificant. Webb-Robertson, B.-J. M., L. A. McCue, and C. E. Lawrence,
2008,PLoSComputBiol4,e1000077.
Mouseexpressiondataby(Suetal.,2004)wasdown-
Zhu,H.,H.Yang,andM.R.Owen,2007,SystSynthBiol1,171.
loadedfromhttp://wombat.gnf.organdthedataselec-
tion and normalization was done as described in (Li
etal.,2007).
Acknowledgments
WewouldliketoacknowledgeEricBonnetforuseful
discussions. AJ is supported by an Early-Stage Marie
CurieFellowship.ThisworkissupportedbyIWT(SBO-
BioFrame)andIUAPP6/25(BioMaGNet).