### 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())