Saving results into a tsv file

29 views
Skip to first unread message

Yujean Kim

unread,
Aug 12, 2021, 6:52:12 AM8/12/21
to slim-discuss
Hi,

Currently I have modeled a gene suppression scenario in snails by using the following code:

//Growth rate parameter:

initialize() {

defineConstant("cut_rate", 0.98); //rate of double strand break occuring

defineConstant("resistance_rate", 0.078); //rate of resistance sequence forming

defineConstant("mutation_rate", 0); //drive specific mutation rate

defineConstant("nonres", 10000); //non-resistant allele (target population)

defineConstant("capacity" , 10000);

defineConstant("drive_size", 1000); //gene drive release (10% of population)

defineConstant("growth_rate", 53);

defineConstant("max_offspring", 14000);

defineConstant("mortality", c(1,0));

defineConstant("selfing_rate", 0.8);

defineConstant("inbreeding_depression", 0.063);

initializeSLiMModelType("nonWF");

//Assume no fitness cost

awt = initializeMutationType("m1", 0.5, "f", 0.0); //wildtype fertility allele

adr = initializeMutationType("m2", 0.5, "f", 0.0); //mutant fertility allele (drive allele)

ar2 = initializeMutationType("m3", 0.5, "f", 0.0); //resistant allele

avariations = c(awt, adr, ar2);

//Defining different alleles into one locus

initializeGenomicElementType("g1", avariations, c(1, 1, 1)); //locus a

//Defining each locus into a chromosome position

initializeGenomicElement(g1, 0, 0); //chromsosome position 0

avariations.mutationStackPolicy = "l"; //makes sure last mutation added is the most recent

avariations.mutationStackGroup = 1; //makes sure there is only one mutation/allele in one locus

initializeMutationRate(0);

initializeRecombinationRate(0.5);



}


//Reproduction rules - ones with no wild type alleles are infertile (doublesex genes are haplosufficient)

reproduction() {

//Calculation for carrying capacity

capacity_fitness_scaling = growth_rate / (((growth_rate-1) * (p1.individualCount / capacity)) + 1);

p = capacity_fitness_scaling * 2 / max_offspring ;

num_offspring = rbinom(1, max_offspring, p); //number of offspring depends on population and carrying capacity

for (i in seqLen(num_offspring)) {

if (subpop.size())

if (runif(1) < selfing_rate)

offspring = subpop.addSelfed(individual); //Add selfing rate

else

offspring = subpop.addCrossed(individual, subpop.sampleIndividuals(1)); //Add outcrossing rate

}

}


//Gene drive implementation in offspring in the germline

//convert to single gene drive


1: modifyChild() {

parent = c(parent1, parent2);

child_chromosome = c(childGenome1, childGenome2);

for (i in c(0,1))

if (sum(parent[i].genomes.countOfMutationsOfType(m2)) & sum(child_chromosome[i].countOfMutationsOfType(m1))) {

if ((runif(1) < cut_rate) & (runif(1) > resistance_rate)){

if (runif(1) < mutation_rate)

child_chromosome[i].addNewDrawnMutation(m3, 0);

else

child_chromosome[i].addNewDrawnMutation(m2, 0);

}

}

if (sum(parent[i].genomes.countOfMutationsOfType(m2)) & sum(child_chromosome[i].countOfMutationsOfType(m1))) {

if ((runif(1) < cut_rate) & (runif(1) < resistance_rate))

childGenome2.addNewDrawnMutation(m3, 0);

}

return T;

if (individual.countOfMutationsOfType(m2) == 2) {

if (parent1 == parent2)

return F;

return T;

}

if (individual.countOfMutationsOfType(m3) == 2) {

if (parent1 == parent2)

return F;

return T;

}

if (individual.countOfMutationsOfType(m1) == 0) {

if (parent1 == parent2)

return F;

return T;

}

all = sim.subpopulations.individuals;

for (ind in all){

if (parent1 == parent2){

if (runif(1) < inbreeding_depression)

ind.fitnessScaling = 1 - inbreeding_depression;

}

}


}


//fitness scaling - discrete generations (non-overlapping generations)


late() {

all = sim.subpopulations.individuals;

for (ind in all) {

age_mortality_rate = mortality[ind.age];

ind.fitnessScaling = 1 - age_mortality_rate;

}

}


