Error in cf[[pvn]] : subscript out of bounds

[nicholas-owen]

Dear @hartleys ,

Thanks you for putting together a great bit of software and continuing to support it.

I have been using JunctionSeq on some RNA-seq data happily for a while but recently have come into a problem with some of my pairwise tests. I have a time series over 3 time points with 2 conditions that I am comparing in a pairwise manner, as well as doing as one.

When I compare the two conditions at a time point JunctionSeq runs without error but when I do across time points with the same condition I am getting a script out of bounds error. Understandably there will be a lot more differences over time compared to between the conditions as the conditions are related here.

I have tried to use method.dispFit='local', and nCores =1 to overcome fitting errors but still get this:

`   rJSA: Estimating effect sizes using effects models... Thu Sep 05 18:25:59 2019.
> Using single-core execution.
Error in cf[[pvn]] : subscript out of bounds
In addition: Warning messages:
1: In MulticoreParam(workers = nCores) :
  MulticoreParam() not supported on Windows, use SnowParam()
2: In MulticoreParam(workers = nCores) :
  MulticoreParam() not supported on Windows, use SnowParam()

`

If you have any suggestions they would be most appreciated. I have 4/5 biological replicates per condition per time point.

Thanks!

Below is the output of the command I am using :
jscs2<-runJunctionSeqAnalyses(sample.files= countFiles,sample.names= decoder$sample.ID, condition=decoder$group.ID,flat.gff.file="./withNovel.forJunctionSeq.gff.gz", nCores= 1,verbose= TRUE,debug.mode= TRUE, method.dispFit='local', optimizeFilteringForAlpha = 0.05)

output:

`> STARTING runJunctionSeqAnalyses (v1.14.0) (Thu Sep 05 18:10:51 2019)
> rJSA: sample.files:  ./1-48-O/QC.spliceJunctionAndExonCounts.forJunctionSeq.txt.gz, ./2-32-O/QC.spliceJunctionAndExonCounts.forJunctionSeq.txt.gz, ./2-48-O/QC.spliceJunctionAndExonCounts.forJunctionSeq.txt.gz, ./3-32-O/QC.spliceJunctionAndExonCounts.forJunctionSeq.txt.gz, ./3-48-O/QC.spliceJunctionAndExonCounts.forJunctionSeq.txt.gz, ./4-32-O/QC.spliceJunctionAndExonCounts.forJunctionSeq.txt.gz, ./4-48-O/QC.spliceJunctionAndExonCounts.forJunctionSeq.txt.gz, ./5-32-O/QC.spliceJunctionAndExonCounts.forJunctionSeq.txt.gz
> rJSA: sample.names:  1-48-O, 2-32-O, 2-48-O, 3-32-O, 3-48-O, 4-32-O, 4-48-O, 5-32-O
> rJSA: condition:  48O, 32O, 48O, 32O, 48O, 32O, 48O, 32O
> 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 Sep 05 18:10:51 2019.
-> STARTING readJunctionSeqCounts (Thu Sep 05 18:10:51 2019)
---> RJSC; (v1.14.0)
---> RJSC: samplenames: 1-48-O,2-32-O,2-48-O,3-32-O,3-48-O,4-32-O,4-48-O,5-32-O
---> RJSC: flat.gff.file: ./withNovel.forJunctionSeq.gff.gz
---> RJSC: use.exons:TRUE
---> RJSC: use.junctions:TRUE
---> RJSC: use.novel.junctions:TRUE
---> File read complete.
---> Extracted counts. Found 615968 features so far.
---> Extracted gene-level counts. Found: 29650 genes and aggregate-genes.
---> Removed gene features. Found: 586318 features to be included so far.
---> Note: 103900 counting bins from overlapping genes
--->          There are 2227 multigene aggregates.
--->          There are 5097 genes that are part of an aggregate.
---> Removed multigene-aggregate features. Found: 482418 features to be included so far.
---> Final feature count: 482418 features to be included in the analysis.
---> Extracted feature counts.
---> counts complete.
-----> reading annotation...
-----> formatting annotation...
-----> initial generation...
-----> creating jscs...
-----> generating count vectors... (Thu Sep 05 18:11:36 2019)
> Using single-core execution.
    getAllJunctionSeqCountVectors: dim(counts) = 482418,8 (Thu Sep 05 18:11:36 2019)
    getAllJunctionSeqCountVectors: dim(gct) = 29650,8
    getAllJunctionSeqCountVectors: out generated. dim = 482418,16 (Thu Sep 05 18:11:51 2019)
-----> count vectors generated (Thu Sep 05 18:11:51 2019)
-----> generating DESeqDataSet... (Thu Sep 05 18:11:51 2019)
  Note: levels of factors in the design contain characters other than
  letters, numbers, '_' and '.'. It is recommended (but not required) to use
  only letters, numbers, and delimiters '_' or '.', as these are safe characters
  for column names in R. [This is a message, not an warning or error]
-----> DESeqDataSet generated (Thu Sep 05 18:11:51 2019)
> rJSA: Count files read. Thu Sep 05 18:11:51 2019.
> rJSA: Estimating Size Factors... Thu Sep 05 18:11:51 2019.
  Note: levels of factors in the design contain characters other than
  letters, numbers, '_' and '.'. It is recommended (but not required) to use
  only letters, numbers, and delimiters '_' or '.', as these are safe characters
  for column names in R. [This is a message, not an warning or error]
> rJSA: Size Factors Done. Size Factors are:.
> rJSA: 1-48-O,2-32-O,2-48-O,3-32-O,3-48-O,4-32-O,4-48-O,5-32-O
> rJSA: 0.950088802644262,0.816513142536775,1.0547054985192,1.13446945547685,1.3975792020582,0.894781568884979,1.08019902113532,0.995497988557683
> rJSA: Estimating Dispersions... Thu Sep 05 18:11:53 2019.
---> STARTING estimateJunctionSeqDispersions: (v1.14.0) (Thu Sep 05 18:11:53 2019)
-----> ejsd: 108422 counting bins are marked 'testable'. across 15344 genes.
             (107183 exonic regions, 1239 known junctions, 0 novel junctions)
---------> Executing DESeq2 call: estimateUnsharedDispersions
  Note: levels of factors in the design contain characters other than
  letters, numbers, '_' and '.'. It is recommended (but not required) to use
  only letters, numbers, and delimiters '_' or '.', as these are safe characters
  for column names in R. [This is a message, not an warning or error]
---------> Finished with DESeq2 call.
-----> ejsd: Dispersion estimation failed for 0 out of 108422 'testable' counting bins. Setting these features to be 'untestable'
---> FINISHED estimateJunctionSeqDispersions (Thu Sep 05 18:22:59 2019)
> rJSA: Dispersions estimated. Thu Sep 05 18:22:59 2019.
> rJSA: Fitting Dispersion Fcn... Thu Sep 05 18:22:59 2019.
> fitDispersionFunction() Starting (Thu Sep 05 18:22:59 2019)
>   (fitType = local)
>   (finalDispersionMethod = shrink)
>   (fitDispersionsForExonsAndJunctionsSeparately = TRUE)
min(means[useForFit], na.rm=T)=0.0894403693299915
>    fdf: Fitting dispersions:
>    fdf: Fitting dispersions of exons and junctions to separate fitted trends.
>    fdf: Fitting exon dispersions:
>    fdf: Fitting splice-junction dispersions:
> fdf(): 'Shrinking' fitted and feature-specific dispersion estimates.
> fdf() Dispersion estimate failed for 0 out of 108422 features.
> fitDispersionFunction() Done. (Thu Sep 05 18:23:49 2019)
> rJSA: Dispersions Fcn Fitted. Thu Sep 05 18:23:49 2019.
> rJSA: Testing for DEU... Thu Sep 05 18:23:49 2019.
> Using single-core execution.
-------> testJunctionsForDiffUsage: Starting hypothesis test iteration. (Thu Sep 05 18:23:49 2019)
-------> testJunctionsForDiffUsage: Finished hypothesis test iteration. (Thu Sep 05 18:25:56 2019)
  Note: levels of factors in the design contain characters other than
  letters, numbers, '_' and '.'. It is recommended (but not required) to use
  only letters, numbers, and delimiters '_' or '.', as these are safe characters
  for column names in R. [This is a message, not an warning or error]
-------> testJunctionsForDiffUsage: Finished compiling hypothesis test results. (Thu Sep 05 18:25:58 2019)
---> tJfDU(): No non-NA maxCooks values. Ignoring cooks.
> Performing final p.adjust filtering.
>      No cook's cutoffs found.
>      Automatically selecting a filtering threshold of 1.26775688708707 to optimize results at the alpha < 0.05 significance level.
>         (Filtering 71212 out of 108422 "testable" features, using baseMean < 1.26775688708707)
>         (Rejected H0 for 1396 out of 37210 features at alpha < 0.05)
> Final p.adjust filtering complete.
> rJSA: DEU tests complete. Thu Sep 05 18:25:59 2019.
> rJSA: Estimating effect sizes using effects models... Thu Sep 05 18:25:59 2019.
> Using single-core execution.
Error in cf[[pvn]] : subscript out of bounds
In addition: Warning messages:
1: In MulticoreParam(workers = nCores) :
  MulticoreParam() not supported on Windows, use SnowParam()
2: In MulticoreParam(workers = nCores) :
  MulticoreParam() not supported on Windows, use SnowParam()`

