Implementing SNP matrix and Gubbins/Iqtree/Raxml changes

config_ali

@@ -118,6 +118,11 @@ openmpi_bin: /openmpi/bin/
 base_cmd: raxmlHPC-HYBRID-SSE3
 parameters: -f a -x 12345 -p 12345 -N autoMRE -m GTRCAT -T 20
 
+[iqtree]
+iqtree_bin: /nfs/esnitkin/bin_group/anaconda3/bin/
+base_cmd: iqtree
+parameters: -nt AUTO -bb 1000 -m GTR+G+ASC
+
 [gubbins]
 # Change this path to wherever gubbins is located/installed. Right now, installed using conda from anaconda3 package installed in bin_group.
 gubbins_bin: /nfs/esnitkin/bin_group/anaconda3/bin/
@@ -964,3 +969,8 @@ Ref_Name: MRSA_Newman.fasta
 # path to the reference genome fasta file.
 Ref_Path: /nfs/esnitkin/bin_group/variant_calling_bin/reference/MRSA_Newman/
 
+[AB030]
+# Name of reference genome fasta file.
+Ref_Name: AB030.fasta
+# path to the reference genome fasta file.
+Ref_Path: /nfs/esnitkin/bin_group/variant_calling_bin/reference/AB030/

modules/gubbins.py

@@ -12,14 +12,15 @@ def gubbins(gubbins_dir, input_fasta, jobrun, logger, Config):
                  logger,
                  'info')
 
-
+    call("module load bioperl python-anaconda2/201607 biopython dendropy reportlab fasttree RAxML fastml/gub gubbins",logger)
+    #os.system("module load bioperl python-anaconda2/201607 biopython dendropy reportlab fasttree RAxML fastml/gub gubbins")
 
     #gubbins_cmd = "%s/%s --prefix %s/%s %s" % (
     # ConfigSectionMap("gubbins", Config)['gubbins_bin'], ConfigSectionMap("gubbins", Config)['base_cmd'], gubbins_dir,
     # (os.path.basename(input_fasta)).replace('.fa', ''), input_fasta)
 
 
