### R code from vignette source 'Glimma.Rnw' ################################################### ### code chunk number 1: style-Sweave (eval = FALSE) ################################################### ## BiocStyle::latex() ################################################### ### code chunk number 2: Rsetup ################################################### set.seed(20161000) options(digits=2) options(width=100) options(browser="false") ################################################### ### code chunk number 3: Rgetstarted ################################################### library(Glimma) library(limma) library(edgeR) data(lymphomaRNAseq) rnaseq <- lymphomaRNAseq rnaseq$samples$group ################################################### ### code chunk number 4: Rpreprocessing ################################################### rnaseq <- rnaseq[rowSums(cpm(rnaseq)>1)>=3,] rnaseq <- calcNormFactors(rnaseq) ################################################### ### code chunk number 5: Rmds ################################################### groups <- rnaseq$samples$group glMDSPlot(rnaseq, groups=groups) ################################################### ### code chunk number 6: Rvoomqw ################################################### design <- model.matrix(~0+groups) contrasts <- cbind(Smchd1null.vs.WT=c(-1,1)) vm <- voomWithQualityWeights(rnaseq, design=design) fit <- lmFit(vm, design=design) fit <- contrasts.fit(fit, contrasts) fit <- eBayes(fit) dt <- decideTests(fit) summary(dt) ################################################### ### code chunk number 7: Rmd ################################################### glMDPlot(fit, status=dt, counts=vm, groups=groups, side.main="Symbols") ################################################### ### code chunk number 8: Rxy ################################################### glXYPlot(x=fit$coef, y=fit$lod, xlab="logFC", ylab="logodds", status=dt, counts=vm, groups=groups, anno=fit$genes) ################################################### ### code chunk number 9: Rmdsmatrix ################################################### lcpm <- cpm(rnaseq, log=TRUE, normalized.lib.sizes=TRUE) glMDSPlot(lcpm, groups=groups) ################################################### ### code chunk number 10: Rmdscheckinput (eval = FALSE) ################################################### ## # Check that MDS plot runs on EList object: PASS ## glMDSPlot(vm, groups=groups) ## # Check that MDS plot runs on DESeqDataSet object: FAIL regarding groups ## library(DESeq2) ## rnaseq.deseq2 <- DESeqDataSetFromMatrix(rnaseq$counts, colData=rnaseq$samples, design=~group) ## glMDSPlot(rnaseq.deseq2, groups=groups) ################################################### ### code chunk number 11: Rmdsgroups ################################################### groups.df <- as.data.frame(cbind( genotype=as.character(groups), lane=c(rep(4,6),3))) groups.df glMDSPlot(lcpm, groups=groups.df) ################################################### ### code chunk number 12: Rmdsimple ################################################### glMDPlot(fit) ################################################### ### code chunk number 13: Rmdage ################################################### groups.age <- runif(ncol(rnaseq), min=3, max=12) groups.age ################################################### ### code chunk number 14: Rmdsamples ################################################### cols <- c("yellow","blue","magenta") sample.cols <- c("limegreen", "purple")[groups] glMDPlot(fit, status=dt, counts=vm, groups=groups.age, sample.cols=sample.cols, cols=cols, side.ylab="logCPM", side.xlab="Age (in months)", side.main="Symbols", main=colnames(dt)) ################################################### ### code chunk number 15: Rmdcheckinput (eval = FALSE) ################################################### ## # limma - PASS ## glMDPlot(fit) ## # edgeR - PASS ## rnaseq.edger <- estimateDisp(rnaseq, design=design) ## fit.edger <- exactTest(rnaseq.edger) ## glMDPlot(fit.edger) ## # DESeq2 - FAIL regarding annotation ## library(DESeq2) ## rnaseq.deseq2 <- DESeqDataSetFromMatrix(rnaseq$counts, colData=rnaseq$samples, design=~group) ## mcols(rnaseq.deseq2) <- DataFrame(mcols(rnaseq.deseq2), rnaseq$genes) ## fit.deseq2 <- DESeq(rnaseq.deseq2) ## glMDPlot(fit.deseq2) ################################################### ### code chunk number 16: Rmdanno ################################################### ID <- paste(fit$genes$Symbols, fit$genes$GeneID) DE <- c("downregulated", "notDE", "upregulated")[as.factor(dt)] anno <- as.data.frame(cbind(ID, DE)) head(anno) glMDPlot(fit, counts=vm, groups=groups, side.main="ID", anno=anno, display.columns=c("ID", "GeneName", "DE")) ################################################### ### code chunk number 17: Rxyvolcanosimple ################################################### glXYPlot(x=fit$coef, y=fit$lod) ################################################### ### code chunk number 18: Rxyvolcanomore ################################################### glXYPlot(x=fit$coef, y=fit$lod, xlab="logFC", ylab="logodds", status=dt, anno=anno, side.main="ID", counts=vm, groups=groups, sample.cols=sample.cols) ################################################### ### code chunk number 19: Rarraydata ################################################### data(arraydata) arrays <- arraydata$arrays targets <- arraydata$targets dim(arrays) targets ################################################### ### code chunk number 20: Rmds.arrays ################################################### glMDSPlot(arrays, groups=targets[,c("Condition", "Chip", "Experiment")]) ################################################### ### code chunk number 21: Rarraysfit ################################################### design <- model.matrix(~0+targets$Condition+as.factor(targets$Experiment)) contrasts <- cbind( DP_Ezh2KO.vs.WT=c(1,-1,0,0,0,0), Lum_Ezh2KO.vs.WT=c(0,0,1,-1,0,0)) fit <- lmFit(arrays, design) fit <- contrasts.fit(fit, contrasts) fit <- eBayes(fit) dt <- decideTests(fit, p.value=0.1) summary(dt) ################################################### ### code chunk number 22: Rarraysmd ################################################### sample.cols <- c("purple", "magenta", "green")[targets$Experiment] for (COEF in 1:2) { glMDPlot(fit, status=dt, coef=COEF, main=colnames(fit)[COEF], counts=arrays, groups=targets$Condition, sample.cols=sample.cols, side.ylab="Log-expression", side.main="ProbeID") } ################################################### ### code chunk number 23: Rarraysxy ################################################### dt2 <- rep(0, nrow(dt)) dt2[rowSums(dt!=0)==1] <- -1 dt2[rowSums(dt!=0)==2] <- 1 table(dt2) cols <- c("black", "grey", "red") glXYPlot(fit$coef[,1], y=fit$coef[,2], xlab="DP", ylab="Lum", status=dt2, cols=cols, anno=fit$genes, side.main="ProbeID", counts=arrays, groups=targets$Condition, sample.cols=sample.cols, side.ylab="Log-expression", main="logFCs") ################################################### ### code chunk number 24: Rmdedgeret (eval = FALSE) ################################################### ## groups <- rnaseq$samples$group ## design <- model.matrix(~groups) ## colnames(design) <- c("WT", "Smchd1null.vs.WT") ## rnaseq.edger <- estimateDisp(rnaseq, design=design) ## fit.edger <- exactTest(rnaseq.edger) ## dt.edger <- decideTestsDGE(fit.edger) ## glMDPlot(fit.edger, status=dt.edger, counts=rnaseq, groups=groups, transform=TRUE) ################################################### ### code chunk number 25: Rmdedgerlrt (eval = FALSE) ################################################### ## fit.edger <- glmFit(rnaseq.edger, design) ## fit.edger <- glmLRT(fit.edger) ## dt.edger <- decideTestsDGE(fit.edger) ## glMDPlot(fit.edger, status=dt.edger, counts=rnaseq, groups=groups, transform=TRUE) ################################################### ### code chunk number 26: Rmddeseq2 (eval = FALSE) ################################################### ## # BUG regarding the scale of sample expression ## library(DESeq2) ## rnaseq.deseq2 <- DESeqDataSetFromMatrix( ## rnaseq$counts, colData=rnaseq$samples, design=~group) ## mcols(rnaseq.deseq2) <- DataFrame(mcols(rnaseq.deseq2), rnaseq$genes) ## rnaseq.deseq2 <- DESeq(rnaseq.deseq2) ## fit.deseq2 <- results(rnaseq.deseq2, contrast=c("group", "Smchd1-null", "WT")) ## dt.deseq2 <- as.numeric(fit.deseq2$padj<0.05) ## glMDPlot(fit.deseq2, status=dt.deseq2, counts=rnaseq, groups=groups, transform=FALSE, ## samples=colnames(rnaseq), anno=rnaseq$genes) ################################################### ### code chunk number 27: Rsessinfo ################################################### sessionInfo()