### R code from vignette source 'Glimma.Rnw'

###################################################
### code chunk number 1: style-Sweave (eval = FALSE)
###################################################
## BiocStyle::latex()


###################################################
### code chunk number 2: Rsetup
###################################################
set.seed(20161000)
options(digits=2)
options(width=100)
options(browser="false")


###################################################
### code chunk number 3: Rgetstarted
###################################################
library(Glimma)
library(limma)
library(edgeR)
data(lymphomaRNAseq)
rnaseq <- lymphomaRNAseq
rnaseq$samples$group


###################################################
### code chunk number 4: Rpreprocessing
###################################################
rnaseq <- rnaseq[rowSums(cpm(rnaseq)>1)>=3,]
rnaseq <- calcNormFactors(rnaseq)


###################################################
### code chunk number 5: Rmds
###################################################
groups <- rnaseq$samples$group
glMDSPlot(rnaseq, groups=groups)


###################################################
### code chunk number 6: Rvoomqw
###################################################
design <- model.matrix(~0+groups)
contrasts <- cbind(Smchd1null.vs.WT=c(-1,1))
vm <- voomWithQualityWeights(rnaseq, design=design)
fit <- lmFit(vm, design=design)
fit <- contrasts.fit(fit, contrasts)
fit <- eBayes(fit)
dt <- decideTests(fit)
summary(dt)


###################################################
### code chunk number 7: Rmd
###################################################
glMDPlot(fit, status=dt, counts=vm, groups=groups, side.main="Symbols")


###################################################
### code chunk number 8: Rxy
###################################################
glXYPlot(x=fit$coef, y=fit$lod, xlab="logFC", ylab="logodds",
  status=dt, counts=vm, groups=groups, anno=fit$genes)


###################################################
### code chunk number 9: Rmdsmatrix
###################################################
lcpm <- cpm(rnaseq, log=TRUE, normalized.lib.sizes=TRUE)
glMDSPlot(lcpm, groups=groups)


###################################################
### code chunk number 10: Rmdscheckinput (eval = FALSE)
###################################################
## # Check that MDS plot runs on EList object: PASS
## glMDSPlot(vm, groups=groups)
## # Check that MDS plot runs on DESeqDataSet object: FAIL regarding groups
## library(DESeq2)
## rnaseq.deseq2 <- DESeqDataSetFromMatrix(rnaseq$counts, colData=rnaseq$samples, design=~group)
## glMDSPlot(rnaseq.deseq2, groups=groups)


###################################################
### code chunk number 11: Rmdsgroups
###################################################
groups.df <- as.data.frame(cbind(
  genotype=as.character(groups),
  lane=c(rep(4,6),3)))
groups.df
glMDSPlot(lcpm, groups=groups.df)


###################################################
### code chunk number 12: Rmdsimple
###################################################
glMDPlot(fit)


###################################################
### code chunk number 13: Rmdage
###################################################
groups.age <- runif(ncol(rnaseq), min=3, max=12)
groups.age


###################################################
### code chunk number 14: Rmdsamples
###################################################
cols <- c("yellow","blue","magenta")
sample.cols <- c("limegreen", "purple")[groups]
glMDPlot(fit, status=dt, counts=vm, groups=groups.age,
  sample.cols=sample.cols, cols=cols,
  side.ylab="logCPM", side.xlab="Age (in months)",
  side.main="Symbols", main=colnames(dt))


###################################################
### code chunk number 15: Rmdcheckinput (eval = FALSE)
###################################################
## # limma - PASS
## glMDPlot(fit)
## # edgeR - PASS
## rnaseq.edger <- estimateDisp(rnaseq, design=design)
## fit.edger <- exactTest(rnaseq.edger)
## glMDPlot(fit.edger)
## # DESeq2 - FAIL regarding annotation
## library(DESeq2)
## rnaseq.deseq2 <- DESeqDataSetFromMatrix(rnaseq$counts, colData=rnaseq$samples, design=~group)
## mcols(rnaseq.deseq2) <- DataFrame(mcols(rnaseq.deseq2), rnaseq$genes)
## fit.deseq2 <- DESeq(rnaseq.deseq2)
## glMDPlot(fit.deseq2)


###################################################
### code chunk number 16: Rmdanno
###################################################
ID <- paste(fit$genes$Symbols, fit$genes$GeneID)
DE <- c("downregulated", "notDE", "upregulated")[as.factor(dt)]
anno <- as.data.frame(cbind(ID, DE))
head(anno)
glMDPlot(fit, counts=vm, groups=groups, side.main="ID",
  anno=anno, display.columns=c("ID", "GeneName", "DE"))


