Using msprime to simulate initial populations with multiple chromosome in SLiM 5.0

[Kamolphat Atsawawaranunt]

Hi Ben,

I am not sure if this is the right place for this question, but I hope that you can point me in the correct direction.

I was hoping to first generate standing genetic diversity using msprime and then use this to start a SLiM simulation in SLiM 5.0 with multiple chromosomes (i.e., use the initializeChromosome() function).

For SLiM 4.3, the simulated chromosomes (one long chromosome with altered recombination rates) from msprime 1.3.3 seems to work. However, how do I achieve this using msprime for SLiM 5.0 if I would like to use initializeChromosome()?

I noticed that there has been some changes in the way treesequence file is stored, so I have tried to mimic this by storing the tree sequence file from msprime for each chromosome (e.g., chromosome_1.trees, chromosome_2.trees) in a folder with a ".trees" suffix (e.g., test.trees). However, when I tried to read these into SLiM, I get the following error:

ERROR (Species::ReadTreeSequenceMetadata): the version of this file appears to be too old to be read, or the file is corrupted; you can try using pyslim to bring an old file version forward to the current version, or generate a new file with the current version of SLiM or pyslim.

I understand that I can still use the old method for SLiM 5.0, but I thought that it would be nice to transition to the newer function. Are there significant computation efficiency gains through the use of initializeChromosome() over the old methods of altering recombination rates?

Yours sincerely,
A

[Ben Haller]

Hi A!

You're on the cutting edge!  Your trees archive, that you've assembled by hand with msprime .trees files, does not have the metadata that is currently required.  You can read about those metadata requirements in chapter 29 of the manual.  You could also generate a trees archive from a SLiM simulation that is configured the way you want, and examine the metadata produced by SLiM from that simulation, and if it matches your msprime simulation closely enough, you could probably even copy that metadata over.  :-O

Normally pyslim would be the tool you'd use to annotate msprime-produced simulations with the appropriate metadata for SLiM.  Unfortunately pyslim is a little behind at the moment, though; the pyslim 1.1b1 version that Peter recently released works with single-chromosome simulations in SLiM 5, but the complete multi-chromosome support needed to do what you're trying to do is not yet there, as I understand it.  Peter is plugging away at it, so it should be along before too terribly long.

So the best thing might just be to be patient and wait for pyslim 1.1, which should make this easy.  But if you're impatient, you can (a) continue with the SLiM 4.3 strategy of modeling "effective chromosomes" with r=0.5 points, and transition over to real chromosomes once the tools are ready (I would expect the performance to be very similar on the SLiM side; on the msprime side, I don't know...), or (b) try to put the correct metadata in yourself, which might just be a matter of copying it over from a saved trees archive from SLiM that is configured identically to your simulation.

I hope this makes sense.  Great to see someone already trying out the new multi-chromosome features; I'm sorry they're not quite there yet.  Happy modeling!

Cheers,
-B.

Benjamin C. Haller
Messer Lab
Cornell University


Kamolphat Atsawawaranunt wrote on 5/18/25 7:22 PM:

A --
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 visit https://groups.google.com/d/msgid/slim-discuss/856e5882-68e6-4818-a257-bfd975756f92n%40googlegroups.com.

[Peter Ralph]

Well gee, I thought it'd be a *little* longer before someone was wanting to do this. I'll be working on it slowly over the next few weeks. 

If you are interested in doing it yourself, here's a very rough outline: the main obstacle to assembling this yourself is that the node tables need to be identical; so, you'd append all the non-sample nodes from one tree sequence to the other one, and re-number references to it in that one, tidying up other bookkeeping. Feel free to open a discussion on pyslim if you want to work through doing it.

-peter

---

From: 'Ben Haller' via slim-discuss <slim-d...@googlegroups.com>
Sent: Sunday, May 18, 2025 4:37 PM
To: Kamolphat Atsawawaranunt <a.kam...@gmail.com>
Cc: slim-discuss <slim-d...@googlegroups.com>
Subject: Re: Using msprime to simulate initial populations with multiple chromosome in SLiM 5.0

 

To view this discussion visit https://groups.google.com/d/msgid/slim-discuss/75f258bf-4734-fcda-e5b6-c5d5ca409fdb%40mac.com.

