# Limma Example / limmaExample.Rnw

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 \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. <>= library(GEOquery) getGEOSuppFiles('GSE16200') untar('GSE16200/GSE16200_RAW.tar') @ <>= 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. <>= 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. <>= 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}