recapitation with recombination map

[Gertjan Bisschop]

Hi Ben, Hi Peter,

Thanks for the great software you have provided us with!

One thing is giving me headaches though.

I've been trying to run a SLiM simulation with a recombination map and recapitate the result using pyslim. However, for some reason the recapitation process finishes without error, but with a certain fraction of trees still having more than one root. The interval of those trees is always extremely small (mean=0.0004) vs 1080 for all other trees (for a particular run of the minimal example I have included). I assume I'm misspecifying something here that trips up pyslim. So I have included both recombination maps that I have used.

`

initialize() {
   
    initializeMutationType("m1",0.5,"f",0.0);
    initializeGenomicElementType("g1", m1, 1);
    initializeGenomicElement(g1, 100,100);
    initializeMutationRate(0.0);

    //read in hapmap
   
    lines = readFile('hapmap_slim_half.txt');
    ends = sapply(strsplit(lines[0],','),'asInteger(applyValue);');
    rates =sapply(strsplit(lines[1],','),'asFloat(applyValue);');

    initializeRecombinationRate(rates,ends);
    initializeTreeSeq(simplificationRatio=NULL, simplificationInterval=1000);

}

// create a population of 500 individuals
1 {
    defineConstant('simID',getSeed());
    defineConstant('m',0.001);
   
    defineConstant("Ne_1", 1000);
    defineConstant("Ne_0", 1000);

    sim.addSubpop('p1',Ne_1);

}



1000 late() {
    sim.addSubpopSplit("p0",Ne_0,p1);
    p1.setMigrationRates(p0,m);
    }

1000:8000 late(){
    if (sim.generation%1000==0){print('sim at time step '+sim.generation);}
     
}

8000 late() {
    sim.treeSeqOutput('min_example.trees');
    sim.simulationFinished();
    }

`

Recapitation

`

#!/usr/bin/env python3

import msprime, pyslim
import numpy as np
import os, glob
import numpy.random

def main():
    directory = os.getcwd()
    treefile='min_example.trees'
   
    with open('hapmap_msprime_half.txt','r') as file:
        positions = [int(element) for element in file.readline().rstrip().split(',')]
        rates = [float(element) for element in file.readline().rstrip().split(',')]
    hapmap = msprime.RecombinationMap(positions,rates)

    mutation_rate = 1.9e-7
    Ne_ancestral = 1000
   
    ts = pyslim.load(os.path.join(directory, treefile))
       
    recap = ts.recapitate(recombination_map=hapmap, Ne=Ne_ancestral)
   
    recap.dump(directory+'/recap_'+treefile)

if __name__ == '__main__':
    main()

`

Thanks,

gertjan

[Peter Ralph]

Hello! Sorry to hear about the headaches. Let's see. After running
your SLiM and python scripts, I get:

ts = pyslim.load("recap_min_example.trees")
max([t.num_roots for t in ts.trees()])
# 1

Ah-ha, but that was with the development version of msprime. if I
instead do this with msprime 0.7.4, I get trees with 2 roots.
So, my guess is that this has something to do with the changes to
recombination that have happened in msprime - I'm not remembering what
might be causing this, though. Is using the development msprime an
option for you?

> --
> SLiM forward genetic simulation: http://messerlab.org/slim/
> ---
> You received this message because you are subscribed to the Google Groups "slim-discuss" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to slim-discuss...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/slim-discuss/7135b3c2-07b7-4f07-a33d-20342a20ec40%40googlegroups.com.

[Gertjan Bisschop]

Thanks Peter for checking this.
I will give it a go with the development msprime.

Kind regards,
gertjan

[Gertjan Bisschop]

Solved.
Thanks Peter!