## ----knitr, echo=FALSE, results="hide"----------------------------------- library("knitr") opts_chunk$set(tidy=FALSE,dev="png",fig.show="hide", fig.width=6.5,fig.height=5.5,fig.keep="high", message=FALSE) ## ----style, eval=TRUE, echo=FALSE, results="asis"------------------------ BiocStyle::latex() ## ----setup,echo=FALSE,results="hide"------------------------------------- options(width=80, signif=3, digits=3, prompt=" ", continue=" ") set.seed(0xdada) suppressWarnings({ library("JunctionSeq") library("JctSeqData") }) ## To create bitmap versions of plots with many dots, circumventing ## Sweave's fig=TRUE mechanism... ## (pdfs are too large) openBitmap = function(nm, rows=1, cols=1, height = 600, width = 800, cex = 1.2) { png(paste("QoRT-", nm, ".png", sep=""), width=width*cols, height=height*rows, pointsize=14) par(mfrow=c(rows, cols), cex=cex) } ## ----biocInstall, eval=FALSE--------------------------------------------- ## source("http://bioconductor.org/biocLite.R") ## biocLite("JunctionSeq") ## ----installer, eval=FALSE----------------------------------------------- ## source("http://hartleys.github.io/JunctionSeq/install/JS.install.R"); ## JS.install(); ## ----installReqPacks, eval=FALSE----------------------------------------- ## #CRAN package dependencies: ## install.packages("statmod"); ## install.packages("plotrix"); ## install.packages("stringr"); ## install.packages("Rcpp"); ## install.packages("RcppArmadillo"); ## install.packages("locfit"); ## install.packages("Hmisc"); ## ## #Bioconductor dependencies: ## source("http://bioconductor.org/biocLite.R"); ## biocLite(); ## biocLite("Biobase"); ## biocLite("BiocGenerics"); ## biocLite("BiocParallel"); ## biocLite("GenomicRanges"); ## biocLite("IRanges"); ## biocLite("S4Vectors"); ## biocLite("genefilter"); ## biocLite("geneplotter"); ## biocLite("SummarizedExperiment"); ## biocLite("DESeq2"); ## ## install.packages("http://hartleys.github.io/JunctionSeq/install/JunctionSeq_LATEST.tar.gz", ## repos = NULL, ## type = "source"); ## ----installExData, eval = FALSE----------------------------------------- ## install.packages("http://hartleys.github.io/JunctionSeq/install/JctSeqData_LATEST.tar.gz", repos=FALSE, type = "source"); ## ----GetExampleDataDirectory--------------------------------------------- decoder.file <- system.file("extdata/annoFiles/decoder.bySample.txt", package="JctSeqData", mustWork=TRUE); decoder <- read.table(decoder.file, header=TRUE, stringsAsFactors=FALSE); gff.file <- system.file( "extdata/cts/withNovel.forJunctionSeq.gff.gz", package="JctSeqData", mustWork=TRUE); ## ----GetCountFiles------------------------------------------------------- countFiles.noNovel <- system.file(paste0("extdata/cts/", decoder$sample.ID, "/QC.spliceJunctionAndExonCounts.forJunctionSeq.txt.gz"), package="JctSeqData", mustWork=TRUE); countFiles <- system.file(paste0("extdata/cts/", decoder$sample.ID, "/QC.spliceJunctionAndExonCounts.withNovel.forJunctionSeq.txt.gz"), package="JctSeqData", mustWork=TRUE); ## ----GetMiniExampleDataDirectory----------------------------------------- gff.file <- system.file( "extdata/tiny/withNovel.forJunctionSeq.gff.gz", package="JctSeqData", mustWork=TRUE); ## ----GetMiniCountFiles--------------------------------------------------- countFiles <- system.file(paste0("extdata/tiny/", decoder$sample.ID, "/QC.spliceJunctionAndExonCounts.withNovel.forJunctionSeq.txt.gz"), package="JctSeqData", mustWork=TRUE); ## ----testForDJU---------------------------------------------------------- jscs <- runJunctionSeqAnalyses(sample.files = countFiles, sample.names = decoder$sample.ID, condition=factor(decoder$group.ID), flat.gff.file = gff.file, nCores = 1, analysis.type = "junctionsAndExons" ); ## ----analysisHelp, eval=FALSE-------------------------------------------- ## help(runJunctionSeqAnalyses); ## ----writeSF------------------------------------------------------------- writeSizeFactors(jscs, file = "sizeFactors.txt"); ## ----testStep0, results="hide", warning=FALSE---------------------------- design <- data.frame(condition = factor(decoder$group.ID)); ## ----testStep0b, results="hide", warning=FALSE--------------------------- geneID.to.symbol.file <- system.file( "extdata/annoFiles/ensid.2.symbol.txt", package="JctSeqData", mustWork=TRUE); ## ----testStep1, results="hide", warning=FALSE---------------------------- jscs = readJunctionSeqCounts(countfiles = countFiles, samplenames = decoder$sample.ID, design = design, flat.gff.file = gff.file, gene.names = geneID.to.symbol.file ); ## ----testStep2, results="hide", warning=FALSE---------------------------- #Generate the size factors and load them into the JunctionSeqCountSet: jscs <- estimateJunctionSeqSizeFactors(jscs); ## ----testStep3, results="hide", warning=FALSE---------------------------- jscs <- estimateJunctionSeqDispersions(jscs, nCores = 1); ## ----testStep4, results="hide", warning=FALSE---------------------------- jscs <- fitJunctionSeqDispersionFunction(jscs); ## ----testStep5, results="hide", warning=FALSE---------------------------- jscs <- testForDiffUsage(jscs, nCores = 1); ## ----testStep6, results="hide", warning=FALSE---------------------------- jscs <- estimateEffectSizes( jscs, nCores = 1); ## ----testStepHelp, eval=FALSE-------------------------------------------- ## help(readJunctionSeqCounts); ## help(estimateJunctionSeqSizeFactors); ## help(estimateJunctionSeqDispersions); ## help(fitJunctionSeqDispersionFunction); ## help(testForDiffUsage); ## help(estimateEffectSizes); ## ----makeRes, results="hide"--------------------------------------------- writeCompleteResults(jscs, outfile.prefix="./test", save.jscs = TRUE ); ## ----makeAllPlots-------------------------------------------------------- buildAllPlots(jscs=jscs, outfile.prefix = "./plots/", use.plotting.device = "png", FDR.threshold = 0.01 ); ## ----summaryplots, results="hide", fig.width=6.5, fig.height=5.5--------- plotDispEsts(jscs); plotMA(jscs, FDR.threshold=0.05); ## ----makeAllPlots2------------------------------------------------------- buildAllPlots(jscs=jscs, outfile.prefix = "./plots2/", use.plotting.device = "png", FDR.threshold = 0.01, expr.plot = FALSE, normCounts.plot = FALSE, rExpr.plot = FALSE, rawCounts.plot = TRUE ); ## ----makeAllPlots3, eval=FALSE------------------------------------------- ## #Make a battery of exon-only plots for one gene only: ## buildAllPlotsForGene(jscs=jscs, geneID = "ENSRNOG00000009281", ## outfile.prefix = "./plots/", ## use.plotting.device = "png", ## colorRed.FDR.threshold = 0.01, ## #Limit plotting to exons only: ## plot.junction.results = FALSE, ## #Change the fill color of significant exons: ## colorList = list(SIG.FEATURE.FILL.COLOR = "red"), ## ); ## ## #Make a set of Junction-Only Plots for a specific list of interesting genes: ## buildAllPlots(jscs=jscs, ## gene.list = c("ENSRNOG00000009281"), ## outfile.prefix = "./plotsJCT/", ## use.plotting.device = "png", ## FDR.threshold = 0.01, ## #Do not graph exonic regions: ## plot.exon.results = FALSE, ## ); ## ----makeAllPlots3BG, echo=FALSE, results="hide"------------------------- #NOTE: This is a version of the above script with additional options to # reduce execution time: buildAllPlotsForGene(jscs=jscs, geneID = "ENSRNOG00000009281", outfile.prefix = "./plots/", use.plotting.device = "png", colorRed.FDR.threshold = 0.01, #Limit plotting to exons only: plot.junction.results = FALSE, #Change the fill color of significant exons: colorList = list(SIG.FEATURE.FILL.COLOR = "red"), #NOTE: The following options reduce the generation of # extra plots not actually needed for this manual: expr.plot = TRUE, normCounts.plot = FALSE, rExpr.plot = FALSE, rawCounts.plot = FALSE, without.TX = FALSE ); #Junction-Only Plots: buildAllPlots(jscs=jscs, gene.list = c("ENSRNOG00000009281"), outfile.prefix = "./plotsJCT/", use.plotting.device = "png", FDR.threshold = 0.01, plot.exon.results = FALSE, #NOTE: The following options reduce the generation of # extra plots not actually needed for this manual: expr.plot = TRUE, normCounts.plot = FALSE, rExpr.plot = FALSE, rawCounts.plot = FALSE, ma.plot = FALSE, variance.plot = FALSE, without.TX = FALSE, writeHTMLresults = FALSE ); ## ----getPlotHelpDocs, eval=FALSE----------------------------------------- ## help(buildAllPlots); ## help(buildAllPlotsForGene); ## help(plotJunctionSeqResultsForGene); ## ----DEXReplicate, eval=FALSE-------------------------------------------- ## jscs.DEX <- runJunctionSeqAnalyses(sample.files = countFiles, ## sample.names = decoder$sample.ID, ## condition = factor(decoder$group.ID), ## flat.gff.file = gff.file, ## nCores = 1, ## analysis.type = "exonsOnly", ## method.countVectors = "sumOfAllBinsForGene", ## method.sizeFactors = "byCountbins", ## method.expressionEstimation = "feature-vs-otherFeatures", ## meanCountTestableThreshold = "auto", ## use.multigene.aggregates = TRUE, ## method.cooksFilter = TRUE, ## optimizeFilteringForAlpha = 0.1 ## ); ## ----DEXReplicateRes, eval=FALSE----------------------------------------- ## writeCompleteResults(jscs.DEX, ## outfile.prefix="./testDEX", ## save.jscs = TRUE ## ); ## buildAllPlots(jscs=jscs.DEX, ## outfile.prefix = "./plotsDEX/", ## use.plotting.device = "png", ## FDR.threshold = 0.01 ## ); ## ----DEXRun, eval=FALSE-------------------------------------------------- ## library("DEXSeq"); ## countFiles.dexseq <- system.file(paste0("extdata/cts/", ## decoder$sample.ID, ## "/QC.exonCounts.formatted.for.DEXSeq.txt.gz"), ## package="JctSeqData", mustWork=TRUE); ## gff.dexseq <- system.file("extdata/annoFiles/DEXSeq.flat.gff.gz", ## package="JctSeqData", mustWork=TRUE); ## dxd <- DEXSeqDataSetFromHTSeq( ## countFiles.dexseq, ## sampleData = design, ## design = ~sample + exon + condition:exon, ## flattenedfile = gff.dexseq ## ); ## dxd <- estimateSizeFactors(dxd); ## dxd <- estimateDispersions(dxd); ## dxd <- testForDEU( dxd); ## dxd <- estimateExonFoldChanges(dxd, fitExpToVar = "condition"); ## dxr <- results(dxd); ## write.table(dxr,file="dxr.out.txt"); ## ----createComplexVar, results="hide", warning=FALSE--------------------- threeLevelVariable <- c("GroupA","GroupA", "GroupB","GroupB", "GroupC","GroupC"); ## ----testForDJUComplexVar, results="hide", warning=FALSE, eval = FALSE---- ## jscs <- runJunctionSeqAnalyses(sample.files = countFiles, ## sample.names = decoder$sample.ID, ## condition=factor(threeLevelVariable), ## flat.gff.file = gff.file, ## nCores = 1, ## analysis.type = "junctionsAndExons" ## ); ## ----createCovar--------------------------------------------------------- #Artificially adding additional samples by using two of the samples twice: # (Note: this is purely for use as an example. Never do this.) countFiles.8 <- c(countFiles, countFiles[3],countFiles[6]); #Make up some sample names, conditions, and covariates for these samples: decoder.8 <- data.frame( sample.names = factor(paste0("SAMP",1:8)), condition = factor(rep(c("CASE","CTRL"),each=4)), smokeStatus = factor(rep(c("Smoker","Nonsmoker"),4)) ) print(decoder.8); ## ----testForDJUCovar, results="hide", warning=FALSE, eval = FALSE-------- ## jscs <- runJunctionSeqAnalyses(sample.files = countFiles.8, ## sample.names = decoder.8$sample.names, ## condition= decoder.8$condition, ## use.covars = decoder.8[,"smokeStatus",drop=F], ## flat.gff.file = gff.file, ## nCores = 1, ## analysis.type = "junctionsAndExons", ## test.formula0 = ~ sample + countbin + smokeStatus : countbin, ## test.formula1 = ~ sample + countbin + smokeStatus : countbin + condition : countbin, ## effect.formula = ~ condition + smokeStatus + countbin + ## smokeStatus : countbin + condition : countbin, ## geneLevel.formula = ~ smokeStatus + condition ## ); ## ----sessionInfo--------------------------------------------------------- sessionInfo()