## ----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()