### R code from vignette source 'vignettes/DEXSeq/inst/doc/DEXSeq.Rnw' ################################################### ### code chunk number 1: style ################################################### BiocStyle::latex() ################################################### ### code chunk number 2: options ################################################### options(digits=2, width=80) ################################################### ### code chunk number 3: systemFile ################################################### system.file( "python_scripts", package="DEXSeq" ) ################################################### ### code chunk number 4: systemFileCheck ################################################### system.file( "python_scripts", package="DEXSeq", mustWork=TRUE ) ################################################### ### code chunk number 5: loadDEXSeq ################################################### library( "DEXSeq" ) ################################################### ### code chunk number 6: sampleTable ################################################### sampleTable <- data.frame( row.names = c( "untreated1", "untreated2", "untreated3", "untreated4", "treated1", "treated2", "treated3" ), countFile = c( "untreated1.counts", "untreated2.counts", "untreated3.counts", "untreated4.counts", "treated1.counts", "treated2.counts","treated3.counts" ), condition = c( "control", "control", "control", "control", "knockdown", "knockdown", "knockdown" ), libType = c( "single-end", "paired-end", "paired-end", "single-end", "single-end", "paired-end", "paired-end" ) ) ################################################### ### code chunk number 7: displaySampleTable ################################################### sampleTable ################################################### ### code chunk number 8: makeecs (eval = FALSE) ################################################### ## ecs <- read.HTSeqCounts( ## sampleTable$countFile, ## sampleTable, ## "Dmel_flattenend.gff" ) ################################################### ### code chunk number 9: start ################################################### library( DEXSeq ) data( "pasillaExons", package="pasilla" ) ################################################### ### code chunk number 10: renameEcs ################################################### ecs <- pasillaExons ecs ################################################### ### code chunk number 11: seeCounts ################################################### head( counts(ecs), 7 ) ################################################### ### code chunk number 12: fData ################################################### head( fData(pasillaExons), 3 ) ################################################### ### code chunk number 13: pData ################################################### design(ecs) ################################################### ### code chunk number 14: sizeFactors1 ################################################### ecs <- estimateSizeFactors( ecs ) ################################################### ### code chunk number 15: inspectSizeFactors ################################################### sizeFactors(ecs) ################################################### ### code chunk number 16: estDisp1 ################################################### ecs <- estimateDispersions( ecs ) ################################################### ### code chunk number 17: fitDisp1 ################################################### ecs <- fitDispersionFunction( ecs ) ################################################### ### code chunk number 18: fitDisp2 ################################################### head( fData(ecs)$dispBeforeSharing ) ecs@dispFitCoefs head( fData(ecs)$dispFitted ) ################################################### ### code chunk number 19: fitdiagnostics ################################################### plotDispEsts( ecs ) ################################################### ### code chunk number 20: testForDEU1 ################################################### ecs <- testForDEU( ecs ) ################################################### ### code chunk number 21: testForDEU2 ################################################### head( fData(ecs)[, c( "pvalue", "padjust" ) ] ) ################################################### ### code chunk number 22: estFC ################################################### ecs <- estimatelog2FoldChanges( ecs ) ################################################### ### code chunk number 23: DEUresults ################################################### res1 <- DEUresultTable(ecs) head( res1 ) ################################################### ### code chunk number 24: tallyExons ################################################### table ( res1$padjust < 0.1 ) ################################################### ### code chunk number 25: tallyGenes ################################################### table ( tapply( res1$padjust < 0.1, geneIDs(pasillaExons), any ) ) ################################################### ### code chunk number 26: MvsA ################################################### plotMA( ecs, FDR=0.1, ylim=c(-4,4), cex=0.8 ) ################################################### ### code chunk number 27: design ################################################### design(pasillaExons) ################################################### ### code chunk number 28: formulas2 ################################################### formulaFullModel <- ~ sample + exon + type:exon + condition:exon formulaReducedModel <- ~ sample + exon + type:exon ################################################### ### code chunk number 29: new_ecs ################################################### ecs2 <- pasillaExons ################################################### ### code chunk number 30: sf_again ################################################### ecs2 <- estimateSizeFactors( ecs2 ) ################################################### ### code chunk number 31: estDisps_again ################################################### ecs2 <- estimateDispersions( ecs2, formula = formulaFullModel ) ################################################### ### code chunk number 32: fit_again ################################################### ecs2 <- fitDispersionFunction( ecs2 ) ################################################### ### code chunk number 33: test_again ################################################### ecs2 <- testForDEU( ecs2, formula0 = formulaReducedModel, formula1 = formulaFullModel ) ################################################### ### code chunk number 34: res_again ################################################### res2 <- DEUresultTable( ecs2 ) head(res2) ################################################### ### code chunk number 35: table2 ################################################### table( res2$padjust < 0.1 ) ################################################### ### code chunk number 36: table3 ################################################### table( before = res1$padjust < 0.1, now = res2$padjust < 0.1 ) ################################################### ### code chunk number 37: oneMore ################################################### stopifnot( sum( res2$padjust<.1, na.rm=TRUE ) == sum( res1$padjust<.1, na.rm=TRUE ) + 1 ) ################################################### ### code chunk number 38: plot1 ################################################### plotDEXSeq( ecs2, "FBgn0010909", legend=TRUE, cex.axis=1.2, cex=1.3, lwd=2 ) ################################################### ### code chunk number 39: checkClaim ################################################### wh = (fData(ecs2)$geneID=="FBgn0010909") stopifnot(sum(fData(ecs2)$padjust[wh] < formals(plotDEXSeq)$FDR)==1) ################################################### ### code chunk number 40: plot2 ################################################### plotDEXSeq( ecs2, "FBgn0010909", displayTranscripts=TRUE, legend=TRUE, cex.axis=1.2, cex=1.3, lwd=2 ) ################################################### ### code chunk number 41: plot3 ################################################### plotDEXSeq( ecs2, "FBgn0010909", expression=FALSE, norCounts=TRUE, legend=TRUE, cex.axis=1.2, cex=1.3, lwd=2 ) ################################################### ### code chunk number 42: plot4 ################################################### plotDEXSeq( ecs2, "FBgn0010909", expression=FALSE, splicing=TRUE, legend=TRUE, cex.axis=1.2, cex=1.3, lwd=2 ) ################################################### ### code chunk number 43: DEXSeqHTML (eval = FALSE) ################################################### ## DEXSeqHTML( pasillaExons, FDR=0.1, color=c("#FF000080", "#0000FF80") ) ################################################### ### code chunk number 44: para1 (eval = FALSE) ################################################### ## data("pasillaExons", package="pasilla") ## library("parallel") ## pasillaExons <- estimateSizeFactors( pasillaExons ) ## pasillaExons <- estimateDispersions( pasillaExons, nCores=16) ## pasillaExons <- fitDispersionFunction( pasillaExons ) ## pasillaExons <- testForDEU( pasillaExons, nCores=16) ################################################### ### code chunk number 45: alldeu ################################################### data( "pasillaExons", package="pasilla" ) pasillaExons <- doCompleteDEUAnalysis( pasillaExons, formula0 = ~ sample + exon + type:exon, formula1 = ~ sample + exon + type:exon + condition:exon, nCores = 1 ) ################################################### ### code chunk number 46: alldeuout ################################################### head( DEUresultTable( pasillaExons ) ) ################################################### ### code chunk number 47: buildExonCountSetLoadPacks (eval = FALSE) ################################################### ## library("GenomicFeatures") ## library("GenomicRanges") ## library("Rsamtools") ################################################### ### code chunk number 48: buildExonCountSetDownloadAnno (eval = FALSE) ################################################### ## hse <- makeTranscriptDbFromBiomart( biomart="ensembl", ## dataset="hsapiens_gene_ensembl" ) ################################################### ### code chunk number 49: buildExonCountSetDisjoin (eval = FALSE) ################################################### ## exonicParts <- disjointExons( hse, aggregateGenes=TRUE ) ################################################### ### code chunk number 50: buildExonCountSet2FindBAMs (eval = FALSE) ################################################### ## bamDir <- system.file( "extdata", package="parathyroidSE", mustWork=TRUE ) ## fls <- list.files( bamDir, pattern="bam$", full=TRUE ) ################################################### ### code chunk number 51: buildExonCountSet2FindBAMs (eval = FALSE) ################################################### ## bamlst <- BamFileList( fls, index=character(), yieldSize=100000, obeyQname=TRUE ) ## SE <- summarizeOverlaps( exonicParts, bamlst, mode="Union", singleEnd=FALSE, ## ignore.strand=TRUE, inter.feature=FALSE, fragments=TRUE ) ################################################### ### code chunk number 52: buildExonCountSet3 (eval = FALSE) ################################################### ## ecs <- buildExonCountSet( SE, c("A", "A", "B"), exonicParts ) ################################################### ### code chunk number 53: countTableForGene ################################################### countTableForGene( pasillaExons, "FBgn0010909" ) ################################################### ### code chunk number 54: countTableForGeneNorm ################################################### head( countTableForGene( pasillaExons, "FBgn0010909", normalized=TRUE ) ) ################################################### ### code chunk number 55: gct ################################################### head( geneCountTable( pasillaExons ) ) ################################################### ### code chunk number 56: modelFrame ################################################### constructModelFrame( pasillaExons ) ################################################### ### code chunk number 57: modelFrame ################################################### mf <- modelFrameForGene( pasillaExons, "FBgn0010909" ) head( mf ) ################################################### ### code chunk number 58: acc ################################################### head( geneIDs(pasillaExons) ) head( exonIDs(pasillaExons) ) ################################################### ### code chunk number 59: sessionInfo ################################################### sessionInfo()