From 407aa9d6ee846a049ec910ef3ccb05bea99b4b93 Mon Sep 17 00:00:00 2001 From: Niema Moshiri Date: Wed, 14 Sep 2022 16:57:04 -0700 Subject: [PATCH] Added models --- .tests/sequence_evolution_GTR-Codon.json | 1 + global.json | 4 ---- plugins/sequence_evolution/__init__.py | 1 + plugins/sequence_evolution/seqgen.py | 10 +++++++--- 4 files changed, 9 insertions(+), 7 deletions(-) create mode 100644 .tests/sequence_evolution_GTR-Codon.json diff --git a/.tests/sequence_evolution_GTR-Codon.json b/.tests/sequence_evolution_GTR-Codon.json new file mode 100644 index 0000000..8c3a45d --- /dev/null +++ b/.tests/sequence_evolution_GTR-Codon.json @@ -0,0 +1 @@ +{"Contact Network": {"model": "Complete", "param": {"n": 100}}, "Transmission Network": {"model": "Susceptible-Infected-Removed (SIR)", "param": {"N_S": 95, "N_I": 5, "N_R": 0, "R_S-I": 0.0, "R_S-I_I": 2.0, "R_I-R": 1.0, "duration": 10.0}}, "Sample Times": {"model": "State Entry (Initial)", "param": {"sampled_states": "R"}}, "Viral Phylogeny (Seeds)": {"model": "Yule", "param": {"height": 0.1}}, "Viral Phylogeny (Transmissions)": {"model": "Transmission Tree", "param": {}}, "Mutation Rates": {"model": "Constant", "param": {"rate": 5.0}}, "Sequence Evolution": {"model": "General Time-Reversible (GTR) + Codon", "param": {"p_A": 0.25, "p_C": 0.25, "p_G": 0.25, "p_T": 0.25, "r_A-C": 1.0, "r_A-G": 1.0, "r_A-T": 1.0, "r_C-G": 1.0, "r_C-T": 1.0, "r_G-T": 1.0, "r_site1": 2.0, "r_site2": 3.0, "r_site3": 4.0}}, "Ancestral Sequence": {"model": "Exact Base Frequencies", "param": {"k": 100, "p_A": 0.1, "p_C": 0.2, "p_G": 0.3, "p_T": 0.4}}} diff --git a/global.json b/global.json index 561e4a1..50b1036 100644 --- a/global.json +++ b/global.json @@ -1370,10 +1370,6 @@ "r_site3": { "TYPE": "non-negative float", "DESC": "Relative rate of substitution at codon site 3" - }, - "prop_invariable": { - "TYPE": "probability", - "DESC": "Proportion of invariable sites" } }, "DESC": "Sequences evolve according to a General Time-Reversible (GTR)\nsubstitution model with codon-specific rate heterogeneity." diff --git a/plugins/sequence_evolution/__init__.py b/plugins/sequence_evolution/__init__.py index 6a4fe35..e279fb6 100644 --- a/plugins/sequence_evolution/__init__.py +++ b/plugins/sequence_evolution/__init__.py @@ -3,6 +3,7 @@ from . import seqgen PLUGIN_FUNCTIONS = { "General Time-Reversible (GTR)": seqgen.seqgen_gtr, + "General Time-Reversible (GTR) + Codon": seqgen.seqgen_gtr_codon, "General Time-Reversible (GTR) + Gamma": seqgen.seqgen_gtr_gamma, "None": DUMMY_PLUGIN_FUNC, } diff --git a/plugins/sequence_evolution/seqgen.py b/plugins/sequence_evolution/seqgen.py index de29e32..23994ec 100644 --- a/plugins/sequence_evolution/seqgen.py +++ b/plugins/sequence_evolution/seqgen.py @@ -20,11 +20,11 @@ def seqgen(mode, params, out_fn, config, verbose=True): seqgen_log_fn = "%s/seqgen.log" % out_fn['intermediate'] f = open(seqgen_tree_fn, 'w'); f.write("1 %d\nROOT %s\n1\n%s" % (len(root_seq),root_seq,treestr)); f.close() command = ['seq-gen', '-of', '-k1'] - if mode in {'GTR', 'GTR+G'}: # add model + if mode in {'GTR', 'GTR+G', 'GTR+Codon'}: # GTR model command += ['-m', 'GTR'] - if mode in {'GTR', 'GTR+G'}: # add base frequencies + if mode in {'GTR', 'GTR+G', 'GTR+Codon'}: # add base frequencies command += (['-f'] + [str(params['p_%s' % n]) for n in 'ACGT']) - if mode in {'GTR', 'GTR+G'}: # add GTR transition rates + if mode in {'GTR', 'GTR+G', 'GTR+Codon'}: # add GTR transition rates command += (['-r'] + [str(params['r_%s-%s' % tuple(pair)]) for pair in ['AC', 'AG', 'AT', 'CG', 'CT', 'GT']]) if mode in {'GTR', 'GTR+G'}: # add proportion invariable command += ['-i', str(params['prop_invariable'])] @@ -32,6 +32,8 @@ def seqgen(mode, params, out_fn, config, verbose=True): command += ['-a', str(params['alpha'])] if params['num_cats'] > 0: command += ['-g', str(params['num_cats'])] + if mode in {'GTR+Codon'}: # add Codon parameters + command += ['-c', str(params['r_site1']), str(params['r_site2']), str(params['r_site3'])] command += [seqgen_tree_fn] if verbose: print_log("Command: %s" % ' '.join(command)) @@ -44,5 +46,7 @@ def seqgen(mode, params, out_fn, config, verbose=True): # model-specific functions def seqgen_gtr(params, out_fn, config, verbose=True): seqgen("GTR", params, out_fn, config, verbose=verbose) +def seqgen_gtr_codon(params, out_fn, config, verbose=True): + seqgen("GTR+Codon", params, out_fn, config, verbose=verbose) def seqgen_gtr_gamma(params, out_fn, config, verbose=True): seqgen("GTR+G", params, out_fn, config, verbose=verbose)