[Kamolphat Atsawawaranunt]

Hi Ben and Peter,

Thank you so much for your reply, and thank you so much for SLiM and pyslim. Peter, I will work on the workarounds so there is no rush or need for this at all. I just thought that it may have been something that was already out there that I am unaware of. I don't think I am skilled enough to manipulate the node tables or the metadata of the tree sequence, so I will patiently wait for future updates or follow the SLiM 4.3 strategy. I thought it was neat to simulate generate standing genetic diversity with msprime, assuring coalescence so I have looked at this. However, my models are actually quite small and I was not intending to run multiple replicates so it is very achievable with just a long(er) SLiM burn-in (if I really want to keep the multiple-chromosome functionality).

Yours sincerely,

A

[Kamolphat Atsawawaranunt]

Hi Ben,

I was not sure if I should start a new thread or not, but I thought that it is somewhat related so I have posted it here.

I was trying the new initializeChromosome function. However, I noticed significantly slower speeds and more memory usage when running the simulation with initializeChromosome() instead of the old strategy of modelling effective chromosomes. For the exact code I ran (see below), using initializeChromosome() was more than three times slower and required more than three times memory than using the old strategy of modelling effective chromosomes. Am I doing something wrong here?

Example code of when I used intializeChromosome() strategy (took me 40.513 CPU s, 103.9 MB memory):

initialize() {

// Initialize nonWF model

initializeSLiMModelType("nonWF"); // Non-WF model

// Define parameters

defineConstant("gen", 500); // number of generations

defineConstant("chrom_len", 5e4);

defineConstant("C_num_QTN", 10); // number of chromosomes with selection (for later part of the simulation)

defineConstant("C_num_Nut", 10); // number of chromosomes with only neutral SNPs

defineConstant("C_num", C_num_QTN + C_num_Nut); // number of chromosomes in total

defineConstant("C_lengths", rep(chrom_len, C_num)); // lengths for each chromosome

defineConstant("R", 1e-5); // recombination rate

defineConstant("MU_base", 1e-7); // base mutation rate overall

defineConstant("N", 4000); // Total size

// Initialize other parameters

initializeSex("A");

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

m1.convertToSubstitution = T;

initializeGenomicElementType("g1", m1, 1.0);

initializeGenomicElementType("g2", m1, 1.0); // Initialise separate genomic element type for now, QTL to be added

// Initialize chromosomes

for (i in seq(1, C_num_QTN)){

initializeChromosome(i, chrom_len);

initializeGenomicElement(g1, 0, chrom_len-1);

initializeMutationRate(MU_base);

initializeRecombinationRate(R);

}

for (i in seq(C_num_QTN + 1, C_num)){

initializeChromosome(i, chrom_len);

initializeGenomicElement(g2, 0, chrom_len-1);

initializeMutationRate(MU_base);

initializeRecombinationRate(R);

}

defineConstant("K", 500); // carrying capacity

}

1:gen reproduction() {

subpops=sim.subpopulations;

for (pop in subpops){

// parents are chosen proportional to fitness

female=which(pop.individuals.sex=="F");

male=which(pop.individuals.sex=="M");

parents1 = sample(pop.individuals[female], N, replace=T);

parents2 = sample(pop.individuals[male], N, replace=T);

for (i in seqLen(N)){

pop.addCrossed(parents1[i], parents2[i]);

}

self.active = 0; // call the reproduction callback only once

}

}

1:gen survival() {

return (individual.age == 0);

}

1 early() {

// Add initial panmictic pops p0

sim.addSubpop("p0", N);

}

Example code of what I effective chromosomes strategy (took me 9.384 CPU s, 31.3 MB memory):

