trying to replicate Fig2 of Schiffels & Durbin (2014) about msmc

[stschiff]

Hi Andrew, it seems you’re running mdmc on only 10Mb of sequence. It is designed to run on the entire human genome, so 200-300 times more data. Also, running for 200 iterations is unnecessary. I have never tried this many and would actually be worried that some artefacts might lead it away from the likelihood maximum (has to do with the sparseness E-M algorithm that we use). So I would suggest to use the default number of iterations, and to run on a lot more data. You can use msprime to simulate 20*100Mb long chromosomes and run the same command line on all of them. Also, I recommend to use msmc2 now (www.github.com/stschiff/msmc2 <http://www.github.com/stschiff/msmc2>). Stephan

…

On 11 Apr 2019, at 20:56, Andrew Kern ***@***.***> wrote: Hey there- I'm trying to replicate Figure 2 from your original paper using the zig-zag simulation that you specify in the supplement of that paper. The particular "problem" that I'm getting is that population size estimates from recent time look way too large. I imagine I'm not performing the inference procedure correctly, but here are my simple steps: #create simulation data, convert to msmc input using msmc-tools ms 4 1 -t 7156.0000000 -r 2000.0000 10000000 -eN 0 5 -eG 0.000582262 1318.18 -eG 0.00232905 -329.546 -eG 0.00931619 82.3865 -eG 0.0372648 -20.5966 -eG 0.149059 5.14916 -eN 0.596236 0.5 -T | ./ms2multihetsep.py 1 10000000 > test.txt #run HMM for 200 iterations msmc --fixedRecombination -o my_msmc_output msmc-tools/test.txt -i 200 I then have plotted the results using the step plot tools from plot_utils.py and get the following result <https://user-images.githubusercontent.com/7674197/55983551-3f42eb00-5c50-11e9-9934-52c754b652b9.png> Clearly I'm doing something wrong here as the most recent population size estimates are orders of magnitude too large. Could you point me in the right direction? thanks, Andy — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#44>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAbQmn9S1qZvAmrj1mluRcG-9UZJYwmYks5vf4VhgaJpZM4cqoxU>.

[andrewkern]

Thanks-that makes sense.

I've gone ahead and done a simulation of 20 100Mb chromosomes using msprime under your Zig-Zag model.

#simulate 20 chromosomes, each 100Mb long
for i in `seq 1 20`; do mspms 4 1 -t 71560.0000000 -r 20000.0000 100000000 -eN 0 5 -eG 0.000582262 1318.18 -eG 0.00232905 -329.546 -eG 0.00931619 82.3865 -eG 0.0372648 -20.5966 -eG 0.149059 5.14916 -eN 0.596236 0.5 -T | ./ms2multihetsep.py 1 100000000 > test_${i}.txt & done

#run msmc on all outputs
msmc --fixedRecombination -o my_msmc_output test_1.txt test_2.txt test_3.txt test_4.txt test_5.txt test_6.txt test_7.txt test_8.txt test_9.txt test_10.txt test_11.txt test_12.txt test_13.txt test_14.txt test_15.txt test_16.txt test_17.txt test_18.txt test_19.txt test_20.txt

Here is the resulting figure

It looks a bit better, but certainly not very similar to your figure. Reading your supplement I was under the assumption that the Zig-Zag model as listed was what you had run to generate Figure 2. Perhaps that is not the case? If not can you please share the correct parameters for those simulations?

thanks again,
Andy

[stschiff]

Hi Andy, hmm, that’s embarrassing, perhaps the Schiffels and Durbin paper Supplementary information on that simulation isn’t correct then. We have just recently rerun this on simulations with msprime and it look totally fine (see our new preprint https://www.biorxiv.org/content/10.1101/585265v1 <https://www.biorxiv.org/content/10.1101/585265v1>). Let me get back to you after having talked to my student Ke. Could you email me directly at [email protected] <mailto:[email protected]>? I’d like to share the simulation script directly for your to test. Happy to post it also in this github thread afterwards. Thanks, Stephan

…

On 11 Apr 2019, at 23:25, Andrew Kern ***@***.***> wrote: Thanks-that makes sense. I've gone ahead and done a simulation of 20 100Mb chromosomes using msprime under your Zig-Zag model. #simulate 20 chromosomes, each 100Mb long for i in `seq 1 20`; do mspms 4 1 -t 71560.0000000 -r 20000.0000 100000000 -eN 0 5 -eG 0.000582262 1318.18 -eG 0.00232905 -329.546 -eG 0.00931619 82.3865 -eG 0.0372648 -20.5966 -eG 0.149059 5.14916 -eN 0.596236 0.5 -T | ./ms2multihetsep.py 1 100000000 > test_${i}.txt & done #run msmc on all outputs msmc --fixedRecombination -o my_msmc_output test_1.txt test_2.txt test_3.txt test_4.txt test_5.txt test_6.txt test_7.txt test_8.txt test_9.txt test_10.txt test_11.txt test_12.txt test_13.txt test_14.txt test_15.txt test_16.txt test_17.txt test_18.txt test_19.txt test_20.txt Here is the resulting figure <https://user-images.githubusercontent.com/7674197/55993950-40324780-5c65-11e9-877e-7ed5c007284a.png> It looks a bit better, but certainly not very similar to your figure. Reading your supplement I was under the assumption that the Zig-Zag model as listed was what you had run to generate Figure 2. Perhaps that is not the case? If not can you please share the correct parameters for those simulations? thanks again, Andy — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#44 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAbQmmJbUT_tdnjjGJ5o6KTHFdQSKt_Pks5vf6higaJpZM4cqoxU>.

[andrewkern]

Sure- will email now.

[andrewkern]

A bit more follow-up. I was concerned that perhaps msmc didn't like the chromosome names that I was handing it, so I reran simulations in the following way and got similar results

for i in `seq 1 20`; do mspms 4 1 -t 71560.0000000 -r 20000.0000 100000000 -eN 0 5 -eG 0.000582262 1318.18 -eG 0.00232905 -329.546 -eG 0.00931619 82.3865 -eG 0.0372648 -20.5966 -eG 0.149059 5.14916 -eN 0.596236 0.5 -T | ./ms2multihetsep.py ${i} 100000000 > test_${i}.txt & done

msmc --fixedRecombination -o my_msmc_output test_1.txt test_2.txt test_3.txt test_4.txt test_5.txt test_6.txt test_7.txt test_8.txt test_9.txt test_10.txt test_11.txt test_12.txt test_13.txt test_14.txt test_15.txt test_16.txt test_17.txt test_18.txt test_19.txt test_20.txt

results below

[stschiff]

continued this discussion offline.