Recapitation memory requirements

[akamolphat]

Hi everyone,

I was planning to utilise recapitation and mutation overlay to add neutral mutations to my treesequence file from a non-neutral SLiM simulation. The non-neutral SLiM simulation was run for 1000 ticks and involved about 2 million individuals. It involved a fixed number of non-neutral mutations (120 mutations) being added at a certain frequency (30%) in the founder population. Due to previous OUT_OF_MEMORY errors during the recapitation stage, I did some subsampling and trimming (splitting one extra-large genome into chromosomes through ts.keep_interval(), ts.rtrim() and ts.ltrim()) of the treesequence file. I end up with a treesequence files of 28229 diploid individuals and the largest chromosome is about 150 Megabases in length. The largest treesequence file is about 800MB and runs out of memory even when using 500 GB RAM on a UNIX cluster. I see that people have been able to recapitate and overlay mutations on larger treesequence files using less memory (128GB memory machine, Code Sample 3 and 4 in the SimHumanity paper). Is my simulated population size too big? Are there any ways around this?

I tried the same subsampling/trimming/recapitation code on a smaller treesequence file with a smaller chromosome (1000 times smaller), same number of non-neutral mutations (120), and a smaller population size (~1000 times smaller) and I was able to recapitate very quickly on my computer so I don't think the issue is with the code.

The input treesequence file is very large so I am not sure if or how I should upload it, but the code I used for the recapitation processs is as follows:

import tskit
import numpy as np
import os
import sys
import warnings
import pyslim
import msprime

intreefile="Chr1.trees"
outrecaptrees="Chr1_recap.trees"
outovrlaytrees="Chr1_recap_overlay.trees"

gt = np.float64(4.0) 
Ne = np.float64(1944663.0)
recomb_rate = 1.01e-8/gt
mut_rate = 4.6e-9/gt
ancestral_Ne = Ne*gt
ts = tskit.load(intreefile)
print("Start recapitation process\n", flush = True)
with warnings.catch_warnings(record=True) as w:
    rts = pyslim.recapitate(ts,
                                recombination_rate=recomb_rate, # rates already divided by gt earlier
                                ancestral_Ne=ancestral_Ne,
                                random_seed=5)
    print("Dump recapitated tree file\n", flush = True)
    rts.dump(outrecaptrees)
    print("Start mutation overlay\n", flush = True)
    mts = msprime.sim_mutations(
                                rts,
                                rate=mut_rate,
                                model=msprime.SLiMMutationModel(type=0),
                                keep=True)
    print("Dump overlaid tree file\n", flush = True)
    mts.dump(outovrlaytrees)

Any help would be greatly appreciated.

A (Kamolphat)

[petrelharp]

Questions:

• which step is using all the memory? (I assume recapitation, but you should check.)
• what model parameters are affecting things? (To find this out, do some experiments; for instance, record max memory usage with ancestral_Ne = 10; and 100, and 500, and 1000, etc; then make a plot to see how it's scaling and when it starts to become a problem.)

[akamolphat]

Hi @petrelharp

Based on the print statements in the *.out file of the slurm script, the crash happens at the pyslim.recapitate() stage.

I will test out the ancestral_Ne parameter as you have suggested and reply to this once I have it. Is it expected that the memory used scales linearly with ancestral_Ne?

A