ms to msprime code conversion

[David-Peede]

Hello All!

I have published ms code compiled from multiple studies and I would like to convert the ms code to msprime format. For simple models, this is relatively straightforward and simple to do by hand, but for more complex models I fear that I may make a mistake in the conversion process. I know there is conversion functionality from msprime pre-version 1.0, but I was wondering if this functionality exists for ms to msprime conversion. If this functionality does not exist what is recommended as "best practices" for handling both ms and msprime models. In an ideal world, I would be able to do everything natively in msprime but it seems like I might have to use mspms plus Seq-Gen.

Thanks in advance for any and all help!
-Dave

[grahamgower]

Hi @David-Peede.

Demes-python has an ms -> Demes converter: see the demes.from_ms() function. The demes model can then be imported to msprime using msprime.Demography.from_demes().

import demes
import msprime

cmd = (
    "-I 3 1 1 1 "
    "-n 1 1.682020 "
    "-n 2 3.736830 "
    "-n 3 7.292050 "
    "-eg 0 2 116.010723 "
    "-eg 0 3 160.246047 "
    "-ma x 0.881098 0.561966 0.881098 x 2.797460 0.561966 2.797460 x "
    "-ej 0.028985 3 2 "
    "-en 0.028985 2 0.287184 "
    "-ema 0.028985 3 x 7.293140 x 7.293140 x x x x x "
    "-ej 0.197963 2 1 "
    "-en 0.303501 1 1"
)
graph = demes.from_ms(cmd, N0=7310, deme_names=["YRI", "CEU", "CHB"])
demography = msprime.Demography.from_demes(graph)

We also have a CLI for this (and a Demes -> ms converter), but these are not yet in a released version of demes-python. Feel free to grab the latest git code if you're curious: https://github.com/popsim-consortium/demes-python/

[David-Peede]

Thanks, @grahamgower I will take a look!

-Dave

[David-Peede]

Actually, on second thought are there instances where it might make more sense to run simulations in ms vs msprime? For example, here is an example of the description of a "simple" published model that I am looking at:

Each simulation consisted of 10,000 non-recombining loci with three ingroup taxa and an outgroup (P4), and the species tree topology (((P1, P2), P3), P4). P1 and P2 split at a time of 1.2N generations; P3 from P1/P2 at 2.4N generations; and the outgroup at 16N generations. Introgression was specified at one of two times (0.2N generations or 0.04N generations), and in one of two directions (P2 → P3 or P3 → P2).

and here an example of the actual ms code used:

ms 4 10000 -T -I 4 1 1 1 1 -ej 4.0 2 1 -ej 0.6 3 2 -ej 0.3 4 3 -es 0.05 3 0.95 -ej 0.05 5 2

However, the N0 value here is 1e-6 which—from my limited past experience simulating—will cause long run times in msprime, however, if I use ms I will be losing a whole host of resources. With that being said the majority of my analyses will require one sample per species/population, are not terribly complex, and I don't necessarily need tree sequences to conduct my analyses—even though I do love the tree sequence ecosystem! Could this be a situation where it might be simpler and faster—but with larger output file sizes—to run the simulations in ms? I have always assumed that msprime is the "gold standard," but I still see people use ms so there must a reason why, right?

Anyways, this is a very long-winded way to ask when and why may someone might use ms over msprime and I apologize upfront if this is a rather silly or simple question that has been answered before—and I greatly appreciate any guidance!

-Dave

[molpopgen]

People use ms mostly out of inertia/historical reasons. That code base is probably longer maintained (Dick retired several years ago), but I haven't asked directly about this. There may be some cases where ms runs faster. We should look at those in msprime with a profiler and see what can be done.

My impression is that most folks still using ms "have always" used it, and change is hard for some. There's probably many PIs still recommending it to their students, etc., but that's because it is hard to replace established tools.

[grahamgower]

There are instances where ms will be faster. E.g. for really small simulations. If ms is faster in other cases, that could indicate a performance bug in msprime and as @molpopgen, it would be good to report this and see if it can be fixed. In your case, I don't think you'll have any problem with speed using msprime.

the N0 value here is 1e-6

I guess that should be 1e6?

time_units: generations
demes:
- name: deme1
  epochs:
  - {end_time: 0, start_size: 1000000.0}
- name: deme2
  start_time: 16000000.0
  ancestors: [deme1]
  epochs:
  - {end_time: 0, start_size: 1000000.0}
- name: deme3
  start_time: 2400000.0
  ancestors: [deme2]
  epochs:
  - {end_time: 0, start_size: 1000000.0}
- name: deme4
  start_time: 1200000.0
  ancestors: [deme3]
  epochs:
  - {end_time: 0, start_size: 1000000.0}
pulses:
- source: deme2
  dest: deme3
  time: 200000.0
  proportion: 0.05

$ time msp ancestry --demography /tmp/a.yaml -L 10000 -k 1 -o /tmp/a.ts deme1:1 deme2:1 deme3:1 deme4:1

real    0m0,315s
user    0m0,418s
sys     0m0,778s

[David-Peede]

Thank you very much @molpopgen and @grahamgower for the very thorough and clear explanations! It seems like my memory and intuition served me wrong about the run times, but I will be sure to report instances where I find ms to be faster in the future!

Thanks again, this was very helpful!
-Dave