Is the parametric bootstrap using simulations mandatory?

[XieHongX]

Hi,
I have been using gIMble to estimate demographic parameters between two populations and then detect possible gene flow barriers. Everything went OK so far and I detected some barrier regions which I think is reasonable according to functional annotations. But I got stuck int the parametric bootstrap step. I don't have a detailed recombination map for this species so I did the simulations using a fixed recombination rate. Firstly, the estimated confidence intervals of the global parameters from parametric bootstraps are close to the parameters from the real data, but they just don't overlap. Secondly, the FPR rate of genomic barriers calculated from the simulations is high, which is not the case in the original article of gIMble. I am confused and don't know what to do. Are the unreasonable results I got from simulations caused by a lack of recombination map, or just some settings of the parameters in my analysis are inappropriate (I am been following the settings in the original gIMble paper with adjustments based on my target species)? In real application of gIMble, is it a must-do to do parametric bootstraps based on a recombination map? I really need some advice from you guys. Thanks!

[GertjanBisschop]

Hi @XieHongX,
As you described yourself, running a complete gIMble analysis requires some well-judged choices. It is quite challenging from my end to determine what might be driving the lack of distribution overlap between the bootstrap and the real data. In case you suspect that recombination might be driving this, you might run the bootstrap without recombination, as well as with an increased recombination rate. This would allow you to gauge the impact of recombination.

[XieHongX]

Hi, indeed I tried both no recombination and a fixed recombination rate. Assuming no recombination resulted in a slightly larger migration rate Me compared to results from real data, while a fixed recombination rate of 1.19 (global estimated recombination rate) resulted in a smaller Me. Maybe using a recombination rate between 0 and 1 will make the CI overlap with the global parameters, but is doing so reasonable?

[GertjanBisschop]

Hmm, I wasn't really proposing to tweak the recombination rate until you get more overlap. I was mainly suggesting that you can use the simulation framework to gauge whether recombination is indeed driving the lack of fit between the optimal model and the bootstrap simulations.

[XieHongX]

Yes， I think recombination rate did matter. Simulating in different mutation rates changes the results. I am thinking maybe a fixed mutation rate does not help much in parametric bootstrap then. But I am not sure about this. The main concern I have is that a great proportion of detected barriers would be rejected based on simulated FPR. Though I already applied a strict definition of barriers

[GertjanBisschop]

This seems like great news. It means the bootstrap helped to quantify the confidence you have in the identified barriers. I don't think you should be meddling with the parameters until you find more barriers.

[XieHongX]

Thank you for the discussion. Then, In this way do you think I should only accept those barrier regions that passed FPR? I identified 1.36% of genomic windows as barriers under strict barrier definition. Then only about 1/5 of these barriers would pass the FPR<0.05 using a fixed recombination rate. Do you think the result is useful?

[KLohse]

Hi @XieHongX,

It's great to hear that you managed to run gIMble on your data successfully! The likelihood calculation implemented in gIMble assumes no recombination within blocks. Recombination affects both the bias and accuracy of estimates. So when running a bootstrap on parameter estimates for data simulated with recombination, you do indeed expect a bias in the direction you mention (e.g. m_e is underestimated). To what extent that bias is acceptable depends on your question and the bias/accuracy trade-off in your data. Given that the global m_e is underestimated, it might make sense to adjust your bootstrap scans for this bias, i.e. when you compute DeltaB for the simulated data it probably makes sense to set the reference m_e (under the null model of no m_e variation) to be value in the m_e grid that the best fits the simulated (rather than the real) data.

[XieHongX]

Hi thank you for your insights.
I did find the estimated m_e of simulated data is lower than the real data. I am just not understanding why the FPR calculated from simulated data is so high for most detected barrier regions. I will double check my codes for calculating FPR and the codes for bootstrap scans. Do you think a detailed recombination map will make big difference compared to a fixed rate?