apeglm code for ASE using ashr / locfdr

apeglm.R

library(apeglm)
suppressPackageStartupMessages(library(SummarizedExperiment))
load("se.oracle.rda")
n <- 10
cts <- assay(wide)[,1:n] + assay(wide)[,(n+1):(2*n)]
dim(cts)
ase_cts <- assay(wide)[,1:n]
mode(ase_cts) <- "integer"

theta_hat <- matrix(nrow=nrow(cts), ncol=niter+1)
theta_hat[,1] <- rep(100,nrow(cts)) # rough initial estimate of dispersion
x <- matrix(rep(1,n),ncol=1)
# 8 seconds
niter <- 2
system.time({
  for (i in 1:niter) {
    param <- cbind(theta_hat[,i], cts)
    fit_mle <- apeglm(Y=ase_cts, x=x, log.lik=NULL, param=param,
                      no.shrink=TRUE, log.link=FALSE, method="betabinC")
    theta_hat[,i+1] <- bbEstDisp(success=ase_cts, size=cts,
                                 x=x, beta=fit_mle$map,
                                 minDisp=.01, maxDisp=500)
  }
})
plot(theta_hat[,2:3], log="xy")

# final parameters
theta_hat_final <- theta_hat[,niter+1]
log10cts <- log10(rowSums(cts) + 1)
#plot(log10cts, theta_hat_final, log="y")

param <- cbind(theta_hat_final, cts)

# 3 seconds
system.time({
  fit_mle <- apeglm(Y=ase_cts, x=x, log.lik=NULL, param=param,
                    no.shrink=TRUE, log.link=FALSE, method="betabinCR")
})

coef <- 1
mle <- cbind(fit_mle$map[,coef], fit_mle$sd[,coef])

### with apeglm... doesn't work well

system.time({
  fit <- apeglm(Y=ase_cts, x=x, log.lik=NULL, param=param, 
                coef=coef, mle=mle, log.link=FALSE, method="betabinCR")
})

qvalue <- fit$svalue 
#cols <- ifelse(qvalue < .01, "dodgerblue", "black")
#plot(log10cts, fit$map[,coef], main="apeglm", ylab="log odds", col=cols, ylim=c(-1,1))
#abline(h=0,col=rgb(1,0,0,.75),lwd=2)
table(qvalue < .01)

### with ashr

library(ashr)
fit <- ash(mle[,1], mle[,2], mixcompdist="uniform", df=n-2, method="fdr")
qvalue <- fit$result[,"qvalue"] # fix this variable name later

### with locfdr

library(locfdr)
fit <- locfdr(mle[,1]/mle[,2], df=n-2, nulltype=1)
qvalue <- fit$fdr # fix this variable name later
  
# eval

abnd <- sapply(strsplit(rownames(cts), split="-"), `[`, 2)
table(abnd == 2)
de <- abnd != 2
target <- c(.01, .05, .1)
names(target) <- target
res <- sapply(target, function(t) {
  sig <- qvalue < t
  c(sens=prop.table(table(de, sig), 1)[2,2],
    fdr=prop.table(table(de, sig), 2)[1,2])
})
#
res