### R code from vignette source 'birte.Rnw' ################################################### ### code chunk number 1: loading library ################################################### rm(list=ls()) library(birte) ################################################### ### code chunk number 2: EColi eSet ################################################### library(Biobase) data(EColiOxygen) EColiOxygen head(exprs(EColiOxygen)) ################################################### ### code chunk number 3: EColi ################################################### # prepare network affinities = list(TF=sapply(names(EColiNetwork$TF), function(tf){w = rep(1, length(EColiNetwork$TF[[tf]])); names(w)= EColiNetwork$TF[[tf]]; w})) affinities = simplify(affinities) affinities$other = proposeInteractions(affinities) # prepare data mydat = exprs(EColiOxygen) colnames(mydat) = make.names(paste(pData(EColiOxygen)$GenotypeVariation, pData(EColiOxygen)$GrowthProtocol, sep=".")) limmamRNA = limmaAnalysis(mydat, design=NULL, "wild.type.anaerobic - wild.type.aerobic") mydat = cbind(mydat[,colnames(mydat) =="wild.type.aerobic"], mydat[,colnames(exprs(EColiOxygen)) == "wild.type.anaerobic"]) ecoli_result = birteLimma(dat.mRNA=mydat, limmamRNA=limmamRNA, affinities=affinities, niter=500, nburnin=5000, thin=1) ################################################### ### code chunk number 4: ecoli log-lik plot dummy (eval = FALSE) ################################################### ## plotConvergence(ecoli_result, title="E. Coli") ################################################### ### code chunk number 5: ecoli log-lik plot ################################################### pdf("loglik_ecoli.pdf") plotConvergence(ecoli_result, title="E. Coli") dev.off() ################################################### ### code chunk number 6: EColi active TFs ################################################### tau = suggestThreshold(ecoli_result$post[,1]) activeTFs = rownames(ecoli_result$post)[ecoli_result$post[,1] > tau] activeTFs ################################################### ### code chunk number 7: ecoli DE genes ################################################### if(length(activeTFs) > 0){ DEgenes = rownames(limmamRNA$pvalue.tab)[limmamRNA$pvalue.tab$adj.P.Val < 0.05 & abs(limmamRNA$pvalue.tab$logFC > 1)] genesetsTF = c(sapply(affinities$TF, names), sapply(affinities$other, names)) DEgenesInTargets = sapply(genesetsTF[intersect(activeTFs, names(genesetsTF))], function(x) c(length(which(x %in% DEgenes)), length(x))) rownames(DEgenesInTargets) = c("#DEgenes", "#targets") DEgenesInTargets[,order(DEgenesInTargets["#targets",], decreasing=TRUE)] } ################################################### ### code chunk number 8: EColi predict expr ################################################### pred = birtePredict(ecoli_result, rownames(mydat)) cor(pred[[1]][[1]]$mean, limmamRNA$pvalue.tab[rownames(mydat), "logFC"]) ################################################### ### code chunk number 9: EColi TF expression ################################################### head(exprs(TFexpr)) ################################################### ### code chunk number 10: EColi TFexpr ################################################### limmaTF = limmamRNA limmaTF$pvalue.tab = limmaTF$pvalue.tab[rownames(limmaTF$pvalue.tab) %in% fData(TFexpr)$Entrez, ] names(limmaTF$lm.fit$sigma) = as.character(fData(EColiOxygen)$symbol[match(names(limmaTF$lm.fit$sigma), fData(EColiOxygen)$Entrez)]) rownames(limmaTF$pvalue.tab) = as.character(fData(EColiOxygen)$symbol[match(rownames(limmaTF$pvalue.tab), fData(EColiOxygen)$Entrez)]) diff.TF = rownames(limmaTF$pvalue.tab)[limmaTF$pvalue.tab$adj.P.Val < 0.05 & abs(limmaTF$pvalue.tab$logFC) > 1] theta.TF = rep(1/length(affinities$TF), length(affinities$TF)) names(theta.TF) = names(affinities$TF) theta.other = rep(1/length(affinities$other), length(affinities$other)) names(theta.other) = names(affinities$other) theta.other[unique(unlist(sapply(diff.TF, function(tf) grep(tf, names(theta.other)))))] = 0.5 # assume an a priori 50% activity probability for differentially expressed TFs init.TF = theta.TF init.TF = (init.TF >= 0.5)*1 init.other = theta.other init.other = (init.other >= 0.5)*1 # note that niter and nburnin are much too small in practice ecoli_TFexpr = birteLimma(dat.mRNA=mydat, data.regulators=list(TF=exprs(TFexpr)), limmamRNA=limmamRNA, limma.regulators=list(TF=limmaTF), theta.regulators=list(TF=theta.TF, other=theta.other), init.regulators=list(TF=init.TF, other=init.other), affinities=affinities, niter=500, nburnin=1000, thin=1, only.diff.TFs=TRUE) ################################################### ### code chunk number 11: EColi TFexpr result ################################################### tau = suggestThreshold(ecoli_TFexpr$post[,1]) activeTFs = ecoli_TFexpr$post[ecoli_TFexpr$post[,1] > tau,1] activeTFs ################################################### ### code chunk number 12: network inference ################################################### DEgenes = rownames(limmamRNA$pvalue.tab)[limmamRNA$pvalue.tab$adj.P.Val < 0.05 & abs(limmamRNA$pvalue.tab$logFC) > 1] net = estimateNetwork(ecoli_TFexpr, thresh=tau, de.genes=DEgenes) library(nem) if(require(Rgraphviz)){ plot(net, transitiveReduction=TRUE) } ################################################### ### code chunk number 13: assignments ################################################### net$mappos ################################################### ### code chunk number 14: network plot ################################################### if(require(Rgraphviz) & require(nem)){ pdf("nemNetwork.pdf") plot(net, transitiveReduction=TRUE) dev.off() } ################################################### ### code chunk number 15: sessionInfo ################################################### toLatex(sessionInfo())