//introduce of drive allele


1 {

sim.addSubpop("p1", nonres);

p1.individuals.genomes.addNewDrawnMutation(m1, 0);

sim.addSubpop("p2", drive_size);

p2.individuals.genomes.addNewDrawnMutation(m2, 0);

drop = p2.individuals;

p1.takeMigrants(drop);


}


late() {

//Sum of each mutation type in population

all = sim.subpopulations.individuals;

num_awt = sum(all.genomes.countOfMutationsOfType(m1) > 0);

num_adr = sum(all.genomes.countOfMutationsOfType(m2) > 0);

num_ar2 = sum(all.genomes.countOfMutationsOfType(m3) > 0);

//Propotion of each mutation type in population

rate_awt = num_awt / (2 * size(all));

rate_adr = num_adr / (2 * size(all));

rate_ar2 = num_ar2 / (2 * size(all));

//output

cat(asInteger(sim.generation) + "\t " + rate_awt + "\t " + rate_adr + "\t " + rate_ar2 + "\t " + size(all) + "\n");

if (rate_awt == 0.0) {

catn("OUTCOME: DRIVE SUCCESS");

sim.simulationFinished();

}

if (rate_adr == 0.0) {

catn("OUTCOME: DRIVE LOST");

sim.simulationFinished();

}

log = sim.createLogFile("~/Desktop/haplosufficient_10.txt", sep="\t");

log.addGeneration();

log.addCustomColumn("wildtype_allele", "rate_awt;");

log.addCustomColumn("mutant_allele", "rate_adr;");

log.addCustomColumn("resistant_allele", "rate_ar2;");

}


//// MAIN END CONDITION.

111 late() {

catn("OUTCOME: DRIVE PERSISTS");

sim.simulationFinished();

}

Currently I want to produce a tsv file which shows the proportion of wildtype alleles, mutant alleles and resistant alleles per generation. I tried to do this by using the following code:

log = sim.createLogFile("~/Desktop/haplosufficient_10.txt", sep="\t");

log.addGeneration();

log.addCustomColumn("wildtype_allele", "rate_awt;");

log.addCustomColumn("mutant_allele", "rate_adr;");

log.addCustomColumn("resistant_allele", "rate_ar2;");


However, this produces a txt file which is empty. I was wondering what the problem was!


Thank you for your help!


Kind regards,

Yujean Kim

Ben Haller

unread,
Aug 12, 2021, 10:07:57 AM8/12/21
to slim-discuss
Hi Yujean.  There are several issues here.  See section 4.2.5 of the SLiM manual for an example of how to set up logging with LogFile.  More details:

First of all, it looks like you're setting up the log file in every generation, in your late() event.  You should set it up just once, at the start of the simulation, and then it will automatically log out its info every generation that it's scheduled to do so.  See section 4.2.5.

Second, you're not actually having the LogFile log anything out.  It can work in two ways: you can set it up with a logging interval, as in section 4.2.5, and then it logs out periodically automatically, or you can call logRow() explicitly to log when you want to.  (These two approaches are not mutually exclusive.)  You're not doing either, so it never actually tries to log anything.

Third, if it *did* try to log, it would fail, because the custom columns you've set up would error out.  You're trying to log out the values of locally defined variables, but by the time the LogFile logs, those variables will no longer exist, and you'll get an "undefined identifier" error.  The code for custom columns gets executed in its own scope, and can't use variables that were defined at the time the log file was configured.  (Since you configure it just once and then it can log in future generations, at different points in the generation cycle, it has to work this way.)  So you could either (a) make the values that your late() event calculates be global variables, using defineGlobal(), so that they exist later on when the log file does its thing, or (b) make user-defined functions that calculate the things you want to log out, and set up your custom columns to call those user-defined functions.  Either approach will work; if you use approach (a) just understand that the values logged out will be the last values set into the global variables by your code, even if that value is now out of date; if you want a new, current value to be calculated dynamically at the moment the log file logs, approach (b) would be better.  But as long as you set new values into the globals in every generation in a late() event, approach (a) ought to work fine.

I hope this helps clarify things!

Cheers,
-B.

Benjamin C. Haller
Messer Lab
Cornell University


Reply all
Reply to author
Forward
0 new messages