## ----setup, echo=FALSE---------------------------------------------------
knitr::opts_chunk$set(cache=TRUE)


## ----deseq---------------------------------------------------------------
require("DESeq2")
require("parathyroidSE")
data("parathyroidGenesSE")
se <- parathyroidGenesSE
colnames(se) <- se$run


## ----deseq-prep----------------------------------------------------------
ddsFull <- DESeqDataSet(se, design = ~ patient + treatment)
ddsCollapsed <-
    collapseReplicates(ddsFull, groupby = ddsFull$sample, 
                       run = ddsFull$run)
dds <- ddsCollapsed[, ddsCollapsed$time == "48h"]
dds$time <- droplevels(dds$time)
dds$treatment <- relevel(dds$treatment, "Control")


## ----deseq-analysis------------------------------------------------------
dds <- DESeq(dds)
res <- results(dds, contrast = c("treatment", "DPN", "Control"))


## ----deseq-interpret-----------------------------------------------------
resSig <- res[which(res$padj < 0.1 ),]
head(resSig[order(resSig$log2FoldChange),])


## ----deseq-MA------------------------------------------------------------
plotMA(res, ylim = c(-1, 1))


## ----deseq-dispest-------------------------------------------------------
plotDispEsts(dds, ylim = c(1e-6, 1e1))


## ----deseq-hist----------------------------------------------------------
hist(res$pvalue, breaks=40, col="grey")