A-J-F-Mackintosh
commented
2 days ago Hi Rowan,
Your question was likely directed to @DRL and @GertjanBisschop , but I will post my answer given that I have been in a similar situation.
It is true that increasing the kmax value would generate a more informative bSFS. However, there should be more than enough information with kmax=2,2,2,2 to fit the IM model (as well as the simpler nested models). Setting higher values is unlikely to lead to very different results (most blocks should be unaffected by the kmax value), but it would increase the computation time substantially. For example, I would expect kmax=4,4,4,4 to be prohibitively slow and introduce more error (blocks with >10 SNPs are likely to be mapping artefacts) than useful information. So, I don't think you need to change this at the moment.
The optimisation behaviour you describe suggests that your bSFS cannot be well described by an IM model with sensible parameter values. There could be a few different reasons for this:
- It could be that the true history is not IM-like (ghost introgression? intra-pop structure?) and so the observed bSFS cannot be well approximated by the models available in gIMble.
- It could be because your blocks are too big. Large blocks with lots of hidden recombination events would lead to a 'underdispersed' bSFS that is hard to fit a sensible model to.
- The divergence time could be too recent. If the split is << Ne generations ago then the bSFS contains little information about the Ne values of the two contemporary populations, and the estimates become biased due to hidden recombination. In particular, intra-block recombination leads to a deficit of monomorphic blocks, and this deficit would translate to massive NA and NB estimates when T is very small.
I encountered the problem of too recent a split time when analysing a pair of subspecies with T ~ 0.1 Ne. There is no real solution other than to sync the three Ne values, or alternatively leave them free but interpret the values with caution.
My suggestion would be to start by reducing the block size (unless it is already set to a small value). There is no standard way to determine the right block size, but when rates of mutation and recombination are similar it is best to set the block size so that the number of segregating sites per-block is something like 1, 1.5 or 2. If you reduce the block size and your parameter estimates change dramatically then this suggests that the problem is (at least partially) due to hidden recombination.
Cheers,
Alex
RSchley
KLohse
First of all, thanks so much for developing this great method.
I am working on an extremely rapid tropical tree radiation, and am aiming to infer barriers to gene flow between several co-occurring species pairs, using whole-genome resequencing data at 20x. These trees are expected to display rampant incomplete lineage sorting because of their large Nc and the fact that they are outcrossing, but they do have relatively short generation times considering they are trees.
I have gone through the gIMble pipeline, but the best-selected models that gIMble optimize returns differs depending on which parameter ranges I specify. Specifically, my initial prior ranges often result in the parameter estimate hitting the upper or lower prior bounds (often making no biological sense), and when I relax or change these priors, the best selected model changes.
I was wondering if this is because I am using the 2,2,2,2 kmax mutation configuration in gIMble tally - in the gIMble paper, it says "... the choice of --kmax involves a trade-off between the information gained by distinguishing the probabilities of rare bSFS configurations and the increase in computational complexity."
If I increased this to, for example, 4,4,4,4, would I capture more information in my sequence data, meaning that I have more power to distinguish different models and make more reliable parameter estimates?
Or is there anything else you might suggest to help me make the most robust inferences I can?
Thanks a lot! Rowan