### R code from vignette source 'vignettes/EDASeq/inst/doc/EDASeq.Rnw' ################################################### ### code chunk number 1: EDASeq.Rnw:38-39 ################################################### options(width=60) ################################################### ### code chunk number 2: data ################################################### require(EDASeq) require(yeastRNASeq) require(leeBamViews) ################################################### ### code chunk number 3: import-raw ################################################### files <- list.files(file.path(system.file(package = "yeastRNASeq"), "reads"), pattern = "fastq", full.names = TRUE) names(files) <- gsub("\\.fastq.*", "", basename(files)) met <- DataFrame(conditions=c(rep("mut",2),rep("wt",2)), row.names=names(files)) fastq <- FastqFileList(files) elementMetadata(fastq) <- met fastq ################################################### ### code chunk number 4: import-aligned ################################################### files <- list.files(file.path(system.file(package = "leeBamViews"), "bam"), pattern = "bam$", full.names = TRUE) names(files) <- gsub("\\.bam", "", basename(files)) gt <- gsub(".*/", "", files) gt <- gsub("_.*", "", gt) lane <- gsub(".*(.)$", "\\1", gt) geno <- gsub(".$", "", gt) pd <- DataFrame(geno=geno, lane=lane, row.names=paste(geno,lane,sep=".")) bfs <- BamFileList(files) elementMetadata(bfs) <- pd bfs ################################################### ### code chunk number 5: plot-total ################################################### colors <- c(rep(rgb(1,0,0,alpha=0.7),2),rep(rgb(0,0,1,alpha=0.7),2), rep(rgb(0,1,0,alpha=0.7),2),rep(rgb(0,1,1,alpha=0.7),2)) barplot(bfs,las=2,col=colors) ################################################### ### code chunk number 6: plot-mean-qual ################################################### plotQuality(bfs,col=colors,lty=1) legend("topright",unique(elementMetadata(bfs)[,1]), fill=unique(colors)) ################################################### ### code chunk number 7: plot-qual-1 ################################################### plotQuality(bfs[[1]],cex.axis=.8) ################################################### ### code chunk number 8: plot-mapped-1 ################################################### barplot(bfs[[1]],las=2) ################################################### ### code chunk number 9: plot-nt-1 ################################################### plotNtFrequency(bfs[[1]]) ################################################### ### code chunk number 10: load-gene-level ################################################### data(geneLevelData) head(geneLevelData) ################################################### ### code chunk number 11: load-lgc ################################################### data(yeastGC) head(yeastGC) data(yeastLength) head(yeastLength) ################################################### ### code chunk number 12: filter ################################################### filter <- apply(geneLevelData,1,function(x) mean(x)>10) table(filter) common <- intersect(names(yeastGC),rownames(geneLevelData[filter,])) length(common) ################################################### ### code chunk number 13: create-object ################################################### feature <- data.frame(gc=yeastGC,length=yeastLength) data <- newSeqExpressionSet(exprs=as.matrix(geneLevelData[common,]), featureData=feature[common,], phenoData=data.frame( conditions=c(rep("mut",2),rep("wt",2)), row.names=colnames(geneLevelData))) data ################################################### ### code chunk number 14: show-data ################################################### head(exprs(data)) pData(data) head(fData(data)) ################################################### ### code chunk number 15: show-offset ################################################### head(offst(data)) ################################################### ### code chunk number 16: boxplot-genelevel ################################################### boxplot(data,col=colors[1:4]) ################################################### ### code chunk number 17: md-plot ################################################### MDPlot(data,c(1,3)) ################################################### ### code chunk number 18: plot-mean-var-mut ################################################### meanVarPlot(data[,1:2],log=T, ylim=c(0,16)) ################################################### ### code chunk number 19: plot-mean-var-all ################################################### meanVarPlot(data,log=T, ylim=c(0,16)) ################################################### ### code chunk number 20: plot-gc ################################################### biasPlot(data,"gc",log=T,ylim=c(1,5)) ################################################### ### code chunk number 21: boxplot-gc ################################################### lfc <- log(exprs(data)[,3]+0.5)-log(exprs(data)[,1]+0.5) biasBoxplot(lfc,fData(data)$gc) ################################################### ### code chunk number 22: normalization ################################################### dataWithin <- withinLaneNormalization(data,"gc",which="full") dataNorm <- betweenLaneNormalization(dataWithin,which="full") ################################################### ### code chunk number 23: plot-gc-norm ################################################### biasPlot(dataNorm,"gc",log=T,ylim=c(1,5)) ################################################### ### code chunk number 24: boxplot-norm ################################################### boxplot(dataNorm,col=colors) ################################################### ### code chunk number 25: norm-offset ################################################### dataOffset <- withinLaneNormalization(data,"gc", which="full",offset=TRUE) dataOffset <- betweenLaneNormalization(dataOffset, which="full",offset=TRUE) ################################################### ### code chunk number 26: edger ################################################### library(edgeR) design <- model.matrix(~conditions, data=pData(dataOffset)) disp <- estimateGLMCommonDisp(exprs(dataOffset), design, offset=-offst(dataOffset)) fit <- glmFit(exprs(dataOffset), design, disp, offset=-offst(dataOffset)) lrt <- glmLRT(fit, coef=2) topTags(lrt) ################################################### ### code chunk number 27: deseq ################################################### library(DESeq) counts <- as(dataNorm,"CountDataSet") sizeFactors(counts) <- rep(1,4) counts <- estimateDispersions(counts) res <- nbinomTest(counts, "wt", "mut" ) head(res) ################################################### ### code chunk number 28: unrounded ################################################### dataNorm <- betweenLaneNormalization(data, round=FALSE) dataOffset <- betweenLaneNormalization(data, offset=TRUE) norm1 <- exprs(dataNorm) norm2 <- exp(log(exprs(dataOffset) + 0.1 ) + offst(dataOffset)) - 0.1 head(norm1 - norm2) ################################################### ### code chunk number 29: rounded ################################################### head(round(exprs(dataNorm)) - round(exprs(dataOffset) * exp(offst(dataOffset)))) ################################################### ### code chunk number 30: sessionInfo ################################################### toLatex(sessionInfo())