Saving results into a tsv file
Yujean Kim
//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
