JiscaH / sequoia

R package for pedigree inference based on SNP data
25 stars 6 forks source link

Error: Sibship number out of bounds #19

Closed OliverPStuart closed 1 year ago

OliverPStuart commented 1 year ago

Hi Jisca,

I'm taking your advice and simulating pedigrees to test the use of sequoia for my intended use case. I'm running many many replicates across some parameters of interest (depth of simulated pedigree, probability of a female choosing to reproduce asexually, tolerance for sib-mating) and I've tested my pedigree simulation function quite well.

I have encountered the following error while trying to compile a markdown file for some results:

Quitting from lines 2031-2148 (analysis.Rmd) 
Error in SeqParSib(ParSib = "sib", FortPARAM, GenoM, LhIN = LifeHistData,  : 
    ERROR! *** setParTmp: Sibship number out of bounds *** Please report this bug.
Calls: <Anonymous> ... suppressMessages -> withCallingHandlers -> sequoia -> SeqParSib
Execution halted

This does not happen every time. In this case, it happened on the 41st replicate scenario. Unfortunately because it occurred while compiling a markdown file, I cannot retrieve the genotype matrix and pedigree than produced the error.

I am not at all sure what this means but it does instuct me to report the bug, so here I am!

Best, Ollie

JiscaH commented 1 year ago

Hi Oliver,

Thanks for the bug report! The ones that only give problems for a few specific parameter combinations, and then only in a small proportion of replicates, are by far the hardest to find & fix.

I'll go over the code to see if I can spot the problem, but since this error is generated by a low-level function that's called nearly a hundred times in the Fortran code, I am not 100% sure that that will work. Could you perhaps email me the parameters / input / markdown chunk you used? I.e. enough info so that I can run it a few hundred times until it hits the error?

Cheers, Jisca

OliverPStuart commented 1 year ago

Hi Jisca,

I reran my code and modified to print out the genotype matrix, life history table and original pedigree for that iteration of the loop. Here is a set of inputs that produces the same error: temp_life_hist.txt temp_test_geno.txt temp_test_ped.txt

The pedigree file is randomly generated with a certain amount of asexual reproduction. As such, there are individuals with only one parent. The founder's genotypes are from real individuals, then I take those and make new offspring genotypes based on the pedigree. The pedigree file has some additional information that is used in making new genotypes. In this case, there are quite a few asexually reproduced individuals in the pedigree and there are 7 generations post-founding. However, the first time I encountered this error, it was a shallow pedigree with only one or two asexually reproduced individuals.

Here is the command that produces the error from these inputs.

sequoia_out <- sequoia(GenoM=test_ped_genotypes,
                                          LifeHistData=life_hist_sim,
                                          Module="ped",
                                          Tfilter = -5,
                                          quiet=T,
                                          args.AP=list(Discrete = TRUE))

Best, Ollie

Edit: here's another set of input files that consistenly produce the error temp_test_geno2.txt temp_life_hist2.txt temp_test_ped2.txt

JiscaH commented 1 year ago

Hi Oliver,

I've put an updated version of the package on here, but it is not a complete fix.

I have been getting the same error occasionally during the simulations I always run before uploading a new version to CRAN. It's the same issue as you had: I got an error but by default don't save the simulated genotypes; I just ran the same parameter combinations 200 times without getting the error again.

I'll let you know when I have a (hopefully) definite complete fix.

Cheers, Jisca

OliverPStuart commented 1 year ago

Thanks for that Jisca.

If it's any help, it appears to be a deterministic problem with specific simulated data. That is to say, re-running the program with the same input data always produces the error for me.

Is the installation method for the recently updated version any different to what is listed on the README file? I am trying remotes::install_github("JiscaH/sequoia", ref="stable") but it is not working:

Error in utils::download.file(url, path, method = method, quiet = quiet,  : 
  cannot open URL 'https://api.github.com/repos/JiscaH/sequoia/tarball/stable'

Best, Ollie

JiscaH commented 1 year ago

Hi Oliver,

Installation without 'ref="stable" ' should work!

And sequoia is fully deterministic - same input always gives same output. But with some bugs (esp. subscript out of array bounds), whether or not it gives an error may very rarely depend on the data you ran before, unless you use special debugging tools.

JiscaH commented 1 year ago

Hi Oliver, I found the issue, and it's definitely fixed now! Cheers, Jisca