Commit 0096d0f1 by Eva Bons

### wrote extension for binary fission

parent f923b725
 ''' Simulation extension to implement reproduction via binary fission In every generation, cells die or reproduce. R0 indicated the relative chance to die or reproduce: an R0 of 1 equals initial equal chances to die or reproduce, R0 of 2 is twice as likely to reproduce as to die. fitness values modify this probability, as f*R0 All other functionalities (calculation of fitness, mutating sequences, etc.) are as in the original simulation code ''' import simulation as sim from simulation import Population, Seq import random class Simulation(sim.Simulation): def survive(self, fitness): '''decides whether a sequence with a certain fitness survives or not''' p = (1-(1.0/(1+fitness*self.settings['R0']))) #probability of survival = f*R0/(1+f*R0) if random.random() self.settings['max_pop']: for delete in range(n_seq-self.settings['max_pop']): to_del = random.randint(0,n_seq-delete-1) self.current_gen.delete_sequence(to_del) elif n_seq == 0: print 'died out' raise NotImplementedError('population died out') self.n_seq = len(self.current_gen) self.average_fitness = sum(all_fitnesses)/len(all_fitnesses) if __name__ == '__main__': from tqdm import tqdm sim_test = Simulation(R0 = 2,max_pop=10000,n_seq_init=1000) for i in tqdm(range(50)): sim_test.new_generation() print sim_test.current_gen
 ... ... @@ -346,7 +346,7 @@ class Simulation(object): (in the current generation) Returns: Nothing (the sequence is edited in `pop`) True/False (did the sequence mutate or not?) """ #get the number of mutations that will take place try: ... ... @@ -378,6 +378,9 @@ class Simulation(object): if (self.settings['ga_increase'] <= 1) or random.random() < (1.0/self.settings['ga_increase'] ): pop.add_change(seq_id_new, where, new_base) success_mut += 1 return True else: return False def new_generation(self,new_gen=None,dieout=False): """ ... ... @@ -390,7 +393,8 @@ class Simulation(object): new_generation if the population died out (False, default) Returns: Nothing. current_gen and other simulation attributes are updated as needed Nothing. Updates current_gen, effective_pop, gen, average_fitness, and n_seq """ self.effective_pop = 0 self.gen += 1 ... ... @@ -526,7 +530,7 @@ class Population(): pos = j[0] string += '{orig}-{pos}-{to}\t{seq}\t{patient}\n'.format(orig=self.sim.sequence[pos], pos=pos, to=j[1], to=self.sim.sequence.translation[j[1]], seq=i, patient=self.sim.settings['name']) return string ... ... @@ -548,7 +552,7 @@ class Population(): Nothing. Output is printed to stdout. ''' if any(np.array(seq_ids)>self.n_seq): raise IndexError, 'seqID out of range' raise IndexError('seqID out of range') string = '#mutID (from-pos-to)\tsequence\tpatient\n' for i in range(self.n_seq): if i in self.changed and i in seq_ids: ... ... @@ -654,10 +658,10 @@ class Population(): Nothing. Population is changed in-place. ''' if pos > len(self.sim.sequence): raise IndexError, 'Pos {} outside sequence length'.format(pos) raise IndexError('Pos {} outside sequence length'.format(pos)) if seq_id > self.n_seq: raise IndexError, 'SeqID {} outside pop size {} {}'.format(seq_id, self.n_seq,self) raise IndexError('SeqID {} outside pop size {} {}'.format(seq_id, self.n_seq,self)) if seq_id in self.changed: ... ... @@ -775,7 +779,7 @@ class Population(): IndexError: when sequence_id is out of bounds ''' if sequence_id > self.n_seq: raise IndexError, 'sequence_id is out of bounds' raise IndexError('sequence_id is out of bounds') elif sequence_id in self.changed: return self.changes[sequence_id] else: ... ...
 from __future__ import print_function import SeqSimEvo #from SeqSimEvo import binary_fission as SeqSimEvo import pytest import numpy as np from copy import deepcopy, copy #tests written for SeqSimEvo.Seq.__init__ # .__str__ ... ... @@ -375,9 +378,17 @@ def test_Simulation_mutate_seq(): assert next_gen.stats()['total_mutations'] == 3 def test_Simulation_new_generation(): seq = SeqSimEvo.Seq(seq='AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA') sim = SeqSimEvo.Simulation(n_seq_init=3,sequence=seq) seq = SeqSimEvo.Seq(seq_len=5000) sim = SeqSimEvo.Simulation(n_seq_init=10000,sequence=seq,mut_rate=0.001) sim.effective_pop = -1 sim.average_fitness = -1 sim.n_seq = -1 old_gen = deepcopy(sim.current_gen.changes) sim.new_generation() assert sim.gen == 1 assert sim.current_gen.changes != old_gen, 'this test can fail by chance. Be worried if it keeps failing.' assert sim.effective_pop != -1 assert sim.average_fitness != -1 assert sim.n_seq != -1
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!