initialize() {

// Initialize nonWF model

initializeSLiMModelType("nonWF"); // Non-WF model

// Define parameters

defineConstant("gen", 500); // number of generations

defineConstant("chrom_len", 5e4);

defineConstant("C_num_QTN", 10); // number of chromosomes with selection (for later part of the simulation)

defineConstant("C_num_Nut", 10); // number of chromosomes with only neutral SNPs

defineConstant("C_num", C_num_QTN + C_num_Nut); // number of chromosomes in total

defineConstant("C_lengths", rep( chrom_len , C_num)); // lengths for each chromosome

defineConstant("R", 1e-5); // recombination rate

defineConstant("MU_base", 1e-7); // base mutation rate overall

defineConstant("N", 4000); // Total size

// Initialize other parameters

initializeMutationRate(MU_base);

initializeSex("A");

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

m1.convertToSubstitution = T;

initializeGenomicElementType("g1", m1, 1.0);

initializeGenomicElementType("g2", m1, 1.0); // Initialise separate genomic element type for now, QTL to be added

// Initialize genome

defineConstant("g2_start", sum(C_lengths[0:(C_num_QTN - 1)]));

defineConstant("gen_len", sum(C_lengths[0:(C_num - 1)]) - 1);

initializeGenomicElement(g1, 0, g2_start-1);

initializeGenomicElement(g2, g2_start, gen_len-1);

// C_num chromosomes of 50 kb each

rates = c(rep(c(R, 0.5), C_num-1), R);

ends = repEach(cumSum(C_lengths), 2);

ends = ends[0:(length(ends) - 2)];

ends = ends - c(rep(c(1,0), C_num-1), 1);

initializeRecombinationRate(rates, ends);

defineConstant("K", 500); // carrying capacity

}

1:gen reproduction() {

subpops=sim.subpopulations;

for (pop in subpops){

// parents are chosen proportional to fitness

female=which(pop.individuals.sex=="F");

male=which(pop.individuals.sex=="M");

parents1 = sample(pop.individuals[female], N, replace=T);

parents2 = sample(pop.individuals[male], N, replace=T);

for (i in seqLen(N)){

pop.addCrossed(parents1[i], parents2[i]);

}

self.active = 0; // call the reproduction callback only once

}

}

1:gen survival() {

return (individual.age == 0);

}

1 early() {

// Add initial panmictic pops p0

sim.addSubpop("p0", N);

}

Yours sincerely,

A

[Ben Haller]

Hi A!

OK, complicated.  I don't think there's anything wrong with your model, at least not that I've noticed so far.  What's happening is, I think, a bit complicated.  Basically, your model is positioned at an unfavorable spot in the parameter space.  SLiM tends to be tuned to make more biologically realistic models run faster, and to make bigger models run faster.  If a model is biologically unrealistic, I don't focus on it as much in my optimization work, because I figure people won't want to run such models anyway; and if a model is small enough to run fairly quickly, I don't worry about it as much since it's already reasonably fast.  I focus primarily on making the really big models run faster, because that's where people really hit the wall in terms of the runtime of the model limiting what they're able to simulate.

The key parameter here, I think, is the chromosome length, 5e4.  I've aimed at making longer chromosome lengths, which are more biologically realistic for most organisms, and which have longer runtimes, run as fast as possible.  If you kick the chromosome length up to 5e6, the multi-chromosome version of the model actually runs *faster*, in my quick tests just now; 304.7 seconds versus 349.5 seconds for the single-chromosome version.  Of course you'd have to run lots of replicates, etc., to get a firm result; I suspect this result would hold when averaged across many runs of the model (I saw similar gains in many tests I ran during development), but I haven't tried it with your model.  The point is, SLiM 5 with multiple chromosomes can be *faster* for large, biologically realistic models than SLiM 4 was, and than a single-chromosome model would be, because that's where I focused my optimization efforts; the underlying architecture does have some significant new optimizations that were not present in SLiM 4.

But interestingly, making the chromosome length *shorter* than your model also makes the multi-chromosome version perform better!  If I make it 5e1 (i.e., 50), the multi-chromosome version runs in 2.53 and the single-chromosome version in 2.14 (for one test run); so the multi-chrom version is slower, but by much less wide of a margin than you observed.  What's going on there?

