lh3 / psmc

Implementation of the Pairwise Sequentially Markovian Coalescent (PSMC) model
Other
146 stars 60 forks source link

Bootstrapping running into erorr 'Use of uninitialized value in division (/)' #43

Open elizeng opened 2 years ago

elizeng commented 2 years ago

Hi all,

I have successfully ran my psmc for a few species, the regular run works out fine with no errors. and now I am running the bootstrapping segment of it but some individuals in some species are running into this error: Use of uninitialized value in division (/) at psmc_plot.pl line 125, <> line 292900.

For species A, running with -N25 -t15 -r5 -p "4+25*2+4+6" did not incur the error, whilst when running it with -N30 -t5 -r5 -p "4+30*2+4+6+10" causes the error to pop up in some but not all samples.

And through looking at the .psmc output, I note that its really just down to just one or two iterations where the RI value is stated as nan, and the QD is excessively huge; an example:

QD -1000000000000000052504760255204420248704468581108159154915854115511802457988908195786371375080447864043704443832883878176942523235360430575644792184786706982848387200926575803737830233794788090059368953234970799945081119038967640880074652742780142494579258788820056842838115669472196386865459400540160.000000 -> -1000000000000000052504760255204420248704468581108159154915854115511802457988908195786371375080447864043704443832883878176942523235360430575644792184786706982848387200926575803737830233794788090059368953234970799945081119038967640880074652742780142494579258788820056842838115669472196386865459400540160.000000

I don't quite get what QD is from the description Q-before-opt Q-after-opt, but its something that has been generated whilst the program runs...

If anyone has any insights to this problem, I would be glad to get some advice.

Thank you!

fbhernandez commented 2 years ago

Hi Elize, in a previous comment you mentioned the new pipeline to get the fastq consensus file. How did you do it? with bcftools? I will appreciate any ideas on how to make it work. Thank you

elizeng commented 2 years ago

Hi Elize, in a previous comment you mentioned the new pipeline to get the fastq consensus file. How did you do it? with bcftools? I will appreciate any ideas on how to make it work. Thank you

Its not possible to do it with the new pipeline, doesn't allow for the output of a fq file, so I just used the pipeline that has been suggested by PSMC. For the old pipeline to work, you'd have to find a version with vcfutils, which I found to be in bcftools1.1.

Good luck!