### R code from vignette source SeattleIntro2012 'RNASeq.Rnw' ################################################### ### code chunk number 1: setup ################################################### options(digits=2) library(SeattleIntro2012) library(edgeR) library(goseq) ################################################### ### code chunk number 2: counts ################################################### ## Read in the Pasilla data set data(counts) dim(counts) grps <- factor(sub("[1-4].*", "", colnames(counts)), levels=c("untreated", "treated")) pairs <- factor(c("single", "paired", "paired", "single", "single", "paired", "paired")) pData <- data.frame(Group=grps, PairType=pairs, row.names=colnames(counts)) ################################################### ### code chunk number 3: DGEList ################################################### library(edgeR) dge <- DGEList(counts, group=pData$Group) dge <- calcNormFactors(dge) ################################################### ### code chunk number 4: DEGList-filter ################################################### m <- sweep(dge$counts, 2, 1e6 / dge$samples$lib.size, `*`) ridx <- rowSums(m > 1) >= 2 table(ridx) # number filtered / retained dge <- dge[ridx,] ################################################### ### code chunk number 5: mds ################################################### plotMDS.DGEList(dge) ################################################### ### code chunk number 6: design ################################################### (design <- model.matrix(~ Group, pData)) ################################################### ### code chunk number 7: common.dispersion ################################################### dge <- estimateTagwiseDisp(dge) mean(sqrt(dge$tagwise.dispersion)) ################################################### ### code chunk number 8: glmFit ################################################### fit <- glmFit(dge, design) ################################################### ### code chunk number 9: lrt ################################################### lrTest <- glmLRT(dge, fit, coef=2) ################################################### ### code chunk number 10: topTags ################################################### tt <- topTags(lrTest, n=10) tt[1:3,] ################################################### ### code chunk number 11: sanity ################################################### sapply(rownames(tt$table)[1:4], function(x) tapply(counts[x,], pData$Group, mean))