BADER.R
In BADER: Bayesian Analysis of Differential Expression in RNA Sequencing Data

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

###################################################
### code chunk number 1: loaddata
###################################################
library(BADER)
datafile <- system.file( "extdata/pasilla_gene_counts.tsv", package="pasilla" )
pasillaCountTable <- read.table( datafile, header=TRUE, row.names=1 )
head(pasillaCountTable)


###################################################
### code chunk number 2: splitup
###################################################
genetable <- pasillaCountTable[,c(3,4,6,7)]
genetable <- genetable[rowSums(genetable) > 0,]
genetable <- genetable[1:500,]
design <- factor(c("A","A","B","B"))


###################################################
### code chunk number 3: runMCMC
###################################################
set.seed(21)
results <- BADER(genetable,design,burn=1e3,reps=1e3,printEvery=1e3)


###################################################
### code chunk number 4: printResults
###################################################
names(results)
rownames(genetable)[order(results$diffProb,decreasing=TRUE)[1:10]]


###################################################
### code chunk number 5: toyDataset
###################################################
set.seed(21)
gam <- c(rnorm(50,0,2),rep(0,450))
muA <- rnorm(500,3,1)
muB <- muA + gam
kA <- t(matrix(rnbinom(2*500,mu=exp(muA),size=exp(1)),nrow=2,byrow=TRUE))
kB <- t(matrix(rnbinom(2*500,mu=exp(muB),size=exp(1)),nrow=2,byrow=TRUE))
genetable <- cbind(kA,kB)
design <- factor(c("A","A","B","B"))

results <- BADER(genetable,design,mode="full",burn=1e3,reps=3e3,printEvery=1e3)


###################################################
### code chunk number 6: postPlot
###################################################
par(mfrow=c(1,2))
temp <- hist(results$logFoldChangeDist[,18],breaks=seq(-10.125,10.125,by=0.25),plot=FALSE)
temp$density <- temp$density/sum(temp$density)
plot(temp, freq=FALSE,ylab="probability",xlab="log fold change",main="gene no. 18",xlim=c(-1,5))
abline(v=gam[18],col="#4DAF4A", lwd=2)
abline(v=results$logFoldChange[18],col="#377EB8",lwd=2)
temp <- hist(results$logFoldChangeDist[,142],breaks=seq(-10.125,10.125,by=0.25),plot=FALSE)
temp$density <- temp$density/sum(temp$density)
plot(temp, freq=FALSE,ylab="probability",xlab="log fold change",main="gene no. 142",xlim=c(-4,1))
abline(v=gam[142],col="#4DAF4A", lwd=2)
abline(v=results$logFoldChange[142],col="#377EB8",lwd=2)


###################################################
### code chunk number 7: gse
###################################################
S <- c(1,2,3,51,52,53)
T <- c(54,55,56)

f <- function(sample,set)
    mean(sample[set] != 0) > mean(sample[-set] != 0)

mean(apply(results$logFoldChangeDist,1,f,set=S))
mean(apply(results$logFoldChangeDist,1,f,set=T))