birte.R
In birte: Bayesian Inference of Regulatory Influence on Expression (biRte)

### 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(mydat) == "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: single sample
###################################################
ss_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, single.sample = TRUE)
ss_ecoli_TFexpr[,colSums(ss_ecoli_TFexpr) > 0] # regulators that are active in at least one sample


###################################################
### code chunk number 16: sessionInfo
###################################################
toLatex(sessionInfo())