Skip to content

Commit

Permalink
further tweaks to push production scale sims through
Browse files Browse the repository at this point in the history
  • Loading branch information
silastittes committed Jun 13, 2024
1 parent 5d5dec4 commit bb13f86
Show file tree
Hide file tree
Showing 4 changed files with 17 additions and 6 deletions.
2 changes: 1 addition & 1 deletion workflows/config/snakemake/production_config_HomSap.yml
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
# snakemake settings
"seed": 12345
"replicates": 3
"output_dir": "results"
"output_dir": "HomSap_results"

# species-specific settings
"species": "HomSap"
Expand Down
8 changes: 5 additions & 3 deletions workflows/dfe.snake
Original file line number Diff line number Diff line change
Expand Up @@ -159,11 +159,11 @@ rule dadi_infer_dfe:
opts = 100,
prefix = output_dir + "/inference/{demog}/dadi/{dfes}/{annots}/{seeds}/pop{ids}/pop{ids}.two_epoch.gamma"
resources:
time=2000
time=4000
threads: 8
shell:
"""
dadi-cli InferDFE --fs {input[0]} --cache1d {input[1]} --demo-popt {input[2]} --output-prefix {params.prefix} --pdf1d {params.dfe} --p0 {params.dfe_p0} --ubounds {params.dfe_ubounds} --lbounds {params.dfe_lbounds} --ratio {params.ratio} --optimizations {params.opts} --cpus {threads} --nomisid --force-convergence 1
dadi-cli InferDFE --fs {input[0]} --cache1d {input[1]} --demo-popt {input[2]} --output-prefix {params.prefix} --pdf1d {params.dfe} --p0 {params.dfe_p0} --ubounds {params.dfe_ubounds} --lbounds {params.dfe_lbounds} --ratio {params.ratio} --optimizations {params.opts} --cpus {threads} --nomisid --force-convergence 0

"""

Expand Down Expand Up @@ -437,7 +437,7 @@ rule plot_results:

if wildcards.demog == 'Constant':
model = stdpopsim.PiecewiseConstantSize(species.population_size)
mutation_rate = 1.29e-08
mutation_rate = species.genome.mean_mutation_rate

dadi_bestfits = [ b for b in dadi_bestfits if 'pop0' in b ]
polydfe_bestfits = [ b for b in polydfe_bestfits if 'pop0' in b ]
Expand All @@ -446,6 +446,8 @@ rule plot_results:
else:
model = species.get_demographic_model(wildcards.demog)
mutation_rate = model.mutation_rate
if mutation_rate is None:
mutation_rate = species.genome.mean_mutation_rate

pop_names = [ model.populations[i].name for i in range(len(model.populations)) ]
plots.plot_all_dfe_results([dadi_bestfits, polydfe_bestfits, dfe_alpha_bestfits, grapes_bestfits], output, mutation_rate, seq_len, 0.7, pop_names)
6 changes: 5 additions & 1 deletion workflows/n_t.snake
Original file line number Diff line number Diff line change
Expand Up @@ -160,6 +160,8 @@ def mutation_rate_helper(demog, species):
mutation_rate = species.genome.mean_mutation_rate
else:
mutation_rate = species.get_demographic_model(demog).mutation_rate
if mutation_rate is None:
mutation_rate = species.genome.mean_mutation_rate
return mutation_rate

rule write_bdd:
Expand Down Expand Up @@ -345,7 +347,7 @@ rule run_msmc:
output_dir + "/inference/{demog}/msmc/{dfes}/{annots}/{seeds}/{pops}/{samps}.trees.multihep.txt.final.txt"
threads: 8
resources:
mem_mb=48000
mem_mb = lambda wildcards, attempt: attempt * 64000
run:
inputs = expand(output_dir + "/inference/{demog}/msmc/{dfes}/{annots}/{seeds}/{pops}/{chrms}.trees.multihep.txt",
demog=wildcards.demog,
Expand Down Expand Up @@ -459,6 +461,8 @@ rule gone_prep_inputs:
output_dir + "/inference/{demog}/gone/{dfes}/{annots}/{seeds}/{pops}/gone.ped",
output_dir + "/inference/{demog}/gone/{dfes}/{annots}/{seeds}/{pops}/gone.map",
threads: 1
resources:
mem_mb=36000
run:
inputs = expand(
output_dir
Expand Down
7 changes: 6 additions & 1 deletion workflows/simulation.snake
Original file line number Diff line number Diff line change
Expand Up @@ -60,14 +60,19 @@ rule simulation:
input:
output:
output_dir + "/simulated_data/{demog}/{dfes}/{annots}/{seeds}/sim_{chrms}.trees"
resources: time=3000, mem_mb=10000
resources:
time = lambda wildcards, attempt: attempt * 3000,
mem_mb = lambda wildcards, attempt: attempt * 50000
run:
if wildcards.demog == 'Constant':
model = stdpopsim.PiecewiseConstantSize(species.population_size)
mutation_rate = species.genome.mean_mutation_rate
else:
model = species.get_demographic_model(wildcards.demog)
mutation_rate = model.mutation_rate
if mutation_rate is None:
mutation_rate = species.genome.mean_mutation_rate

samples = {f"{model.populations[i].name}": m for i, m in enumerate(demo_sample_size_dict[wildcards.demog])}
genetic_map_id = config["genetic_map"]
contig = species.get_contig(wildcards.chrms, genetic_map=genetic_map_id)
Expand Down

0 comments on commit bb13f86

Please sign in to comment.