Enable seed specification

python/gubbins/common.py

@@ -744,13 +744,13 @@ def return_algorithm_choices(args,i):
 def return_algorithm(algorithm_choice, model, input_args, node_labels = None, extra = None):
     initialised_algorithm = None
     if algorithm_choice == "fasttree":
-        initialised_algorithm = FastTree(threads = input_args.threads, model = model, bootstrap = input_args.bootstrap, verbose = input_args.verbose, additional_args = extra)
+        initialised_algorithm = FastTree(threads = input_args.threads, model = model, seed = input_args.seed, bootstrap = input_args.bootstrap, verbose = input_args.verbose, additional_args = extra)
     elif algorithm_choice == "raxml":
-        initialised_algorithm = RAxML(threads = input_args.threads, model = model, bootstrap = input_args.bootstrap, internal_node_prefix = node_labels, verbose = input_args.verbose, additional_args = extra)
+        initialised_algorithm = RAxML(threads = input_args.threads, model = model, seed = input_args.seed, bootstrap = input_args.bootstrap, internal_node_prefix = node_labels, verbose = input_args.verbose, additional_args = extra)
     elif algorithm_choice == "raxmlng":
-        initialised_algorithm = RAxMLNG(threads = input_args.threads, model = model, bootstrap = input_args.bootstrap, internal_node_prefix = node_labels, verbose = input_args.verbose, additional_args = extra)
+        initialised_algorithm = RAxMLNG(threads = input_args.threads, model = model, seed = input_args.seed, bootstrap = input_args.bootstrap, internal_node_prefix = node_labels, verbose = input_args.verbose, additional_args = extra)
     elif algorithm_choice == "iqtree":
-        initialised_algorithm = IQTree(threads = input_args.threads, model = model, bootstrap = input_args.bootstrap, internal_node_prefix = node_labels, verbose = input_args.verbose, use_best = (model is None and input_args.best_model), additional_args = extra)
+        initialised_algorithm = IQTree(threads = input_args.threads, model = model, seed = input_args.seed, bootstrap = input_args.bootstrap, internal_node_prefix = node_labels, verbose = input_args.verbose, use_best = (model is None and input_args.best_model), additional_args = extra)
     elif algorithm_choice == "rapidnj":
         initialised_algorithm = RapidNJ(threads = input_args.threads, model = model, bootstrap = input_args.bootstrap, verbose = input_args.verbose, additional_args = extra)
     elif algorithm_choice == "star":

python/gubbins/run_gubbins.py

@@ -76,6 +76,8 @@ def parse_input_args():
                                                       default = False, action = 'store_true')
     treeGroup.add_argument('--sh-test',               help='Perform an SH test of node likelihoods', default = False,
                                                       action = 'store_true')
