birta.R
In birta: Bayesian Inference of Regulation of Transcriptional Activity

### R code from vignette source 'birta.Rnw'

###################################################
### code chunk number 1: loading library
###################################################
library(birta)


###################################################
### code chunk number 2: load simulated data
###################################################
data(humanSim)
str(head(genesets$TF))
str(head(genesets$miRNA))
head(sim$dat.mRNA)


###################################################
### code chunk number 3: simulation limma
###################################################
design = model.matrix(~0+factor(c(rep("control", 5), rep("treated", 5))))
colnames(design) = c("control", "treated")
contrasts = "treated - control"
limmamRNA = limmaAnalysis(sim$dat.mRNA, design, contrasts)
limmamiRNA = limmaAnalysis(sim$dat.miRNA, design, contrasts)


###################################################
### code chunk number 4: simulation MCMC
###################################################
sim_result = birta(sim$dat.mRNA, sim$dat.miRNA, limmamRNA=limmamRNA, 
limmamiRNA=limmamiRNA, nrep=c(5,5,5,5), genesets=genesets, 
model="all-plug-in", niter=50000, nburnin=10000, 
sample.weights=FALSE, potential_swaps=potential_swaps)


###################################################
### code chunk number 5: simulation log-lik plot dummy (eval = FALSE)
###################################################
## plotConvergence(sim_result, nburnin=10000, title="simulation")


###################################################
### code chunk number 6: simulation log-lik plot
###################################################
pdf("loglik_sim.pdf")
plotConvergence(sim_result, nburnin=10000, title="simulation")
dev.off()	


###################################################
### code chunk number 7: simulation result
###################################################
sim$TFstates[sim$TFstates == 1]
sim_result$TFActivitySwitch[sim_result$TFActivitySwitch > 0]
sim$miRNAstates[sim$miRNAstates == 1]
sim_result$miRNAstates1[sim_result$miRNAstates1 > 0]
sim_result$miRNAstates2[sim_result$miRNAstates2 > 0]


###################################################
### code chunk number 8: Only miRNA-target graph
###################################################
genesetsmiRNA = genesets["miRNA"]
result_miRonly = birta(sim$dat.mRNA, sim$dat.miRNA, limmamRNA=limmamRNA, 
limmamiRNA=limmamiRNA, nrep=c(5,5,5,5), genesets=genesetsmiRNA, 
model="all-plug-in", niter=50000, nburnin=10000, 
sample.weights=FALSE, potential_swaps=potential_swaps)


###################################################
### code chunk number 9: Only miRNA-target graph result
###################################################
result_miRonly$miRNAstates1[result_miRonly$miRNAstates1 > 0]
result_miRonly$miRNAstates2[result_miRonly$miRNAstates2 > 0]


###################################################
### code chunk number 10: EColi eSet
###################################################
data(EColiOxygen)
EColiOxygen
head(exprs(EColiOxygen))


###################################################
### code chunk number 11: EColi
###################################################
design = model.matrix(~0+factor(pData(EColiOxygen)$GrowthProtocol))
colnames(design) = c("aerobic.growth", "anaerobic.growth")
contrasts = "anaerobic.growth - aerobic.growth"
limmamRNA = limmaAnalysis(EColiOxygen, design, contrasts)
ecoli_result = birta(EColiOxygen, nrep=c(0, 0, 3, 4), 
genesets=EColiNetwork, limmamRNA=limmamRNA, 
model="all-plug-in", niter=50000, nburnin=10000, 
sample.weights=FALSE)


###################################################
### code chunk number 12: ecoli log-lik plot dummy (eval = FALSE)
###################################################
## plotConvergence(ecoli_result, nburnin=10000, title="E. Coli")


###################################################
### code chunk number 13: ecoli log-lik plot
###################################################
pdf("loglik_ecoli.pdf")
plotConvergence(ecoli_result, nburnin=10000, title="E. Coli")
dev.off()	


###################################################
### code chunk number 14: EColi active TFs
###################################################
activeTFs = ecoli_result$TFActivitySwitch[ecoli_result$TFActivitySwitch > 0.9]
sort(activeTFs)


###################################################
### code chunk number 15: ecoli DE genes
###################################################
DEgenes = limmamRNA$pvalue.tab$ID[limmamRNA$pvalue.tab$adj.P.Val < 0.05]
DEgenesInTargets = sapply(ecoli_result$genesetsTF[names(activeTFs)], 
function(x) c(length(which(x %in% DEgenes)), length(x)))
rownames(DEgenesInTargets) = c("#DEgenes", "#targets")
DEgenesInTargets[,order(DEgenesInTargets["#targets",], decreasing=T)]


###################################################
### code chunk number 16: EColi TF expression
###################################################
head(exprs(TFexpr))


###################################################
### code chunk number 17: EColi TFexpr
###################################################
limmaTF = limmaAnalysis(TFexpr, design, contrasts)
ecoli_TFexpr = birta(EColiOxygen, nrep=c(0, 0, 3, 4), 
genesets=EColiNetwork, TFexpr=TFexpr, limmaTF=limmaTF, 
limmamRNA=limmamRNA, model="all-plug-in", niter=50000, 
nburnin=10000, sample.weights=FALSE)


###################################################
### code chunk number 18: EColi TFexpr result
###################################################
sort(ecoli_TFexpr$TFstates1[ecoli_TFexpr$TFstates1 > 0.9])
sort(ecoli_TFexpr$TFstates2[ecoli_TFexpr$TFstates2 > 0.9])
sort(ecoli_TFexpr$TFActivitySwitch[ecoli_TFexpr$TFActivitySwitch > 0.9])


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