samseq_remove_batch_hack.R

library(samr)
library(limma)
library(DESeq2)

dds <- makeExampleDESeqDataSet(m=100, betaSD=rep(c(0,.5),c(900,100)))

# remove batch effect and re-exponentiate
design_matrix <- model.matrix(design(dds), colData(dds))
batch <- runif(100)
mat <- log2(counts(dds)+1)
mat_post <- removeBatchEffect(mat, covariates=cbind(batch), design=design_matrix)
x <- round(2^mat_post)

plot(rowMeans(counts(dds)), rowMeans(x), log="xy"); abline(0,1,col=2)

# SAMseq
condition <- dds$condition
samfit <- SAMseq(x, condition, resp.type="Two class unpaired", fdr.output=1)

# make container for results
res <- data.frame(gene=rownames(x), padj=rep(1,nrow(x)), row.names=rownames(x))

# pull out results from the `up` and `lo` tables
for (type in c("up","lo")) {
  num <- paste0("ngenes.",type)
  if (samfit$siggenes.table[[num]] > 0) {
    tab <- paste0("genes.",type)
    dat <- samfit$siggenes.table[[tab]]
    if (!is.matrix(dat)) { 
      dat <- matrix(dat, nrow=1)
    }
    idx <- as.numeric(dat[,"Gene Name"])
    res$padj[idx] <- 1/100 * as.numeric(dat[,"q-value(%)"])
  }
}

table(res$padj < .1)
plot(-log10(res$padj + 1e-3))