Sampling individuals and calculating microsat expected Heterozygosity

[Ravi Vishwakarma]

Hello All, 

This is regarding SLiM recipe 14.11 (Modelling microsatellite). I modified the script to calculate the expected heterozygosity at each locus and average it across all loci. Here is the working part of the code:

10000 late() {

// print frequency information for each microsatellite site

all_msats = sim.mutationsOfType(m2);

i = 0.;

He_sum = 0.;

for (pos in sort(msatPositions))

{

catn("Microsatellite at " + pos + ":");

msatsAtPos = all_msats[all_msats.position == pos];

for (msat in sortBy(msatsAtPos, "tag")){

catn(" variant with " + msat.tag + " repeats: " +

sim.mutationFrequencies(NULL, msat));

He_locus = 1. - sum(sim.mutationFrequencies(NULL, msatsAtPos)^2) ;

}

i = i+1.;

print(He_locus);

He_sum = He_sum + He_locus;

}

print(He_sum + "_");

}

The above code works well. The problem is when I want to sample some individuals from the population and calculate heterozygosity just for the sampled individuals.

Here is what I was trying, which gives an error.

10000 late() {

// print frequency information for each microsatellite site

sampled_ind = p1.sampleIndividuals(5);

all_msats = sampled_ind.genomes.mutationsOfType(m2);

i = 0.;

He_sum = 0.;

for (pos in sort(msatPositions))

{

catn("Microsatellite at " + pos + ":");

msatsAtPos = all_msats[all_msats.position == pos];

for (msat in sortBy(msatsAtPos, "tag")){

catn(" variant with " + msat.tag + " repeats: " +

sim.mutationFrequencies(NULL, msat));

He_locus = 1. - sum(sim.mutationFrequencies(sampled_ind, msatsAtPos)^2) ;

}

i = i+1.;

print(He_locus);

He_sum = He_sum + He_locus;

}

print(He_sum + "_");

}

I understand the error but have not been able to resolve it. Any suggestions on how to sample individuals and calculate heterozygosity would be greatly appreciated!

Thank you in advance for help!

Regards,

Ravi

[Ben Haller]

Hi Ravi!  Please supply a full model that can be run in SLiMgui (as minimal as possible; just a skeleton within which the code you posted can run).  And what error message are you getting?  Hard to guess what the problem might be without that information.

Cheers,
-B.

Benjamin C. Haller
Messer Lab
Cornell University


Ravi Vishwakarma wrote on 6/13/24 12:33 PM:

--
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/50549a4d-ee35-4782-932d-6e857354809bn%40googlegroups.com.

[Ravi Vishwakarma]

Hi Ben! 

Here is the full model that runs in SLiMgui.

// Keywords:

initialize() {

defineConstant("L", 1e6); // chromosome length

defineConstant("msatCount", 10); // number of microsats

defineConstant("msatMu", 0.0001); // mutation rate per microsat

defineConstant("msatUnique", T); // T = unique msats, F = lineages

initializeMutationRate(1e-7);

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

initializeGenomicElementType("g1", m1, 1.0);

initializeGenomicElement(g1, 0, L-1);

initializeRecombinationRate(1e-8);

// microsatellite mutation type; also neutral, but magenta

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

m2.convertToSubstitution = F;

m2.color = "#900090";

}

1 late() {

sim.addSubpop("p1", 500);

// create some microsatellites at random positions

genomes = sim.subpopulations.genomes; //take all the genomes

positions = rdunif(msatCount, 0, L-1); //msatCount random positions across the genome

repeats = rpois(msatCount, 20) + 5; // generate msatCount numbers with poisson dist with mean 20

for (msatIndex in 0:(msatCount-1)) //going over each msatCount (each microsat locus)

{

pos = positions[msatIndex]; //pos taking the position value

mut = genomes.addNewDrawnMutation(m2, pos); //add new mutation to the target genome(s) with specified mutation type

mut.tag = repeats[msatIndex];//

}

// remember the microsat positions for later

defineConstant("msatPositions", positions);

}

modifyChild() {

// mutate microsatellites with rate msatMu

for (genome in child.genomes)

{

mutCount = rpois(1, msatMu * msatCount);

if (mutCount)

{

mutSites = sample(msatPositions, mutCount);

msats = genome.mutationsOfType(m2);

for (mutSite in mutSites)

{

msat = msats[msats.position == mutSite];

repeats = msat.tag;

// modify the number of repeats by adding -1 or +1

repeats = repeats + (rdunif(1, 0, 1) * 2 - 1);

if (repeats < 5)

next;

// if we're uniquing microsats, do so now

if (msatUnique)

{

all_msats = sim.mutationsOfType(m2);

msatsAtSite = all_msats[all_msats.position == mutSite];

matchingMut = msatsAtSite[msatsAtSite.tag == repeats];

if (matchingMut.size() == 1)

{

genome.removeMutations(msat);

genome.addMutations(matchingMut);

next;

}

}

// make a new mutation with the new repeat count

genome.removeMutations(msat);

msat = genome.addNewDrawnMutation(m2, mutSite);

msat.tag = repeats;

}

}

}

return T;

}

// taking whole population and calculating the Heterozygosity

