diff --git a/example.config b/example.config index d222ce5..2a3f692 100644 --- a/example.config +++ b/example.config @@ -30,3 +30,5 @@ read_length = 150 #maximum value with example error profile is 150, use a defaul fragment_size = 380 stdev_frag_size = 120 +random_seed = 12345 + diff --git a/treetoreads.py b/treetoreads.py index a047f2a..18365b7 100755 --- a/treetoreads.py +++ b/treetoreads.py @@ -5,6 +5,7 @@ import os import sys import argparse +import subprocess from subprocess import call import dendropy @@ -25,9 +26,6 @@ class TreeToReads(object): _vargen = 0 def __init__(self, configfi, run=1, main=None): """initialized object, most attributes generated through self._check_args using config file.""" - self.seed = random.randint(0, sys.maxsize) - sys.stdout.write("Random seed is {}\n".format(self.seed)) - random.seed(self.seed) self.configfi = configfi self.run = run self.main = main @@ -42,6 +40,13 @@ def __init__(self, configfi, run=1, main=None): self.prefix = self.config['prefix'] else: self.prefix = '' + if 'random_seed' in self.config: + self.seed = self.config['random_seed'] + random.seed(self.seed) + else: + self.seed = random.randint(0, sys.maxsize) + sys.stdout.write("Random seed is {}\n".format(self.seed)) + random.seed(self.seed) if self.run: if (self.no_art == 0) and (self.config.get('coverage')): self.run_art() @@ -111,7 +116,8 @@ def read_args(self): 'error_model1', 'error_model2', 'indel_model', - 'indel_rate'] + 'indel_rate', + 'random_seed'] for lin in config: lii = lin.split('=') self.config[lii[0].strip()] = lii[-1].split('#')[0].strip() @@ -677,6 +683,7 @@ def mut_genomes_indels(self):#TODO does not account for SNPs in indsertions def run_art(self, coverage=None): """Runs ART to simulate reads from the simulated genomes""" + if not self._genmut: if self.get_arg("indel_model"): if call(['which', 'indelible'], stdout=open('/dev/null', 'w')) == 1: @@ -687,6 +694,9 @@ def run_art(self, coverage=None): self.mut_genomes_indels() else: self.mut_genomes_no_indels() + if call(['which', 'art_illumina'], stdout=open('/dev/null', 'w')) == 1: + sys.stderr.write('''ERROR: Art not found. install and add to path. Exiting\n''') + sys.exit() if coverage is None: coverarg = self.get_arg('cov') else: @@ -740,6 +750,7 @@ def run_art(self, coverage=None): '-2', self.get_arg('error_model2'), '-na', #Don't output alignment file '-p', #for paired end reads + '-rs', str(self.seed), '-i', '{}/fasta_files/{}{}.fasta'.format(self.outd, self.prefix, seq), '-l', '{}'.format(read_length), '-f', str(cov[seq]), @@ -756,6 +767,7 @@ def run_art(self, coverage=None): '-na', #Don't output alignment file '-i', '{}/fasta_files/{}{}.fasta'.format(self.outd, self.prefix, seq), '-l', '{}'.format(read_length), + '-rs', str(self.seed), '-f', str(cov[seq]), '-m', '{}'.format(fragment_size), '-s', '{}'.format(stdev_frag_size), @@ -764,7 +776,8 @@ def run_art(self, coverage=None): seq, self.prefix, seq)] - call(artparam, stdout=open('{}/art_log'.format(self.outd), 'w'), stderr=open('{}/art_log'.format(self.outd), 'a')) + #sys.stdout.write("Running: "+" ".join(artparam)) + process = subprocess.check_call(artparam, stdout=open('{}/art_log'.format(self.outd), 'w'), stderr=open('{}/art_log'.format(self.outd), 'a')) # print("called {}".format(" ".join(artparam))) assert os.path.exists('{}/fastq/{}{}/{}{}_1.fq'.format(self.outd, self.prefix, seq, self.prefix, seq)) gzippar = ['gzip',