That's less a result of deliberate optimization, and more just luck.  I *think* what is happening with your model has to do with "mutation runs", a genetic structure used under the hood to capture haplotype structure and leverage it to make runs faster – basically, SLiM can detect that haplotype structure exists, and reduce the amount of computation needed by doing some calculations just once, for a given haplotype, and reusing that computation for all occurrences of the haplotype across the model.  I *think* the reason that your chromosome length of 5e4 performs so much better as a single-chromosome model is perhaps that, at that length scale, with the particular mutation rate and recombination rate you're using, haplotype structure actually extends across more than one chromosome to a useful degree – calculations can be shared, not just within a single chromosome, but across multiple "effective chromosomes" in the single-chromosome version of the model.  That is not possible with SLiM's architecture, in the multi-chromosome version of the model; each chromosome is tracking haplotype structure independently.  So the single-chromosome version of the model gets a much bigger boost from that optimization.  And that happens particularly around the length that you chose, of 5e4; the effect of it is much less pronounced both at longer chromosome lengths (where the odds of a haplotype usefully spanning more than one chromosome diminish) and at shorter lengths (where the complexity of those haplotype-based calculations probably makes SLiM a little bit slower than it would be if it just ignored haplotypes and claculated everything by brute force).

That's just a guess.  The internals of SLiM are complex enough that things like this are difficult to analyze, and depend upon parameter values as well as your particular hardware architecture, operating system platform, etc.  Things like the size of the CPU cache can make a big difference to how well a particular model performs.  It's really hard to make a solid prediction.  But the take-home point, I guess, is: if performance is really important to you, then try writing your model in several different ways, and see which performs the best.  But you have to be careful – one version might perform the best on your own work machine, but a different version might perform the best on the computing cluster you use for your production runs!  Also, you have to be careful – one version might perform the best with a chromosome length of 5e4, but a different version might perform the best with a length of 5e6!  There's not really anything to do but experiment and measure, and you have to do that with the final parameterization of the model, on the final hardware you intend to run on, otherwise it might not be worth much.  Performance optimization is a complex topic.  :->

Cheers,
-B.

Benjamin C. Haller
Messer Lab
Cornell University

Kamolphat Atsawawaranunt wrote on 5/21/25 7:49 PM:

To view this discussion visit https://groups.google.com/d/msgid/slim-discuss/96111e59-79b6-4af3-b849-405842ffd8a3n%40googlegroups.com.

[Kamolphat Atsawawaranunt]

Starting a multi-chromosome SLiM simulation using treesequence files from msprime
Hi everyone,

Recapitation in my scenario seems to take a long time (probably due to sequence length and the number of individuals involved), so I decided to try to start a SLiM simulation from treesequence files generated from msprime. The simulation in msprime is under WF, and my SLiM simulation is nonWF, but I am hoping to use the files generated from msprime to kick start the SLiM simulation and save some time.

I decided to post the question here in a thread that I have asked a similar question before, but I am not sure if this is the appropriate place. 

In summary, I get the following error every time I read in treesequence files created by msprime into SLiM when it had nine or more chromosomes:

Simulation Runtime Error
ERROR (Species::__TabulateSubpopulationsFromTreeSequence): unexpected node metadata length; this file cannot be read.
This error has invalidated the simulation; it cannot be run further. Once the script is fixed, you can recycle the simulation and try again.

I tried simulating many different lengths of chromosomes, and also with varying lengths, and all it seems work as long as it has fewer than nine chromosomes.

To avoid the problem with node table not being identical across different chromosome treeseq files, I have gone down the following workflow:
1. Use msprime and recombination map to create a pseudo-chromosome treeseq file
2. Split pseudo-chromosome treeseq file into separate treeseq file for each chromosome
3. Update the metadata of the split/separated treeseq file.
4. Output split treeseq file into a folder with a .trees suffix.
The code for generating tree sequence files for multichromosome SLiM models are as follows (please alter the chr_no variable to 9 to recreate the error):

import msprime, pyslim, tskit, random
import numpy as np
import pandas as pd
import os

# 0.a. Define input parameters
chr_no = 8 # outputs from chr_no = 9 throws error in SLiM
chr_len = 5000
recomb = 1e-8
nsamp = 2000
ne = nsamp
mu_rate_scaled = 1e-8
outcoalesce = "coalesce.trees"
outoverlaid = "mut_overlaid.trees"
outfoldname = "test_Multichrom.trees/"
output_recomb = "test_recomb_as_simulated.csv"

# Create example recombination file
dt = pd.DataFrame({'Chromosome_ID': np.arange(1, chr_no + 1), 'chr_len': np.repeat(chr_len, chr_no), 'predicted_recomb_R': np.repeat(recomb, chr_no)})

