### R code from vignette source 'vignettes/DEXSeq/inst/doc/DEXSeq.Rnw' ################################################### ### code chunk number 1: options ################################################### options(digits=3, width=106) ################################################### ### code chunk number 2: start ################################################### library("DEXSeq") library("pasilla") data("pasillaExons", package="pasilla") ################################################### ### code chunk number 3: pData ################################################### design(pasillaExons) ################################################### ### code chunk number 4: counts ################################################### head( counts(pasillaExons), 20 ) ################################################### ### code chunk number 5: options ################################################### options(width=96) ################################################### ### code chunk number 6: fData ################################################### head(fData(pasillaExons)[,c(1,2,9:12)]) ################################################### ### code chunk number 7: tabtab1 ################################################### tabGeneIDsPasillaExons = table(geneIDs(pasillaExons)) ngenesPasillaExons = sum(tabGeneIDsPasillaExons>0) tabtabGeneIDsPasillaExons = table(tabGeneIDsPasillaExons) stopifnot(tabtabGeneIDsPasillaExons["36"]==1, tabtabGeneIDsPasillaExons["16"]==3) ################################################### ### code chunk number 8: tabtab2 ################################################### table(table(geneIDs(pasillaExons))) ################################################### ### code chunk number 9: sizeFactors1 ################################################### pasillaExons <- estimateSizeFactors(pasillaExons) sizeFactors(pasillaExons) ################################################### ### code chunk number 10: helperfunctions ################################################### plotDispEsts = function( cds, ymin, linecol="#ff000080", xlab = "mean of normalized counts", ylab = "dispersion", log = "xy", cex = 0.45, ... ) { px = rowMeans( counts( cds, normalized=TRUE ) ) sel = (px>0) px = px[sel] py = fData(cds)$dispBeforeSharing[sel] if(missing(ymin)) ymin = 10^floor(log10(min(py[py>0], na.rm=TRUE))-0.1) plot(px, pmax(py, ymin), xlab=xlab, ylab=ylab, log=log, pch=ifelse(py=0.1, "gray32", "red3"), linecol = "#ff000080", xlab = "mean of normalized counts", ylab = expression(log[2]~fold~change), log = "x", cex=0.45, ...) { if (!(is.data.frame(x) && all(c("baseMean", "log2FoldChange") %in% colnames(x)))) stop("'x' must be a data frame with columns named 'baseMean', 'log2FoldChange'.") x = subset(x, baseMean!=0) py = x$log2FoldChange if(missing(ylim)) ylim = c(-1,1) * quantile(abs(py[is.finite(py)]), probs=0.99) * 1.1 plot(x$baseMean, pmax(ylim[1], pmin(ylim[2], py)), log=log, pch=ifelse(pyylim[2], 2, 16)), cex=cex, col=col, xlab=xlab, ylab=ylab, ylim=ylim, ...) abline(h=0, lwd=4, col=linecol) } ################################################### ### code chunk number 11: estDisp1 ################################################### pasillaExons <- estimateDispersions(pasillaExons) ################################################### ### code chunk number 12: fitDisp1 ################################################### pasillaExons <- fitDispersionFunction(pasillaExons) ################################################### ### code chunk number 13: fitDisp2 ################################################### head(fData(pasillaExons)$dispBeforeSharing) pasillaExons@dispFitCoefs head(fData(pasillaExons)$dispFitted) ################################################### ### code chunk number 14: fitdiagnostics ################################################### plotDispEsts( pasillaExons ) ################################################### ### code chunk number 15: testForDEU1 ################################################### pasillaExons <- testForDEU( pasillaExons ) ################################################### ### code chunk number 16: testForDEU2 ################################################### head( fData(pasillaExons)[, c( "pvalue", "padjust" ) ] ) ################################################### ### code chunk number 17: estFC ################################################### pasillaExons <- estimatelog2FoldChanges( pasillaExons ) ################################################### ### code chunk number 18: DEUresults ################################################### res1 <- DEUresultTable(pasillaExons) head( res1 ) ################################################### ### code chunk number 19: tallyExons ################################################### table ( res1$padjust < 0.1 ) ################################################### ### code chunk number 20: tallyGenes ################################################### table ( tapply( res1$padjust < 0.1, geneIDs(pasillaExons), any ) ) ################################################### ### code chunk number 21: MvsA ################################################### plotMA(with(res1, data.frame(baseMean = meanBase, log2FoldChange = `log2fold(untreated/treated)`, padj = padjust)), ylim=c(-4,4), cex=0.8) ################################################### ### code chunk number 22: design ################################################### design(pasillaExons) ################################################### ### code chunk number 23: formuladispersion ################################################### formuladispersion <- count ~ sample + ( condition + type ) * exon pasillaExons <- estimateDispersions( pasillaExons, formula = formuladispersion ) pasillaExons <- fitDispersionFunction(pasillaExons) ################################################### ### code chunk number 24: formula1 ################################################### formula0 <- count ~ sample + type * exon + condition formula1 <- count ~ sample + type * exon + condition * I(exon == exonID) ################################################### ### code chunk number 25: formula2 ################################################### pasillaExons <- testForDEU( pasillaExons, formula0=formula0, formula1=formula1 ) res2 <- DEUresultTable( pasillaExons ) ################################################### ### code chunk number 26: formula3 ################################################### table( res2$padjust < 0.1 ) table(res1$padjust < 0.1, res2$padjust < 0.1) ################################################### ### code chunk number 27: figcomparep ################################################### bottom = function(x) pmax(x, 1e-6) plot(bottom(res1$padjust), bottom(res2$padjust), log="xy", pch=20) abline(a=0,b=1, col="red3") abline(v=0.1, h=0.1, col="blue") ################################################### ### code chunk number 28: plot1 ################################################### plotDEXSeq(pasillaExons, "FBgn0010909", cex.axis=1.2, cex=1.3, lwd=2, legend=TRUE) ################################################### ### code chunk number 29: checkClaim ################################################### wh = (fData(pasillaExons)$geneID=="FBgn0010909") stopifnot(sum(fData(pasillaExons)$padjust[wh] < formals(plotDEXSeq)$FDR)==1) ################################################### ### code chunk number 30: plot2 ################################################### plotDEXSeq(pasillaExons, "FBgn0010909", displayTranscripts=TRUE, cex.axis=1.2, cex=1.3, lwd=2, legend=TRUE) ################################################### ### code chunk number 31: plot3 ################################################### plotDEXSeq(pasillaExons, "FBgn0010909", expression=FALSE, norCounts=TRUE, cex.axis=1.2, cex=1.3, lwd=2, legend=TRUE) ################################################### ### code chunk number 32: plot4 ################################################### plotDEXSeq(pasillaExons, "FBgn0010909", expression=FALSE, splicing=TRUE, cex.axis=1.2, cex=1.3, lwd=2, legend=TRUE ) ################################################### ### code chunk number 33: DEXSeqHTML (eval = FALSE) ################################################### ## DEXSeqHTML( pasillaExons, FDR=0.1, color=c("#FF000080", "#0000FF80") ) ################################################### ### code chunk number 34: para1 (eval = FALSE) ################################################### ## data("pasillaExons", package="pasilla") ## library(parallel) ## pasillaExons <- estimateSizeFactors( pasillaExons ) ## pasillaExons <- estimateDispersions( pasillaExons, nCores=3, quiet=TRUE) ## pasillaExons <- fitDispersionFunction( pasillaExons ) ## pasillaExons <- testForDEU( pasillaExons, nCores=3) ################################################### ### code chunk number 35: alldeu ################################################### data("pasillaExons", package="pasilla") pasillaExons <- makeCompleteDEUAnalysis( pasillaExons, formulaDispersion=formuladispersion, formula0=formula0, formula1=formula1, nCores=1) ################################################### ### code chunk number 36: pkgRoot ################################################### pkgDir = system.file(package="DEXSeq") pkgDir list.files(pkgDir) list.files(file.path(pkgDir, "python_scripts")) ################################################### ### code chunk number 37: dirpasilla ################################################### dir(system.file("extdata",package="pasilla")) ################################################### ### code chunk number 38: ecswithout ################################################### bare <- newExonCountSet( countData = counts(pasillaExons), design=design(pasillaExons), geneIDs=geneIDs(pasillaExons), exonIDs=exonIDs(pasillaExons)) ################################################### ### code chunk number 39: buildExonCountSet (eval = FALSE) ################################################### ## ## library(GenomicFeatures) ## hse <- makeTranscriptDbFromBiomart(biomart="ensembl", dataset="hsapiens_gene_ensembl") ## exonicParts <- prepareAnnotationForDEXSeq( hse, aggregateGenes=TRUE ) ################################################### ### code chunk number 40: buildExonCountSet2 (eval = FALSE) ################################################### ## ## bamDir <- system.file("extdata",package="parathyroidSE",mustWork=TRUE) ## fls <- list.files(bamDir, pattern="bam$",full=TRUE) ## bamlst <- BamFileList(fls) ## SE <- countReadsForDEXSeq( exonicParts, bamlst ) ################################################### ### code chunk number 41: buildExonCountSet3 (eval = FALSE) ################################################### ## ## ecs <- buildExonCountSet( SE, c("A", "A", "B"), exonicParts ) ################################################### ### code chunk number 42: countTableForGene ################################################### head(countTableForGene(pasillaExons,"FBgn0010909")) ################################################### ### code chunk number 43: countTableForGeneNorm ################################################### head(countTableForGene(pasillaExons,"FBgn0010909", normalized=TRUE)) ################################################### ### code chunk number 44: modelFrame ################################################### mf <- modelFrameForGene( pasillaExons, "FBgn0010909" ) head( mf ) ################################################### ### code chunk number 45: mffg1 ################################################### mf <- estimateExonDispersionsForModelFrame( modelFrameForGene( pasillaExons, "FBgn0010909" ) ) ################################################### ### code chunk number 46: mffg2 ################################################### testGeneForDEU( pasillaExons, "FBgn0010909" ) ################################################### ### code chunk number 47: gct ################################################### head(geneCountTable(pasillaExons)) ################################################### ### code chunk number 48: acc ################################################### head(geneIDs(pasillaExons)) head(exonIDs(pasillaExons)) ################################################### ### code chunk number 49: sessionInfo ################################################### sessionInfo()