PoMo errors: ascertainment bias, polymorphism, and scaffold numbering
adtw...@gmail.com
Hi all,
I’ve come up against a series of errors when doing a polymorphism aware analysis using SNP data. These come from low-coverage GBS SNPs (i.e. lots of missing data) called relative to a highly fragmented reference genome (1000s of scaffolds).
1. When using ascertainment bias correction (+ASC) it says this isn’t yet implemented in PoMo. Do you know if there are plans to do this (or if there’s a beta version with this implemented)?
2. If I then try it without ASC, I get this error “Provided alignment is exclusively used by PoMo but model string does not contain, e.g., "+P".” What is “+P”? I can’t find information on it anywhere.
3. When I then run it with “+P” it runs, but with an error “WARNING: Numerical underflow for lh-derivative”, it then runs for a long time producing the standard output, but then fails, with the following error:
WARNING: Numerical underflow for lh-derivative
ERROR: iqtree.cpp:3149: virtual pair<int, int> IQTree::optimizeNNI(bool): Assertion `curScore > appliedNNIs.at(0).newloglh - params->loglh_epsilon' failed.
ERROR: STACK TRACE FOR DEBUGGING:
ERROR:
ERROR: *** IQ-TREE CRASHES WITH SIGNAL ABORTED
ERROR: *** For bug report please send to developers:
ERROR: *** Log file: /localdisk/home/atwyford/Programs/cflib/scripts/Euphrasia_UK_maxmiss020318.cf.log
ERROR: *** Alignment files (if possible)
Aborted
I’ve attached the logfile (relevant part line 1255 onward) and the first 10 lines of the alignment. I note that the scaffold is missing in this file (reported as NA)—is this part of the problem? I formatted the data with cflib, but I might need to go back and check the format for scaffold numbering in the genome file.
In the meantime, I'll generate a sequence file with invariant sites, in the hope this gets round some of these issues.
Many thanks in advance for your help.
Alex
Dominik Schrempf
I apologize for not coming back to you earlier. Thank you very much for using
PoMo with your data set, which seems to be a perfect use case!
As you already observed, ascertainment bias correction is not implemented for
PoMo. At the moment, we do not plan to add this feature. One reason is that the
level of polymorphism (or heterozygosity) cannot be estimated when constant
sites are not present in the alignment. A possibility would be to let the user
define a heterozygosity (which is in fact already possible), if it is known and
constant sites are not available at all. You mentioned that you could prepare
data including constant sites, so I would opt for this path at the moment.
Second, the +P is the +polymorphism (PoMo) flag. We decided to force an explicit
statement of +P so that the used model is completely defined when looking at the
model string only (and not only by the type of data). The complete model string
will be printed in the .iqtree file when the run has completed. You can check
the possible models with `iqtree -h`. There is a section about POLYMORPHISM
AWARE MODELS. For instance, you might want to decrease N to 7 or 5, this also
greatly increases numerical stability.
Third, I observed that you started your runs from a previous checkpoint. This is
not recommended when the command line has changed (from GTR+ASC to GTR+P in your
case). Please use the `-redo` flag to force a complete rerun.
Finally, I believe that because of not including constant sites, the estimated
level of polymorphism (0.35 as far as I can see) is extremely large. That is why
the Eigendecomposition fails and the error is reported. I expect that this will
be not the case when you run IQ-TREE with PoMo on your full data, but I am very
happy to help in case you encounter problems.
Thanks again for your interest,
Dominik
> You received this message because you are subscribed to the Google Groups
> "IQ-TREE" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to iqtree+un...@googlegroups.com.
> To post to this group, send email to iqt...@googlegroups.com.
> Visit this group at https://groups.google.com/group/iqtree.
> For more options, visit https://groups.google.com/d/optout.
adtw...@gmail.com
Dominik Schrempf
I am glad that it worked! You are right. You can compare the "PoMo branch
lengths" to branch lengths of standard substitution models. The conversion
factor is 1/(N**2), in words, one over N squared (if you did not specify N, it
will be 9). I.e., PoMo branch lengths are longer because they also include
frequency shifts (i.e., genetic drift, and not only substitutions).
However, it is expected that the converted branch lengths will differ (because
PoMo accounts for polymorphisms and thereby, e.g., terminal branch lengths are
expected to be much shorter).
I added some code to the `pomo_latest` branch, so that trees with both types of
branch lengths are printed. At the moment, you need to compile IQ-TREE yourself,
if you want to use this code. I am also not sure, if you want to rerun your
analysis :), so a manual conversion might be better for you.
You can find the derivation of the conversion factor in our publication, in
section 2.8.
Schrempf, D., Minh, B. Q. B. Q., De Maio, N., von Haeseler, A., & Kosiol, C.
(2016). Reversible polymorphism-aware phylogenetic models and their application
to tree inference. Journal of Theoretical Biology, 407, 362–370.
http://doi.org/10.1016/j.jtbi.2016.07.042
I hope that I could help!
Kind regards,
Dominik