Question: Replacement error in PermTest
gravatar for AlejandroMA
4 months ago by
AlejandroMA0 wrote:


I was trying to run a permutation test with the regioneR package these days but I keep getting the same error:

Error in if (alternative == "less" & orig.ev > mean(rand.ev, na.rm = TRUE)) message("Alternative is less and the observed statistic is greater than the permuted statistic mean. Maybe you want to use recomputePermTest to change the alternative hypothesis.") : 
  missing value where TRUE/FALSE needed
In addition: Warning messages:
1: In mclapply(c(1:ntimes), randomize_and_evaluate, ...) :
  all scheduled cores encountered errors in user code
2: In mean.default(rand.ev, na.rm = TRUE) :
  argument is not numeric or logical: returning NA
3: In mean.default(rand.ev, na.rm = TRUE) :
  argument is not numeric or logical: returning NA

So I understand from a previous post here ( that I have some character in my tables and that's why the rand.ev function is not "numeric or logical". When I follow the suggestions of the post and set force.parallel to False I get this:

Error in `$<`(`*tmp*`, "start", value = numeric(0)) : 
  replacement has 0 rows, data has 51

I've checked the bed files but no extra character had slipped so I don't get why rand.ev is NA then. Any Idea of what might be driving the error?

Thank you in advance,

Data Both my bed files have this structure:

chr1 1179385 1179385 chr1 11858584 11858584 chr1 11859390 11859390


one <- import("one.bed")
two <- import("two.bed")

pt1 <- permTest(A=one, B=two, ntimes=10,
                randomize.function=circularRandomizeRegions, alternative = "less",
                evaluate.function=numOverlaps, force.parallel = FALSE)


R version 3.6.1 (2019-07-05)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/

 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=es_ES.UTF-8       
 [7] LC_PAPER=es_ES.UTF-8       LC_NAME=C                  LC_ADDRESS=C              

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

other attached packages:
 [1] BSgenome.Hsapiens.UCSC.hg19.masked_1.3.99 BSgenome.Hsapiens.UCSC.hg19_1.4.0        
 [3] BSgenome_1.52.0                           Biostrings_2.52.0                        
 [5] XVector_0.24.0                            rtracklayer_1.44.0                       
 [7] regioneR_1.16.2                           GenomicRanges_1.36.0                     
 [9] GenomeInfoDb_1.20.0                       IRanges_2.18.1                           
[11] S4Vectors_0.22.0                          BiocGenerics_0.30.0                      
[13] readxl_1.3.1                              readr_1.3.1                              

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.1                  cellranger_1.1.0            pillar_1.4.2               
 [4] compiler_3.6.1              bitops_1.0-6                tools_3.6.1                
 [7] zlibbioc_1.30.0             digest_0.6.20               zeallot_0.1.0              
[10] memoise_1.1.0               tibble_2.1.3                lattice_0.20-38            
[13] pkgconfig_2.0.2             rlang_0.4.0                 Matrix_1.2-17              
[16] DelayedArray_0.10.0         GenomeInfoDbData_1.2.1      vctrs_0.2.0                
[19] hms_0.5.0                   grid_3.6.1                  Biobase_2.44.0             
[22] R6_2.4.0                    XML_3.98-1.20               BiocParallel_1.18.0        
[25] backports_1.1.4             Rsamtools_2.0.0             matrixStats_0.54.0         
[28] GenomicAlignments_1.20.1    SummarizedExperiment_1.14.0 RCurl_1.95-4.12            
[31] crayon_1.3.4    
regioner • 83 views
ADD COMMENTlink written 4 months ago by AlejandroMA0

Hi Alejandro,

Can you check the size of the segments in "one"? I assume you are trying to represent segments of length 1, but BED is a 0-based format and you are representing segments of width 0. What's the output of width(one)?

I think this might be at least part of the problem


ADD REPLYlink written 4 months ago by bernatgel120


Thank you for the reply. Managed to make it work, it was indeed a problem with the SNPs bed file (one) being 1-based even if it was a bed file. Took me some time to realize because I was convinced I had processed it properly before in the work pipeline, but apparently I missed it. So I just changed import and used this instead:

one <- read_table2("one")
one <-
one <- toGRanges(one)

I should have noticed before, silly error. Thank you for the help!


ADD REPLYlink written 4 months ago by AlejandroMA0

Great :) Glad it worked in the end. I've been bitten by this more than once!

ADD REPLYlink written 4 months ago by bernatgel120
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 244 users visited in the last hour