To receive notifications about scheduled maintenance, please subscribe to the mailing-list gitlab-operations@sympa.ethz.ch. You can subscribe to the mailing-list at https://sympa.ethz.ch

Commit 656ec2a7 authored by Eva Bons's avatar Eva Bons

removed outdated scenarios

parent 5a25a06a
import SeqSimEvo as seq_sim
import numpy as np
from collections import Counter
from copy import deepcopy
from tqdm import tqdm
class Simulation(seq_sim.Simulation):
def __init__(self,interactionTable=None,interactionFraction=0.3, ratio=0.5,**kwargs):
seq_sim.Simulation.__init__(self,**kwargs)
if interactionTable is None:
self.interactionTable = self.generateInteractions(fraction=interactionFraction, ratio=ratio)
else:
self.interactionTable = interactionTable
def transform_fitnessTable(self,fraction=0.3,transform=['inc',0.05]):
n_to_change = int(len(self.sequence)*3*fraction)
#print self.sequence.sequence
options = [(i,x) for i in range(len(self.sequence)) for x in range(4) if x != self.sequence[i]]
indx = np.random.choice(range(len(options)),n_to_change,replace=False)
for i in indx:
old_value = self.fitness_table[options[i][1],options[i][0]]
if transform[0] == 'product':
self.fitness_table[options[i][1],options[i][0]] = old_value * transform[1]
elif transform[0] == 'add':
self.fitness_table[options[i][1],options[i][0]] = old_value + transform[1]
if self.fitness_table[options[i][1],options[i][0]] < 0:
self.fitness_table[options[i][1],options[i][0]] = 0
elif transform[0] == 'set':
self.fitness_table[options[i][1],options[i][0]] = transform[1]
elif transform[0] == 'scramble':
self.fitness_table[options[i][1],options[i][0]] = np.random.choice(self.fitness_table.flatten())
else:
raise ValueError('transform[0] must be one of "product","add","set","scramble"')
def generateInteractions(self,fraction,ratio=0.5):
n_interactions = int(len(self.sequence)*(len(self.sequence)-1)*fraction)
combis = {}
for i in tqdm(range(n_interactions)):
while True:
combi = [[np.random.randint(len(self.sequence)),np.random.randint(4)] for i in range(2)]
if combi[0][0] != combi[1][0]:
break
if combi[0][0]>combi[1][0]:
combi = [combi[1],combi[0]]
if np.random.random() < ratio:
#combi.append(1)
fit1 = self.get_fitness_effect(combi[0][0],combi[0][1])
fit2 = self.get_fitness_effect(combi[1][0],combi[1][1])
cur = fit1*fit2
new = 1+np.random.exponential(self.settings['parameters']['lb'])
if new>cur:
combi.append(new)
else:
combi.append(cur)
else:
combi.append(0)
try:
combis[tuple(combi[0])].append(combi[1:])
except KeyError:
combis[tuple(combi[0])] = [combi[1:]]
return combis
def get_nr_offspring(self, sequence_id, return_fitness=False):
"""returns the number of offspring of a sequence according to the fitness
of that sequence"""
R0 = self.settings['R0']
changes = self.current_gen.get_seq(sequence_id)
fitness = 1
if changes is not None:
#for pos, base in zip(changes[:, 0], changes[:, 1]):
# fitness *= (self.fitness_table[int(base), int(pos)])
sort = np.argsort(changes[:,0])
already_done = []
for i,change1 in enumerate(changes[sort]):
if tuple(change1) in self.interactionTable:
ch1 = [tuple(c[0]) for c in self.interactionTable[tuple(change1)]]
chs2 = set([tuple(c) for c in changes[sort][i+1:]])
interactions = set(ch1).intersection(chs2)
for interaction in interactions:
already_done += [tuple(change1),interaction]
value = [self.interactionTable[tuple(change1)][i][-1] for i in range(len(ch1)) if ch1[i]==interaction]
fitness*=np.product(value)
done=False
if len(already_done)>0:
if tuple(change1) in already_done:
done=True
if not done:
fitness*=(self.fitness_table[int(change1[1]), int(change1[0])])
if return_fitness:
return np.random.poisson(fitness*R0), fitness
return np.random.poisson(fitness*R0)
def copy(self,name):
return Simulation(simulation_settings=deepcopy(self.settings), sequence = self.sequence,
fitness_table=deepcopy(self.fitness_table),name= name,
interactionTable=deepcopy(self.interactionTable))
class Population(seq_sim.Population):
def __init__(self,**kwargs):
seq_sim.Population.__init__(self,**kwargs)
class Seq(seq_sim.Seq):
def __init__(self,**kwargs):
seq_sim.Seq.__init__(self,**kwargs)
def analyze_passage(sim,f,passage):
try:
changes = np.vstack(list(sim.current_gen.changes.values()))
except ValueError:
return 0
n_seq = float(sim.n_seq)
changes = [tuple(i) for i in changes]
counts = Counter(changes)
for i in counts:
f.write("{}\t{}\t{}\t{}\t{}\n".format(i[0], i[1], counts[i]/n_seq,
passage,sim.settings['name']))
def run_custom():
interactionfile = sys.argv[1]
resultfile = sys.argv[2]
timingsfile = sys.argv[3]
interactionFraction = float(sys.argv[4])
ratio = float(sys.argv[5])
transformType = sys.argv[6]
transformValue = float(sys.argv[7])
fraction = float(sys.argv[8])
env1_1 = Simulation(name='env1_1',interactionFraction=interactionFraction,
ratio=ratio,seq_len=500)
sims = [env1_1]+[env1_1.copy('env1_'+str(i)) for i in range(2,6)]
env2_1 = env1_1.copy('env2_1')
env2_1.transform_fitnessTable(transform=[transformType,transformValue],fraction=fraction)
sims += [env2_1]+[env2_1.copy('env2_'+str(i)) for i in range(2,6)]
with open(interactionfile,'w') as f:
pickle.dump(env1_1.interactionTable,f)
f = open(resultfile,'w',buffering=1)
t = open(timingsfile,'w',buffering=1)
for i in range(201):
for sim in sims:
if i%2 == 0:
sim.settings['max_pop'] = 100
sim.new_generation()
sim.settings['max_pop'] = 100000
if i%10 == 0:
t.write("{} {}\n".format(time.time(), i))
for sim in sims:
analyze_passage(sim, f,i)
f.close()
t.close()
if __name__ == '__main__':
import time
import pickle
import sys
interactionfile = sys.argv[1]
resultfile = sys.argv[2]
timingsfile = sys.argv[3]
env1_1 = Simulation(name='env1_1',interactionFraction=0.05,
ratio=0,seq_len=500)
sims = [env1_1]+[env1_1.copy('env1_'+str(i)) for i in range(2,6)]
env2_1 = env1_1.copy('env2_1')
#nv2_1.generateInteractions(float(sys.argv[4]),1)
#env2_1.transform_fitnessTable(transform=['scramble',0],fraction=0.3)
#env2_1.transform_fitnessTable(transform=['product',1.1],fraction=1)
#env2_1.transform_fitnessTable(transform=['set',0],fraction=0.3)
#env2_1.transform_fitnessTable(transform=['set',1],fraction=0.3)
sims += [env2_1]+[env2_1.copy('env2_'+str(i)) for i in range(2,6)]
with open(interactionfile,'wb') as f:
pickle.dump(env1_1.interactionTable,f)
f = open(resultfile,'w',buffering=1)
t = open(timingsfile,'w',buffering=1)
for i in range(10):
for sim in sims:
if i%2 == 0:
if 'env2' in sim.settings['name']:
sim.settings['max_pop'] = int(sys.argv[4])
else:
sim.settings['max_pop'] = 1000
sim.new_generation()
sim.settings['max_pop'] = 10000
if i%10 == 0:
t.write("{} {}\n".format(time.time(), i))
for sim in sims:
analyze_passage(sim, f,i)
f.close()
t.close()
from __future__ import print_function
import SeqSimEvo as seq_sim
import numpy as np
from collections import Counter
from copy import deepcopy
from tqdm import tqdm
import sys
class Simulation(seq_sim.Simulation):
def __init__(self,epistasis={}, ratio=0.5,fitness_table=None,**kwargs):
seq_sim.Simulation.__init__(self,**kwargs)
if fitness_table is None:
self.fitness_table = self.__getfitness_table()
else:
self.fitness_table = fitness_table
self.settings['epistasis'] = epistasis
def __getfitness_table(self):
return {}
def assign_value(self,changes):
if len(changes) == 1: #single mutation: draw from MFED
if self.settings['model'] == 'exponential':
pars = self.settings['parameters']
random_nr = np.random.random()
if random_nr < pars['fl']: #lethal mutation
return 0
elif random_nr < (pars['fl']+pars['fd']):
fit = 1-np.random.exponential(pars['ld'])
if fit<0:
fit=0
return fit
elif random_nr > (pars['fl']+pars['fd']+pars['fn']):
return 1+np.random.exponential(pars['lb'])
if self.settings['model'] == 'neutral':
return 1
if self.settings['model'] == 'lognormal':
random_nr = np.random.random()
pars = self.settings['parameters']
if random_nr < pars['fl']: #lethal mutation
return 0
else:
return np.random.lognormal(mean=pars['mu'],sigma=pars['sigma'])
else:
raise('model {} not implemented'.format(self.settings['model']))
else:
epi = self.settings['epistasis']
random_nr = np.random.random()
if epi['model'] == 'exponential':
random_nr = np.random.random()
if random_nr < epi['par']['fl']:
e = 0
elif random_nr < epi['par']['fl']+epi['par']['fd']:
e = np.exp(-np.random.exponential(epi['par']['ld']))
else:
e = np.exp(np.random.exponential(epi['par']['lb']))
if epi['model'] == 'none':
e = 1
else:
print('model {} not implemented'.format(epi['model']))
raise()
if len(changes) == 2:
str1 = '{}-{}'.format(changes[0][0],changes[0][1])
str2 = '{}-{}'.format(changes[1][0],changes[1][1])
try:
f1 = self.fitness_table[str1]
except KeyError:
self.fitness_table[str1] = self.assign_value([changes[0]])
f1 = self.fitness_table[str1]
try:
f2 = self.fitness_table[str2]
except KeyError:
self.fitness_table[str2] = self.assign_value([changes[1]])
f2 = self.fitness_table[str2]
else:
found=False
for c in range(len(changes)):
subset = np.vstack([changes[:c],changes[c+1:]])
changes_str = ['{}-{}'.format(i[0],i[1]) for i in subset]
changes_str = ' '.join(sorted(changes_str))
if changes_str in self.fitness_table.keys():
f1 = self.fitness_table[changes_str]
found=True
break
if not found:
self.fitness_table[changes_str] = self.assign_value(subset)
f1 = self.fitness_table[changes_str]
try:
f2 = self.fitness_table['{}-{}'.format(changes[c][0],changes[c][1])]
except KeyError:
self.fitness_table['{}-{}'.format(changes[c][0],changes[c][1])] = self.assign_value([changes[c]])
f2 = self.fitness_table['{}-{}'.format(changes[c][0],changes[c][1])]
return e*f1*f2
def get_nr_offspring(self, sequence_id, return_fitness=False):
"""returns the number of offspring of a sequence according to the fitness
of that sequence"""
R0 = self.settings['R0']
changes = self.current_gen.get_seq(sequence_id)
if changes is None:
fitness = 1
else:
changes_str = ['{}-{}'.format(i[0],i[1]) for i in changes]
changes_str = ' '.join(sorted(changes_str))
try:
fitness = self.fitness_table[changes_str]
except KeyError:
self.fitness_table[changes_str] = self.assign_value(changes)
fitness = self.fitness_table[changes_str]
if return_fitness:
return np.random.poisson(R0*fitness), fitness
else:
return np.random.poisson(R0*fitness)
def copy(self,name):
return Simulation(simulation_settings=deepcopy(self.settings),
epistasis=self.settings['epistasis'],
sequence = self.sequence,
fitness_table=deepcopy(self.fitness_table),
name= name)
class Population(seq_sim.Population):
def __init__(self,**kwargs):
seq_sim.Population.__init__(self,**kwargs)
class Seq(seq_sim.Seq):
def __init__(self,**kwargs):
seq_sim.Seq.__init__(self,**kwargs)
if __name__ == '__main__':
n_gen = 200
pars = {'fl': 0.1,'mu': -0.25,'sigma': 0.15}
epistasis = {'same':1, 'other':0,
'par':{'fb':0,'lb':0,'fl':0,'fd':0,'ld':0,'fn':1},
'model':'none'}
sim = Simulation(model='lognormal',
parameters = pars,
epistasis=epistasis,
n_seq_init=1000,
max_pop=10000,
mut_rate=1e-5,
seq_len=1000)
for gen in range(n_gen):
sim.new_generation()
if gen%10 == 0:
print('mean: {}'.format(sim.average_fitness))
from __future__ import print_function
import SeqSimEvo as sim
import numpy as np
master_sim = sim.Simulation(simulation_settings='HIV')
n_gen = 10
#settings of the different compartments
compartment_names = ['sanctuary', 'liver','blood']
compartment_sizes = [1e3, 1e5, 1e2]
initial_sizes = [0, 0, 10]
migration_rates = np.array([[0, 0, 1e-3], #row=from, col = to
[0, 0, 1e-2],
[1e-3, 1e-2, 0]])
replication = [0.2, 1, 0.166]
mut_rate = [0.2, 1, 0]
fasta = ''
#setup different compartments:
compartments = []
for i in range(len(compartment_names)):
compartments.append(master_sim.copy(compartment_names[i], n_seq=initial_sizes[i]))
compartments[-1].settings['max_pop'] = compartment_sizes[i]
compartments[-1].settings['mut_rate'] = master_sim.settings['mut_rate']*mut_rate[i]
compartments[-1].settings['R0'] = master_sim.settings['R0']*replication[i]
#do simulations
for gen in range(n_gen):
#new generation
for comp in compartments:
comp.new_generation(dieout=True)
#do migration
additions = [[] for i in compartment_names]
for comp1 in range(len(compartments)):
for comp2 in range(len(compartments)):
migration = migration_rates[comp1, comp2]
if migration>0:
n_migrate = np.random.poisson(migration*len(compartments[comp1].current_gen))
for i in range(n_migrate):
migrant = np.random.randint(len(compartments[comp1].current_gen))
additions[comp2].append(compartments[comp1].current_gen.get_seq(migrant))
compartments[comp1].current_gen.delete_sequence(migrant)
for i, add in enumerate(additions):
for seq in add:
compartments[i].current_gen.add_sequence(changes = seq)
#liver transplant
# if gen == 50:
# compartments[1].settings['max_pop'] = 0
#
# if gen == 51:
# compartments[1].settings['max_pop'] = compartment_sizes[1]
#sampling
if gen in [4,20,50,99]:#55,70,99]:
fasta += compartments[0].current_gen.to_fasta(n_seq=10,description='-1-{}'.format(gen))
fasta += compartments[1].current_gen.to_fasta(n_seq=10,description='-2-{}'.format(gen))
#fasta += compartments[2].current_gen.to_fasta(n_seq=10,description='-3-{}'.format(gen))
#print some stats
print('gen {}: nseq'.format(gen),end=' ')
for comp in compartments:
print(comp.current_gen.stats()['n_seq'],end=' ')
print('')
with open('3comp.fasta', 'w') as f:
f.write(fasta)
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment