## ----style, echo = FALSE, results = 'asis'---------------------------------
knitr::opts_chunk$set(
    eval=as.logical(Sys.getenv("KNITR_EVAL", "TRUE")),
    cache=as.logical(Sys.getenv("KNITR_CACHE", "TRUE")))

## ----setup, echo=FALSE-----------------------------------------------------
suppressPackageStartupMessages({
    library(DESeq2)
    library(gplots)
    library(ggplot2)
    library(RColorBrewer)
    library(org.Hs.eg.db)
})

## --------------------------------------------------------------------------
library(airway)
data("airway")
se <- airway

## --------------------------------------------------------------------------
head(assay(se))

## --------------------------------------------------------------------------
colSums(assay(se))

## --------------------------------------------------------------------------
colData(se)

## ----rowRanges-------------------------------------------------------------
rowRanges(se)

## --------------------------------------------------------------------------
library("DESeq2")
dds <- DESeqDataSet(se, design = ~ cell + dex)

## --------------------------------------------------------------------------
dds$dex <- relevel(dds$dex, "untrt")

## --------------------------------------------------------------------------
dds <- DESeq(dds)

## --------------------------------------------------------------------------
(res <- results(dds))

## --------------------------------------------------------------------------
mcols(res, use.names=TRUE)

## --------------------------------------------------------------------------
summary(res)

## --------------------------------------------------------------------------
sum(res$pvalue < 0.05, na.rm=TRUE)
sum(!is.na(res$pvalue))

## --------------------------------------------------------------------------
sum(res$padj < 0.1, na.rm=TRUE)

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

## --------------------------------------------------------------------------
head(resSig[ order( -resSig$log2FoldChange ), ])

## ----plotcounts, fig.width=5, fig.height=5---------------------------------
topGene <- rownames(res)[which.min(res$padj)]
data <- plotCounts(dds, gene=topGene, intgroup=c("dex"), returnData=TRUE)

## ----ggplotcountsdot, fig.height=5-----------------------------------------
library(ggplot2)
ggplot(data, aes(x=dex, y=count, fill=dex)) +
  scale_y_log10() + 
  geom_dotplot(binaxis="y", stackdir="center")

## ----plotma, eval=FALSE----------------------------------------------------
#  plotMA(res, ylim=c(-5,5))

## ----plotma2, eval=FALSE---------------------------------------------------
#  plotMA(res, ylim=c(-5,5))
#  with(res[topGene, ], {
#    points(baseMean, log2FoldChange, col="dodgerblue", cex=2, lwd=2)
#    text(baseMean, log2FoldChange, topGene, pos=2, col="dodgerblue")
#  })

## ----plotdispests----------------------------------------------------------
plotDispEsts(dds)

## ----histpvalue------------------------------------------------------------
hist(res$pvalue, breaks=20, col="grey50", border="white")

## ----histpvalue2-----------------------------------------------------------
hist(res$pvalue[res$baseMean > 1], breaks=20, col="grey50", border="white")

## ----sensitivityovermean, fig.height=4-------------------------------------
# create bins using the quantile function
qs <- c(0, quantile(res$baseMean[res$baseMean > 0], 0:7/7))
# cut the genes into the bins
bins <- cut(res$baseMean, qs)
# rename the levels of the bins using the middle point
levels(bins) <- paste0("~",round(.5*qs[-1] + .5*qs[-length(qs)]))
# calculate the ratio of $p$ values less than .01 for each bin
ratios <- tapply(res$pvalue, bins, function(p) mean(p < .01, na.rm=TRUE))
# plot these ratios
barplot(ratios, xlab="mean normalized count", ylab="ratio of small p values")

## --------------------------------------------------------------------------
library(org.Hs.eg.db)

## --------------------------------------------------------------------------
columns(org.Hs.eg.db)
res$hgnc_symbol <- 
    unname(mapIds(org.Hs.eg.db, rownames(res), "SYMBOL", "ENSEMBL"))
res$entrezgene <- 
    unname(mapIds(org.Hs.eg.db, rownames(res), "ENTREZID", "ENSEMBL"))

## --------------------------------------------------------------------------
resOrdered <- res[order(res$pvalue),]
head(resOrdered)

## ----eval=FALSE------------------------------------------------------------
#  write.csv(as.data.frame(resOrdered), file="results.csv")

## --------------------------------------------------------------------------
sessionInfo()