Commits

Sean Davis committed a53a2d4

Example RNW and pdf file

Comments (0)

Files changed (2)

+\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}

limmaExample.pdf

Binary file added.