## ----echo=FALSE, include=FALSE--------------------------------------------- knitr::opts_chunk$set(tidy = FALSE, cache = TRUE, dev = "png", message = FALSE, error = FALSE, warning = TRUE) ## ----systemFile------------------------------------------------------------ pythonScriptsDir = system.file( "python_scripts", package="DEXSeq" ) list.files(pythonScriptsDir) system.file( "python_scripts", package="DEXSeq", mustWork=TRUE ) ## ----loadDEXSeq------------------------------------------------------------ inDir = system.file("extdata", package="pasilla") countFiles = list.files(inDir, pattern="fb.txt$", full.names=TRUE) basename(countFiles) flattenedFile = list.files(inDir, pattern="gff$", full.names=TRUE) basename(flattenedFile) ## ----sampleTable----------------------------------------------------------- sampleTable = data.frame( row.names = c( "treated1", "treated2", "treated3", "untreated1", "untreated2", "untreated3", "untreated4" ), condition = c("knockdown", "knockdown", "knockdown", "control", "control", "control", "control" ), libType = c( "single-end", "paired-end", "paired-end", "single-end", "single-end", "paired-end", "paired-end" ) ) ## ----check----------------------------------------------------------------- sampleTable ## ----message=FALSE--------------------------------------------------------- library( "DEXSeq" ) dxd = DEXSeqDataSetFromHTSeq( countFiles, sampleData=sampleTable, design= ~ sample + exon + condition:exon, flattenedfile=flattenedFile ) ## ----startVignette--------------------------------------------------------- genesForSubset = read.table( file.path(inDir, "geneIDsinsubset.txt"), stringsAsFactors=FALSE)[[1]] dxd = dxd[geneIDs( dxd ) %in% genesForSubset,] ## -------------------------------------------------------------------------- colData(dxd) ## -------------------------------------------------------------------------- head( counts(dxd), 5 ) ## -------------------------------------------------------------------------- split( seq_len(ncol(dxd)), colData(dxd)$exon ) ## -------------------------------------------------------------------------- head( featureCounts(dxd), 5 ) ## -------------------------------------------------------------------------- head( rowRanges(dxd), 3 ) ## -------------------------------------------------------------------------- sampleAnnotation( dxd ) ## -------------------------------------------------------------------------- dxd = estimateSizeFactors( dxd ) ## -------------------------------------------------------------------------- dxd = estimateDispersions( dxd ) ## ----fitdiagnostics, fig.small=TRUE, fig.cap="Fit Diagnostics. The initial per-exon dispersion estimates (shown by black points), the fitted mean-dispersion values function (red line), and the shrinked values in blue"---- plotDispEsts( dxd ) ## ----testForDEU1----------------------------------------------------------- dxd = testForDEU( dxd ) ## ----estFC----------------------------------------------------------------- dxd = estimateExonFoldChanges( dxd, fitExpToVar="condition") ## ----results1-------------------------------------------------------------- dxr1 = DEXSeqResults( dxd ) dxr1 ## ----results2-------------------------------------------------------------- mcols(dxr1)$description ## ----tallyExons------------------------------------------------------------ table ( dxr1$padj < 0.1 ) ## ----tallyGenes------------------------------------------------------------ table ( tapply( dxr1$padj < 0.1, dxr1$groupID, any ) ) ## ----MvsA, fig.small=TRUE, fig.cap="MA plot. Mean expression versus $log{2}$ fold change plot. Significant hits at an $FDR=0.1$ are coloured in red."---- plotMA( dxr1, cex=0.8 ) ## ----design---------------------------------------------------------------- sampleAnnotation(dxd) ## ----formulas2------------------------------------------------------------- formulaFullModel = ~ sample + exon + libType:exon + condition:exon formulaReducedModel = ~ sample + exon + libType:exon ## ----estDisps_again-------------------------------------------------------- dxd = estimateDispersions( dxd, formula = formulaFullModel ) ## ----test_again------------------------------------------------------------ dxd = testForDEU( dxd, reducedModel = formulaReducedModel, fullModel = formulaFullModel ) ## ----res_again------------------------------------------------------------- dxr2 = DEXSeqResults( dxd ) ## ----table2---------------------------------------------------------------- table( dxr2$padj < 0.1 ) ## ----table3---------------------------------------------------------------- table( before = dxr1$padj < 0.1, now = dxr2$padj < 0.1 ) ## ----plot1, fig.height=8, fig.width=12, fig.cap="Fitted expression. The plot represents the expression estimates from a call to `testForDEU()`. Shown in red is the exon that showed significant differential exon usage."---- plotDEXSeq( dxr2, "FBgn0010909", legend=TRUE, cex.axis=1.2, cex=1.3, lwd=2 ) ## ----checkClaim------------------------------------------------------------ wh = (dxr2$groupID=="FBgn0010909") stopifnot(sum(dxr2$padj[wh] < formals(plotDEXSeq)$FDR)==1) ## ----plot2, fig.height=8, fig.width=12, fig.cap="Transcripts. As in the previous plots, but including the annotated transcript models"---- plotDEXSeq( dxr2, "FBgn0010909", displayTranscripts=TRUE, legend=TRUE, cex.axis=1.2, cex=1.3, lwd=2 ) ## ----plot3, fig.height=8, fig.width=12, fig.cap="Normalized counts. As in the previous plots, with normalized count values of each exon in each of the samples."---- plotDEXSeq( dxr2, "FBgn0010909", expression=FALSE, norCounts=TRUE, legend=TRUE, cex.axis=1.2, cex=1.3, lwd=2 ) ## ----plot4, fig.height=8, fig.width=12, fig.cap="Fitted splicing. As in the previous plots, but after subtraction of overall changes in gene expression."---- plotDEXSeq( dxr2, "FBgn0010909", expression=FALSE, splicing=TRUE, legend=TRUE, cex.axis=1.2, cex=1.3, lwd=2 ) ## ----DEXSeqHTML, cache=TRUE, eval=FALSE------------------------------------ # DEXSeqHTML( dxr2, FDR=0.1, color=c("#FF000080", "#0000FF80") ) ## ----para1,cache=TRUE,results='hide', eval=FALSE--------------------------- # BPPARAM = MultiCoreParam(4) # dxd = estimateSizeFactors( dxd ) # dxd = estimateDispersions( dxd, BPPARAM=BPPARAM) # dxd = testForDEU( dxd, BPPARAM=BPPARAM) # dxd = estimateExonFoldChanges(dxd, BPPARAM=BPPARAM) ## ----alldeu---------------------------------------------------------------- dxr = DEXSeq(dxd) class(dxr) ## ----pergeneqval----------------------------------------------------------- numbOfGenes <- sum( perGeneQValue(dxr) < 0.1) numbOfGenes ## ----buildExonCountSetLoadPacks, eval=FALSE-------------------------------- # library(GenomicRanges) # library(GenomicFeatures) # library(GenomicAlignments) ## ----buildExonCountSetDownloadAnno, eval=FALSE----------------------------- # hse = makeTxDbFromBiomart( biomart="ensembl", dataset="hsapiens_gene_ensembl", host="grch37.ensembl.org" ) ## ----buildExonCountSetDisjoin,cache=TRUE, eval=FALSE----------------------- # exonicParts = exonicParts( hse, linked.to.single.gene.only = TRUE ) ## ----buildExonCountSet2FindBAMs, cache=TRUE, eval=FALSE-------------------- # bamDir = system.file( "extdata", package="parathyroidSE", mustWork=TRUE ) # fls = list.files( bamDir, pattern="bam$", full=TRUE ) ## ----buildExonCountSet2FindBAMs2, 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 ) ## ----buildExonCountSet3, eval=FALSE---------------------------------------- # colData(SE)$condition = c("A", "A", "B") # DEXSeqDataSetFromSE( SE, design= ~ sample + exon + condition:exon ) ## ----acc------------------------------------------------------------------- head( geneIDs(dxd) ) head( exonIDs(dxd) ) ## ----grmethods------------------------------------------------------------- interestingRegion = GRanges( "chr2L", IRanges(start=3872658, end=3875302) ) subsetByOverlaps( x=dxr, ranges=interestingRegion ) findOverlaps( query=dxr, subject=interestingRegion ) ## ----sessionInfo----------------------------------------------------------- sessionInfo()