search

wangke16 / MSMC-IM

A tool for estimating time-dependent migration rates based on cross-/within-population coalescent rates from MSMC
16 stars 2 forks source link

Differences in split time estimates between MSMC2 and MSMC-IM #12

Open rancilhac opened 7 months ago

rancilhac commented 7 months ago

Hello, I am using MSMC2 and MSMC-IM to estimate the timing of splits across four lineages. I have noticed that the split times from both approaches (time at which the CCR crosses the 0.5 threshold in MSMC2 and time at which M(t) reaches a plateau around 1 in MSMC-IM) are always different, with the latter being more recent. Here is an example:

MSMC_plots_1

With MSMC2, the split time is estimated at 414,867ya (left plot). However, MSMC-IM (center and right plots) suggests a split 50,000ya. I am wondering how to interpret this difference. The MSMC2 estimates are also in contradiction with the order of the splits according to phylogenetic analyses, hence I am not sure how much they can be trusted.

Another odd result that I got is that in analyses involving a lineage with a large Ne (>1.5e+6), the split time inferred by MSMC-IM is extremely recent. Here is an example, where MSMC2 estimates a split at 835,209ya and MSMC-IM around 2,000ya:

plots_msmc_2

Here are examples of commands I used to run MSMC2: msmc2 -t 32 -i 20 -s -p 1*2+13*1+1*2+1*3 -I 0-4,0-5,0-6,0-7,1-4,1-5,1-6,1-7,2-4,2-5,2-6,2-7,3-4,3-5,3-6,3-7 -o ./pusillus_extoni_Run1_allopatrics /path/to/multihetsep/*_MSMC_input.txt and MSMC-IM: python3 MSMC_IM.py -mu 4.007e-9 -o $(pwd)/pusillus_extoni_MSMC-IM -N1 87789 -N2 214368 -p 1*2+13*1+1*2+1*3 -b 1e-8,1e-6 --printfittingdetails --plotfittingdetails --xlog –ylog ./pusillus_extoni_Run1_allopatrics_combined.msmc.final.txt

For MSMC-IM, I have used the Ne estimates for the most recent time segment in the MSMC2 results.

I would be thankful for any comment that could help me understand these results. Best, Loïs

wangke16 commented 5 months ago

hi Loïs, in both cases, the rCCR plot and m(t) results seem to suggest gene flow in recent times (i.e. the first peak in rCCR/mt(t)). depending on the extensity of gene flow, it likely reflects either post-split gene flow event or population split event. To differentiate two cases, it might be useful to check the mt(t) plot with M(t) quantile plotted in shades (similar to Fig2 in MSMC-IM paper).

best, Ke