plots.Rmd

# Bottomly plots

```{r echo=FALSE, message=FALSE}
# data and scripts come from the Bottomly analysis in the DESeq2 paper.
#
# paper:
# https://genomebiology.biomedcentral.com/articles/10.1186/s13059-014-0550-8
# data and scripts:
# http://www-huber.embl.de/DESeq2paper/
#
# data downloaded:
# random_subsets.txt
# bottomly_sumexp.RData
#
# scripts downloaded:
# diffExpr.R 
# runScripts.R (edgeR code updated for 2016)
#
# quick look at the random subsets:
randomSubsets <- read.table("random_subsets.txt",strings=FALSE)
randomSubsets <- as.matrix(randomSubsets)
library(Biobase)
library(SummarizedExperiment)
load("bottomly_sumexp.RData")
bottomly <- updateObject(bottomly)
strain <- colData(bottomly)[,"strain",drop=FALSE]
exper <- colData(bottomly)[,"experiment.number",drop=FALSE]
exper[,1] <- factor(exper[,1])
```

### Condition and batch in test and heldout

```{r cond_batch}
# plot the test and heldout sets coloring condition and batch
library(rafalib)
bigpar(2,2)
cols <- c("orange","purple","dodgerblue")
image(sapply(1:30, function(i) as.integer(strain[randomSubsets[i,1:6],])),
      col=cols, main="test cond")
image(sapply(1:30, function(i) as.integer(exper[randomSubsets[i,1:6],])),
      col=cols, main="test batch")
image(sapply(1:30, function(i) as.integer(strain[randomSubsets[i,7:21],])),
      col=cols, main="out cond")
image(sapply(1:30, function(i) as.integer(exper[randomSubsets[i,7:21],])),
      col=cols, main="out batch")
```

```{r echo=FALSE}
library(rafalib)
# first run diffExpr.R and save result
# load the data from DE calling
load("sensFDR.rda")
# define some functions for compiling results
getTestCalls <- function(alpha) {
  t(sapply(1:nreps, function(i) sapply(namesAlgos, function(algo) {
    sum((resTest[[i]][[algo]] < alpha))
  })))
}
getHeldoutCalls <- function(alpha) {
  t(sapply(1:nreps, function(i) sapply(namesAlgos, function(algo) {
    sum((resHeldout[[i]][[algo]] < alpha))
  })))
}
getSensitivity <- function(alpha) {
  t(sapply(1:nreps, function(i) sapply(namesAlgos, function(algo) {
    sigHeldout <- resHeldout[[i]][[algo]] < alpha
    if (sum(sigHeldout) == 0) return(0)
    mean((resTest[[i]][[algo]] < alpha)[sigHeldout])
  })))
}
getFDR <- function(alpha) {
  t(sapply(1:nreps, function(i) sapply(namesAlgos, function(algo) {
    sigTest <- resTest[[i]][[algo]] < alpha
    if (sum(sigTest) == 0) return(0)
    mean((resHeldout[[i]][[algo]] > alpha)[sigTest])
  })))
}
nreps <- length(resTest)
test <- getTestCalls(.1)
held <- getHeldoutCalls(.1)
sens <- getSensitivity(.1)
fdr <- getFDR(.1)
```

### Number of calls

```{r num_calls}
bigpar(1,2,mar=c(10,5,3,1))
boxplot(test, las=2, ylim=c(0,5000), main="n=3 #pos")
boxplot(held, las=2, ylim=c(0,5000), main="n=7/8 #pos")
```

Mean of calls in test set (n=3) and heldout set (n=7/8)

```{r}
round(colMeans(test))
round(colMeans(held))
```

### FDR and sensitivity against heldout

```{r fdr_sens}
lines <- function() abline(h=0:10/10,col=rgb(0,0,0,.1))
bigpar(1,2,mar=c(10,5,3,1))
boxplot(fdr, las=2, ylim=c(0,0.5), main="rough est. FDR")
lines()
boxplot(sens, las=2, ylim=c(0,0.5), main="sensitivity")
lines()
```

