reading rec map from file: slim complaining that ends not in ascending order - but they are?

24 views
Skip to first unread message

Jeffrey Groh

unread,
Mar 30, 2022, 8:25:47 PM3/30/22
to slim-discuss
Hi Slim team,

I expect that I must be missing something glaringly obvious - but here goes.

I'm getting the error message:

(SLiMSim::ExecuteContextFunction_initializeRecombinationRate): initializeRecombinationRate() requires ends to be in strictly ascending order." 

when trying to set up the simulation below by using the example code in the manual to read in a recombination map from a file (attached). I'm confused because I've convinced myself the ends are in fact sorted in ascending order, e.g. when I run `any(ends != sort(ends))` the output is `F` as desired. Is there something simple I'm missing?

initialize() {

   initializeTreeSeq();

   initializeMutationRate(0);

   initializeMutationType("m1", 0.5, "f", 0.0);

   initializeGenomicElementType("g1", m1, 1.0);

   // number of windows present in rec map

   defineConstant("L", 23600);

   initializeGenomicElement(g1, 0, L-1);

   // read recombination map for hg1

   lines = readFile("/Users/jeff/workspace/selection-against-introgression/datasets/human/Kong_etal_recmaps/hg1_recmap.txt");

   rates = NULL;

   ends = NULL;

   for (line in lines)

   {    

       components = strsplit(line, "\t");

       ends = c(ends, asInteger(components[0]));

       rates = c(rates, asFloat(components[1]));

   }

   ends = ends/10000;

   ends = ends - min(ends) + 1;

   ends = asInteger(c(ends[1:(size(ends)-1)] - 2, L-1)); 

 

   rates = rates*0.01;

   initializeRecombinationRate(rates, ends);

}

1 {

   sim.addSubpop("p1", N);

}

2 late() {

   sim.simulationFinished();

}

hg1_recmap.txt

Ben Haller

unread,
Mar 30, 2022, 10:08:32 PM3/30/22
to slim-discuss
Hi!  "Strictly ascending" means that each successive value must be *greater* than the previous value; being equal does not cut it.  This is because if two successive values are equal, then you are in effect trying to set the recombination rate at a given position twice.  The `ends` vector produced by your input file and your code has two successive values of 1446 in it, violating this rule, as revealed by a bit of poking around in the Eidos console.  Fix that, and you should be good to go.  :->  Happy modeling!

Cheers,
-B.

Benjamin C. Haller
Messer Lab
Cornell University

Jeff

unread,
Mar 31, 2022, 1:53:29 PM3/31/22
to slim-discuss
Ah, OK this makes sense. Many thanks for spotting this which I should have seen myself!

Cheers,
Jeff
Reply all
Reply to author
Forward
0 new messages