Problems running JunctionSeq
0
0
Entering edit mode
jbono ▴ 10
@jbono-7682
Last seen 25 days ago
United States

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 

junctionseq • 507 views
ADD COMMENT

Login before adding your answer.

Traffic: 199 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6