Mean of the rough estimate of FDR (%) using the 7 vs 8 heldout set as ground truth:

```{r}
100 * round(colMeans(fdr),3)
```

```{r echo=FALSE}
# examine overlap for methods for test and heldout sets
getOverlap <- function(a, b, res, alpha) {
  out <- sapply(1:nreps, function(i) {
    a.padj <- res[[i]][[a]]
    b.padj <- res[[i]][[b]]
    over <- sum(a.padj < alpha & b.padj < alpha)
    c(over/sum(a.padj < alpha), over/sum(b.padj < alpha))
  })
  out <- t(out)
  colnames(out) <- c(a,b)
  out
}
```

### Overlap of method pairs in test

```{r over_test}
bigpar(3,2,mar=c(5,5,1,1))
ylims <- c(0, 1)
boxplot(getOverlap("DESeq2","edgeR",resTest,.1), ylim=ylims)
boxplot(getOverlap("DESeq2","edgeRQL",resTest,.1), ylim=ylims)
boxplot(getOverlap("DESeq2","limma.voom",resTest,.1), ylim=ylims)
boxplot(getOverlap("edgeR","edgeRQL",resTest,.1), ylim=ylims)
boxplot(getOverlap("edgeR","limma.voom",resTest,.1), ylim=ylims)
boxplot(getOverlap("edgeRQL","limma.voom",resTest,.1), ylim=ylims)
```

### Overlap of method pairs in heldout

```{r over_heldout}
bigpar(3,2,mar=c(5,5,1,1))
ylims <- c(0.6, 1)
boxplot(getOverlap("DESeq2","edgeR",resHeldout,.1), ylim=ylims)
boxplot(getOverlap("DESeq2","edgeRQL",resHeldout,.1), ylim=ylims)
boxplot(getOverlap("DESeq2","limma.voom",resHeldout,.1), ylim=ylims)
boxplot(getOverlap("edgeR","edgeRQL",resHeldout,.1), ylim=ylims)
boxplot(getOverlap("edgeR","limma.voom",resHeldout,.1), ylim=ylims)
boxplot(getOverlap("edgeRQL","limma.voom",resHeldout,.1), ylim=ylims)
```

Which genes called by DESeq2 but not be edgeR in the heldout set? 
Maybe we can characterize these. We'll just look at the first iteration,
and look at the top genes called exclusively by DESeq2.

```{r message=FALSE}
deseq2.padj <- resHeldout[[1]][["DESeq2"]]
edger.padj <- resHeldout[[1]][["edgeR"]]
deseq2.exclusive <- deseq2.padj < 0.1 & edger.padj > 0.1
table(deseq2.exclusive)
names(deseq2.padj) <- rownames(bottomly)
names(edger.padj) <- rownames(bottomly)
exc.padj <- deseq2.padj[deseq2.exclusive]
top.padj <- head(sort(exc.padj),4)
edger.padj[names(top.padj)] # top were likely filtered by edgeR
table(edger.padj[names(exc.padj)] == 1) # how many likely filtered?
library(DESeq2)
dds <- DESeqDataSet(bottomly, ~strain)
dds <- estimateSizeFactors(dds)
heldout.idx <- randomSubsets[1,7:21]
```

It appears these didn't pass the CPM filter we used for edgeR 
of 10/L in 3 or more samples,
where L is the number of millions of reads in the smallest library.
Top row is the heldout set on which the samples were compared,
bottom row is adding in all the samples.

```{r, deseq2.counts}
par(mfcol=c(2,4))
for (i in 1:4) {
  gene <- names(top.padj)[i]
  cts <- counts(dds, normalized=TRUE)[gene,]
  stripchart(cts[heldout.idx] ~ strain[heldout.idx,], 
             method="jitter", vertical=TRUE, ylab="", 
             main=substr(gene,14,18))
  stripchart(cts ~ dds$strain, 
             method="jitter", vertical=TRUE, ylab="", 
             main=substr(gene,14,18))
}
```