Hi,
I am trying to run the JunctionSeq package, but am running into a problem with the runJunctionsSeqAnalyses function that I can't seem to solve. I entered the following:
> jscs <- runJunctionSeqAnalyses(
+ sample.files = countFiles,
+ sample.names = decoder$sample.ID,
+ condition = decoder$group.ID,
+ flat.gff.file = "outputData/countTables/withNovel.forJunctionSeq.gff.gz",
+ nCores = 1,
+ verbose=TRUE,
+ debug.mode = TRUE
It appears there might be an issue with DESeq (see below). Furthermore,the function does not seem to produce the jscs object because when I try the writeCompleteResults function it says it cannot find the jscs object. Any help on what might be causing the problem would be greatly appreciated.
This is what it returns from runJunctionSeqAnalyses (DEseq part is at the end):
> STARTING runJunctionSeqAnalyses (v1.3.5) (Thu Nov 17 14:47:55 2016)
> rJSA: sample.files: outputData/countTables/C1/QC.spliceJunctionAndExonCounts.withNovel.forJunctionSeq.txt.gz, outputData/countTables/C2/QC.spliceJunctionAndExonCounts.withNovel.forJunctionSeq.txt.gz, outputData/countTables/C3/QC.spliceJunctionAndExonCounts.withNovel.forJunctionSeq.txt.gz, outputData/countTables/O1/QC.spliceJunctionAndExonCounts.withNovel.forJunctionSeq.txt.gz, outputData/countTables/O2/QC.spliceJunctionAndExonCounts.withNovel.forJunctionSeq.txt.gz, outputData/countTables/O3/QC.spliceJunctionAndExonCounts.withNovel.forJunctionSeq.txt.gz
> rJSA: sample.names: C1, C2, C3, O1, O2, O3
> rJSA: condition: caper, caper, caper, control, control, control
> rJSA: analysis.type: junctionsAndExons
> rJSA: use.junctions: TRUE
> rJSA: use.novel.junctions: TRUE
> rJSA: use.exons: TRUE
> rJSA: nCores: 1
> rJSA: use.covars:
> rJSA: test.formula0: ~ sample + countbin
> rJSA: test.formula1: ~ sample + countbin + condition:countbin
> rJSA: use.multigene.aggregates: FALSE
> rJSA: Reading Count files... Thu Nov 17 14:47:55 2016.
-> STARTING readJunctionSeqCounts (Thu Nov 17 14:47:55 2016)
---> RJSC; (v1.3.5)
---> RJSC: samplenames: C1,C2,C3,O1,O2,O3
---> RJSC: flat.gff.file: outputData/countTables/withNovel.forJunctionSeq.gff.gz
---> RJSC: use.exons:TRUE
---> RJSC: use.junctions:TRUE
---> RJSC: use.novel.junctions:TRUE
---> File read complete.
---> Extracted counts. Found 160953 features so far.
---> Extracted gene-level counts. Found: 16101 genes and aggregate-genes.
---> Removed gene features. Found: 144852 features to be included so far.
---> Note: 108353 counting bins from overlapping genes
---> There are 7248 multigene aggregates.
---> There are 25080 genes that are part of an aggregate.
---> Removed multigene-aggregate features. Found: 36499 features to be included so far.
---> Final feature count: 36499 features to be included in the analysis.
---> Extracted feature counts.
---> counts complete.
-----> reading annotation...
-----> formatting annotation...
-----> initial generation...
-----> creating jscs...
-----> generating count vectors... (Thu Nov 17 14:48:10 2016)
> Using single-core execution.
getAllJunctionSeqCountVectors: dim(counts) = 36499,6 (Thu Nov 17 14:48:10 2016)
getAllJunctionSeqCountVectors: dim(gct) = 16101,6
getAllJunctionSeqCountVectors: out generated. dim = 36499,12 (Thu Nov 17 14:48:10 2016)
-----> count vectors generated (Thu Nov 17 14:48:10 2016)
-----> generating DESeqDataSet... (Thu Nov 17 14:48:10 2016)
Error in DESeqDataSet(se, design = design, ignoreRank) :
all samples have 0 counts for all genes. check the counting script.
This is returned from the writeCompleteResults function:
> STARTING writeCompleteResults (v1.3.5) (Thu Nov 17 14:57:36 2016)
> wcr: outfile.prefix: outputData/junctionseq
> wcr: FDR.threshold: 0.01
> wcr: save.allGenes: TRUE
> wcr: save.sigGenes: TRUE
> wcr: save.fit: FALSE
> wcr: save.VST: FALSE
> wcr: bedtrack.format: BED
Error in save(jscs, file = paste0(outfile.prefix, "jscs.RData")) :
object 'jscs' not found
Thank you, I appreciate it!
Jeremy
The output of the sessioninfo command is below:
R version 3.2.3 (2015-12-10)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.11.3 (El Capitan)
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats4 parallel stats graphics grDevices utils datasets methods base
other attached packages:
[1] geneplotter_1.48.0 annotate_1.48.0 XML_3.98-1.5 AnnotationDbi_1.32.3
[5] genefilter_1.52.1 BiocParallel_1.4.3 Hmisc_4.0-0 ggplot2_2.2.0
[9] Formula_1.2-1 survival_2.40-1 lattice_0.20-34 locfit_1.5-9.1
[13] stringr_1.1.0 plotrix_3.6-3 statmod_1.4.26 DESeq2_1.10.1
[17] RcppArmadillo_0.7.500.0.0 Rcpp_0.12.7 JunctionSeq_1.3.5 SummarizedExperiment_1.0.2
[21] Biobase_2.30.0 GenomicRanges_1.22.4 GenomeInfoDb_1.6.3 IRanges_2.4.8
[25] S4Vectors_0.8.11 BiocGenerics_0.16.1 BiocInstaller_1.20.3
loaded via a namespace (and not attached):
[1] splines_3.2.3 colorspace_1.3-0 htmltools_0.3.5 chron_2.3-47 DBI_0.5-1
[6] foreign_0.8-67 RColorBrewer_1.1-2 lambda.r_1.1.9 plyr_1.8.4 zlibbioc_1.16.0
[11] munsell_0.4.3 gtable_0.2.0 futile.logger_1.4.3 latticeExtra_0.6-28 knitr_1.15
[16] htmlTable_1.7 acepack_1.4.1 xtable_1.8-2 scales_0.4.1 XVector_0.10.0
[21] gridExtra_2.2.1 digest_0.6.10 stringi_1.1.2 grid_3.2.3 tools_3.2.3
[26] magrittr_1.5 RSQLite_1.0.0 lazyeval_0.2.0 tibble_1.2 cluster_2.0.5
[31] futile.options_1.0.0 Matrix_1.2-7.1 data.table_1.9.6 assertthat_0.1 rpart_4.1-10
[36] nnet_7.3-12
