##############################################################################
##############################################################################
###
### Running command:
###
### C:\Users\biocbuild\bbs-3.11-bioc\R\bin\R.exe CMD check --force-multiarch --install=check:JunctionSeq.install-out.txt --library=C:\Users\biocbuild\bbs-3.11-bioc\R\library --no-vignettes --timings JunctionSeq_1.18.0.tar.gz
###
##############################################################################
##############################################################################
* using log directory 'C:/Users/biocbuild/bbs-3.11-bioc/meat/JunctionSeq.Rcheck'
* using R version 4.0.3 (2020-10-10)
* using platform: x86_64-w64-mingw32 (64-bit)
* using session charset: ISO8859-1
* using option '--no-vignettes'
* checking for file 'JunctionSeq/DESCRIPTION' ... OK
* this is package 'JunctionSeq' version '1.18.0'
* checking package namespace information ... OK
* checking package dependencies ... OK
* checking if this is a source package ... OK
* checking if there is a namespace ... OK
* checking for hidden files and directories ... OK
* checking for portable file names ... OK
* checking whether package 'JunctionSeq' can be installed ... OK
* checking installed package size ... OK
* checking package directory ... OK
* checking 'build' directory ... OK
* checking DESCRIPTION meta-information ... OK
* checking top-level files ... OK
* checking for left-over files ... OK
* checking index information ... OK
* checking package subdirectories ... OK
* checking R files for non-ASCII characters ... OK
* checking R files for syntax errors ... OK
* loading checks for arch 'i386'
** checking whether the package can be loaded ... OK
** checking whether the package can be loaded with stated dependencies ... OK
** checking whether the package can be unloaded cleanly ... OK
** checking whether the namespace can be loaded with stated dependencies ... OK
** checking whether the namespace can be unloaded cleanly ... OK
* loading checks for arch 'x64'
** checking whether the package can be loaded ... OK
** checking whether the package can be loaded with stated dependencies ... OK
** checking whether the package can be unloaded cleanly ... OK
** checking whether the namespace can be loaded with stated dependencies ... OK
** checking whether the namespace can be unloaded cleanly ... OK
* checking dependencies in R code ... NOTE
Packages in Depends field not imported from:
'Rcpp' 'RcppArmadillo'
These packages need to be imported from (in the NAMESPACE file)
for when this namespace is loaded but not attached.
* checking S3 generic/method consistency ... OK
* checking replacement functions ... OK
* checking foreign function calls ... OK
* checking R code for possible problems ... OK
* checking Rd files ... NOTE
prepare_Rd: setJunctionSeqCompiledSourcePackage.Rd:23-25: Dropping empty section \examples
* checking Rd metadata ... OK
* checking Rd cross-references ... OK
* checking for missing documentation entries ... OK
* checking for code/documentation mismatches ... OK
* checking Rd \usage sections ... WARNING
Documented arguments not in \usage in documentation object 'estimateJunctionSeqSizeFactors':
'...'
Documented arguments not in \usage in documentation object 'fitJunctionSeqDispersionFunction':
'...'
Functions with \usage entries need to have the appropriate \alias
entries, and all their arguments documented.
The \usage entries must correspond to syntactically valid R code.
See chapter 'Writing R documentation files' in the 'Writing R
Extensions' manual.
* checking Rd contents ... OK
* checking for unstated dependencies in examples ... OK
* checking line endings in C/C++/Fortran sources/headers ... OK
* checking line endings in Makefiles ... OK
* checking compilation flags in Makevars ... OK
* checking for GNU extensions in Makefiles ... OK
* checking for portable use of $(BLAS_LIBS) and $(LAPACK_LIBS) ... OK
* checking use of PKG_*FLAGS in Makefiles ... OK
* checking compiled code ... NOTE
Note: information on .o files for i386 is not available
Note: information on .o files for x64 is not available
File 'C:/Users/biocbuild/bbs-3.11-bioc/R/library/JunctionSeq/libs/i386/JunctionSeq.dll':
Found 'abort', possibly from 'abort' (C), 'runtime' (Fortran)
Found 'exit', possibly from 'exit' (C), 'stop' (Fortran)
Found 'printf', possibly from 'printf' (C)
File 'C:/Users/biocbuild/bbs-3.11-bioc/R/library/JunctionSeq/libs/x64/JunctionSeq.dll':
Found 'abort', possibly from 'abort' (C), 'runtime' (Fortran)
Found 'exit', possibly from 'exit' (C), 'stop' (Fortran)
Found 'printf', possibly from 'printf' (C)
Compiled code should not call entry points which might terminate R nor
write to stdout/stderr instead of to the console, nor use Fortran I/O
nor system RNGs. The detected symbols are linked into the code but
might come from libraries and not actually be called.
See 'Writing portable packages' in the 'Writing R Extensions' manual.
* checking installed files from 'inst/doc' ... OK
* checking files in 'vignettes' ... OK
* checking examples ...
** running examples for arch 'i386' ... ERROR
Running examples in 'JunctionSeq-Ex.R' failed
The error most likely occurred in:
> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: buildAllPlots
> ### Title: Create and save a full battery of JunctionSeq expression plots.
> ### Aliases: buildAllPlots
>
> ### ** Examples
>
>
> data(exampleDataSet,package="JctSeqData");
> buildAllPlots(jscs);
> buildAllPlots: Found 5 genes with at least one significant exon, at adjusted-p-value threshold 0.01
> buildAllPlots: Starting plotting...
Warning in dir.create(outfile.prefix) : '.' already exists
> buildAllPlots: Generating Dispersion Plot
Starting device: png (./dispersion-plot.png)
abundance ranges from 0.146009972023612 to 5452.18841073905
dispersion ranges from 1.94121697570967e-06 to 2.08746163174998
Plotting dispersions from 2.71367228072472e-05 to 2.08746163174998
> buildAllPlots: Generating MA-Plot (log2FC(CTRL/CASE))
Starting device: png (./ma-plot-log2FC(CTRLvsCASE).png)
> buildAllPlots: Writing HTML results index.
Copying default css stylesheet.
Writing html index. Sat Oct 17 05:30:35 2020
Writing formula data. Sat Oct 17 05:30:35 2020
Writing methods data. Sat Oct 17 05:30:35 2020
Writing sample data. Sat Oct 17 05:30:35 2020
Writing dispersion data. Sat Oct 17 05:30:35 2020
Writing summary plots. Sat Oct 17 05:30:35 2020
Compiling data table. Sat Oct 17 05:30:35 2020
Writing data table. Sat Oct 17 05:30:35 2020
Html index complete. Sat Oct 17 05:30:35 2020
Writing pages. Sat Oct 17 05:30:35 2020
Finished all html files. Sat Oct 17 05:30:36 2020
> buildAllPlots: Finished writing HTML results index.
> buildAllPlots: starting geneID: ENSRNOG00000008079 (1 of 5)
starting buildAllPlotsForGene() for geneID: ENSRNOG00000008079
Starting nested heights...
Starting device: png (.//exprTX/1-expr-TX.png)
> pJSRfG(): ENSRNOG00000008079, plot.type: expr (with TX)
Starting device: png (.//expr/1-expr.png)
> pJSRfG(): ENSRNOG00000008079, plot.type: expr
Starting device: png (.//normCountsTX/1-normCts-TX.png)
> pJSRfG(): ENSRNOG00000008079, plot.type: normCounts (with TX)
Starting device: png (.//normCounts/1-normCts.png)
> pJSRfG(): ENSRNOG00000008079, plot.type: normCounts
Starting device: png (.//rExprTX/1-rExpr-TX.png)
> pJSRfG(): ENSRNOG00000008079, plot.type: rExpr (with TX)
Starting device: png (.//rExpr/1-rExpr.png)
> pJSRfG(): ENSRNOG00000008079, plot.type: rExpr
> buildAllPlots: starting geneID: ENSRNOG00000009281 (2 of 5)
starting buildAllPlotsForGene() for geneID: ENSRNOG00000009281
Starting nested heights...
Starting device: png (.//exprTX/2-expr-TX.png)
> pJSRfG(): ENSRNOG00000009281, plot.type: expr (with TX)
----------- FAILURE REPORT --------------
--- failure: length > 1 in coercion to logical ---
--- srcref ---
:
--- package (from environment) ---
JunctionSeq
--- call from context ---
doTryCatch(return(expr), name, parentenv, handler)
--- call from argument ---
any(c(pw[1] - 1, pw[1] + 1) %in% pwo) && pw%%2 == 0
--- R stacktrace ---
where 1: doTryCatch(return(expr), name, parentenv, handler)
where 2: tryCatchOne(expr, names, parentenv, handlers[[1L]])
where 3: tryCatchList(expr, classes, parentenv, handlers)
where 4: tryCatch({
if (debug.mode)
message("> Step 6.9b")
pval.drop <- rep(-pval.height, length(p.values))
pw <- which(p.values != "")
pwo <- pw
while (length(pw) > 0) {
if (any(c(pw[1] - 1, pw[1] + 1) %in% pwo) && pw%%2 ==
0) {
pval.drop[pw[1]] <- pval.drop[pw[1]] + pval.height
}
else {
}
pw <- pw[-1]
}
pval.y <- pval.ceiling + pval.drop
}, error = function(e) {
message("WARNING: staggering p-values failed. Falling back to simpler method. This warning should never appear. If you see this, reporting this to the developer would be appreciated.")
})
where 5: drawPlot(matr = count, ylimn, jscs, intervals, rango, textAxis = y.axis.title,
geneLevelAxisTitle = y.axis.title.right, rt = rt, color.count = color.count,
colorlines = vertline.col, countbinIDs = merged.data$countbinID[rt],
use.vst = use.vst, use.log = use.log, plot.type = plot.type,
main.title = main.title, draw.legend = draw.legend, color.key = color,
condition.names = condition.names, p.values = p.values.labels,
draw.p.values = label.p.vals, plot.lwd = plot.lwd, axes.lwd = axes.lwd,
anno.lwd = anno.lwd, par.cex = par.cex, anno.cex.text = anno.cex.text,
anno.cex.axis = anno.cex.axis, anno.cex.main = anno.cex.main,
fit.countbin.names = fit.countbin.names, debug.mode = debug.mode,
plot.gene.level.expression = plot.gene.level.expression,
geneCount = geneCount, color.geneCount = color.geneCount,
yAxisLabels.inExponentialForm = yAxisLabels.inExponentialForm,
italicize.label = italicize.label, condition.legend.text = condition.legend.text,
annolink.col = annolink.col, exonlty = exonlty, graph.margins = graph.margins,
plotWindowXmax = plotWindowXmax, fit.labels = fit.labels,
...)
where 6: doTryCatch(return(expr), name, parentenv, handler)
where 7: tryCatchOne(expr, names, parentenv, handlers[[1L]])
where 8: tryCatchList(expr, classes, parentenv, handlers)
where 9: tryCatch({
GENE.annotation.height <- GENE.annotation.relative.height *
10
TX.annotation.height <- TX.annotation.relative.height * 10
CONNECTIONS.height <- CONNECTIONS.relative.height * 10
SPLICE.annotation.height <- SPLICE.annotation.relative.height *
10
if (!any(geneID == fData(jscs)$geneID)) {
stop(paste0("FATAL ERROR: GeneID \"", geneID, "\" not found in dataset!"))
}
flat.gff.data <- jscs@flatGffData
merged.data <- fData(jscs)
condition <- jscs@phenoData$condition
exonRescaleFunction <- match.arg(exonRescaleFunction)
sequencing.type <- match.arg(sequencing.type)
splice.junction.drawing.style <- match.arg(splice.junction.drawing.style)
plot.type <- match.arg(plot.type)
truncateBelowOne <- TRUE
if (is.null(plot.exon.results)) {
plot.exon.results <- any(merged.data$featureType == "exonic_part")
}
if (is.null(plot.junction.results)) {
plot.junction.results <- any(merged.data$featureType ==
"splice_site" | merged.data$featureType == "novel_splice_site")
}
if (is.null(plot.novel.junction.results)) {
if (plot.junction.results) {
plot.novel.junction.results <- any(merged.data$featureType ==
"novel_splice_site")
}
else {
plot.novel.junction.results <- FALSE
}
}
flip.splicing <- if (plot.junction.results) {
FALSE
}
else {
TRUE
}
geneName <- jscs@flatGffGeneData$gene_name[jscs@flatGffGeneData$geneID ==
geneID]
if (is.null(condition.legend.text)) {
condition.legend.text <- levels(jscs@phenoData$condition)
names(condition.legend.text) <- condition.legend.text
}
else {
if (is.null(names(condition.legend.text))) {
warning("names(condition.legend.text is NULL! condition.legend.text mis-formatted. Must be a list or character vector, with element names equal to the levels of pData(jscs)$condition. Falling back.")
condition.legend.text <- levels(jscs@phenoData$condition)
names(condition.legend.text) <- condition.legend.text
}
else if (!all(levels(jscs@phenoData$condition) %in% names(condition.legend.text))) {
warning("Not all levels contained in names(condition.legend.text)! condition.legend.text mis-formatted. Must be a list or character vector, with element names equal to the levels of pData(jscs)$condition. Falling back.")
condition.legend.text <- levels(jscs@phenoData$condition)
names(condition.legend.text) <- condition.legend.text
}
}
if (!is.null(jscs@flatGffGeneData[["aggregateGeneStrand"]])) {
geneStrand <- as.character(jscs@flatGffGeneData[["aggregateGeneStrand"]][jscs@flatGffGeneData[["geneID"]] ==
geneID])
txSetString <- as.character(jscs@flatGffGeneData[["tx_set"]][jscs@flatGffGeneData[["geneID"]] ==
geneID])
txStrandString <- as.character(jscs@flatGffGeneData[["tx_strands"]][jscs@flatGffGeneData[["geneID"]] ==
geneID])
txSet <- strsplit(as.character(txSetString), "+", fixed = TRUE)[[1]]
txStrand <- strsplit(as.character(txStrandString), ",",
fixed = TRUE)[[1]]
txStrandMap <- as.list(txStrand)
names(txStrandMap) <- txSet
}
else {
txStrandMap <- list()
geneStrand <- "."
}
final.color.list <- overmerge.list(JUNCTIONSEQ.DEFAULT.COLOR.LIST,
colorList)
SIG.VERTLINE.COLOR = final.color.list[["SIG.VERTLINE.COLOR"]]
NOSIG.VERTLINE.COLOR = final.color.list[["NOSIG.VERTLINE.COLOR"]]
UNTESTABLE.VERTLINE.COLOR = final.color.list[["UNTESTABLE.VERTLINE.COLOR"]]
SIG.FEATURE.COLOR = final.color.list[["SIG.FEATURE.COLOR"]]
NOSIG.FEATURE.COLOR = final.color.list[["NOSIG.FEATURE.COLOR"]]
UNTESTABLE.FEATURE.COLOR = final.color.list[["UNTESTABLE.FEATURE.COLOR"]]
EXCLUDED.FEATURE.COLOR = final.color.list[["EXCLUDED.FEATURE.COLOR"]]
SIG.FEATURE.BORDER.COLOR = final.color.list[["SIG.FEATURE.BORDER.COLOR"]]
NOSIG.FEATURE.BORDER.COLOR = final.color.list[["NOSIG.FEATURE.BORDER.COLOR"]]
UNTESTABLE.FEATURE.BORDER.COLOR = final.color.list[["UNTESTABLE.FEATURE.BORDER.COLOR"]]
EXCLUDED.FEATURE.BORDER.COLOR = final.color.list[["EXCLUDED.FEATURE.BORDER.COLOR"]]
SIG.FEATURE.FILL.COLOR = final.color.list[["SIG.FEATURE.FILL.COLOR"]]
NOSIG.FEATURE.FILL.COLOR = final.color.list[["NOSIG.FEATURE.FILL.COLOR"]]
UNTESTABLE.FEATURE.FILL.COLOR = final.color.list[["UNTESTABLE.FEATURE.FILL.COLOR"]]
EXCLUDED.FEATURE.FILL.COLOR = final.color.list[["EXCLUDED.FEATURE.FILL.COLOR"]]
PLOTTING.LINE.COLORS = final.color.list[["PLOTTING.LINE.COLORS"]]
gene.level.buffer <- 0.5
if (is.null(plot.gene.level.expression)) {
if (use.vst) {
plot.gene.level.expression <- FALSE
}
else {
plot.gene.level.expression <- TRUE
}
}
if (plot.gene.level.expression && use.vst) {
warning("WARNING: plotting of gene-level expression is not supported for vst-transformed plots. Errors are likely to follow.")
}
FDR <- colorRed.FDR.threshold
if (verbose) {
displayTXstring <- if (displayTranscripts)
" (with TX)"
else ""
message("> pJSRfG(): ", geneID, ", plot.type: ", plot.type,
displayTXstring)
}
chrom.label <- as.character(merged.data$chr[merged.data$geneID ==
geneID][1])
rt <- merged.data$geneID == geneID
if (!plot.exon.results) {
if (debug.mode)
message("> Removing ", sum(rt & merged.data$featureType ==
"exonic_part"), " exonic_part features. ", sum(rt &
merged.data$featureType != "exonic_part"), " features remaining")
rt <- rt & merged.data$featureType != "exonic_part"
}
if (!plot.junction.results) {
if (debug.mode)
message("> Removing ", sum(rt & merged.data$featureType ==
"splice_site"), " splice_site features. ", sum(rt &
merged.data$featureType != "splice_site"), " features remaining")
rt <- rt & merged.data$featureType != "splice_site"
}
if (!plot.novel.junction.results) {
if (debug.mode)
message("> Removing ", sum(rt & merged.data$featureType ==
"novel_splice_site"), " novel_splice_site features. ",
sum(rt & merged.data$featureType != "novel_splice_site"),
" features remaining")
rt <- rt & merged.data$featureType != "novel_splice_site"
}
untestable.rt <- which(rt & (!merged.data$testable))
if (!plot.untestable.results) {
if (debug.mode)
message("> Removing ", sum(rt & !merged.data$testable),
" untestable features. ", sum(rt & merged.data$testable),
" features remaining")
rt <- rt & merged.data$testable
}
if (debug.mode)
message("> Plotting ", sum(rt), " features.")
rt <- which(rt)
if (length(rt) == 0) {
message("NO FEATURES TO PLOT!")
}
else {
rt.allExon <- which(flat.gff.data$gene_id == geneID &
flat.gff.data$featureType == "exonic_part")
rango.allExon <- 1:length(rt.allExon)
rt.allJunction <- which(flat.gff.data$gene_id == geneID &
(flat.gff.data$featureType == "splice_site" | flat.gff.data$featureType ==
"novel_splice_site"))
rango.allJunction <- 1:length(rt.allJunction)
rescale.iv <- generate.interval.scale(data.frame(start = c(flat.gff.data$start[rt.allExon],
flat.gff.data$start[rt.allJunction]), end = c(flat.gff.data$end[rt.allExon],
flat.gff.data$end[rt.allJunction]), is.exon = c(rep(TRUE,
length(rt.allExon)), rep(FALSE, length(rt.allJunction)))),
exon.rescale.factor, exonRescaleFunction, debug.mode = debug.mode)
rel <- data.frame(start = rescale.coords(merged.data$start[rt],
rescale.iv), end = rescale.coords(merged.data$end[rt],
rescale.iv))
if (sort.features) {
rt <- rt[order(((rel$end - rel$start)/2) + rel$start)]
}
rango <- 1:length(rt)
draw.legend <- TRUE
condition.names <- levels(condition)
sample.names <- sampleNames(jscs@phenoData)
colorcode.title <- TRUE
numcond <- length(condition.names)
if (debug.mode && verbose)
message("> pJSRforGene(): ", "Reached step 2.")
default.color.list <- PLOTTING.LINE.COLORS
if (numcond > length(default.color.list)) {
message("Too many condition values, the default color selection may not look good! Set your own colors by setting the \"color\" parameter.")
color <- rgb(colorRamp(c("#D7191C", "#FFFFBF", "#2B83BA"))(seq(0,
1, length.out = numcond)), maxColorValue = 255,
alpha = 175)
}
else {
color <- color2transparentVector(default.color.list[1:numcond],
t = 175)
}
if (debug.mode && verbose)
message("> pJSRforGene(): ", "color = ", paste0(color,
collapse = ","))
if (debug.mode && verbose)
message("> pJSRforGene(): ", "length(color) = ",
length(color))
if (debug.mode && verbose)
message("> pJSRforGene(): ", "length(condition.names) = ",
length(condition.names))
if (debug.mode && verbose)
message("> pJSRforGene(): ", "condition.names = ",
paste0(condition.names, collapse = ","))
names(color) <- condition.names
y.axis.title <- ""
main.title <- ""
if (debug.mode && verbose)
message("> pJSRforGene(): ", "Reached step 3.")
if (truncateBelowOne) {
convertY <- function(y) {
ifelse(y < 0, ifelse(is.infinite(y), INTERNAL.NINF.VALUE,
(1 - exp(y)) * INTERNAL.NINF.VALUE), y)
}
}
else {
convertY <- function(y) {
y
}
}
if (plot.type == "rExpr") {
count <- if (use.vst) {
vst(jscs@plottingEstimates[["relExprEstimate"]][rt,
, drop = FALSE], jscs)
}
else if (use.log) {
apply(log10(jscs@plottingEstimates[["relExprEstimate"]][rt,
, drop = FALSE]), c(1, 2), FUN = convertY)
}
else {
jscs@plottingEstimates[["relExprEstimate"]][rt,
, drop = FALSE]
}
color.count <- rep(color[condition.names], each = nrow(count))
y.axis.title <- "Relative Coverage"
y.axis.title.right <- "Gene-Level Mean Normalized Counts"
main.title <- paste0("Relative Coverage (", geneName,
")")
if (plot.gene.level.expression) {
geneCount <- if (use.vst) {
plot.gene.level.expression <- FALSE
NULL
}
else if (use.log) {
sapply(log10(jscs@geneLevelPlottingEstimates[["geneLevelEstModeled"]][rownames(jscs@geneLevelPlottingEstimates[["geneLevelEstModeled"]]) ==
geneID, ]), FUN = convertY)
}
else {
jscs@geneLevelPlottingEstimates[["geneLevelEstModeled"]][rownames(jscs@geneLevelPlottingEstimates[["geneLevelEstModeled"]]) ==
geneID, ]
}
color.geneCount <- color[condition.names]
}
}
else if (plot.type == "normCounts") {
count <- if (use.vst) {
vst(jscs@plottingEstimates[["normCounts"]][rt,
1:length(sample.names), drop = FALSE], jscs)
}
else if (use.log) {
apply(log10(jscs@plottingEstimates[["normCounts"]][rt,
1:length(sample.names), drop = FALSE]), c(1,
2), FUN = convertY)
}
else {
jscs@plottingEstimates[["normCounts"]][rt, 1:length(sample.names),
drop = FALSE]
}
if (plot.gene.level.expression) {
geneCount <- if (use.vst) {
plot.gene.level.expression <- FALSE
NULL
}
else if (use.log) {
sapply(log10(jscs@geneCountData[rownames(jscs@geneCountData) ==
geneID, ]/sizeFactors(jscs)), FUN = convertY)
}
else {
jscs@geneCountData[rownames(jscs@geneCountData) ==
geneID, ]/sizeFactors(jscs)
}
color.geneCount <- color[as.character(condition)]
}
color.count <- rep(color[as.character(condition)],
each = nrow(count))
y.axis.title <- "Normalized Counts"
y.axis.title.right <- "Gene-Level Normalized Counts"
main.title <- paste0("Normalized Counts (", geneName,
")")
}
else if (plot.type == "rawCounts") {
count <- if (use.vst) {
vst(jscs@countVectors[rt, 1:length(sample.names),
drop = FALSE], jscs)
}
else if (use.log) {
apply(log10(jscs@countVectors[rt, 1:length(sample.names),
drop = FALSE]), c(1, 2), FUN = function(y) {
max(INTERNAL.NINF.VALUE, y)
})
}
else {
jscs@countVectors[rt, 1:length(sample.names),
drop = FALSE]
}
if (plot.gene.level.expression) {
geneCount <- if (use.vst) {
plot.gene.level.expression <- FALSE
NULL
}
else if (use.log) {
sapply(log10(jscs@geneCountData[rownames(jscs@geneCountData) ==
geneID, ]), FUN = convertY)
}
else {
jscs@geneCountData[rownames(jscs@geneCountData) ==
geneID, ]
}
color.geneCount <- color[as.character(condition)]
}
color.count <- rep(color[as.character(condition)],
each = nrow(count))
y.axis.title <- "Raw Counts"
y.axis.title.right <- "Raw Gene-Level Counts"
main.title <- paste0("Raw Counts (", geneName, ")")
}
else if (plot.type == "expr") {
count <- if (use.vst) {
vst(jscs@plottingEstimates[["exprEstimate"]][rt,
, drop = FALSE], jscs)
}
else if (use.log) {
apply(log10(jscs@plottingEstimates[["exprEstimate"]][rt,
, drop = FALSE]), c(1, 2), FUN = convertY)
}
else {
jscs@plottingEstimates[["exprEstimate"]][rt,
, drop = FALSE]
}
if (plot.gene.level.expression) {
geneCount <- if (use.vst) {
plot.gene.level.expression <- FALSE
NULL
}
else if (use.log) {
sapply(log10(jscs@geneLevelPlottingEstimates[["geneLevelEstModeled"]][rownames(jscs@geneLevelPlottingEstimates[["geneLevelEstModeled"]]) ==
geneID, ]), FUN = convertY)
}
else {
jscs@geneLevelPlottingEstimates[["geneLevelEstModeled"]][rownames(jscs@geneLevelPlottingEstimates[["geneLevelEstModeled"]]) ==
geneID, ]
}
color.geneCount <- color[condition.names]
}
color.count <- rep(color[condition.names], each = nrow(count))
y.axis.title <- "Mean Normalized Counts"
y.axis.title.right <- "Gene-Level Mean Normalized Counts"
main.title <- paste0("Mean Normalized Coverage (",
geneName, ")")
}
else {
stop(paste0("FATAL ERROR: Unknown plot type! plot.type = \"",
plot.type, "\""))
}
fragment.label <- if (sequencing.type == "paired-end")
"Read-Pair"
else "Read"
if (plot.type == "rExpr") {
main.title <- paste0("Relative Coverage (", geneName,
")")
y.axis.title <- paste0("Relative Coverage")
y.axis.title.right <- paste0(fragment.label, "s per Sample, Gene-Level")
}
else if (plot.type == "normCounts") {
main.title <- paste0("Normalized Counts (", geneName,
")")
y.axis.title <- paste0("Normalized Counts")
y.axis.title.right <- paste0("Normalized Counts, Gene-Level")
}
else if (plot.type == "rawCounts") {
main.title <- paste0("Raw Counts (", geneName, ")")
y.axis.title <- paste0("Raw ", fragment.label, " Counts")
y.axis.title.right <- paste0("Raw ", fragment.label,
" Counts, Gene-Level")
}
else if (plot.type == "expr") {
main.title <- paste0("Mean Normalized Coverage (",
geneName, ")")
y.axis.title <- paste0(fragment.label, "s per Sample")
y.axis.title.right <- paste0(fragment.label, "s per Sample, Gene-Level")
}
else {
stop(paste0("FATAL ERROR: Unknown plot type! plot.type = \"",
plot.type, "\""))
}
if (!is.null(title.main))
main.title <- title.main
if (!is.null(title.ylab))
y.axis.title <- title.ylab
if (!is.null(title.ylab.right))
y.axis.title.right <- title.ylab.right
if (!plot.gene.level.expression) {
geneCount <- NULL
color.geneCount <- NULL
}
count <- as.matrix(count)
if (debug.mode && verbose)
message("> pJSRforGene(): ", "Reached step 4.")
intervals <- (0:nrow(count))/nrow(count)
numexons <- nrow(count)
each <- merged.data$padjust[rt]
vertline.col <- ifelse(merged.data$testable[rt], ifelse(f.na(merged.data$padjust[rt] <=
FDR), SIG.VERTLINE.COLOR, NOSIG.VERTLINE.COLOR),
UNTESTABLE.VERTLINE.COLOR)
annolink.col <- ifelse(merged.data$testable[rt], ifelse(f.na(merged.data$padjust[rt] <=
FDR), SIG.FEATURE.COLOR, NOSIG.FEATURE.COLOR), UNTESTABLE.FEATURE.COLOR)
exonlty <- rep(final.color.list[["EXON.CONNECTION.LTY"]],
length(vertline.col))
exonlty[as.character(merged.data$featureType[rt]) ==
"novel_splice_site"] <- final.color.list[["NOVEL.SPLICE.CONNECTION.LTY"]]
exonlty[as.character(merged.data$featureType[rt]) ==
"splice_site"] <- final.color.list[["KNOWN.SPLICE.CONNECTION.LTY"]]
is.sig.feature <- f.na(each <= FDR)
sig.feature <- which(is.sig.feature)
if (debug.mode && verbose)
message("> pJSRforGene(): ", "Reached step 5.")
sub <- data.frame(start = merged.data$start[rt], end = merged.data$end[rt],
chr = merged.data$chr[rt], strand = merged.data$strand[rt],
is.exon = (merged.data$featureType[rt] == "exonic_part"),
is.testable = merged.data$testable[rt], is.sig = is.sig.feature,
col = vertline.col, featureType = merged.data$featureType[rt],
featureID = rownames(merged.data)[rt], stringsAsFactors = FALSE)
sub.allExon <- data.frame(start = flat.gff.data$start[rt.allExon],
end = flat.gff.data$end[rt.allExon], chr = flat.gff.data$chrom[rt.allExon],
strand = flat.gff.data$strand[rt.allExon], is.exon = (flat.gff.data$featureType[rt.allExon] ==
"exonic_part"), featureID = as.character(flat.gff.data$featureName[rt.allExon]),
stringsAsFactors = FALSE)
sub.allJunction <- data.frame(start = flat.gff.data$start[rt.allJunction],
end = flat.gff.data$end[rt.allJunction], chr = flat.gff.data$chrom[rt.allJunction],
strand = flat.gff.data$strand[rt.allJunction], is.exon = (flat.gff.data$featureType[rt.allJunction] ==
"exonic_part"), feature.type = flat.gff.data$featureType[rt.allJunction],
is.novel = (flat.gff.data$featureType[rt.allJunction] ==
"novel_splice_site"), featureID = as.character(flat.gff.data$featureName[rt.allJunction]),
stringsAsFactors = FALSE)
sub.allJunction$is.plotted <- sub.allJunction$featureID %in%
sub$featureID
testable.featureIDs <- sub$featureID[sub$is.testable]
sig.featureIDs <- sub$featureID[sub$is.sig]
untestable.featureIDs <- rownames(merged.data)[untestable.rt]
if (debug.mode && verbose) {
message("testable.featureIDs: ", paste0(testable.featureIDs,
collapse = ","))
}
if (debug.mode && verbose) {
message("sig.featureIDs: ", paste0(sig.featureIDs,
collapse = ","))
}
if (debug.mode && verbose) {
message("untestable.featureIDs: ", paste0(untestable.featureIDs,
collapse = ","))
}
sub.allExon$is.testable <- sub.allExon$featureID %in%
testable.featureIDs
sub.allExon$is.sig <- sub.allExon$featureID %in% sig.featureIDs
sub.allExon$is.untestable <- sub.allExon$featureID %in%
untestable.featureIDs
sub.allJunction$is.testable <- sub.allJunction$featureID %in%
testable.featureIDs
sub.allJunction$is.sig <- sub.allJunction$featureID %in%
sig.featureIDs
sub.allJunction$is.untestable <- sub.allJunction$featureID %in%
untestable.featureIDs
sub.allExon$lineColor <- ifelse(sub.allExon$is.testable,
ifelse(sub.allExon$is.sig, SIG.FEATURE.COLOR, NOSIG.FEATURE.COLOR),
ifelse(sub.allExon$is.untestable, UNTESTABLE.FEATURE.COLOR,
EXCLUDED.FEATURE.COLOR))
sub.allExon$fillColor <- ifelse(sub.allExon$is.testable,
ifelse(sub.allExon$is.sig, SIG.FEATURE.FILL.COLOR,
NOSIG.FEATURE.FILL.COLOR), ifelse(sub.allExon$is.untestable,
UNTESTABLE.FEATURE.FILL.COLOR, EXCLUDED.FEATURE.FILL.COLOR))
sub.allExon$borderColor <- ifelse(sub.allExon$is.testable,
ifelse(sub.allExon$is.sig, SIG.FEATURE.BORDER.COLOR,
NOSIG.FEATURE.COLOR), ifelse(sub.allExon$is.untestable,
UNTESTABLE.FEATURE.BORDER.COLOR, EXCLUDED.FEATURE.BORDER.COLOR))
sub.allJunction$lineColor <- ifelse(sub.allJunction$is.testable,
ifelse(sub.allJunction$is.sig, SIG.FEATURE.COLOR,
NOSIG.FEATURE.COLOR), ifelse(sub.allJunction$is.untestable,
UNTESTABLE.FEATURE.COLOR, EXCLUDED.FEATURE.COLOR))
sub.allJunction$lty <- ifelse(sub.allJunction$is.novel,
final.color.list[["NOVEL.SPLICE.LTY"]], final.color.list[["KNOWN.SPLICE.LTY"]])
sig.feature.names <- sub$featureID[is.sig.feature & sub$is.exon]
allExon.isSig <- sub.allExon$featureID %in% sig.featureIDs
allExon.exonCol <- ifelse(allExon.isSig, "#F219ED", "#CCCCCC")
if (debug.mode == 2 && verbose) {
message("> Debugging Info:")
message(" Exons:")
message(" ", paste0(names(sub.allExon),
collapse = "\t"))
for (i in 1:length(sub.allExon$start)) {
message(" ", paste0(sub.allExon[i,
], collapse = "\t"))
}
message("> dim(count) = ", paste0(dim(count),
collapse = ","))
message("> rt: ", paste0(rt, collapse = ","))
message("> rango: ", paste0(rango, collapse = ","))
message("> ncol(count): ", ncol(count))
}
sub.sig <- sub[sig.feature, , drop = FALSE]
rel.calc.min <- min(sub.allJunction$start, sub.allExon$start)
rel.calc.max <- max(sub.allJunction$end, sub.allExon$end)
transcripts <- sapply(sapply(flat.gff.data$transcripts[rt.allExon],
toString), function(x) {
strsplit(x, "+", fixed = TRUE)
})
trans <- Reduce(union, transcripts)
if (displayTranscripts == TRUE) {
mat <- 1:4
hei <- c(10, CONNECTIONS.height, GENE.annotation.height +
SPLICE.annotation.height, TX.annotation.height *
(length(trans) + TX.margins[1] + TX.margins[2]))
}
else {
mat <- 1:4
hei <- c(10, CONNECTIONS.height, GENE.annotation.height +
SPLICE.annotation.height, TX.annotation.height *
TX.margins[2])
}
layout(matrix(mat), heights = hei)
if (debug.mode && verbose) {
message("> length(trans) = ", length(trans))
message("> FINAL LAYOUT:")
message("> heights = [", paste0(hei, collapse = ","),
"]")
}
if (debug.mode && verbose)
message("> pJSRforGene(): ", "Reached step 6.")
ylimn <- c(min(min(count, na.rm = TRUE), 0), max(count,
na.rm = TRUE))
if ((!use.vst) && use.log)
ylimn[1] <- INTERNAL.NINF.VALUE
p.values.labels <- ifelse(f.na(each <= FDR), format(each,
digits = 3), "")
if (any(sub$is.exon)) {
italicize.label <- !sub$is.exon
}
else {
italicize.label <- NULL
}
plotWindowXmax <- if (plot.gene.level.expression) {
(length(intervals) + 1)/length(intervals)
}
else {
1
}
intervals <- drawPlot(matr = count, ylimn, jscs, intervals,
rango, textAxis = y.axis.title, geneLevelAxisTitle = y.axis.title.right,
rt = rt, color.count = color.count, colorlines = vertline.col,
countbinIDs = merged.data$countbinID[rt], use.vst = use.vst,
use.log = use.log, plot.type = plot.type, main.title = main.title,
draw.legend = draw.legend, color.key = color, condition.names = condition.names,
p.values = p.values.labels, draw.p.values = label.p.vals,
plot.lwd = plot.lwd, axes.lwd = axes.lwd, anno.lwd = anno.lwd,
par.cex = par.cex, anno.cex.text = anno.cex.text,
anno.cex.axis = anno.cex.axis, anno.cex.main = anno.cex.main,
fit.countbin.names = fit.countbin.names, debug.mode = debug.mode,
plot.gene.level.expression = plot.gene.level.expression,
geneCount = geneCount, color.geneCount = color.geneCount,
yAxisLabels.inExponentialForm = yAxisLabels.inExponentialForm,
italicize.label = italicize.label, condition.legend.text = condition.legend.text,
annolink.col = annolink.col, exonlty = exonlty, graph.margins = graph.margins,
plotWindowXmax = plotWindowXmax, fit.labels = fit.labels,
...)
if (debug.mode && verbose)
message("> pJSRforGene(): ", "Reached end of step 6.")
if (debug.mode && verbose)
message("> pJSRforGene(): ", "Reached step 7.")
if (USE.MARGIN.MEX) {
par(mar = c(0, graph.margins[2], 0, graph.margins[4]),
cex = par.cex, mex = anno.cex.text)
}
else {
par(mar = c(0, graph.margins[2], 0, graph.margins[4]),
cex = par.cex)
}
plot.new()
plot.window(xlim = c(0, plotWindowXmax), ylim = c(0,
1), xaxs = "i")
rel <- data.frame(start = rescale.coords(sub$start, rescale.iv),
end = rescale.coords(sub$end, rescale.iv))
connection.lines.bottom <- apply((rbind(rel[rango, 2],
rel[rango, 1])), 2, median) * plotWindowXmax
connection.lines.top <- apply(rbind(intervals[rango],
intervals[rango + 1] - ((intervals[rango + 1] - intervals[rango]) *
0.2)), 2, median)
segments(connection.lines.bottom, 0, connection.lines.top,
1, col = annolink.col, lty = exonlty, lwd = gene.lwd,
cex = anno.cex.text, cex.axis = anno.cex.main, cex.main = anno.cex.main,
xpd = NA, ...)
par(mar = c(1.5, graph.margins[2], 0, graph.margins[4]),
cex = par.cex)
startSites <- c()
endSites <- c()
for (i in 1:length(trans)) {
logicexons <- sapply(transcripts, function(x) {
any(x == trans[i])
})
startSites <- c(startSites, sub.allExon$start[logicexons][1])
endSites <- c(endSites, sub.allExon$end[logicexons][sum(logicexons)])
}
startSites <- unique(startSites)
endSites <- unique(endSites)
drawGene(rel.calc.min, rel.calc.max, tr = sub, tr.allExon = sub.allExon,
tr.allJunction = sub.allJunction, rango, rescale.iv = rescale.iv,
exoncol = annolink.col, allExon.exonCol = allExon.exonCol,
names, trName = "Gene model", anno.cex.text = anno.cex.text,
par.cex = par.cex, exonlty = exonlty, plot.lwd = gene.lwd,
anno.lwd = anno.lwd, show.strand.arrows = show.strand.arrows,
geneStrand = geneStrand, cex.axis = anno.cex.axis,
cex.main = anno.cex.main, draw.untestable.annotation = draw.untestable.annotation,
draw.start.end.sites = draw.start.end.sites, startSites = startSites,
endSites = endSites, cex.arrows = cex.arrows, chrom.label = chrom.label,
label.chromosome = label.chromosome, splice.junction.drawing.style = splice.junction.drawing.style,
draw.nested.SJ = draw.nested.SJ, merge.exon.parts = merge.exon.parts,
plot.untestable.results = plot.untestable.results,
exon.height = GENE.annotation.height/(SPLICE.annotation.height +
GENE.annotation.height), INTERNAL.VARS = INTERNAL.VARS,
flip.splicing = flip.splicing, ...)
include.endpoints.on.coordinates <- FALSE
num.coord.miniticks.per.tick <- 10
if (drawCoordinates) {
if (!is.null(rescale.iv)) {
pretty.x <- pretty(c(rel.calc.min, rel.calc.max),
n = 5)
pretty.interval <- pretty.x[2] - pretty.x[1]
pretty.x <- pretty.x[pretty.x > rel.calc.min &
pretty.x < rel.calc.max]
rescaled.pretty.x <- rescale.coords(pretty.x,
rescale.iv)
if (num.coord.miniticks.per.tick > 0) {
rel.coord.miniticks <- pretty.interval * (1:(num.coord.miniticks.per.tick -
1))/num.coord.miniticks.per.tick
unscaled.coord.miniticks <- unlist(lapply(pretty.x,
function(a) {
a + rel.coord.miniticks
}))
unscaled.coord.miniticks <- c(pretty.x[1] -
rel.coord.miniticks, unscaled.coord.miniticks)
rescaled.coord.miniticks <- rescale.coords(unscaled.coord.miniticks,
rescale.iv) * (rel.calc.max - rel.calc.min) +
rel.calc.min
}
else {
rescaled.coord.miniticks <- FALSE
}
if (include.endpoints.on.coordinates) {
if (min(rescaled.pretty.x) > 0.05) {
pretty.x <- c(rel.calc.min, pretty.x)
rescaled.pretty.x <- c(0, rescaled.pretty.x)
}
if (max(rescaled.pretty.x) < 0.95) {
pretty.x <- c(pretty.x, rel.calc.max)
rescaled.pretty.x <- c(rescaled.pretty.x,
1)
}
}
rescaled.pretty.x <- rescaled.pretty.x * (rel.calc.max -
rel.calc.min) + rel.calc.min
}
else {
pretty.x <- pretty(c(rel.calc.min, rel.calc.max),
n = 5)
coord.miniticks <- FALSE
}
usr <- par("usr")
cxy <- par("cxy")
pretty.x <- sprintf("%0.f", pretty.x)
smallest.width.coordAxis <- min(abs(rescaled.pretty.x[-1] -
rescaled.pretty.x[-length(rescaled.pretty.x)]))
if (fit.genomic.axis) {
anno.cex.coordAxis <- shrink.character.vector(paste0(pretty.x,
"0"), curr.cex = anno.cex.axis, max.width = smallest.width.coordAxis)
if (anno.cex.coordAxis * 2 < anno.cex.axis) {
anno.cex.coordAxis <- anno.cex.axis/2
coordAxis.widths <- abs(rescaled.pretty.x[-1] -
rescaled.pretty.x[-length(rescaled.pretty.x)])
for (i in 2:length(pretty.x)) {
if ((strwidth(pretty.x[i - 1], cex = anno.cex.coordAxis)/2) +
(strwidth(pretty.x[i], cex = anno.cex.coordAxis)/2) >
abs(rescaled.pretty.x[i - 1] - rescaled.pretty.x[i])) {
pretty.x[i] <- ""
}
}
}
}
else {
anno.cex.coordAxis <- anno.cex.axis
}
devlim <- device.limits()
coord.ticks.top <- usr[3]
coord.mainTicks.bottom <- usr[3] - (cxy[2]/2)
coord.text.top <- usr[3] - (cxy[2] * (3/4))
coord.miniTicks.bottom <- coord.ticks.top - abs(coord.mainTicks.bottom -
coord.ticks.top)/2
segments(x0 = rescaled.pretty.x, y0 = usr[3], x1 = rescaled.pretty.x,
y1 = usr[3] - (cxy[2]/2), xpd = NA, lwd = anno.lwd,
...)
lines(c(rel.calc.min, rel.calc.max), c(par("usr")[3],
par("usr")[3]), lwd = axes.lwd, xpd = NA, ...)
segments(x0 = rescaled.coord.miniticks, y0 = usr[3],
x1 = rescaled.coord.miniticks, y1 = usr[3] -
(cxy[2]/4), xpd = NA, lwd = anno.lwd, ...)
text(rescaled.pretty.x, usr[3] - (cxy[2] * (3/4)),
pretty.x, cex = anno.cex.coordAxis, xpd = NA,
adj = c(0.5, 1), ...)
}
if (label.chromosome) {
chrom.label.width.max <- abs(par("usr")[1] - device.limits()[1]) *
0.9
chrom.label.cex <- shrink.character.vector(chrom.label,
curr.cex = anno.cex.text, max.width = chrom.label.width.max)
text(par("usr")[1], par("usr")[3], chrom.label, cex = chrom.label.cex,
adj = c(1.1, 0.5), xpd = NA, font = 2, ...)
}
if (displayTranscripts) {
if (USE.MARGIN.MEX) {
par(cex = par.cex, mar = c(0, graph.margins[2],
0, graph.margins[4]), mex = anno.cex.text)
}
else {
par(cex = par.cex, mar = c(0, graph.margins[2],
0, graph.margins[4]))
}
plot.new()
plot.window(xlim = c(rel.calc.min, rel.calc.max),
ylim = c(-TX.margins[2], length(trans) + TX.margins[1]),
xaxs = "i", yaxs = "i")
for (i in 1:length(trans)) {
ymin <- length(trans) - i
if (include.TX.names) {
trName = trans[i]
}
else {
trName = NULL
}
logicexons <- sapply(transcripts, function(x) {
any(x == trans[i])
})
tr <- sub.allExon[logicexons, ]
curr.exoncol <- ifelse(allExon.isSig[logicexons],
"#F219ED", "#CCCCCC")
drawTranscript(rel.calc.min, rel.calc.max, ymin = ymin,
tr = tr, tr.allJunction = sub.allJunction,
rango = 1:nrow(tr), rescale.iv = rescale.iv,
names = c(), trName = trName, trStrand = txStrandMap[[trName]],
draw.strand = geneStrand == ".", par.cex = par.cex,
anno.cex.text = anno.cex.text, sub.sig = sub.sig,
anno.lwd = gene.lwd, cex.axis = anno.cex.axis,
cex.main = anno.cex.main, cex.arrows = cex.arrows,
...)
}
}
else {
par(mar = c(0, 0, 0, 0))
plot.new()
}
par(mar = c(0, 0, 0, 0))
if (debug.mode && verbose)
message("> pJSRfG(): ", " Done.")
}
}, error = function(e) {
message("Error caught while attempting plotJunctionSeqResultsForGene")
message("---------------------")
message(" Error text:")
message(" ", e)
message("")
message("---------------------")
message(" Input parameters:")
message(" geneID = ", geneID)
message(" plot.type = ", plot.type)
message(" displayTranscripts = ", displayTranscripts)
message("---------------------")
warning("Error caught while attempting plotJunctionSeqResultsForGene")
warning("Error text:")
warning(e)
})
where 10: plotJunctionSeqResultsForGene(geneID, jscs, colorRed.FDR.threshold = FDR,
colorList = colorList, plot.type = plot.type, use.vst = use.vst,
use.log = use.log, exon.rescale.factor = exon.rescale.factor,
displayTranscripts = TRUE, plot.lwd = plot.lwd, axes.lwd = axes.lwd,
anno.lwd = anno.lwd, par.cex = par.cex, anno.cex.text = anno.cex.text,
plot.exon.results = plot.exon.results, plot.junction.results = plot.junction.results,
plot.novel.junction.results = plot.novel.junction.results,
plot.untestable.results = plot.untestable.results, anno.cex.axis = anno.cex.axis,
anno.cex.main = anno.cex.main, drawCoordinates = drawCoordinates,
graph.margins = graph.margins, yAxisLabels.inExponentialForm = yAxisLabels.inExponentialForm,
plot.gene.level.expression = plot.gene.level.expression,
condition.legend.text = condition.legend.text, include.TX.names = include.TX.names,
GENE.annotation.relative.height = GENE.annotation.relative.height,
TX.annotation.relative.height = TX.annotation.relative.height,
draw.start.end.sites = draw.start.end.sites, show.strand.arrows = show.strand.arrows,
debug.mode = debug.mode, sequencing.type = sequencing.type,
gene.lwd = gene.lwd, CONNECTIONS.relative.height = CONNECTIONS.relative.height,
TX.margins = TX.margins, SPLICE.annotation.relative.height = SPLICE.annotation.relative.height,
draw.nested.SJ = draw.nested.SJ, INTERNAL.VARS = INTERNAL.VARS,
...)
where 11: buildAllPlotsForGene(geneID = geneID, jscs = jscs, outfile.prefix = outfile.prefixes,
use.plotting.device = use.plotting.device, use.vst = use.vst,
use.log = use.log, exon.rescale.factor = exon.rescale.factor,
plot.gene.level.expression = plot.gene.level.expression,
with.TX = with.TX, without.TX = without.TX, expr.plot = expr.plot,
normCounts.plot = normCounts.plot, rExpr.plot = rExpr.plot,
rawCounts.plot = rawCounts.plot, colorRed.FDR.threshold = colorRed.FDR.threshold,
colorList = colorList, plot.exon.results = plot.exon.results,
plot.junction.results = plot.junction.results, plot.novel.junction.results = plot.novel.junction.results,
plot.untestable.results = plot.untestable.results, plot.lwd = plot.lwd,
axes.lwd = axes.lwd, anno.lwd = anno.lwd, gene.lwd = gene.lwd,
drawCoordinates = drawCoordinates, show.strand.arrows = show.strand.arrows,
graph.margins = graph.margins, par.cex = par.cex, anno.cex.text = anno.cex.text,
anno.cex.axis = anno.cex.axis, anno.cex.main = anno.cex.main,
base.plot.height = base.plot.height, base.plot.width = base.plot.width,
base.plot.units = base.plot.units, GENE.annotation.relative.height = GENE.annotation.relative.height,
TX.annotation.relative.height = TX.annotation.relative.height,
CONNECTIONS.relative.height = CONNECTIONS.relative.height,
TX.margins = TX.margins, SPLICE.annotation.relative.height = SPLICE.annotation.relative.height,
autoscale.height.to.fit.TX.annotation = autoscale.height.to.fit.TX.annotation,
autoscale.width.to.fit.bins = autoscale.width.to.fit.bins,
name.files.with.geneID = name.files.with.geneID, plotting.device.params = plotting.device.params,
yAxisLabels.inExponentialForm = yAxisLabels.inExponentialForm,
condition.legend.text = condition.legend.text, include.TX.names = include.TX.names,
draw.start.end.sites = draw.start.end.sites, verbose = verbose,
debug.mode = debug.mode, sequencing.type = sequencing.type,
minimalImageFilenames = minimalImageFilenames, INTERNAL.VARS = INTERNAL.VARS,
...)
where 12: buildAllPlots(jscs)
--- value of length: 3 type: logical ---
[1] FALSE TRUE FALSE
--- function from context ---
function (expr, name, parentenv, handler)
{
.Internal(.addCondHands(name, list(handler), parentenv, environment(),
FALSE))
expr
}
<bytecode: 0x01ee70e0>
<environment: 0x2100cc60>
--- function search by body ---
----------- END OF FAILURE REPORT --------------
Fatal error: length > 1 in coercion to logical
** running examples for arch 'x64' ... ERROR
Running examples in 'JunctionSeq-Ex.R' failed
The error most likely occurred in:
> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: buildAllPlots
> ### Title: Create and save a full battery of JunctionSeq expression plots.
> ### Aliases: buildAllPlots
>
> ### ** Examples
>
>
> data(exampleDataSet,package="JctSeqData");
> buildAllPlots(jscs);
> buildAllPlots: Found 5 genes with at least one significant exon, at adjusted-p-value threshold 0.01
> buildAllPlots: Starting plotting...
Warning in dir.create(outfile.prefix) : '.' already exists
> buildAllPlots: Generating Dispersion Plot
Starting device: png (./dispersion-plot.png)
abundance ranges from 0.146009972023612 to 5452.18841073905
dispersion ranges from 1.94121697570967e-06 to 2.08746163174998
Plotting dispersions from 2.71367228072472e-05 to 2.08746163174998
> buildAllPlots: Generating MA-Plot (log2FC(CTRL/CASE))
Starting device: png (./ma-plot-log2FC(CTRLvsCASE).png)
> buildAllPlots: Writing HTML results index.
Copying default css stylesheet.
Writing html index. Sat Oct 17 05:31:11 2020
Writing formula data. Sat Oct 17 05:31:11 2020
Writing methods data. Sat Oct 17 05:31:11 2020
Writing sample data. Sat Oct 17 05:31:11 2020
Writing dispersion data. Sat Oct 17 05:31:11 2020
Writing summary plots. Sat Oct 17 05:31:11 2020
Compiling data table. Sat Oct 17 05:31:11 2020
Writing data table. Sat Oct 17 05:31:11 2020
Html index complete. Sat Oct 17 05:31:11 2020
Writing pages. Sat Oct 17 05:31:11 2020
Finished all html files. Sat Oct 17 05:31:11 2020
> buildAllPlots: Finished writing HTML results index.
> buildAllPlots: starting geneID: ENSRNOG00000008079 (1 of 5)
starting buildAllPlotsForGene() for geneID: ENSRNOG00000008079
Starting nested heights...
Starting device: png (.//exprTX/1-expr-TX.png)
> pJSRfG(): ENSRNOG00000008079, plot.type: expr (with TX)
Starting device: png (.//expr/1-expr.png)
> pJSRfG(): ENSRNOG00000008079, plot.type: expr
Starting device: png (.//normCountsTX/1-normCts-TX.png)
> pJSRfG(): ENSRNOG00000008079, plot.type: normCounts (with TX)
Starting device: png (.//normCounts/1-normCts.png)
> pJSRfG(): ENSRNOG00000008079, plot.type: normCounts
Starting device: png (.//rExprTX/1-rExpr-TX.png)
> pJSRfG(): ENSRNOG00000008079, plot.type: rExpr (with TX)
Starting device: png (.//rExpr/1-rExpr.png)
> pJSRfG(): ENSRNOG00000008079, plot.type: rExpr
> buildAllPlots: starting geneID: ENSRNOG00000009281 (2 of 5)
starting buildAllPlotsForGene() for geneID: ENSRNOG00000009281
Starting nested heights...
Starting device: png (.//exprTX/2-expr-TX.png)
> pJSRfG(): ENSRNOG00000009281, plot.type: expr (with TX)
----------- FAILURE REPORT --------------
--- failure: length > 1 in coercion to logical ---
--- srcref ---
:
--- package (from environment) ---
JunctionSeq
--- call from context ---
doTryCatch(return(expr), name, parentenv, handler)
--- call from argument ---
any(c(pw[1] - 1, pw[1] + 1) %in% pwo) && pw%%2 == 0
--- R stacktrace ---
where 1: doTryCatch(return(expr), name, parentenv, handler)
where 2: tryCatchOne(expr, names, parentenv, handlers[[1L]])
where 3: tryCatchList(expr, classes, parentenv, handlers)
where 4: tryCatch({
if (debug.mode)
message("> Step 6.9b")
pval.drop <- rep(-pval.height, length(p.values))
pw <- which(p.values != "")
pwo <- pw
while (length(pw) > 0) {
if (any(c(pw[1] - 1, pw[1] + 1) %in% pwo) && pw%%2 ==
0) {
pval.drop[pw[1]] <- pval.drop[pw[1]] + pval.height
}
else {
}
pw <- pw[-1]
}
pval.y <- pval.ceiling + pval.drop
}, error = function(e) {
message("WARNING: staggering p-values failed. Falling back to simpler method. This warning should never appear. If you see this, reporting this to the developer would be appreciated.")
})
where 5: drawPlot(matr = count, ylimn, jscs, intervals, rango, textAxis = y.axis.title,
geneLevelAxisTitle = y.axis.title.right, rt = rt, color.count = color.count,
colorlines = vertline.col, countbinIDs = merged.data$countbinID[rt],
use.vst = use.vst, use.log = use.log, plot.type = plot.type,
main.title = main.title, draw.legend = draw.legend, color.key = color,
condition.names = condition.names, p.values = p.values.labels,
draw.p.values = label.p.vals, plot.lwd = plot.lwd, axes.lwd = axes.lwd,
anno.lwd = anno.lwd, par.cex = par.cex, anno.cex.text = anno.cex.text,
anno.cex.axis = anno.cex.axis, anno.cex.main = anno.cex.main,
fit.countbin.names = fit.countbin.names, debug.mode = debug.mode,
plot.gene.level.expression = plot.gene.level.expression,
geneCount = geneCount, color.geneCount = color.geneCount,
yAxisLabels.inExponentialForm = yAxisLabels.inExponentialForm,
italicize.label = italicize.label, condition.legend.text = condition.legend.text,
annolink.col = annolink.col, exonlty = exonlty, graph.margins = graph.margins,
plotWindowXmax = plotWindowXmax, fit.labels = fit.labels,
...)
where 6: doTryCatch(return(expr), name, parentenv, handler)
where 7: tryCatchOne(expr, names, parentenv, handlers[[1L]])
where 8: tryCatchList(expr, classes, parentenv, handlers)
where 9: tryCatch({
GENE.annotation.height <- GENE.annotation.relative.height *
10
TX.annotation.height <- TX.annotation.relative.height * 10
CONNECTIONS.height <- CONNECTIONS.relative.height * 10
SPLICE.annotation.height <- SPLICE.annotation.relative.height *
10
if (!any(geneID == fData(jscs)$geneID)) {
stop(paste0("FATAL ERROR: GeneID \"", geneID, "\" not found in dataset!"))
}
flat.gff.data <- jscs@flatGffData
merged.data <- fData(jscs)
condition <- jscs@phenoData$condition
exonRescaleFunction <- match.arg(exonRescaleFunction)
sequencing.type <- match.arg(sequencing.type)
splice.junction.drawing.style <- match.arg(splice.junction.drawing.style)
plot.type <- match.arg(plot.type)
truncateBelowOne <- TRUE
if (is.null(plot.exon.results)) {
plot.exon.results <- any(merged.data$featureType == "exonic_part")
}
if (is.null(plot.junction.results)) {
plot.junction.results <- any(merged.data$featureType ==
"splice_site" | merged.data$featureType == "novel_splice_site")
}
if (is.null(plot.novel.junction.results)) {
if (plot.junction.results) {
plot.novel.junction.results <- any(merged.data$featureType ==
"novel_splice_site")
}
else {
plot.novel.junction.results <- FALSE
}
}
flip.splicing <- if (plot.junction.results) {
FALSE
}
else {
TRUE
}
geneName <- jscs@flatGffGeneData$gene_name[jscs@flatGffGeneData$geneID ==
geneID]
if (is.null(condition.legend.text)) {
condition.legend.text <- levels(jscs@phenoData$condition)
names(condition.legend.text) <- condition.legend.text
}
else {
if (is.null(names(condition.legend.text))) {
warning("names(condition.legend.text is NULL! condition.legend.text mis-formatted. Must be a list or character vector, with element names equal to the levels of pData(jscs)$condition. Falling back.")
condition.legend.text <- levels(jscs@phenoData$condition)
names(condition.legend.text) <- condition.legend.text
}
else if (!all(levels(jscs@phenoData$condition) %in% names(condition.legend.text))) {
warning("Not all levels contained in names(condition.legend.text)! condition.legend.text mis-formatted. Must be a list or character vector, with element names equal to the levels of pData(jscs)$condition. Falling back.")
condition.legend.text <- levels(jscs@phenoData$condition)
names(condition.legend.text) <- condition.legend.text
}
}
if (!is.null(jscs@flatGffGeneData[["aggregateGeneStrand"]])) {
geneStrand <- as.character(jscs@flatGffGeneData[["aggregateGeneStrand"]][jscs@flatGffGeneData[["geneID"]] ==
geneID])
txSetString <- as.character(jscs@flatGffGeneData[["tx_set"]][jscs@flatGffGeneData[["geneID"]] ==
geneID])
txStrandString <- as.character(jscs@flatGffGeneData[["tx_strands"]][jscs@flatGffGeneData[["geneID"]] ==
geneID])
txSet <- strsplit(as.character(txSetString), "+", fixed = TRUE)[[1]]
txStrand <- strsplit(as.character(txStrandString), ",",
fixed = TRUE)[[1]]
txStrandMap <- as.list(txStrand)
names(txStrandMap) <- txSet
}
else {
txStrandMap <- list()
geneStrand <- "."
}
final.color.list <- overmerge.list(JUNCTIONSEQ.DEFAULT.COLOR.LIST,
colorList)
SIG.VERTLINE.COLOR = final.color.list[["SIG.VERTLINE.COLOR"]]
NOSIG.VERTLINE.COLOR = final.color.list[["NOSIG.VERTLINE.COLOR"]]
UNTESTABLE.VERTLINE.COLOR = final.color.list[["UNTESTABLE.VERTLINE.COLOR"]]
SIG.FEATURE.COLOR = final.color.list[["SIG.FEATURE.COLOR"]]
NOSIG.FEATURE.COLOR = final.color.list[["NOSIG.FEATURE.COLOR"]]
UNTESTABLE.FEATURE.COLOR = final.color.list[["UNTESTABLE.FEATURE.COLOR"]]
EXCLUDED.FEATURE.COLOR = final.color.list[["EXCLUDED.FEATURE.COLOR"]]
SIG.FEATURE.BORDER.COLOR = final.color.list[["SIG.FEATURE.BORDER.COLOR"]]
NOSIG.FEATURE.BORDER.COLOR = final.color.list[["NOSIG.FEATURE.BORDER.COLOR"]]
UNTESTABLE.FEATURE.BORDER.COLOR = final.color.list[["UNTESTABLE.FEATURE.BORDER.COLOR"]]
EXCLUDED.FEATURE.BORDER.COLOR = final.color.list[["EXCLUDED.FEATURE.BORDER.COLOR"]]
SIG.FEATURE.FILL.COLOR = final.color.list[["SIG.FEATURE.FILL.COLOR"]]
NOSIG.FEATURE.FILL.COLOR = final.color.list[["NOSIG.FEATURE.FILL.COLOR"]]
UNTESTABLE.FEATURE.FILL.COLOR = final.color.list[["UNTESTABLE.FEATURE.FILL.COLOR"]]
EXCLUDED.FEATURE.FILL.COLOR = final.color.list[["EXCLUDED.FEATURE.FILL.COLOR"]]
PLOTTING.LINE.COLORS = final.color.list[["PLOTTING.LINE.COLORS"]]
gene.level.buffer <- 0.5
if (is.null(plot.gene.level.expression)) {
if (use.vst) {
plot.gene.level.expression <- FALSE
}
else {
plot.gene.level.expression <- TRUE
}
}
if (plot.gene.level.expression && use.vst) {
warning("WARNING: plotting of gene-level expression is not supported for vst-transformed plots. Errors are likely to follow.")
}
FDR <- colorRed.FDR.threshold
if (verbose) {
displayTXstring <- if (displayTranscripts)
" (with TX)"
else ""
message("> pJSRfG(): ", geneID, ", plot.type: ", plot.type,
displayTXstring)
}
chrom.label <- as.character(merged.data$chr[merged.data$geneID ==
geneID][1])
rt <- merged.data$geneID == geneID
if (!plot.exon.results) {
if (debug.mode)
message("> Removing ", sum(rt & merged.data$featureType ==
"exonic_part"), " exonic_part features. ", sum(rt &
merged.data$featureType != "exonic_part"), " features remaining")
rt <- rt & merged.data$featureType != "exonic_part"
}
if (!plot.junction.results) {
if (debug.mode)
message("> Removing ", sum(rt & merged.data$featureType ==
"splice_site"), " splice_site features. ", sum(rt &
merged.data$featureType != "splice_site"), " features remaining")
rt <- rt & merged.data$featureType != "splice_site"
}
if (!plot.novel.junction.results) {
if (debug.mode)
message("> Removing ", sum(rt & merged.data$featureType ==
"novel_splice_site"), " novel_splice_site features. ",
sum(rt & merged.data$featureType != "novel_splice_site"),
" features remaining")
rt <- rt & merged.data$featureType != "novel_splice_site"
}
untestable.rt <- which(rt & (!merged.data$testable))
if (!plot.untestable.results) {
if (debug.mode)
message("> Removing ", sum(rt & !merged.data$testable),
" untestable features. ", sum(rt & merged.data$testable),
" features remaining")
rt <- rt & merged.data$testable
}
if (debug.mode)
message("> Plotting ", sum(rt), " features.")
rt <- which(rt)
if (length(rt) == 0) {
message("NO FEATURES TO PLOT!")
}
else {
rt.allExon <- which(flat.gff.data$gene_id == geneID &
flat.gff.data$featureType == "exonic_part")
rango.allExon <- 1:length(rt.allExon)
rt.allJunction <- which(flat.gff.data$gene_id == geneID &
(flat.gff.data$featureType == "splice_site" | flat.gff.data$featureType ==
"novel_splice_site"))
rango.allJunction <- 1:length(rt.allJunction)
rescale.iv <- generate.interval.scale(data.frame(start = c(flat.gff.data$start[rt.allExon],
flat.gff.data$start[rt.allJunction]), end = c(flat.gff.data$end[rt.allExon],
flat.gff.data$end[rt.allJunction]), is.exon = c(rep(TRUE,
length(rt.allExon)), rep(FALSE, length(rt.allJunction)))),
exon.rescale.factor, exonRescaleFunction, debug.mode = debug.mode)
rel <- data.frame(start = rescale.coords(merged.data$start[rt],
rescale.iv), end = rescale.coords(merged.data$end[rt],
rescale.iv))
if (sort.features) {
rt <- rt[order(((rel$end - rel$start)/2) + rel$start)]
}
rango <- 1:length(rt)
draw.legend <- TRUE
condition.names <- levels(condition)
sample.names <- sampleNames(jscs@phenoData)
colorcode.title <- TRUE
numcond <- length(condition.names)
if (debug.mode && verbose)
message("> pJSRforGene(): ", "Reached step 2.")
default.color.list <- PLOTTING.LINE.COLORS
if (numcond > length(default.color.list)) {
message("Too many condition values, the default color selection may not look good! Set your own colors by setting the \"color\" parameter.")
color <- rgb(colorRamp(c("#D7191C", "#FFFFBF", "#2B83BA"))(seq(0,
1, length.out = numcond)), maxColorValue = 255,
alpha = 175)
}
else {
color <- color2transparentVector(default.color.list[1:numcond],
t = 175)
}
if (debug.mode && verbose)
message("> pJSRforGene(): ", "color = ", paste0(color,
collapse = ","))
if (debug.mode && verbose)
message("> pJSRforGene(): ", "length(color) = ",
length(color))
if (debug.mode && verbose)
message("> pJSRforGene(): ", "length(condition.names) = ",
length(condition.names))
if (debug.mode && verbose)
message("> pJSRforGene(): ", "condition.names = ",
paste0(condition.names, collapse = ","))
names(color) <- condition.names
y.axis.title <- ""
main.title <- ""
if (debug.mode && verbose)
message("> pJSRforGene(): ", "Reached step 3.")
if (truncateBelowOne) {
convertY <- function(y) {
ifelse(y < 0, ifelse(is.infinite(y), INTERNAL.NINF.VALUE,
(1 - exp(y)) * INTERNAL.NINF.VALUE), y)
}
}
else {
convertY <- function(y) {
y
}
}
if (plot.type == "rExpr") {
count <- if (use.vst) {
vst(jscs@plottingEstimates[["relExprEstimate"]][rt,
, drop = FALSE], jscs)
}
else if (use.log) {
apply(log10(jscs@plottingEstimates[["relExprEstimate"]][rt,
, drop = FALSE]), c(1, 2), FUN = convertY)
}
else {
jscs@plottingEstimates[["relExprEstimate"]][rt,
, drop = FALSE]
}
color.count <- rep(color[condition.names], each = nrow(count))
y.axis.title <- "Relative Coverage"
y.axis.title.right <- "Gene-Level Mean Normalized Counts"
main.title <- paste0("Relative Coverage (", geneName,
")")
if (plot.gene.level.expression) {
geneCount <- if (use.vst) {
plot.gene.level.expression <- FALSE
NULL
}
else if (use.log) {
sapply(log10(jscs@geneLevelPlottingEstimates[["geneLevelEstModeled"]][rownames(jscs@geneLevelPlottingEstimates[["geneLevelEstModeled"]]) ==
geneID, ]), FUN = convertY)
}
else {
jscs@geneLevelPlottingEstimates[["geneLevelEstModeled"]][rownames(jscs@geneLevelPlottingEstimates[["geneLevelEstModeled"]]) ==
geneID, ]
}
color.geneCount <- color[condition.names]
}
}
else if (plot.type == "normCounts") {
count <- if (use.vst) {
vst(jscs@plottingEstimates[["normCounts"]][rt,
1:length(sample.names), drop = FALSE], jscs)
}
else if (use.log) {
apply(log10(jscs@plottingEstimates[["normCounts"]][rt,
1:length(sample.names), drop = FALSE]), c(1,
2), FUN = convertY)
}
else {
jscs@plottingEstimates[["normCounts"]][rt, 1:length(sample.names),
drop = FALSE]
}
if (plot.gene.level.expression) {
geneCount <- if (use.vst) {
plot.gene.level.expression <- FALSE
NULL
}
else if (use.log) {
sapply(log10(jscs@geneCountData[rownames(jscs@geneCountData) ==
geneID, ]/sizeFactors(jscs)), FUN = convertY)
}
else {
jscs@geneCountData[rownames(jscs@geneCountData) ==
geneID, ]/sizeFactors(jscs)
}
color.geneCount <- color[as.character(condition)]
}
color.count <- rep(color[as.character(condition)],
each = nrow(count))
y.axis.title <- "Normalized Counts"
y.axis.title.right <- "Gene-Level Normalized Counts"
main.title <- paste0("Normalized Counts (", geneName,
")")
}
else if (plot.type == "rawCounts") {
count <- if (use.vst) {
vst(jscs@countVectors[rt, 1:length(sample.names),
drop = FALSE], jscs)
}
else if (use.log) {
apply(log10(jscs@countVectors[rt, 1:length(sample.names),
drop = FALSE]), c(1, 2), FUN = function(y) {
max(INTERNAL.NINF.VALUE, y)
})
}
else {
jscs@countVectors[rt, 1:length(sample.names),
drop = FALSE]
}
if (plot.gene.level.expression) {
geneCount <- if (use.vst) {
plot.gene.level.expression <- FALSE
NULL
}
else if (use.log) {
sapply(log10(jscs@geneCountData[rownames(jscs@geneCountData) ==
geneID, ]), FUN = convertY)
}
else {
jscs@geneCountData[rownames(jscs@geneCountData) ==
geneID, ]
}
color.geneCount <- color[as.character(condition)]
}
color.count <- rep(color[as.character(condition)],
each = nrow(count))
y.axis.title <- "Raw Counts"
y.axis.title.right <- "Raw Gene-Level Counts"
main.title <- paste0("Raw Counts (", geneName, ")")
}
else if (plot.type == "expr") {
count <- if (use.vst) {
vst(jscs@plottingEstimates[["exprEstimate"]][rt,
, drop = FALSE], jscs)
}
else if (use.log) {
apply(log10(jscs@plottingEstimates[["exprEstimate"]][rt,
, drop = FALSE]), c(1, 2), FUN = convertY)
}
else {
jscs@plottingEstimates[["exprEstimate"]][rt,
, drop = FALSE]
}
if (plot.gene.level.expression) {
geneCount <- if (use.vst) {
plot.gene.level.expression <- FALSE
NULL
}
else if (use.log) {
sapply(log10(jscs@geneLevelPlottingEstimates[["geneLevelEstModeled"]][rownames(jscs@geneLevelPlottingEstimates[["geneLevelEstModeled"]]) ==
geneID, ]), FUN = convertY)
}
else {
jscs@geneLevelPlottingEstimates[["geneLevelEstModeled"]][rownames(jscs@geneLevelPlottingEstimates[["geneLevelEstModeled"]]) ==
geneID, ]
}
color.geneCount <- color[condition.names]
}
color.count <- rep(color[condition.names], each = nrow(count))
y.axis.title <- "Mean Normalized Counts"
y.axis.title.right <- "Gene-Level Mean Normalized Counts"
main.title <- paste0("Mean Normalized Coverage (",
geneName, ")")
}
else {
stop(paste0("FATAL ERROR: Unknown plot type! plot.type = \"",
plot.type, "\""))
}
fragment.label <- if (sequencing.type == "paired-end")
"Read-Pair"
else "Read"
if (plot.type == "rExpr") {
main.title <- paste0("Relative Coverage (", geneName,
")")
y.axis.title <- paste0("Relative Coverage")
y.axis.title.right <- paste0(fragment.label, "s per Sample, Gene-Level")
}
else if (plot.type == "normCounts") {
main.title <- paste0("Normalized Counts (", geneName,
")")
y.axis.title <- paste0("Normalized Counts")
y.axis.title.right <- paste0("Normalized Counts, Gene-Level")
}
else if (plot.type == "rawCounts") {
main.title <- paste0("Raw Counts (", geneName, ")")
y.axis.title <- paste0("Raw ", fragment.label, " Counts")
y.axis.title.right <- paste0("Raw ", fragment.label,
" Counts, Gene-Level")
}
else if (plot.type == "expr") {
main.title <- paste0("Mean Normalized Coverage (",
geneName, ")")
y.axis.title <- paste0(fragment.label, "s per Sample")
y.axis.title.right <- paste0(fragment.label, "s per Sample, Gene-Level")
}
else {
stop(paste0("FATAL ERROR: Unknown plot type! plot.type = \"",
plot.type, "\""))
}
if (!is.null(title.main))
main.title <- title.main
if (!is.null(title.ylab))
y.axis.title <- title.ylab
if (!is.null(title.ylab.right))
y.axis.title.right <- title.ylab.right
if (!plot.gene.level.expression) {
geneCount <- NULL
color.geneCount <- NULL
}
count <- as.matrix(count)
if (debug.mode && verbose)
message("> pJSRforGene(): ", "Reached step 4.")
intervals <- (0:nrow(count))/nrow(count)
numexons <- nrow(count)
each <- merged.data$padjust[rt]
vertline.col <- ifelse(merged.data$testable[rt], ifelse(f.na(merged.data$padjust[rt] <=
FDR), SIG.VERTLINE.COLOR, NOSIG.VERTLINE.COLOR),
UNTESTABLE.VERTLINE.COLOR)
annolink.col <- ifelse(merged.data$testable[rt], ifelse(f.na(merged.data$padjust[rt] <=
FDR), SIG.FEATURE.COLOR, NOSIG.FEATURE.COLOR), UNTESTABLE.FEATURE.COLOR)
exonlty <- rep(final.color.list[["EXON.CONNECTION.LTY"]],
length(vertline.col))
exonlty[as.character(merged.data$featureType[rt]) ==
"novel_splice_site"] <- final.color.list[["NOVEL.SPLICE.CONNECTION.LTY"]]
exonlty[as.character(merged.data$featureType[rt]) ==
"splice_site"] <- final.color.list[["KNOWN.SPLICE.CONNECTION.LTY"]]
is.sig.feature <- f.na(each <= FDR)
sig.feature <- which(is.sig.feature)
if (debug.mode && verbose)
message("> pJSRforGene(): ", "Reached step 5.")
sub <- data.frame(start = merged.data$start[rt], end = merged.data$end[rt],
chr = merged.data$chr[rt], strand = merged.data$strand[rt],
is.exon = (merged.data$featureType[rt] == "exonic_part"),
is.testable = merged.data$testable[rt], is.sig = is.sig.feature,
col = vertline.col, featureType = merged.data$featureType[rt],
featureID = rownames(merged.data)[rt], stringsAsFactors = FALSE)
sub.allExon <- data.frame(start = flat.gff.data$start[rt.allExon],
end = flat.gff.data$end[rt.allExon], chr = flat.gff.data$chrom[rt.allExon],
strand = flat.gff.data$strand[rt.allExon], is.exon = (flat.gff.data$featureType[rt.allExon] ==
"exonic_part"), featureID = as.character(flat.gff.data$featureName[rt.allExon]),
stringsAsFactors = FALSE)
sub.allJunction <- data.frame(start = flat.gff.data$start[rt.allJunction],
end = flat.gff.data$end[rt.allJunction], chr = flat.gff.data$chrom[rt.allJunction],
strand = flat.gff.data$strand[rt.allJunction], is.exon = (flat.gff.data$featureType[rt.allJunction] ==
"exonic_part"), feature.type = flat.gff.data$featureType[rt.allJunction],
is.novel = (flat.gff.data$featureType[rt.allJunction] ==
"novel_splice_site"), featureID = as.character(flat.gff.data$featureName[rt.allJunction]),
stringsAsFactors = FALSE)
sub.allJunction$is.plotted <- sub.allJunction$featureID %in%
sub$featureID
testable.featureIDs <- sub$featureID[sub$is.testable]
sig.featureIDs <- sub$featureID[sub$is.sig]
untestable.featureIDs <- rownames(merged.data)[untestable.rt]
if (debug.mode && verbose) {
message("testable.featureIDs: ", paste0(testable.featureIDs,
collapse = ","))
}
if (debug.mode && verbose) {
message("sig.featureIDs: ", paste0(sig.featureIDs,
collapse = ","))
}
if (debug.mode && verbose) {
message("untestable.featureIDs: ", paste0(untestable.featureIDs,
collapse = ","))
}
sub.allExon$is.testable <- sub.allExon$featureID %in%
testable.featureIDs
sub.allExon$is.sig <- sub.allExon$featureID %in% sig.featureIDs
sub.allExon$is.untestable <- sub.allExon$featureID %in%
untestable.featureIDs
sub.allJunction$is.testable <- sub.allJunction$featureID %in%
testable.featureIDs
sub.allJunction$is.sig <- sub.allJunction$featureID %in%
sig.featureIDs
sub.allJunction$is.untestable <- sub.allJunction$featureID %in%
untestable.featureIDs
sub.allExon$lineColor <- ifelse(sub.allExon$is.testable,
ifelse(sub.allExon$is.sig, SIG.FEATURE.COLOR, NOSIG.FEATURE.COLOR),
ifelse(sub.allExon$is.untestable, UNTESTABLE.FEATURE.COLOR,
EXCLUDED.FEATURE.COLOR))
sub.allExon$fillColor <- ifelse(sub.allExon$is.testable,
ifelse(sub.allExon$is.sig, SIG.FEATURE.FILL.COLOR,
NOSIG.FEATURE.FILL.COLOR), ifelse(sub.allExon$is.untestable,
UNTESTABLE.FEATURE.FILL.COLOR, EXCLUDED.FEATURE.FILL.COLOR))
sub.allExon$borderColor <- ifelse(sub.allExon$is.testable,
ifelse(sub.allExon$is.sig, SIG.FEATURE.BORDER.COLOR,
NOSIG.FEATURE.COLOR), ifelse(sub.allExon$is.untestable,
UNTESTABLE.FEATURE.BORDER.COLOR, EXCLUDED.FEATURE.BORDER.COLOR))
sub.allJunction$lineColor <- ifelse(sub.allJunction$is.testable,
ifelse(sub.allJunction$is.sig, SIG.FEATURE.COLOR,
NOSIG.FEATURE.COLOR), ifelse(sub.allJunction$is.untestable,
UNTESTABLE.FEATURE.COLOR, EXCLUDED.FEATURE.COLOR))
sub.allJunction$lty <- ifelse(sub.allJunction$is.novel,
final.color.list[["NOVEL.SPLICE.LTY"]], final.color.list[["KNOWN.SPLICE.LTY"]])
sig.feature.names <- sub$featureID[is.sig.feature & sub$is.exon]
allExon.isSig <- sub.allExon$featureID %in% sig.featureIDs
allExon.exonCol <- ifelse(allExon.isSig, "#F219ED", "#CCCCCC")
if (debug.mode == 2 && verbose) {
message("> Debugging Info:")
message(" Exons:")
message(" ", paste0(names(sub.allExon),
collapse = "\t"))
for (i in 1:length(sub.allExon$start)) {
message(" ", paste0(sub.allExon[i,
], collapse = "\t"))
}
message("> dim(count) = ", paste0(dim(count),
collapse = ","))
message("> rt: ", paste0(rt, collapse = ","))
message("> rango: ", paste0(rango, collapse = ","))
message("> ncol(count): ", ncol(count))
}
sub.sig <- sub[sig.feature, , drop = FALSE]
rel.calc.min <- min(sub.allJunction$start, sub.allExon$start)
rel.calc.max <- max(sub.allJunction$end, sub.allExon$end)
transcripts <- sapply(sapply(flat.gff.data$transcripts[rt.allExon],
toString), function(x) {
strsplit(x, "+", fixed = TRUE)
})
trans <- Reduce(union, transcripts)
if (displayTranscripts == TRUE) {
mat <- 1:4
hei <- c(10, CONNECTIONS.height, GENE.annotation.height +
SPLICE.annotation.height, TX.annotation.height *
(length(trans) + TX.margins[1] + TX.margins[2]))
}
else {
mat <- 1:4
hei <- c(10, CONNECTIONS.height, GENE.annotation.height +
SPLICE.annotation.height, TX.annotation.height *
TX.margins[2])
}
layout(matrix(mat), heights = hei)
if (debug.mode && verbose) {
message("> length(trans) = ", length(trans))
message("> FINAL LAYOUT:")
message("> heights = [", paste0(hei, collapse = ","),
"]")
}
if (debug.mode && verbose)
message("> pJSRforGene(): ", "Reached step 6.")
ylimn <- c(min(min(count, na.rm = TRUE), 0), max(count,
na.rm = TRUE))
if ((!use.vst) && use.log)
ylimn[1] <- INTERNAL.NINF.VALUE
p.values.labels <- ifelse(f.na(each <= FDR), format(each,
digits = 3), "")
if (any(sub$is.exon)) {
italicize.label <- !sub$is.exon
}
else {
italicize.label <- NULL
}
plotWindowXmax <- if (plot.gene.level.expression) {
(length(intervals) + 1)/length(intervals)
}
else {
1
}
intervals <- drawPlot(matr = count, ylimn, jscs, intervals,
rango, textAxis = y.axis.title, geneLevelAxisTitle = y.axis.title.right,
rt = rt, color.count = color.count, colorlines = vertline.col,
countbinIDs = merged.data$countbinID[rt], use.vst = use.vst,
use.log = use.log, plot.type = plot.type, main.title = main.title,
draw.legend = draw.legend, color.key = color, condition.names = condition.names,
p.values = p.values.labels, draw.p.values = label.p.vals,
plot.lwd = plot.lwd, axes.lwd = axes.lwd, anno.lwd = anno.lwd,
par.cex = par.cex, anno.cex.text = anno.cex.text,
anno.cex.axis = anno.cex.axis, anno.cex.main = anno.cex.main,
fit.countbin.names = fit.countbin.names, debug.mode = debug.mode,
plot.gene.level.expression = plot.gene.level.expression,
geneCount = geneCount, color.geneCount = color.geneCount,
yAxisLabels.inExponentialForm = yAxisLabels.inExponentialForm,
italicize.label = italicize.label, condition.legend.text = condition.legend.text,
annolink.col = annolink.col, exonlty = exonlty, graph.margins = graph.margins,
plotWindowXmax = plotWindowXmax, fit.labels = fit.labels,
...)
if (debug.mode && verbose)
message("> pJSRforGene(): ", "Reached end of step 6.")
if (debug.mode && verbose)
message("> pJSRforGene(): ", "Reached step 7.")
if (USE.MARGIN.MEX) {
par(mar = c(0, graph.margins[2], 0, graph.margins[4]),
cex = par.cex, mex = anno.cex.text)
}
else {
par(mar = c(0, graph.margins[2], 0, graph.margins[4]),
cex = par.cex)
}
plot.new()
plot.window(xlim = c(0, plotWindowXmax), ylim = c(0,
1), xaxs = "i")
rel <- data.frame(start = rescale.coords(sub$start, rescale.iv),
end = rescale.coords(sub$end, rescale.iv))
connection.lines.bottom <- apply((rbind(rel[rango, 2],
rel[rango, 1])), 2, median) * plotWindowXmax
connection.lines.top <- apply(rbind(intervals[rango],
intervals[rango + 1] - ((intervals[rango + 1] - intervals[rango]) *
0.2)), 2, median)
segments(connection.lines.bottom, 0, connection.lines.top,
1, col = annolink.col, lty = exonlty, lwd = gene.lwd,
cex = anno.cex.text, cex.axis = anno.cex.main, cex.main = anno.cex.main,
xpd = NA, ...)
par(mar = c(1.5, graph.margins[2], 0, graph.margins[4]),
cex = par.cex)
startSites <- c()
endSites <- c()
for (i in 1:length(trans)) {
logicexons <- sapply(transcripts, function(x) {
any(x == trans[i])
})
startSites <- c(startSites, sub.allExon$start[logicexons][1])
endSites <- c(endSites, sub.allExon$end[logicexons][sum(logicexons)])
}
startSites <- unique(startSites)
endSites <- unique(endSites)
drawGene(rel.calc.min, rel.calc.max, tr = sub, tr.allExon = sub.allExon,
tr.allJunction = sub.allJunction, rango, rescale.iv = rescale.iv,
exoncol = annolink.col, allExon.exonCol = allExon.exonCol,
names, trName = "Gene model", anno.cex.text = anno.cex.text,
par.cex = par.cex, exonlty = exonlty, plot.lwd = gene.lwd,
anno.lwd = anno.lwd, show.strand.arrows = show.strand.arrows,
geneStrand = geneStrand, cex.axis = anno.cex.axis,
cex.main = anno.cex.main, draw.untestable.annotation = draw.untestable.annotation,
draw.start.end.sites = draw.start.end.sites, startSites = startSites,
endSites = endSites, cex.arrows = cex.arrows, chrom.label = chrom.label,
label.chromosome = label.chromosome, splice.junction.drawing.style = splice.junction.drawing.style,
draw.nested.SJ = draw.nested.SJ, merge.exon.parts = merge.exon.parts,
plot.untestable.results = plot.untestable.results,
exon.height = GENE.annotation.height/(SPLICE.annotation.height +
GENE.annotation.height), INTERNAL.VARS = INTERNAL.VARS,
flip.splicing = flip.splicing, ...)
include.endpoints.on.coordinates <- FALSE
num.coord.miniticks.per.tick <- 10
if (drawCoordinates) {
if (!is.null(rescale.iv)) {
pretty.x <- pretty(c(rel.calc.min, rel.calc.max),
n = 5)
pretty.interval <- pretty.x[2] - pretty.x[1]
pretty.x <- pretty.x[pretty.x > rel.calc.min &
pretty.x < rel.calc.max]
rescaled.pretty.x <- rescale.coords(pretty.x,
rescale.iv)
if (num.coord.miniticks.per.tick > 0) {
rel.coord.miniticks <- pretty.interval * (1:(num.coord.miniticks.per.tick -
1))/num.coord.miniticks.per.tick
unscaled.coord.miniticks <- unlist(lapply(pretty.x,
function(a) {
a + rel.coord.miniticks
}))
unscaled.coord.miniticks <- c(pretty.x[1] -
rel.coord.miniticks, unscaled.coord.miniticks)
rescaled.coord.miniticks <- rescale.coords(unscaled.coord.miniticks,
rescale.iv) * (rel.calc.max - rel.calc.min) +
rel.calc.min
}
else {
rescaled.coord.miniticks <- FALSE
}
if (include.endpoints.on.coordinates) {
if (min(rescaled.pretty.x) > 0.05) {
pretty.x <- c(rel.calc.min, pretty.x)
rescaled.pretty.x <- c(0, rescaled.pretty.x)
}
if (max(rescaled.pretty.x) < 0.95) {
pretty.x <- c(pretty.x, rel.calc.max)
rescaled.pretty.x <- c(rescaled.pretty.x,
1)
}
}
rescaled.pretty.x <- rescaled.pretty.x * (rel.calc.max -
rel.calc.min) + rel.calc.min
}
else {
pretty.x <- pretty(c(rel.calc.min, rel.calc.max),
n = 5)
coord.miniticks <- FALSE
}
usr <- par("usr")
cxy <- par("cxy")
pretty.x <- sprintf("%0.f", pretty.x)
smallest.width.coordAxis <- min(abs(rescaled.pretty.x[-1] -
rescaled.pretty.x[-length(rescaled.pretty.x)]))
if (fit.genomic.axis) {
anno.cex.coordAxis <- shrink.character.vector(paste0(pretty.x,
"0"), curr.cex = anno.cex.axis, max.width = smallest.width.coordAxis)
if (anno.cex.coordAxis * 2 < anno.cex.axis) {
anno.cex.coordAxis <- anno.cex.axis/2
coordAxis.widths <- abs(rescaled.pretty.x[-1] -
rescaled.pretty.x[-length(rescaled.pretty.x)])
for (i in 2:length(pretty.x)) {
if ((strwidth(pretty.x[i - 1], cex = anno.cex.coordAxis)/2) +
(strwidth(pretty.x[i], cex = anno.cex.coordAxis)/2) >
abs(rescaled.pretty.x[i - 1] - rescaled.pretty.x[i])) {
pretty.x[i] <- ""
}
}
}
}
else {
anno.cex.coordAxis <- anno.cex.axis
}
devlim <- device.limits()
coord.ticks.top <- usr[3]
coord.mainTicks.bottom <- usr[3] - (cxy[2]/2)
coord.text.top <- usr[3] - (cxy[2] * (3/4))
coord.miniTicks.bottom <- coord.ticks.top - abs(coord.mainTicks.bottom -
coord.ticks.top)/2
segments(x0 = rescaled.pretty.x, y0 = usr[3], x1 = rescaled.pretty.x,
y1 = usr[3] - (cxy[2]/2), xpd = NA, lwd = anno.lwd,
...)
lines(c(rel.calc.min, rel.calc.max), c(par("usr")[3],
par("usr")[3]), lwd = axes.lwd, xpd = NA, ...)
segments(x0 = rescaled.coord.miniticks, y0 = usr[3],
x1 = rescaled.coord.miniticks, y1 = usr[3] -
(cxy[2]/4), xpd = NA, lwd = anno.lwd, ...)
text(rescaled.pretty.x, usr[3] - (cxy[2] * (3/4)),
pretty.x, cex = anno.cex.coordAxis, xpd = NA,
adj = c(0.5, 1), ...)
}
if (label.chromosome) {
chrom.label.width.max <- abs(par("usr")[1] - device.limits()[1]) *
0.9
chrom.label.cex <- shrink.character.vector(chrom.label,
curr.cex = anno.cex.text, max.width = chrom.label.width.max)
text(par("usr")[1], par("usr")[3], chrom.label, cex = chrom.label.cex,
adj = c(1.1, 0.5), xpd = NA, font = 2, ...)
}
if (displayTranscripts) {
if (USE.MARGIN.MEX) {
par(cex = par.cex, mar = c(0, graph.margins[2],
0, graph.margins[4]), mex = anno.cex.text)
}
else {
par(cex = par.cex, mar = c(0, graph.margins[2],
0, graph.margins[4]))
}
plot.new()
plot.window(xlim = c(rel.calc.min, rel.calc.max),
ylim = c(-TX.margins[2], length(trans) + TX.margins[1]),
xaxs = "i", yaxs = "i")
for (i in 1:length(trans)) {
ymin <- length(trans) - i
if (include.TX.names) {
trName = trans[i]
}
else {
trName = NULL
}
logicexons <- sapply(transcripts, function(x) {
any(x == trans[i])
})
tr <- sub.allExon[logicexons, ]
curr.exoncol <- ifelse(allExon.isSig[logicexons],
"#F219ED", "#CCCCCC")
drawTranscript(rel.calc.min, rel.calc.max, ymin = ymin,
tr = tr, tr.allJunction = sub.allJunction,
rango = 1:nrow(tr), rescale.iv = rescale.iv,
names = c(), trName = trName, trStrand = txStrandMap[[trName]],
draw.strand = geneStrand == ".", par.cex = par.cex,
anno.cex.text = anno.cex.text, sub.sig = sub.sig,
anno.lwd = gene.lwd, cex.axis = anno.cex.axis,
cex.main = anno.cex.main, cex.arrows = cex.arrows,
...)
}
}
else {
par(mar = c(0, 0, 0, 0))
plot.new()
}
par(mar = c(0, 0, 0, 0))
if (debug.mode && verbose)
message("> pJSRfG(): ", " Done.")
}
}, error = function(e) {
message("Error caught while attempting plotJunctionSeqResultsForGene")
message("---------------------")
message(" Error text:")
message(" ", e)
message("")
message("---------------------")
message(" Input parameters:")
message(" geneID = ", geneID)
message(" plot.type = ", plot.type)
message(" displayTranscripts = ", displayTranscripts)
message("---------------------")
warning("Error caught while attempting plotJunctionSeqResultsForGene")
warning("Error text:")
warning(e)
})
where 10: plotJunctionSeqResultsForGene(geneID, jscs, colorRed.FDR.threshold = FDR,
colorList = colorList, plot.type = plot.type, use.vst = use.vst,
use.log = use.log, exon.rescale.factor = exon.rescale.factor,
displayTranscripts = TRUE, plot.lwd = plot.lwd, axes.lwd = axes.lwd,
anno.lwd = anno.lwd, par.cex = par.cex, anno.cex.text = anno.cex.text,
plot.exon.results = plot.exon.results, plot.junction.results = plot.junction.results,
plot.novel.junction.results = plot.novel.junction.results,
plot.untestable.results = plot.untestable.results, anno.cex.axis = anno.cex.axis,
anno.cex.main = anno.cex.main, drawCoordinates = drawCoordinates,
graph.margins = graph.margins, yAxisLabels.inExponentialForm = yAxisLabels.inExponentialForm,
plot.gene.level.expression = plot.gene.level.expression,
condition.legend.text = condition.legend.text, include.TX.names = include.TX.names,
GENE.annotation.relative.height = GENE.annotation.relative.height,
TX.annotation.relative.height = TX.annotation.relative.height,
draw.start.end.sites = draw.start.end.sites, show.strand.arrows = show.strand.arrows,
debug.mode = debug.mode, sequencing.type = sequencing.type,
gene.lwd = gene.lwd, CONNECTIONS.relative.height = CONNECTIONS.relative.height,
TX.margins = TX.margins, SPLICE.annotation.relative.height = SPLICE.annotation.relative.height,
draw.nested.SJ = draw.nested.SJ, INTERNAL.VARS = INTERNAL.VARS,
...)
where 11: buildAllPlotsForGene(geneID = geneID, jscs = jscs, outfile.prefix = outfile.prefixes,
use.plotting.device = use.plotting.device, use.vst = use.vst,
use.log = use.log, exon.rescale.factor = exon.rescale.factor,
plot.gene.level.expression = plot.gene.level.expression,
with.TX = with.TX, without.TX = without.TX, expr.plot = expr.plot,
normCounts.plot = normCounts.plot, rExpr.plot = rExpr.plot,
rawCounts.plot = rawCounts.plot, colorRed.FDR.threshold = colorRed.FDR.threshold,
colorList = colorList, plot.exon.results = plot.exon.results,
plot.junction.results = plot.junction.results, plot.novel.junction.results = plot.novel.junction.results,
plot.untestable.results = plot.untestable.results, plot.lwd = plot.lwd,
axes.lwd = axes.lwd, anno.lwd = anno.lwd, gene.lwd = gene.lwd,
drawCoordinates = drawCoordinates, show.strand.arrows = show.strand.arrows,
graph.margins = graph.margins, par.cex = par.cex, anno.cex.text = anno.cex.text,
anno.cex.axis = anno.cex.axis, anno.cex.main = anno.cex.main,
base.plot.height = base.plot.height, base.plot.width = base.plot.width,
base.plot.units = base.plot.units, GENE.annotation.relative.height = GENE.annotation.relative.height,
TX.annotation.relative.height = TX.annotation.relative.height,
CONNECTIONS.relative.height = CONNECTIONS.relative.height,
TX.margins = TX.margins, SPLICE.annotation.relative.height = SPLICE.annotation.relative.height,
autoscale.height.to.fit.TX.annotation = autoscale.height.to.fit.TX.annotation,
autoscale.width.to.fit.bins = autoscale.width.to.fit.bins,
name.files.with.geneID = name.files.with.geneID, plotting.device.params = plotting.device.params,
yAxisLabels.inExponentialForm = yAxisLabels.inExponentialForm,
condition.legend.text = condition.legend.text, include.TX.names = include.TX.names,
draw.start.end.sites = draw.start.end.sites, verbose = verbose,
debug.mode = debug.mode, sequencing.type = sequencing.type,
minimalImageFilenames = minimalImageFilenames, INTERNAL.VARS = INTERNAL.VARS,
...)
where 12: buildAllPlots(jscs)
--- value of length: 3 type: logical ---
[1] FALSE TRUE FALSE
--- function from context ---
function (expr, name, parentenv, handler)
{
.Internal(.addCondHands(name, list(handler), parentenv, environment(),
FALSE))
expr
}
<bytecode: 0x0000000005071700>
<environment: 0x00000000317393c0>
--- function search by body ---
----------- END OF FAILURE REPORT --------------
Fatal error: length > 1 in coercion to logical
* checking for unstated dependencies in vignettes ... OK
* checking package vignettes in 'inst/doc' ... OK
* checking running R code from vignettes ... SKIPPED
* checking re-building of vignette outputs ... SKIPPED
* checking PDF version of manual ... OK
* DONE
Status: 2 ERRORs, 1 WARNING, 3 NOTEs
See
'C:/Users/biocbuild/bbs-3.11-bioc/meat/JunctionSeq.Rcheck/00check.log'
for details.