# 0.b.iv. Use result of 0.b.ii. to create cumulative length of chromosome end points minus 1 (msprime is zero-based) anduse value + 1 for chromosome start (excluding the last chromosome)
chr_cumsum = dt['chr_len'].cumsum() - 1
chr_cumsum_plus_one = chr_cumsum + 1
chr_cumsum_plus_one = chr_cumsum_plus_one[:-1]
chr_len_ls = np.zeros(((len(chr_cumsum)*2)), dtype = int)
chr_len_ls[1::2] = chr_cumsum
chr_len_ls[2::2] = chr_cumsum_plus_one
chr_len_ls[-1] += 1

# 0.b.v. Use result from 0.b.iii to create recombination rates matching the cumulative length from 0.b.iv.
chr_R = dt['predicted_recomb_R']
chr_brk_val = 0.5
chr_R_ls = np.full(((len(chr_R)*2) - 1), chr_brk_val)
chr_R_ls[::2] = chr_R

# 0.b.vi. Use result from 0.b.iv. and 0.b.v. to create an msprime rate map
recomb_map = msprime.RateMap(position=chr_len_ls, rate=chr_R_ls)
# Also create a start_stop matrix to be added to dt and output to csv
nchrom = int(len(chr_len_ls)/2)
start_stop = chr_len_ls.reshape(nchrom,2)
dt['start'] = start_stop[:,0]
dt['end'] = start_stop[:,1]
dt['chr_size'] = dt['end'] - dt['start']
# 0.b.vii. Output rate map to CSV for SLiM simulation. We can use the same input file here, but I have noticed potential rounding problems in SLiM
dt.to_csv(output_recomb, index = False)
# 1. Use msprime and recombination map to create a pseudo-chromosome treeseq file
seq_len = chr_len_ls[-1]
# 1.a. simulate ancestry
ts = msprime.sim_ancestry(samples=nsamp, population_size=ne, sequence_length=seq_len, recombination_rate=recomb_map)
# 1.b. dump tables to allow altering metadata
tables = ts.dump_tables()
# nonWF model, start at tick = 1
pyslim.annotate_tables(tables, model_type="nonWF", tick=1)
# Loop through each individual to assign sex
individual_metadata = [ind.metadata for ind in tables.individuals]

for md in individual_metadata:
md["sex"] = random.choice([pyslim.INDIVIDUAL_TYPE_FEMALE, pyslim.INDIVIDUAL_TYPE_MALE])
md["age"] = 0 # I only assigned values of zero for simplicity, this value will get re-written in SLiM anyways

# Assign new values to SLiM metadata
ims = tables.individuals.metadata_schema

tables.individuals.packset_metadata([ims.validate_and_encode_row(md) for md in individual_metadata])
# Dump out the treesequence file as checkpoint
print("Dump empty treesequence file with coalescence information")
slim_ts = tables.tree_sequence()
slim_ts.dump(outcoalesce)

#slim_ts = tskit.load(outcoalesce)
print("Overlay treesequence file with neutral mutations")
mut_model = msprime.SLiMMutationModel(type=1)
ms_ts = msprime.sim_mutations(
            slim_ts,
            rate=mu_rate_scaled,
            model=mut_model,
            keep=True,
            random_seed=12)

print("Dump overlaid treeseq files")
ms_ts.dump(outoverlaid)

# ms_ts = tskit.load("overlaid.trees")
# Loop through each chromosome, split the treesequence file accordingly, update the metadata accordingly
# First, there is a header that is shared across all treesequence file
ls = []
for i in np.arange(1, (dt.shape[0] + 1)):
print(i)
ls.append({'id': int(i), 'index': int(i-1), 'symbol': str(i), 'type': 'A'})

