Friday, July 18, 2014

Expression analysis with R: package edgeR

Sample code for edgeR and DESeq for differential expressed gene

Best step-by-step code

Comparision of edgeR and DESeq

# edgeR 
## Make design matrix
condition < relevel(factor(meta$condition), ref="untreated")
libType < factor(meta$libType)
edesign < model.matrix(~libType+condition)
## Make new DGEList, normalize by library size, and estimate dispersion allowing
 possible trend with average count size
e < DGEList(counts=counttable)
e < calcNormFactors(e)
e < estimateGLMCommonDisp(e, edesign)
e < estimateGLMTrendedDisp(e, edesign)
e < estimateGLMTagwiseDisp(e, edesign)
## MDS Plot
plotMDS(e, main="edgeR MDS Plot")
## Biological coefficient of variation plot
plotBCV(e, cex=0.4, main="edgeR: Biological coefficient of variation (BCV) vs abundance")
## Fit the model, testing the coefficient for the treated vs untreated comparison
efit < glmFit(e, edesign)
efit < glmLRT(efit, coef="conditiontreated")
## Make a table of results
etable < topTags(efit, n=nrow(e))$table
etable < etable[order(etable$FDR), ]
## ~MA Plot
with(etable, plot(logCPM, logFC, pch=20, main="edgeR: Fold change vs abundance"))
with(subset(etable, FDR<0.05), points(logCPM, logFC, pch=20, col="red"))
abline(h=c(-1,1), col="blue")