In aghozlane/DESeq2shaman: Differential gene expression analysis based on the negative binomial distribution

Test suite

medianPercentShrinkage <- function(res,...) {
  idx <- res$baseMean > 10
  baseMean <- res$baseMean[idx]
  qs <- quantile(baseMean, 0:10/10)
  nms <- unname( round( (qs[-1] + qs[-length(qs)])/2 ) )
  f <- cut(baseMean, qs)
  delta <- 100 * ( (res$lfcMLE - res$log2FoldChange) / res$lfcMLE )[ idx ]
  barplot(tapply(delta, f, median, na.rm=TRUE), las=2,
          names=nms, xlab="mean expression", ylab="median percent LFC shrinkage",
          ...)
}

summarizeDESeqRun <- function(x,Wald=TRUE) {
  name <- x$name
  time <- x$time
  res <- x$res
  dds <- x$dds
  cat(name,"\n")
  cat(paste(as.character(design(dds)),collapse=""),"\n")
  cat(paste(nrow(dds),"genes",ncol(dds),"samples \n"))
  cat(paste(round(unname(time[3])),"seconds \n"))
  summary(res)
  if (Wald) {
    par(mfrow=c(2,2))
    yext <- max(abs(res$log2FoldChange),na.rm=TRUE)
    plotMA(res,ylim=c(-yext,yext),main=name)
    medianPercentShrinkage(res,main=name)
  } else {
    par(mfrow=c(1,2))
  }
  plotDispEsts(dds,main=name)
  hist(res$pvalue[res$baseMean > 10],col="grey",main="p-values | base-mean > 10",xlab="")
}

library("GenomicRanges")

recount2SE <- function(name) {
  filename <- paste0(name,"_eset.RData")
  if (!file.exists(filename)) download.file(paste0(
    "http://bowtie-bio.sourceforge.net/recount/ExpressionSets/",
    filename),filename)
  load(filename)
  e <- get(paste0(name,".eset"))
  se <- SummarizedExperiment(SimpleList(counts=exprs(e)),
                             colData=DataFrame(pData(e)))
  se                   
}

library("airway")
data(airway)
dds <- DESeqDataSet(airway, ~ cell + dex)
time <- system.time({ dds <- DESeq(dds) })
res <- results(dds,addMLE=TRUE)
airwayRes <- list(name="airway", time=time, res=res, dds=dds)
rm(time, res, dds)

library("pasilla")
library("Biobase")
data("pasillaGenes")
countData <- counts(pasillaGenes)
colData <- pData(pasillaGenes)[,c("condition","type")]
dds <- DESeqDataSetFromMatrix(countData = countData,
                              colData = colData,
                              design = ~ type + condition)
dds$condition <- relevel(dds$condition, "untreated","treated")
time <- system.time({ dds <- DESeq(dds) })
res <- results(dds,addMLE=TRUE)
pasillaRes <- list(name="pasilla", time=time, res=res, dds=dds)
rm(time, res, dds)

se <- recount2SE("hammer")
se$Time[4] <- "2 months"
se$Time <- droplevels(se$Time)
dds <- DESeqDataSet(se, ~ Time + protocol)
dds$protocol <- relevel(dds$protocol, "control")
time <- system.time({ dds <- DESeq(dds) })
res <- results(dds,addMLE=TRUE)
hammerRes <- list(name="hammer", time=time, res=res, dds=dds)
rm(time, res, dds)

se <- recount2SE("bottomly")
dds <- DESeqDataSet(se, ~ strain)
time <- system.time({ dds <- DESeq(dds) })
res <- results(dds,addMLE=TRUE)
bottomlyRes <- list(name="bottomly", time=time, res=res, dds=dds)
rm(time, res, dds)

library("DESeq2")
library("parathyroidSE")
data(parathyroidGenesSE)
se <- parathyroidGenesSE
dds0 <- DESeqDataSet(se, ~ patient + treatment)
dds0 <- dds0[,dds0$treatment != "OHT" & dds0$time == "48h"]
dds <- collapseReplicates(dds0, groupby = dds0$sample, run = dds0$run)
dds$treatment <- factor(dds$treatment, levels=c("Control","DPN"))
time <- system.time({ dds <- DESeq(dds) })
res <- results(dds,addMLE=TRUE)
parathyroidRes <- list(name="parathyroid", time=time, res=res, dds=dds)
rm(time, res, dds)

library("fission")
data(fission)
dds <- DESeqDataSet(fission, ~ strain + minute + strain:minute)
time <- system.time({ dds <- DESeq(dds, test="LRT", reduced= ~ strain + minute) })
res <- results(dds)
fissionRes <- list(name="fission", time=time, res=res, dds=dds)
rm(time, res, dds)

summarizeDESeqRun(airwayRes)

summarizeDESeqRun(pasillaRes)

summarizeDESeqRun(hammerRes)

summarizeDESeqRun(bottomlyRes)

summarizeDESeqRun(parathyroidRes)

summarizeDESeqRun(fissionRes, Wald=FALSE)

gene <- rownames(fissionRes$res)[which.min(fissionRes$res$pvalue)]
data <- plotCounts(fissionRes$dds, gene, intgroup=c("minute","strain"),
                   returnData=TRUE, transform=TRUE)
library("ggplot2")
ggplot(data, aes(minute, count, color=strain, group=strain)) +
  ylab("log2 count") + geom_point() + geom_smooth(se=FALSE,method="loess") +
  ggtitle(gene)

sapply(list(airway=airwayRes, pasilla=pasillaRes, hammer=hammerRes,
            bottomly=bottomlyRes, parathyroid=parathyroidRes, fission=fissionRes),
       function(z) unname(z$time[3]))

sessionInfo()

aghozlane/DESeq2shaman documentation built on Nov. 1, 2019, 9:01 p.m.