###################################################
### code chunk number 17: Rxyvolcanosimple
###################################################
glXYPlot(x=fit$coef, y=fit$lod)


###################################################
### code chunk number 18: Rxyvolcanomore
###################################################
glXYPlot(x=fit$coef, y=fit$lod, xlab="logFC", ylab="logodds",
  status=dt, anno=anno, side.main="ID",
  counts=vm, groups=groups, sample.cols=sample.cols)


###################################################
### code chunk number 19: Rarraydata
###################################################
data(arraydata)
arrays <- arraydata$arrays
targets <- arraydata$targets
dim(arrays)
targets


###################################################
### code chunk number 20: Rmds.arrays
###################################################
glMDSPlot(arrays, groups=targets[,c("Condition", "Chip", "Experiment")])


###################################################
### code chunk number 21: Rarraysfit
###################################################
design <- model.matrix(~0+targets$Condition+as.factor(targets$Experiment))
contrasts <- cbind(
  DP_Ezh2KO.vs.WT=c(1,-1,0,0,0,0),
  Lum_Ezh2KO.vs.WT=c(0,0,1,-1,0,0))
fit <- lmFit(arrays, design)
fit <- contrasts.fit(fit, contrasts)
fit <- eBayes(fit)
dt <- decideTests(fit, p.value=0.1)
summary(dt)


###################################################
### code chunk number 22: Rarraysmd
###################################################
sample.cols <- c("purple", "magenta", "green")[targets$Experiment]
for (COEF in 1:2) {
  glMDPlot(fit, status=dt, coef=COEF, main=colnames(fit)[COEF],
    counts=arrays, groups=targets$Condition, sample.cols=sample.cols,
    side.ylab="Log-expression", side.main="ProbeID")
}


###################################################
### code chunk number 23: Rarraysxy
###################################################
dt2 <- rep(0, nrow(dt))
dt2[rowSums(dt!=0)==1] <- -1
dt2[rowSums(dt!=0)==2] <- 1
table(dt2)
cols <- c("black", "grey", "red")
glXYPlot(fit$coef[,1], y=fit$coef[,2], xlab="DP", ylab="Lum",
  status=dt2, cols=cols, anno=fit$genes, side.main="ProbeID",
  counts=arrays, groups=targets$Condition, sample.cols=sample.cols,
  side.ylab="Log-expression", main="logFCs")


###################################################
### code chunk number 24: Rmdedgeret (eval = FALSE)
###################################################
## groups <- rnaseq$samples$group
## design <- model.matrix(~groups)
## colnames(design) <- c("WT", "Smchd1null.vs.WT")
## rnaseq.edger <- estimateDisp(rnaseq, design=design)
## fit.edger <- exactTest(rnaseq.edger)
## dt.edger <- decideTestsDGE(fit.edger)
## glMDPlot(fit.edger, status=dt.edger, counts=rnaseq, groups=groups, transform=TRUE)


###################################################
### code chunk number 25: Rmdedgerlrt (eval = FALSE)
###################################################
## fit.edger <- glmFit(rnaseq.edger, design)
## fit.edger <- glmLRT(fit.edger)
## dt.edger <- decideTestsDGE(fit.edger)
## glMDPlot(fit.edger, status=dt.edger, counts=rnaseq, groups=groups, transform=TRUE)


###################################################
### code chunk number 26: Rmddeseq2 (eval = FALSE)
###################################################
## # BUG regarding the scale of sample expression
## library(DESeq2)
## rnaseq.deseq2 <- DESeqDataSetFromMatrix(
##   rnaseq$counts, colData=rnaseq$samples, design=~group)
## mcols(rnaseq.deseq2) <- DataFrame(mcols(rnaseq.deseq2), rnaseq$genes)
## rnaseq.deseq2 <- DESeq(rnaseq.deseq2)
## fit.deseq2 <- results(rnaseq.deseq2, contrast=c("group", "Smchd1-null", "WT"))
## dt.deseq2 <- as.numeric(fit.deseq2$padj<0.05)
## glMDPlot(fit.deseq2, status=dt.deseq2, counts=rnaseq, groups=groups, transform=FALSE,
##   samples=colnames(rnaseq), anno=rnaseq$genes)


###################################################
### code chunk number 27: Rsessinfo
###################################################
sessionInfo()