convert_msmc_to_demes.py

#!/usr/bin/env python3

import argparse
import plot_utils


parser = argparse.ArgumentParser(description="Convert MSMC final result files to a Demes model (see https://popsim-consortium.github.io/demes-spec-docs)")
parser.add_argument("msmc_final_file", metavar="FILE", type=str, help="final output file from MSMC or MSMC2 run.")
parser.add_argument("mu", metavar="MUT", type=float, help="per-generation mutation rate")

args = parser.parse_args()

msmc_obj = plot_utils.MSMCresult(args.msmc_final_file)

n = len(msmc_obj.lambdas)

print("# This Demes model was generated by convert_msmc_to_demes.py from https://github.com/stschiff/msmc-tools")
print("description:\n  Piecewise-constant population size model as inferred by MSMC/MSMC2")
print("time_units: generations")
print("demes:\n  - name: generic\n    epochs:")
for t, la in zip(msmc_obj.times_left[::-1], msmc_obj.lambdas[0][::-1]):
    t_scaled = round(t / args.mu)
    popsize = (1.0 / la) / (2.0 * args.mu)
    print(f"    - {{end_time: {t_scaled}, start_size: {popsize}}}")