sessionInfo:

`R version 3.6.0 (2019-04-26)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 17134)

Matrix products: default

Random number generation:
 RNG:     Mersenne-Twister 
 Normal:  Inversion 
 Sample:  Rounding 
 
locale:
[1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United Kingdom.1252    LC_MONETARY=English_United Kingdom.1252
[4] LC_NUMERIC=C                            LC_TIME=English_United Kingdom.1252    

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] QoRTs_1.3.6                 JunctionSeq_1.14.0          RcppArmadillo_0.9.600.4.0   Rcpp_1.0.2                  edgeR_3.26.6               
 [6] limma_3.40.6                DEXSeq_1.30.0               RColorBrewer_1.1-2          AnnotationDbi_1.46.0        DESeq2_1.24.0              
[11] SummarizedExperiment_1.14.1 DelayedArray_0.10.0         BiocParallel_1.18.0         matrixStats_0.54.0          Biobase_2.44.0             
[16] GenomicRanges_1.36.0        GenomeInfoDb_1.20.0         IRanges_2.18.1              S4Vectors_0.22.0            BiocGenerics_0.30.0        
[21] ballgown_2.16.0            

loaded via a namespace (and not attached):
 [1] nlme_3.1-141             bitops_1.0-6             bit64_0.9-7              progress_1.2.2           httr_1.4.1               tools_3.6.0             
 [7] backports_1.1.4          R6_2.4.0                 rpart_4.1-15             Hmisc_4.2-0              DBI_1.0.0                lazyeval_0.2.2          
[13] mgcv_1.8-28              colorspace_1.4-1         nnet_7.3-12              tidyselect_0.2.5         gridExtra_2.3            prettyunits_1.0.2       
[19] bit_1.1-14               compiler_3.6.0           htmlTable_1.13.1         rtracklayer_1.44.2       scales_1.0.0             checkmate_1.9.4         
[25] genefilter_1.66.0        stringr_1.4.0            digest_0.6.20            Rsamtools_2.0.0          foreign_0.8-71           XVector_0.24.0          
[31] base64enc_0.1-3          pkgconfig_2.0.2          htmltools_0.3.6          plotrix_3.7-6            htmlwidgets_1.3          rlang_0.4.0             
[37] rstudioapi_0.10          RSQLite_2.1.2            hwriter_1.3.2            acepack_1.4.1            dplyr_0.8.3              RCurl_1.95-4.12         
[43] magrittr_1.5             GenomeInfoDbData_1.2.1   Formula_1.2-3            Matrix_1.2-17            munsell_0.5.0            stringi_1.4.3           
[49] zlibbioc_1.30.0          grid_3.6.0               blob_1.2.0               crayon_1.3.4             lattice_0.20-38          Biostrings_2.52.0       
[55] splines_3.6.0            annotate_1.62.0          hms_0.5.0                locfit_1.5-9.1           zeallot_0.1.0            knitr_1.23              
[61] pillar_1.4.2             geneplotter_1.62.0       biomaRt_2.40.3           XML_3.98-1.20            glue_1.3.1               latticeExtra_0.6-28     
[67] BiocManager_1.30.4       data.table_1.12.2        vctrs_0.2.0              gtable_0.3.0             purrr_0.3.2              assertthat_0.2.1        
[73] ggplot2_3.2.0            xfun_0.8                 xtable_1.8-4             survival_2.44-1.1        tibble_2.1.3             GenomicAlignments_1.20.1
[79] memoise_1.1.0            cluster_2.0.8            sva_3.32.1               statmod_1.4.32 `