# Loop through the start_stop values
os.makedirs(outfoldname, exist_ok = True)
for i in range(0, nchrom):
    # 2. Split pseudo-chromosome treeseq file into separate treeseq file for each chromosome
    this_chrom = i + 1
    # Start and end position
    start = start_stop[i, 0]
    end = start_stop[i, 1]
    # The keep_intervals method handles the complex edge manipulations
    # I think we might need to simplify = F because we want to recapitate (and this will likely keep the root nodes).
    split_ts = ms_ts.keep_intervals([[start, end]], simplify=False).ltrim().rtrim()
    # 3. Update the metadata of the split/separated treeseq file.
    ms_tab = split_ts.dump_tables()
    ts_metadata = ms_tab.metadata
    ts_metadata["SLiM"]["chromosomes"] = ls
    this_chrom_dict= ls[i]
    ts_metadata["SLiM"]["this_chromosome"] = this_chrom_dict
    ms_tab.metadata = ts_metadata
    ms_tab.sequence_length = end - start
    split_ts = ms_tab.tree_sequence() # overwrite file
    # Save the split tree sequence
    split_ts.dump(outfoldname + "/" +"chromosome_"+ str(this_chrom) + ".trees")

The SLiM code for reading in these outputs is as follows:

initialize() {
if (exists("slimgui")){
defineConstant("in_tree", "test_Multichrom.trees");
defineConstant("in_file", "test_recomb_as_simulated.csv"); // Table of genome sizes and recombination rates
}
initializeSLiMModelType("nonWF");
initializeSex();

initializeMutationType("m1", 0.5, "f", 0.0);
m1.convertToSubstitution=T;

initializeMutationType("m2", 0.5, "f", 0.1);
initializeGenomicElementType("g1", c(m1,m2), c(1.0, 0.001));
// Initialize genome
dt = readCSV(in_file);
defineGlobal("ChromID", dt.getValue("Chromosome_ID"));
defineGlobal("ChromSize", asInteger(dt.getValue("chr_size")));
defineGlobal("ChromRecomb", dt.getValue("predicted_recomb_R")); // 1000 times higher recombination rate
C_num = length(ChromSize);
for (i in seqLen(C_num)){
// Note that Chromosome ID starts at 1 and not 0
initializeChromosome(ChromID[i], ChromSize[i]);
initializeGenomicElement(g1, 0, ChromSize[i]-1);
initializeMutationRate(0.0);
initializeRecombinationRate(ChromRecomb[i]);
}
}

1 early() {
catn("read treeseq");
sim.readFromPopulationFile(in_tree);
catn("Finished reading treeseq");
}

reproduction(NULL, "F") {
subpop.addCrossed(individual, subpop.sampleIndividuals(1, sex="M"));
}

early() {
p0.fitnessScaling = 2000 / p0.individualCount;
}

2000 early() { sim.simulationFinished(); }

Yours sincerely,

Kamolphat (A)

[Ben Haller]

Hi A,

This is presumably a bug, although the issues surrounding it are complex and will require investigation.  I'd suggest that you open a new issue (on SLiM) and Peter and I will look into it when we have time (which is not right now :->).  Thanks!

Cheers,
-B.

Benjamin C. Haller
Messer Lab
Cornell University

To view this discussion visit https://groups.google.com/d/msgid/slim-discuss/bf523ea5-d9f5-4c13-a487-b55f691ed4fan%40googlegroups.com.

[Peter Ralph]

Thanks for the report! Amazingly, given the complexity of this example, I know what the problem is: "works below 8 chromosomes" gives it away. I've written a brief explanation and (possibly, only very lightly tested) workaround here: https://github.com/tskit-dev/pyslim/issues/410

If you do use that code *please* make sure it's doinig the right thing, and let me know if you hit any issues.

thanks,

--peter

To view this discussion visit https://groups.google.com/d/msgid/slim-discuss/a1f97403-9220-4f3f-90a6-fc7a8f23a9ee%40mac.com.

[Ben Haller]

Hey, awesome, thanks Peter!

Cheers,
-B.

To view this discussion visit https://groups.google.com/d/msgid/slim-discuss/CAO0xWMBUwaLoXGS8%2BZ6HHCoQUhrhmC4ePu07QeLQ-boxrOHvvw%40mail.gmail.com.

[Kamolphat Atsawawaranunt]

Hi Ben and Peter,

Thank you so much for your replies, and for these programs and providing such good support on it.

Peter, thank you so much for presenting the workaround. I have gone down a slightly different pathway so I have yet to test out the codes. However, when I do test out the codes to see if it is doing the right thing, I will let you know of any issues on the pyslim GitHub page.

Yours sincerely,
A