Core dumped

[jshleap]

Hi! I am trying to run accuMUlate in my 64 samples, but it is core-sumping with the following error:

$ accuMUlate -c params.ini -b CAN.bam -r ../reference_clean.fasta > putative_mutations.tsv
terminate called after throwing an instance of 'boost::exception_detail::clone_impl<boost::exception_detail::error_info_injector<boost::program_options::required_option> >'
  what():  the option '--mu' is required but missing
Aborted (core dumped)

The CAN.bam was prepared following your instructions, as well as the atached params.init.
Am I doing something wrong?
params.txt

PS: The ancestor has not been sequenced, so is only in label

[dwinter]

Hi @jshleap , the error message suggests the parser is not finding a value for "mu", the approximate mutation rate used for the mutation-calling.

Is one being passed from the command line (there isn't one in the params file)

David

[jshleap]

Thanks, I assume it has to do something about that, however, it is not clear in the docs (at least not to me). Also, could you point me to some sources on how do you use this mu as a prior? I am assuming is a starting point of some sort since one of the goal of accumulate is to estimate the rate, am I wrong? I cant find the estimation method in your paper or references (at least not one that require a mu) but i am sure I am missing something important. How much bias this approximate mu has in the estimation of the mu by accumulate and denominate?
Thanks!

[jshleap]

Hi! I tried to play around with different values of theta, phi, and seq-error, however, it always throws the following dump:

terminate called after throwing an instance of 'std::out_of_range'
what(): basic_string::substr: __pos (which is 18446744073709551615) > this->size() (which is 52)
Aborted (core dumped)

The mu has previously been described for the species. I would like to know if you have a systematic way to estimate the overdispersion (phi) and theta (which I assume is the measure of genetic diversity?)

I am attaching thei params.ini I used last.

Best!

[dwinter]

Hi @jshleap,

Sorry about the delay in getting to you, we tend to have a long summer/end of year holiday down here in NZ.

This error message suggests an integer is overflowing at some stage, but there doesn't seem to be any issue with the params file you've included. If you can send me a subsample of your aligned reads I might be able to make some progress in debugging it. (something like

samtools view -h [large.bam] | head -n [nrows, incl header] | samtools view -Sb > sample.bam

Theta is expected heterozygosity. It would be possible to estimate this using the framework in accumulate but it's not easy to do at the moment. You really just want some _a prior _ estimate of how many hets you will see (maybe a very low number of a lab-cultured organism?)

Do get a number of phi will likely require some exploratory analysis. You can get an initial call-set with a value like 0.01 (which I found worked for illuimina sequencing on the Arabidopsis genome). You can then re-run accumulate on just the putative mutations while changing the value of phi. You can do this informally (and see how dependant on this value the set of called mutations are) or mroe formally keeping track of the over-all likelihood of the run from each value of phi (the tenth coloum in the output)

[jshleap]

Hi @dwinter , thanks for the reply, and don't worry about the delay. I am still not able to run accumulate on my data. I am attaching a 200 line bam file, and the current param file. I used angsd to compute the average theta, which I have included in the param file.

I hope this can help you debug, or point to my error. Let me know if you need anything else.

Best,

[jshleap]

Hi @dwinter, re-estimated the expected heterozygocity (which is higher of what I sent before), but I stull get the same error:

terminate called after throwing an instance of 'std::out_of_range'
  what():  basic_string::substr: __pos (which is 18446744073709551615) > this->size() (which is 52)
Aborted (core dumped)

Any ideas?