1000 late() {

// print frequency information for each microsatellite site

all_msats = sim.mutationsOfType(m2);

i = 0.;

He_sum = 0.;

for (pos in sort(msatPositions))

{

catn("Microsatellite at " + pos + ":");

msatsAtPos = all_msats[all_msats.position == pos];

for (msat in sortBy(msatsAtPos, "tag")){

catn(" variant with " + msat.tag + " repeats: " +

sim.mutationFrequencies(NULL, msat));

He_locus = 1. - sum(sim.mutationFrequencies(NULL, msatsAtPos)^2) ;

}

i = i+1.;

print(He_locus);

He_sum = He_sum + He_locus;

}

print(He_sum + "_");

}

// taking a sample and calculating the heterozygosity

1001 late() {

// print frequency information for each microsatellite site

sampled_ind = p1.sampleIndividuals(5);

all_msats = sampled_ind.genomes.mutationsOfType(m2);

i = 0.;

He_sum = 0.;

for (pos in sort(msatPositions))

{

catn("Microsatellite at " + pos + ":");

msatsAtPos = all_msats[all_msats.position == pos];

for (msat in sortBy(msatsAtPos, "tag")){

catn(" variant with " + msat.tag + " repeats: " +

sim.mutationFrequencies(NULL, msat));

He_locus = 1. - sum(sim.mutationFrequencies(sampled_ind, msatsAtPos)^2) ;

}

i = i+1.;

print(He_locus);

He_sum = He_sum + He_locus;

}

print(He_sum + "_");

}

Regards,

Ravi

[Ravi Vishwakarma]

Hi Ben!

My apologies! Here is the code that works in SLiMgui

// Keywords:

initialize() {

defineConstant("L", 1e6); // chromosome length

defineConstant("msatCount", 10); // number of microsats

defineConstant("msatMu", 0.0001); // mutation rate per microsat

defineConstant("msatUnique", T); // T = unique msats, F = lineages

initializeMutationRate(1e-7);

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

initializeGenomicElementType("g1", m1, 1.0);

initializeGenomicElement(g1, 0, L-1);

initializeRecombinationRate(1e-8);

// microsatellite mutation type; also neutral, but magenta

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

m2.convertToSubstitution = F;

m2.color = "#900090";

}

1 late() {

sim.addSubpop("p1", 500);

// create some microsatellites at random positions

genomes = sim.subpopulations.genomes; 

positions = rdunif(msatCount, 0, L-1); 

repeats = rpois(msatCount, 20) + 5; 

for (msatIndex in 0:(msatCount-1)) 

{

pos = positions[msatIndex]; 

mut = genomes.addNewDrawnMutation(m2, pos); 

1000 late() {

// print frequency information for each microsatellite site

all_msats = sim.mutationsOfType(m2);

i = 0.;

He_sum = 0.;

for (pos in sort(msatPositions))

{

catn("Microsatellite at " + pos + ":");

msatsAtPos = all_msats[all_msats.position == pos];

for (msat in sortBy(msatsAtPos, "tag")){

catn(" variant with " + msat.tag + " repeats: " +

sim.mutationFrequencies(NULL, msat));

He_locus = 1. - sum(sim.mutationFrequencies(NULL, msatsAtPos)^2) ;

}

i = i+1.;

print(He_locus);

He_sum = He_sum + He_locus;

}

print(He_sum + "_");

}

// taking a sample and calculating the heterozygosity

1001 late() {

// print frequency information for each microsatellite site

sampled_ind = p1.sampleIndividuals(5);

all_msats = sampled_ind.genomes.mutationsOfType(m2);

i = 0.;

He_sum = 0.;

for (pos in sort(msatPositions))

{

catn("Microsatellite at " + pos + ":");

msatsAtPos = all_msats[all_msats.position == pos];

for (msat in sortBy(msatsAtPos, "tag")){

catn(" variant with " + msat.tag + " repeats: " +

sim.mutationFrequencies(NULL, msat));

He_locus = 1. - sum(sim.mutationFrequencies(sampled_ind, msatsAtPos)^2) ;

}

i = i+1.;

print(He_locus);

He_sum = He_sum + He_locus;

}

print(He_sum + "_");

}

Regards,

Ravi

On Thursday 13 June 2024 at 12:44:53 UTC+1 Ben Haller wrote:

[Ben Haller]

So, Ravi is having trouble posting to the group (if you have difficulties, I recommend using the web interface for slim-discuss rather than posting by email; but Google Groups is just really flaky these days :->).  So I'm now posting the resolution of the thread, which happened off-list.  The error message he was getting was:

"ERROR (EidosCallSignature::CheckArgument): argument 1 cannot be object element type Individual for method mutationFrequencies(); expected object element type Subpopulation."

And I wrote:

Hi Ravi!  OK, just from the error message I think I know what the problem is.  As the error says, mutationFrequencies() expects a subpopulation as its first argument; it provides the frequency of a mutation or mutations within a given set of subpopulations (see the doc for more info).  If you want the frequency within a set of genomes instead – like the two genomes of an individual, to get a frequency of 0.0, 0.5, or 1.0 – use the Genome method mutationFrequenciesInGenomes() instead.  (The doc for mutationFrequencies() points you to that, in fact; it's always a good idea to check the doc when you're puzzled by something!  :->)

Just for posterity, in case others are reading this thread.  :->

Cheers,
-B.

Benjamin C. Haller
Messer Lab
Cornell University