Source

Limma Example / limmaExample.Rnw

Full commit
\documentclass{article}

\begin{document}
\SweaveOpts{concordance=TRUE}

\section{Get Data from NCBI GEO}

I simply need an example dataset similar to that you have.  Just a quick look brings up GSE16200 as something that might be simple to play with.  Normally, you would not need these lines and would start with the Analysis section.

<<loadData>>=
library(GEOquery)
getGEOSuppFiles('GSE16200')
untar('GSE16200/GSE16200_RAW.tar')
@

<<loadAnnotation>>=
gse = getGEO('GSE16200')[[1]]
pd = pData(gse)
knockdown = pd$characteristics_ch1!="sirna: control"
collagen  = pd$characteristics_ch1.1!="plate: plastic"
newpd = data.frame(knockdown,collagen)
rownames(newpd) = rownames(pd)
head(newpd)
@

\section{Analysis}

Use the affy package to load and normalize the data.

<<normalizeCelFiles>>=
library(affy)
eset = justRMA(filenames=paste0(rownames(pd),'.CEL.gz'),compress=TRUE)
pData(eset)=newpd
sampleNames(eset) = rownames(pd)
@

The 'eset' object is an ExpressionSet and contains the normalized gene expression measures as well as the sample
annotation.  They cannot get out-of-sync, now.

Use the limma package for analysis.  I am keeping it simple here, but adding multiple factors is just a little
more complex.  In any case, the next 10 lines do the entire analysis and produce a heatmap.

<<>>=
library(limma)
dm = model.matrix( ~ 0 + knockdown,pData(eset))
cm = makeContrasts(knockdownTRUE-knockdownFALSE,levels=dm)
fit = lmFit(eset,dm)
fit2 = contrasts.fit(fit,cm)
fit3 = eBayes(fit2)
topTable(fit3)
tt = topTable(fit3,number=30)
@

Now, make a heatmap of the top genes.  Dumping out data for excel, making online reports, or pdfs like
this one are straightforward.  I often make large PDFs of all DE genes for collaborators to look at on
their own.  You can also export in a format that TreeView understands so that collaborators can explore
on their own.

<<fig=TRUE>>=
library(gplots)
heatmap.2(exprs(eset[tt$ID,]),trace='none',col=greenred(59),
          scale='row',labCol=as.character(pData(eset)$knockdown))
@


\section{sessionInfo}

<<>>=
sessionInfo()
@


\end{document}