ts: Add a runMacsTs() version of runMacs()

[gregorgorjanc]

Current runMacs() does:

1. Simulates local trees
2. Simulates mutation events
3. Drops mutations down the local trees to generate haplotypes (this is a memory bottleneck with the currently-implemented high mutation rates ~1e-8)
4. Downsamples sites so we don't ran out of memory in downstream work with AlphaSimR

With @LynxJinyangii and @bryo-han wediscussed today that we would like to:

• Avoid the memory bottleneck in generating haplotypes (and possibly also downsampling sites)
  • To address the above we could just implement lower mutation rate in the current species models
  • And we set it via a new mutRate argument to runMacs() so users can manipulate it as they need
• runMacs() to return haplotypes - as MapPop, so that we can continue with the current standard AlphaSimR genome workflow: founderGenomes <- runMacs(...); SP <- SimParam$new(founderGenomes)
• runMacs() to possibly also return tree sequence - as RcppTskit::TreeSequence and we then do ts <- runMacs(...); founderGenomes <- asMapPop(ts, segSites=n); SP <- SimParam$new(founderGenomes)
  • convert MaCS local trees into RcppTskit::TreeSequence (is this easy or very hard?)
  • add returnTreeSeq argument to runMacs() to control whether we return RcppTskit::TreeSequence or MapPop
• asMapPop(ts, segSites=n) could be used to downsample sites if needed
• get genetic and physical maps from MaCS so we have both for downstream work

Are we missing anything in the above summary @LynxJinyangii?

[LynxJinyangii]

It looks good to me. A few minor points: 1. Should we first create a parallel function runMacsTs()? 2. I was thinking that if we eventually implement the conversion of local trees to ts and can add mutations, then no matter if the final output is haplotypes or ts, I think we can do it in this way. For haplotypes, we can estimate the number of SNPs from branch length and mutation rate, the new argument could be about scaling factor or approximate segSites (I'm not sure if mutRate is intuitive enough for sampling...) 3. Does runMACs allow for custom recombination map? I can't seem to find this parameter (I'm wondering if it's supposed to be the recombination rate r or something else?)

[gregorgorjanc]

@LynxJinyangii I agree that we want go away from working with haplotype matrices and I like your idea of estimating how many SNP we could get from given local trees / ts (via mutRate or segSites or some function thereof), BUT we need a PATHWAY to get to this future situation - we can's juts rip out current runMacs/haplotypes and replace it with ts because that will have all sorts of downstream implications (whole of AlphaSimR is based on haplotypes in pop@geno so we need the haplotypes for now). Hence we the need to keep the old functionality and build the ts functionality in parallel. We keep getting back to this point. It is critical to have this parallelism or we will break the current functionality, which is not what we want.

To avoid this issue, we discussed that you would work on the parallel runMacsTs().

[LynxJinyangii]

I understand your concerns, but what I mean is not to discard the haplotypes output, but to extract the genotype matrix (haplotypes) from ts after adding mutations on ts (if both 1. record tree from runMacs 2. adding mutations can work), and the downstream would keep using haplotypes.

[gregorgorjanc]

OK, then we are on the same page! When I said let's work on parallel runMacsTs() was exactly this, so that we get the new functionality and we still retain the old. Once both work, we can think of the next steps in managing their lifetime in the package.