adding limma lab

course4/using_limma.Rmd

@@ -0,0 +1,45 @@
+# Using limma for microarray analysis
+
+```{r}
+biocLite("SpikeInSubset")
+library(SpikeInSubset)
+data(rma95)
+library(genefilter)
+fac <- factor(rep(1:2,each=3))
+tt <- rowttests(exprs(rma95),fac)
+mask <- with(tt, abs(dm) < .2 & p.value < .01)
+spike <- rownames(rma95) %in% colnames(pData(rma95))
+cols <- ifelse(mask,"red",ifelse(spike,"dodgerblue","black"))
+
+with(tt, plot(-dm, -log10(p.value), cex=.8, pch=16,
+     xlim=c(-1,1), ylim=c(0,5),
+     xlab="difference in means",
+     col=cols))
+abline(h=2,v=c(-.2,.2), lty=2)
+
+tt$s <- apply(exprs(rma95), 1, function(row) sqrt(.5 * (var(row[1:3]) + var(row[4:6]))))
+with(tt, plot(s, -log10(p.value), cex=.8, pch=16,
+              log="x",xlab="estimate of standard deviation",
+              col=cols))
+
+library(limma)
+fit <- lmFit(rma95, model.matrix(~ fac))
+ebfit <- ebayes(fit)
+limmares <- data.frame(dm=coef(fit)[,"fac2"], p.value=ebfit$p.value[,"fac2"])
+with(limmares, plot(dm, -log10(p.value),cex=.8, pch=16,
+     col=cols,xlab="difference in means",
+     xlim=c(-1,1), ylim=c(0,5)))
+abline(h=2,v=c(-.2,.2), lty=2)
+
+
+n <- 40
+qs <- seq(from=0,to=.2,length=n)
+idx <- sapply(seq_len(n),function(i) which(as.integer(cut(tt$s^2,qs)) == i)[1])
+idx <- idx[!is.na(idx)]
+par(mar=c(5,5,2,2))
+plot(1,1,xlim=c(0,.21),ylim=c(0,1),type="n",
+     xlab="variance estimates",ylab="",yaxt="n")
+axis(2,at=c(.1,.9),c("before","after"),las=2)
+segments((tt$s^2)[idx],rep(.1,n),
+         ebfit$s2.post[idx],rep(.9,n))
+```