-    gubbins_cmd = "%s --prefix %s/%s %s" % (
+    gubbins_cmd = "%s --threads 6 --prefix %s/%s %s" % (
     ConfigSectionMap("gubbins", Config)['base_cmd'], gubbins_dir,
     (os.path.basename(input_fasta)).replace('.fa', ''), input_fasta)
     keep_logging('\nRunning Gubbins on: %s' % input_fasta, '\nRunning Gubbins: %s\n' % input_fasta,
@@ -36,9 +37,9 @@ def gubbins(gubbins_dir, input_fasta, jobrun, logger, Config):
     elif jobrun == "parallel-cluster":
         job_file_name = "%s/gubbins_%s.pbs" % (gubbins_dir, os.path.basename(input_fasta))
         job_name = os.path.basename(job_file_name)
-        job_print_string = "#PBS -N %s\n#PBS -M %s\n#PBS -m %s\n#PBS -V\n#PBS -l nodes=1:ppn=4,mem=47000mb,walltime=76:00:00\n#PBS -q %s\n#PBS -A %s\n#PBS -l qos=flux\ncd %s\n%s" % (job_name, ConfigSectionMap("scheduler", Config)['email'], ConfigSectionMap("scheduler", Config)['notification'], ConfigSectionMap("scheduler", Config)['queue'], ConfigSectionMap("scheduler", Config)['flux_account'], gubbins_dir, gubbins_cmd)
+        job_print_string = "#PBS -N %s\n#PBS -M %s\n#PBS -m %s\n#PBS -V\n#PBS -l nodes=1:ppn=4,mem=47000mb,walltime=250:00:00\n#PBS -q %s\n#PBS -A %s\n#PBS -l qos=flux\ncd %s\n%s" % (job_name, ConfigSectionMap("scheduler", Config)['email'], ConfigSectionMap("scheduler", Config)['notification'], ConfigSectionMap("scheduler", Config)['queue'], ConfigSectionMap("scheduler", Config)['flux_account'], gubbins_dir, gubbins_cmd)
         f1=open(job_file_name, 'w+')
         f1.write(job_print_string)
         f1.close()
         #os.system("qsub %s" % job_file_name)
-        call("qsub %s" % job_file_name, logger)
+        call("qsub %s" % job_file_name, logger)

modules/iqtree.py

@@ -0,0 +1,35 @@
+__author__ = 'alipirani'
+import os
+from config_settings import ConfigSectionMap
+from logging_subprocess import *
+from log_modules import *
+
+
+def iqtree(tree_dir, input_fasta, jobrun, logger, Config):
+    keep_logging('Running iqtree on input: %s' % input_fasta, 'Running iqtree on input: %s' % input_fasta, logger, 'info')
+    iqtree_cmd = "%s/%s -s %s %s -pre %s" % (
+        ConfigSectionMap("iqtree", Config)['iqtree_bin'],
+        ConfigSectionMap("iqtree", Config)['base_cmd'],
+        input_fasta, ConfigSectionMap("iqtree", Config)['parameters'],
+        (os.path.basename(input_fasta)).replace('.fa', ''))
+
+    keep_logging('%s' % iqtree_cmd, '%s' % iqtree_cmd, logger, 'info')
+    if jobrun == "parallel-local" or jobrun == "local":
+        call("cd %s" % tree_dir, logger)
+        call(iqtree_cmd, logger)
+    elif jobrun == "cluster":
+        call("cd %s" % tree_dir, logger)
+        call(iqtree_cmd, logger)
+    elif jobrun == "parallel-cluster":
+        job_file_name = "%s/iqtree_%s.pbs" % (tree_dir, os.path.basename(input_fasta))
+        job_name = os.path.basename(job_file_name)
+        job_print_string = "#PBS -N %s\n#PBS -M %s\n#PBS -m %s\n#PBS -V\n#PBS -l nodes=1:ppn=12,mem=47000mb,walltime=250:00:00\n#PBS -q %s\n#PBS -A %s\n#PBS -l qos=flux\ncd %s\n%s" % (
+        job_name, ConfigSectionMap("scheduler", Config)['email'], ConfigSectionMap("scheduler", Config)['notification'],
+        ConfigSectionMap("scheduler", Config)['queue'], ConfigSectionMap("scheduler", Config)['flux_account'], tree_dir,
+        iqtree_cmd)
+        f1 = open(job_file_name, 'w+')
+        f1.write(job_print_string)
+        f1.close()
+        # os.system("qsub %s" % job_file_name)
+        keep_logging('qsub %s' % job_file_name, 'qsub %s' % job_file_name, logger, 'info')
+	#call("qsub %s" % job_file_name, logger)

modules/raxml.py

@@ -24,7 +24,7 @@ def raxml(tree_dir, input_fasta, jobrun, logger, Config):
     elif jobrun == "parallel-cluster":
         job_file_name = "%s/raxml_%s.pbs" % (tree_dir, os.path.basename(input_fasta))
         job_name = os.path.basename(job_file_name)
-        job_print_string = "#PBS -N %s\n#PBS -M %s\n#PBS -m %s\n#PBS -V\n#PBS -l nodes=1:ppn=4,mem=47000mb,walltime=76:00:00\n#PBS -q %s\n#PBS -A %s\n#PBS -l qos=flux\ncd %s\n%s" % (job_name, ConfigSectionMap("scheduler", Config)['email'], ConfigSectionMap("scheduler", Config)['notification'], ConfigSectionMap("scheduler", Config)['queue'], ConfigSectionMap("scheduler", Config)['flux_account'], tree_dir, raxml_cmd)
+        job_print_string = "#PBS -N %s\n#PBS -M %s\n#PBS -m %s\n#PBS -V\n#PBS -l nodes=1:ppn=12,mem=47000mb,walltime=250:00:00\n#PBS -q %s\n#PBS -A %s\n#PBS -l qos=flux\ncd %s\n%s" % (job_name, ConfigSectionMap("scheduler", Config)['email'], ConfigSectionMap("scheduler", Config)['notification'], ConfigSectionMap("scheduler", Config)['queue'], ConfigSectionMap("scheduler", Config)['flux_account'], tree_dir, raxml_cmd)
         f1=open(job_file_name, 'w+')
         f1.write(job_print_string)
         f1.close()