+    treeGroup.add_argument('--seed',                  help='Set seed for reproducibility of analysis',
+                                                      default = None, type = int)
 
     modelGroup = parser.add_argument_group('Nucleotide substitution model options')
     modelGroup.add_argument('--model',          '-M', help='Nucleotide substitution model (not all available for all '

python/gubbins/tests/test_dependencies.py

@@ -79,6 +79,18 @@ def test_fasttree(self):
         self.cleanup('multiple_recombinations')
         assert exit_code == 0
 
+    def test_fasttree_seed(self):
+        exit_code = 1
+        parser = run_gubbins.parse_input_args()
+        common.parse_and_run(parser.parse_args(["--tree-builder", "fasttree",
+                                                    "--verbose", "--iterations", "3",
+                                                    "--seed","42",
+                                                    "--threads", "1",
+                                                    os.path.join(data_dir, 'multiple_recombinations.aln')]))
+        exit_code = self.check_for_output_files('multiple_recombinations')
+        self.cleanup('multiple_recombinations')
+        assert exit_code == 0
+
     # Test resuming a default analysis
     def test_fasttree_resume(self):
         exit_code = 1
@@ -122,6 +134,18 @@ def test_iqtree(self):
         self.cleanup('multiple_recombinations')
         assert exit_code == 0
 
+    def test_iqtree_seed(self):
+        exit_code = 1
+        parser = run_gubbins.parse_input_args()
+        common.parse_and_run(parser.parse_args(["--tree-builder", "iqtree",
+                                                "--verbose", "--iterations", "3",
+                                                "--seed","42",
+                                                "--threads", "1",
+                                                os.path.join(data_dir, 'multiple_recombinations.aln')]))
+        exit_code = self.check_for_output_files('multiple_recombinations')
+        self.cleanup('multiple_recombinations')
+        assert exit_code == 0
+
     def test_iqtree_custom_model(self):
         exit_code = 1
         parser = run_gubbins.parse_input_args()
@@ -246,6 +270,18 @@ def test_raxml(self):
         self.cleanup('multiple_recombinations')
         assert exit_code == 0
 
+    def test_raxml_seed(self):
+        exit_code = 1
+        parser = run_gubbins.parse_input_args()
+        common.parse_and_run(parser.parse_args(["--tree-builder", "raxml",
+                                                    "--verbose", "--iterations", "3",
+                                                    "--seed","42",
+                                                    "--threads", "1",
+                                                    os.path.join(data_dir, 'multiple_recombinations.aln')]))
+        exit_code = self.check_for_output_files('multiple_recombinations')
+        self.cleanup('multiple_recombinations')
+        assert exit_code == 0
+
     def test_raxml_custom_model(self):
         exit_code = 1
         parser = run_gubbins.parse_input_args()
@@ -301,6 +337,19 @@ def test_raxmlng(self):
         self.cleanup('multiple_recombinations')
         assert exit_code == 0
 
+    def test_raxmlng_seed(self):
+        exit_code = 1
+        parser = run_gubbins.parse_input_args()
+        common.parse_and_run(parser.parse_args(["--tree-builder", "raxmlng",
+                                                    "--model","GTR",
+                                                    "--verbose", "--iterations", "3",
+                                                    "--seed","42",
+                                                    "--threads", "1",
+                                                    os.path.join(data_dir, 'multiple_recombinations.aln')]))
+        exit_code = self.check_for_output_files('multiple_recombinations')
+        self.cleanup('multiple_recombinations')
+        assert exit_code == 0
+
     def test_raxmlng_custom_model(self):
         exit_code = 1
         parser = run_gubbins.parse_input_args()

python/gubbins/tests/test_utils.py

@@ -107,3 +107,9 @@ def test_verbose_printer(self):
             printer.print(["AAA", "BBB"])
             printed = f.getvalue()
         assert printed == "AAA-BBB\n"
+
+    def test_seed(self):
+        set_seed_val = utils.set_seed(42)
+        assert set_seed_val == "42"
+        random_seed_val = utils.set_seed(None)
+        assert(int(random_seed_val) < 10001)

python/gubbins/treebuilders.py

@@ -20,7 +20,6 @@
 import sys
 import os
 import subprocess
-from random import randint
 
 from Bio import SeqIO
 
@@ -129,7 +128,7 @@ def bootstrapping_command(self, alignment_filename: str, input_tree: str, basena
 class FastTree:
     """Class for operations with the FastTree executable"""
 
-    def __init__(self, threads: int, bootstrap = 0, model='GTRCAT', verbose=False, additional_args = None):
+    def __init__(self, threads: int, bootstrap = 0, model='GTRCAT', seed = None, verbose=False, additional_args = None):
         """Initialises the object"""
         self.verbose = verbose
         self.threads = threads
@@ -139,6 +138,7 @@ def __init__(self, threads: int, bootstrap = 0, model='GTRCAT', verbose=False, a
         self.alignment_suffix = ".snp_sites.aln"
         self.bootstrap = bootstrap
         self.additional_args = additional_args
+        self.seed = utils.set_seed(seed)
 
         # Identify executable
         self.potential_executables = ["FastTree", "fasttree"]
@@ -164,6 +164,7 @@ def __init__(self, threads: int, bootstrap = 0, model='GTRCAT', verbose=False, a
             command.extend(["-gtr"])
         else:
             command.extend([self.model])
+        command.extend(["-seed",self.seed])
         # Additional arguments
         if self.additional_args is not None:
             command.extend([self.additional_args])
@@ -256,7 +257,7 @@ def get_bootstrapped_trees_file(self, tmp: str, basename: str) -> str:
 class IQTree:
     """Class for operations with the IQTree executable"""
 
-    def __init__(self, threads: 1, model: str, bootstrap = 0, internal_node_prefix="", verbose=False, use_best=False, additional_args = None):
+    def __init__(self, threads: 1, model: str, bootstrap = 0, seed = None, internal_node_prefix="", verbose=False, use_best=False, additional_args = None):
         """Initialises the object"""
         self.verbose = verbose
         self.threads = threads
@@ -271,6 +272,7 @@ def __init__(self, threads: 1, model: str, bootstrap = 0, internal_node_prefix="
         self.internal_node_prefix = internal_node_prefix
         self.bootstrap = bootstrap
         self.use_best = use_best
+        self.seed = utils.set_seed(seed)
         self.additional_args = additional_args
 
         # Construct base command
@@ -303,6 +305,7 @@ def __init__(self, threads: 1, model: str, bootstrap = 0, internal_node_prefix="
             command.extend(["-m","GTR+G4"])
         else:
             command.extend(["-m",self.model])
+        command.extend(["-seed",self.seed])
         # Additional arguments
         if self.additional_args is not None:
             command.extend([self.additional_args])
@@ -432,7 +435,7 @@ def run_model_comparison(self, alignment_filename: str, basename: str) -> str:
 class RAxML:
     """Class for operations with the RAxML executable"""
 
-    def __init__(self, threads: 1, model='GTRCAT', bootstrap = 0, internal_node_prefix="", verbose=False, additional_args = None):
+    def __init__(self, threads: 1, model='GTRCAT', bootstrap = 0, seed = None, internal_node_prefix="", verbose=False, additional_args = None):
         """Initialises the object"""
         self.verbose = verbose
         self.threads = threads
@@ -446,6 +449,7 @@ def __init__(self, threads: 1, model='GTRCAT', bootstrap = 0, internal_node_pref
         self.alignment_suffix = ".phylip"
         self.internal_node_prefix = internal_node_prefix
         self.bootstrap = bootstrap
+        self.seed = utils.set_seed(seed)
         self.additional_args = additional_args
 
         self.single_threaded_executables = ['raxmlHPC-AVX2', 'raxmlHPC-AVX', 'raxmlHPC-SSE3', 'raxmlHPC']
@@ -465,9 +469,6 @@ def __init__(self, threads: 1, model='GTRCAT', bootstrap = 0, internal_node_pref
         if self.threads > 1:
             command.extend(["-T", str(self.threads)])
 
-        # Set a seed
-        command.extend(["-p",str(randint(0, 10000))])
-
         # Add flags
         command.extend(["-safe"])
         if self.model == 'JC':
@@ -482,6 +483,7 @@ def __init__(self, threads: 1, model='GTRCAT', bootstrap = 0, internal_node_pref
             command.extend(["-m","GTRGAMMA"])
         else:
             command.extend(["-m", self.model])
+        command.extend(["-p",self.seed])
         # Additional arguments
         if self.additional_args is not None:
             command.extend([self.additional_args])
@@ -579,9 +581,7 @@ def bootstrapping_command(self, alignment_filename: str, input_tree: str, basena
         command = self.base_command.copy()
         command.extend(["-s", alignment_filename, "-n", basename + ".bootstrapped_trees"])
         command.extend(["-w",tmp])
-        p_seed = str(randint(0, 10000))
-        command.extend(["-p",p_seed])
-        command.extend(["-x",p_seed])
+        command.extend(["-x",self.seed])
         command.extend(["-#",str(self.bootstrap)])
         # Output
         if not self.verbose:
@@ -592,8 +592,6 @@ def bootstrapping_command(self, alignment_filename: str, input_tree: str, basena
     def sh_test(self, alignment_filename: str, input_tree: str, basename: str, tmp: str) -> str:
         """Runs a single branch support test"""
         command = self.base_command.copy()
-        p_seed = str(randint(0, 10000))
-        command.extend(["-p",p_seed])
         command.extend(["-f", "J"])
         command.extend(["-s", alignment_filename, "-n", input_tree + ".sh_support"])
         command.extend(["-t", input_tree])
@@ -610,7 +608,7 @@ def get_bootstrapped_trees_file(self, tmp: str, basename: str) -> str:
 class RAxMLNG:
     """Class for operations with the RAxML executable"""
 
-    def __init__(self, threads: 1, model: str, bootstrap = 0, internal_node_prefix = "", verbose = False, additional_args = None):
+    def __init__(self, threads: 1, model: str, bootstrap = 0, seed = None, internal_node_prefix = "", verbose = False, additional_args = None):
         """Initialises the object"""
         self.verbose = verbose
         self.threads = threads
@@ -624,6 +622,7 @@ def __init__(self, threads: 1, model: str, bootstrap = 0, internal_node_prefix =
         self.alignment_suffix = ".phylip"
         self.internal_node_prefix = internal_node_prefix
         self.bootstrap = bootstrap
+        self.seed = utils.set_seed(seed)
         self.additional_args = additional_args
 
         self.single_threaded_executables = ['raxml-ng']
@@ -655,6 +654,7 @@ def __init__(self, threads: 1, model: str, bootstrap = 0, internal_node_prefix =
             command.extend(["GTR+G"])
         else:
             command.extend([self.model])
+        command.extend(["--seed",self.seed])
         # Additional arguments
         if self.additional_args is not None:
             command.extend([self.additional_args])

python/gubbins/utils.py

@@ -23,6 +23,7 @@
 import re
 import numpy as np
 import collections
+from random import randint
 try:
     from multiprocessing.managers import SharedMemoryManager
     NumpyShared = collections.namedtuple('NumpyShared', ('name', 'shape', 'dtype'))
@@ -195,3 +196,11 @@ def extend_args(var,add):
         var.extend([add])
         var = " ".join(var)
     return var
+
+def set_seed(seed):
+    """Set seed when specified"""
+    if seed is None:
+        seed = str(randint(0, 10000))
+    else:
+        seed = str(seed)
+    return seed