search

cqgd / pky

Other
3 stars 3 forks source link

Error in if (dv > olddv && itn >= 2 && auto == TRUE) { : missing value where TRUE/FALSE needed #9

Open karolnowosad opened 1 year ago

karolnowosad commented 1 year ago

Hi, I try to use the Peaky with my dataset but it generates the error message followed by warning message:

"Error in if (dv > olddv && itn >= 2 && auto == TRUE) { : 
  missing value where TRUE/FALSE needed
  In addition: Warning message:
  In pnbinom(q, size = 1/sigma, mu = mu, lower.tail = lower.tail,  :
  NaNs produced"

I use the Peaky on a large subset, however, I set the max distance 1Mb to obtain high variance of N. 

 BI_fly = bin_interactions(
 interactions_fly_2,
 fragments_fly_2,
 bins=5,
 min_dist=2000, 
 max_dist=1e6
 )

I tried to build the model using all interactions (it also ends with the error), but for purpose of visualisation I set the subsample_size=1000.

I run separately the model for each bin as multiple times Peaky output the error while computing the model for different bins using the code: "models = by(BI_fly$interactions, BI_fly$interactions$dist.bin, model_bin, subsample_size=1000)".

As the bin No 1 has the highest variance of N (please see the output below) I decided to run the model multiple times on the bin No 1 with different results. Either it fails (generating the abovementioned error) or finish with success. It is understandable that this depends on the selected subsed but what leads to the generation of the error?

Could you please help me to solve this issue?

BI_fly = readRDS(file = file_path)

#extract BIN and interaction table from BI_fly object
interactions <- BI_fly$interactions
bins <- BI_fly$bins

Check bins from BI object

bins

##    dist.bin dist.abs.min dist.abs.max interactions
## 1:        1         2000        71986      2851024
## 2:        2        71987       179570      2850992
## 3:        3       179571       343639      2850987
## 4:        4       343640       601404      2851002
## 5:        5       601405      1000000      2850993
## 6:     <NA>           NA           NA     43326975

Check interaction table from BI object

head(interactions)

##    baitID preyID  N b.chr b.mid b.length p.chr p.mid p.length   dist b.trans
## 1:     83      2  2 chr2L 22222      255 chr2L   153      148 -22069    6790
## 2:     83      4 14 chr2L 22222      255 chr2L   452      167 -21770    6790
## 3:     83      5  1 chr2L 22222      255 chr2L   611      148 -21611    6790
## 4:     83      6  1 chr2L 22222      255 chr2L   756      140 -21466    6790
## 5:     83      7  4 chr2L 22222      255 chr2L   910      167 -21312    6790
## 6:     83      9  2 chr2L 22222      255 chr2L  1214      140 -21008    6790
##    b.trans_res p.trans_res dist.bin
## 1:   0.7830051           0        1
## 2:   0.7830051           0        1
## 3:   0.7830051           0        1
## 4:   0.7830051           0        1
## 5:   0.7830051           0        1
## 6:   0.7830051           0        1

check the variance of the reads in each BIN

variance_N <- function(data, bin){
  variance_N_bin <- data$interactions %>%
  subset(dist.bin == bin) %>%
  select(N) %>%
  table %>%
  t %>%
  as.data.frame()
variance_N_bin$Var1 <- rep(paste0("BIN_",bin), nrow(variance_N_bin))
variance_N_bin
}

#test 
bin_No <- 1:5

First 100 rows are displayed 

#BIN_1
variance_N_bin_1 <- variance_N(data = BI_fly, bin = bin_No[1])
variance_N_bin_1

##      Var1     N   Freq
## 1   BIN_1     1 558184
## 2   BIN_1     2 443991
## 3   BIN_1     3 340851
## 4   BIN_1     4 260714
## 5   BIN_1     5 201910
## 6   BIN_1     6 159406
## 7   BIN_1     7 126499
## 8   BIN_1     8 102263
## 9   BIN_1     9  84611
## 10  BIN_1    10  69901
## 11  BIN_1    11  57920
## 12  BIN_1    12  49296
## 13  BIN_1    13  42116
## 14  BIN_1    14  36020
## 15  BIN_1    15  31224
## 16  BIN_1    16  27137
## 17  BIN_1    17  23456
## 18  BIN_1    18  20587
## 19  BIN_1    19  18422
## 20  BIN_1    20  16247
## 21  BIN_1    21  14181
## 22  BIN_1    22  12936
## 23  BIN_1    23  11542
## 24  BIN_1    24  10530
## 25  BIN_1    25   9471
## 26  BIN_1    26   8460
## 27  BIN_1    27   7856
## 28  BIN_1    28   7049
## 29  BIN_1    29   6391
## 30  BIN_1    30   5849
## 31  BIN_1    31   5484
## 32  BIN_1    32   5073
## 33  BIN_1    33   4622
## 34  BIN_1    34   4193
## 35  BIN_1    35   3888
## 36  BIN_1    36   3587
## 37  BIN_1    37   3348
## 38  BIN_1    38   2983
## 39  BIN_1    39   2833
## 40  BIN_1    40   2675
## 41  BIN_1    41   2407
## 42  BIN_1    42   2360
## 43  BIN_1    43   2095
## 44  BIN_1    44   2034
## 45  BIN_1    45   1884
## 46  BIN_1    46   1745
## 47  BIN_1    47   1683
## 48  BIN_1    48   1513
## 49  BIN_1    49   1474
## 50  BIN_1    50   1377
## 51  BIN_1    51   1278
## 52  BIN_1    52   1142
## 53  BIN_1    53   1185
## 54  BIN_1    54   1063
## 55  BIN_1    55   1023
## 56  BIN_1    56    985
## 57  BIN_1    57    849
## 58  BIN_1    58    821
## 59  BIN_1    59    803
## 60  BIN_1    60    791
## 61  BIN_1    61    704
## 62  BIN_1    62    675
## 63  BIN_1    63    606
## 64  BIN_1    64    613
## 65  BIN_1    65    588
## 66  BIN_1    66    595
## 67  BIN_1    67    522
## 68  BIN_1    68    505
## 69  BIN_1    69    509
## 70  BIN_1    70    487
## 71  BIN_1    71    440
## 72  BIN_1    72    461
## 73  BIN_1    73    396
## 74  BIN_1    74    377
## 75  BIN_1    75    392
## 76  BIN_1    76    355
## 77  BIN_1    77    325
## 78  BIN_1    78    322
## 79  BIN_1    79    343
## 80  BIN_1    80    317
## 81  BIN_1    81    297
## 82  BIN_1    82    263
## 83  BIN_1    83    256
## 84  BIN_1    84    226
## 85  BIN_1    85    249
## 86  BIN_1    86    232
## 87  BIN_1    87    228
## 88  BIN_1    88    216
## 89  BIN_1    89    191
## 90  BIN_1    90    209
## 91  BIN_1    91    205
## 92  BIN_1    92    190
## 93  BIN_1    93    185
## 94  BIN_1    94    159
## 95  BIN_1    95    161
## 96  BIN_1    96    172
## 97  BIN_1    97    153
## 98  BIN_1    98    133
## 99  BIN_1    99    143
## 100 BIN_1   100    127

BIN_2

variance_N_bin_2 <- variance_N(data = BI_fly, bin = bin_No[2])
variance_N_bin_2

##      Var1   N    Freq
## 1   BIN_2   1 1257772
## 2   BIN_2   2  660499
## 3   BIN_2   3  357558
## 4   BIN_2   4  202612
## 5   BIN_2   5  121474
## 6   BIN_2   6   76449
## 7   BIN_2   7   49570
## 8   BIN_2   8   33335
## 9   BIN_2   9   23141
## 10  BIN_2  10   16204
## 11  BIN_2  11   11654
## 12  BIN_2  12    8664
## 13  BIN_2  13    6502
## 14  BIN_2  14    5000
## 15  BIN_2  15    3704
## 16  BIN_2  16    3038
## 17  BIN_2  17    2320
## 18  BIN_2  18    1918
## 19  BIN_2  19    1494
## 20  BIN_2  20    1283
## 21  BIN_2  21     996
## 22  BIN_2  22     860
## 23  BIN_2  23     704
## 24  BIN_2  24     578
## 25  BIN_2  25     458
## 26  BIN_2  26     394
## 27  BIN_2  27     366
## 28  BIN_2  28     314
## 29  BIN_2  29     225
## 30  BIN_2  30     206
## 31  BIN_2  31     152
## 32  BIN_2  32     149
## 33  BIN_2  33     140
## 34  BIN_2  34     131
## 35  BIN_2  35     106
## 36  BIN_2  36      88
## 37  BIN_2  37      92
## 38  BIN_2  38      72
## 39  BIN_2  39      80
## 40  BIN_2  40      58
## 41  BIN_2  41      45
## 42  BIN_2  42      42
## 43  BIN_2  43      28
## 44  BIN_2  44      37
## 45  BIN_2  45      36
## 46  BIN_2  46      42
## 47  BIN_2  47      39
## 48  BIN_2  48      29
## 49  BIN_2  49      28
## 50  BIN_2  50      17
## 51  BIN_2  51      19
## 52  BIN_2  52      13
## 53  BIN_2  53      18
## 54  BIN_2  54      14
## 55  BIN_2  55      12
## 56  BIN_2  56      14
## 57  BIN_2  57      11
## 58  BIN_2  58       6
## 59  BIN_2  59       9
## 60  BIN_2  60      14
## 61  BIN_2  61      10
## 62  BIN_2  62      12
## 63  BIN_2  63      13
## 64  BIN_2  64       7
## 65  BIN_2  65       3
## 66  BIN_2  66       5
## 67  BIN_2  67       9
## 68  BIN_2  68       8
## 69  BIN_2  69       4
## 70  BIN_2  70       3
## 71  BIN_2  71       6
## 72  BIN_2  72       2
## 73  BIN_2  73       4
## 74  BIN_2  74       1
## 75  BIN_2  75       1
## 76  BIN_2  76       3
## 77  BIN_2  77       2
## 78  BIN_2  78       3
## 79  BIN_2  79       3
## 80  BIN_2  80       3
## 81  BIN_2  81       3
## 82  BIN_2  82       4
## 83  BIN_2  85       4
## 84  BIN_2  86       1
## 85  BIN_2  88       1
## 86  BIN_2  89       1
## 87  BIN_2  90       1
## 88  BIN_2  91       4
## 89  BIN_2  93       2
## 90  BIN_2  95       1
## 91  BIN_2  96       3
## 92  BIN_2  97       2
## 93  BIN_2  98       1
## 94  BIN_2  99       2
## 95  BIN_2 100       1
## 96  BIN_2 101       1
## 97  BIN_2 102       1
## 98  BIN_2 103       2
## 99  BIN_2 106       3
## 100 BIN_2 114       1

BIN_3

variance_N_bin_3 <- variance_N(data = BI_fly, bin = bin_No[3])
variance_N_bin_3

##     Var1   N    Freq
## 1  BIN_3   1 1731385
## 2  BIN_3   2  627390
## 3  BIN_3   3  251117
## 4  BIN_3   4  112801
## 5  BIN_3   5   54832
## 6  BIN_3   6   29205
## 7  BIN_3   7   16416
## 8  BIN_3   8    9670
## 9  BIN_3   9    5914
## 10 BIN_3  10    3747
## 11 BIN_3  11    2416
## 12 BIN_3  12    1643
## 13 BIN_3  13    1195
## 14 BIN_3  14     810
## 15 BIN_3  15     567
## 16 BIN_3  16     406
## 17 BIN_3  17     318
## 18 BIN_3  18     202
## 19 BIN_3  19     178
## 20 BIN_3  20     139
## 21 BIN_3  21     119
## 22 BIN_3  22      87
## 23 BIN_3  23      63
## 24 BIN_3  24      54
## 25 BIN_3  25      43
## 26 BIN_3  26      40
## 27 BIN_3  27      26
## 28 BIN_3  28      29
## 29 BIN_3  29      21
## 30 BIN_3  30      11
## 31 BIN_3  31      12
## 32 BIN_3  32      18
## 33 BIN_3  33       9
## 34 BIN_3  34      17
## 35 BIN_3  35       9
## 36 BIN_3  36      11
## 37 BIN_3  37       6
## 38 BIN_3  38       5
## 39 BIN_3  39       5
## 40 BIN_3  40       6
## 41 BIN_3  42       3
## 42 BIN_3  43       1
## 43 BIN_3  44       2
## 44 BIN_3  45       3
## 45 BIN_3  46       3
## 46 BIN_3  47       2
## 47 BIN_3  48       5
## 48 BIN_3  49       4
## 49 BIN_3  50       1
## 50 BIN_3  52       2
## 51 BIN_3  53       2
## 52 BIN_3  55       1
## 53 BIN_3  62       1
## 54 BIN_3  65       1
## 55 BIN_3  70       1
## 56 BIN_3  71       1
## 57 BIN_3  72       1
## 58 BIN_3  75       1
## 59 BIN_3  79       1
## 60 BIN_3  80       1
## 61 BIN_3  82       1
## 62 BIN_3  83       1
## 63 BIN_3  85       1
## 64 BIN_3  96       1
## 65 BIN_3  98       1
## 66 BIN_3 160       1
## 67 BIN_3 167       1
## 68 BIN_3 427       1

BIN_4

variance_N_bin_4 <- variance_N(data = BI_fly, bin = bin_No[4])
variance_N_bin_4

##     Var1   N    Freq
## 1  BIN_4   1 2081584
## 2  BIN_4   2  520863
## 3  BIN_4   3  153707
## 4  BIN_4   4   53607
## 5  BIN_4   5   21406
## 6  BIN_4   6    9501
## 7  BIN_4   7    4474
## 8  BIN_4   8    2359
## 9  BIN_4   9    1291
## 10 BIN_4  10     783
## 11 BIN_4  11     457
## 12 BIN_4  12     291
## 13 BIN_4  13     218
## 14 BIN_4  14     126
## 15 BIN_4  15      83
## 16 BIN_4  16      53
## 17 BIN_4  17      42
## 18 BIN_4  18      23
## 19 BIN_4  19      25
## 20 BIN_4  20      18
## 21 BIN_4  21      18
## 22 BIN_4  22       9
## 23 BIN_4  23      10
## 24 BIN_4  24       8
## 25 BIN_4  25      10
## 26 BIN_4  26       5
## 27 BIN_4  27       1
## 28 BIN_4  28       4
## 29 BIN_4  30       5
## 30 BIN_4  31       1
## 31 BIN_4  32       2
## 32 BIN_4  33       5
## 33 BIN_4  35       1
## 34 BIN_4  37       1
## 35 BIN_4  39       2
## 36 BIN_4  43       1
## 37 BIN_4  45       1
## 38 BIN_4  46       2
## 39 BIN_4  48       1
## 40 BIN_4  61       1
## 41 BIN_4  70       1
## 42 BIN_4  72       1
## 43 BIN_4 162       1

BIN_5

variance_N_bin_5 <- variance_N(data = BI_fly, bin = bin_No[5]) 
variance_N_bin_5

##     Var1   N    Freq
## 1  BIN_5   1 2321560
## 2  BIN_5   2  404334
## 3  BIN_5   3   87821
## 4  BIN_5   4   24002
## 5  BIN_5   5    7654
## 6  BIN_5   6    2919
## 7  BIN_5   7    1217
## 8  BIN_5   8     603
## 9  BIN_5   9     318
## 10 BIN_5  10     173
## 11 BIN_5  11      98
## 12 BIN_5  12      83
## 13 BIN_5  13      42
## 14 BIN_5  14      31
## 15 BIN_5  15      30
## 16 BIN_5  16      13
## 17 BIN_5  17      11
## 18 BIN_5  18      15
## 19 BIN_5  19       8
## 20 BIN_5  20       9
## 21 BIN_5  21       9
## 22 BIN_5  22       7
## 23 BIN_5  23       9
## 24 BIN_5  25       4
## 25 BIN_5  26       3
## 26 BIN_5  30       1
## 27 BIN_5  31       1
## 28 BIN_5  33       1
## 29 BIN_5  36       1
## 30 BIN_5  38       1
## 31 BIN_5  39       1
## 32 BIN_5  41       1
## 33 BIN_5  48       1
## 34 BIN_5  50       2
## 35 BIN_5  51       1
## 36 BIN_5  64       1
## 37 BIN_5  80       1
## 38 BIN_5  98       1
## 39 BIN_5 110       1
## 40 BIN_5 125       1
## 41 BIN_5 129       1
## 42 BIN_5 157       1
## 43 BIN_5 306       1
## 44 BIN_5 351       1

Create the model for 1st bin

#extract bin_No_1

bin_1 <- BI_fly$interactions[dist.bin==1,]
bin_1$dist <- abs(bin_1$dist)

message("Build model for bin: 1")

## Build model for bin: 1

models = model_bin(
     bin_1,
     subsample_size=1000
     )

## A truncated family of distributions from NBI has been generated 
##  and saved under the names:  
##  dNBI.0tr pNBI.0tr qNBI.0tr rNBI.0tr NBI.0tr 
## The type of truncation is left 
##  and the truncation parameter is 0  
## 18-01-2023 16:31:50
## Subsampling 1000/2851024 interactions for the null model regression...
## 18-01-2023 16:31:50
## Fitting with a maximum of 200 iterations...
## Using formula:
## N ~ log(abs(dist)) + b.trans_res  + p.trans_res + sqrt(b.length) + sqrt(p.length)
## GAMLSS-RS iteration 1: Global Deviance = 5206.992 
## GAMLSS-RS iteration 2: Global Deviance = 5176.838 
## GAMLSS-RS iteration 3: Global Deviance = 5169.337 
## GAMLSS-RS iteration 4: Global Deviance = 5167.785 
## GAMLSS-RS iteration 5: Global Deviance = 5167.566 
## GAMLSS-RS iteration 6: Global Deviance = 5167.539 
## 18-01-2023 16:31:51
## Converged: YES
## Iterations: 6
## 
## Coefficients:
## (Intercept)  8.67743046425092
##  log(abs(dist))  -0.725532500725227
##  b.trans_res 0.584344117331707
##  p.trans_res NA
##  sqrt(b.length)  0.000720678557023248
##  sqrt(p.length)  0.00727315178719011
## 18-01-2023 16:31:51
## Obtaining normalized randomized quantile residuals for the full dataset...

message("Finished!")

## Finished!

Session info

sessionInfo()

## R version 4.2.0 (2022-04-22)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Rocky Linux 8.7 (Green Obsidian)
## 
## Matrix products: default
## BLAS/LAPACK: /g/easybuild/x86_64/Rocky/8/haswell/software/FlexiBLAS/3.0.4-GCC-11.2.0/lib64/libflexiblas.so.3.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  splines   stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] forcats_0.5.1           stringr_1.4.0           dplyr_1.0.8            
##  [4] purrr_0.3.4             readr_2.1.2             tidyr_1.2.0            
##  [7] tibble_3.1.8            ggplot2_3.4.0           tidyverse_1.3.1        
## [10] peaky_0.0.0.9002        R2BGLiMS_0.1-26-03-2017 gamlss.tr_5.1-7        
## [13] gamlss_5.4-10           nlme_3.1-157            gamlss.data_6.0-2      
## [16] gamlss.dist_6.0-5       MASS_7.3-57            
## 
## loaded via a namespace (and not attached):
##  [1] lubridate_1.8.0   lattice_0.20-45   assertthat_0.2.1  digest_0.6.29    
##  [5] utf8_1.2.2        R6_2.5.1          cellranger_1.1.0  backports_1.4.1  
##  [9] reprex_2.0.1      evaluate_0.15     httr_1.4.2        pillar_1.8.1     
## [13] rlang_1.0.6       readxl_1.4.0      rstudioapi_0.13   data.table_1.14.6
## [17] Matrix_1.4-1      rmarkdown_2.14    munsell_0.5.0     broom_0.8.0      
## [21] compiler_4.2.0    modelr_0.1.8      xfun_0.30         pkgconfig_2.0.3  
## [25] htmltools_0.5.2   tidyselect_1.1.2  fansi_1.0.3       crayon_1.5.1     
## [29] withr_2.5.0       tzdb_0.3.0        dbplyr_2.1.1      grid_4.2.0       
## [33] jsonlite_1.8.0    gtable_0.3.1      lifecycle_1.0.3   DBI_1.1.2        
## [37] magrittr_2.0.3    scales_1.2.1      cli_3.6.0         stringi_1.7.6    
## [41] fs_1.5.2          xml2_1.3.3        ellipsis_0.3.2    generics_0.1.2   
## [45] vctrs_0.5.1       cowplot_1.1.1     tools_4.2.0       glue_1.6.2       
## [49] hms_1.1.1         fastmap_1.1.0     survival_3.3-1    yaml_2.3.5       
## [53] colorspace_2.0-3  rvest_1.0.2       knitr_1.38        haven_2.5.0
lovegrovec10 commented 2 months ago

Was this ever resolved?