Question: Replacement error in PermTest
gravatar for AlejandroMA
13 days 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 • 32 views
ADD COMMENTlink written 13 days 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 13 days ago by bernatgel100


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 12 days ago by AlejandroMA0

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

ADD REPLYlink written 12 days ago by bernatgel100
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: 333 users visited in the last hour