phylociraptor

#!/usr/bin/env python3

import sys, os, io
sys.path.insert(0, os.getcwd()+"/bin")
from phylociraptor_modules.usageinfo import *
from phylociraptor_modules.outfiles import *
import time
import subprocess
import argparse
import glob

if sys.version_info[0] < 3:
	raise Exception("Must be using Python 3")
	exit(1)

njobs = "10000"
latency_wait = "10" #in seconds
singularity_bindpoints = "-B $(pwd)/.usr_tmp/:/usertmp"
debug = False #turn debugging mode on (True) and off (False)
cluster_config_defaults= {"slurm": "data/cluster-config-SLURM.yaml.template", "sge":"data/cluster-config-SGE.yaml.template", "torque":"data/cluster-config-TORQUE.yaml.template"}


def now():
	return time.strftime("%Y-%m-%d %H:%M") + " -"

def progressbar(it, progress, prefix="", size=60, file=sys.stdout):
	count = len(it)
	def show(j):
		x = int(size*j/count)
		return "%s[%s%s] %i/%i\r" % (prefix, "#"*x, "."*(size-x), j, count)
	return show(progress)

def help_message(mes):
	return mes

def determine_submission_mode(flag):
	cmd = []
	if "serial" in flag:
		return ["--cores" + flag.replace("serial","")]
	elif "sge" in flag: # for SGE the dependencies need to be under quotes, because there are () characters in the dependency string.
		return ["--cluster", "bin/immediate-submit/immediate_submit.py '{dependencies}' %s" % flag, "--immediate-submit", "--jobs", njobs, "--notemp"]
	else:
		return ["--cluster", 'bin/immediate-submit/immediate_submit.py {dependencies} %s' % flag, "--immediate-submit", "--jobs", njobs, "--notemp"]

def get_flags(flags):
	mapdict ={
	#"t": '--cluster "bin/immediate_submit.py {dependencies} ', "cluster": '--cluster "bin/immediate_submit.py {dependencies} ',
	"c": "--cluster-config", "cluster_config": "--cluster-config",
	"FORCE": "-F",
	"force": "-f",
	"dry": "-n"
	}
	cmd = []
	if debug:
		print(now(), "DEBUG: flags: ", flags)
	if flags["cluster_config"] == None and not "serial" in flags["cluster"]: #in case no cluster config file was specified, get default value
		ccf = cluster_config_defaults[flags["cluster"]]
		flags["cluster_config"] = ccf
		print(now(), "INFO: No cluster config file specified. Will try to use default file: %s" % ccf)	
	else:
		ccf = flags["cluster_config"]
		if ccf != None:
			if not os.path.isfile(ccf):
				print(now(), "ERROR: Specified cluster config file:", ccf, "not found.")
				sys.exit(1)
	for flag in flags.keys():
		if flag in mapdict.keys() and flags[flag] != None:
			if flag == "t" or flag == "cluster": #handle cluster specification
				arg = mapdict[flag]
				arg = arg + " "+flags[flag]+'"'
				cmd.append(arg) 
			if flag == "c" or flag == "cluster_config": #handle cluster config file
				#print("here")
				#arg = mapdict[flag]
				#arg = arg + " " + flags[flag]
				cmd.append(mapdict[flag])
				cmd.append(ccf)
			else:
				if flags[flag]:
					cmd.append(mapdict[flag])
	
	return cmd

def check_required_files(runmode):
	for f in outfile_dict[runmode]:
		if not os.path.isfile(f):
			if debug:
				print(now(),"DEBUG, check_required_files: File not found:", f)
			return f
	return

def check_directories(mode):
	for directory in outdir_dict[mode]:
		if not os.path.isdir(directory):
			if debug:
				print(now(), "DEBUG, check_directories: Directory not found:", directory)
			return directory
	return

def check_hasrun(mode):
	f = checkpoint_file_dict[mode]
	if not os.path.isfile(f):
		if debug:
			print(now(), "DEBUG, check_hasrun: Not found:", f)
		return f
	else:
		return

def check_config_file(f):
	if os.path.isfile(f):
		return ["--configfile", f]
	else:
		print(now(), "ERROR: Config file specified with --config-file not found:", f)
		sys.exit(1)

def check_for_errors(result):
	if result.startswith("WorkflowError"):
		return now()+" ERROR: phylociraptor encountered an error. You may try to run with --verbose to diagnose what has gone wrong.\n"
	if result.startswith("The singularity command"):
		return result
	if result.startswith("Error") or result.startswith("error"):
		return now()+" ERROR: There was an error. Maybe run with --verbose to diagnose. The error occurred here: %s " % result	
	if result.startswith("Directory cannot be locked."):
		return now()+" ERROR: "+result	
	if result.startswith("IncompleteFilesException"):
		return now()+" ERROR: There seems to be a problem with incomplete output files from a previous run.\nYou can run --verbose to see which files are incomplete and how to resolve the problem.\n"
	if result.startswith("KeyError"):
		return now()+" ERROR: "+result
	return ""


def check_required_software():	
	try:
		singularity_version = str(subprocess.check_output(['singularity', 'version']).decode('ascii').strip())
		print(now(), "Singularity version:", singularity_version)
	except FileNotFoundError:
		print(now(), "ERROR: Singularity command not found.")
		sys.exit(1)
	try:
		snakemake_version = str(subprocess.check_output(['snakemake', '-v']).decode('ascii').strip())
		print(now(), "Snakemake version:", snakemake_version)
	except FileNotFoundError:
		print(now(),"ERROR: Snakemake command not found.")
		sys.exit(1)
	
	
def execute_command(cmd, verbose):
	# this should also correctly parse and display:
	# IncompleteFilesException
	# when singularity command is not available
	popen = subprocess.Popen(cmd, stdout=subprocess.PIPE, stderr=subprocess.STDOUT)
	print(now(), "INFO: If you execute this command for the first time, singularity containers need to be downloaded which can take several minutes.")
	line = ""
	jobcounts = False
	njobs = 0
	char = ""
	nline = 0
	curr_task = 0
	if "-n" in cmd:
		print(now(), "INFO: --dry specified, will only perform a dry run of the analysis")
		for line in io.TextIOWrapper(popen.stdout, encoding="utf-8"):
			nline += 1
			if verbose:
				yield line
			else:
				result = line
				if check_for_errors(result):
					yield check_for_errors(result)
				if result.startswith("Nothing"):
					yield now() + " There is nothing to do. If you want to force a rerun of this step try: -f or --force\n"
				if result.startswith("Job counts"):
					jobcounts = True
				elif jobcounts and result != "" and nline <= 30: #only keep jobcount info from the beginning of the output
					if len(result.split("\t")) == 2: # the line with the total number of jobs has two elements
						yield now()+" Total number of tasks to run/submit: %s\n" % result.split("\t")[1].strip() 
				elif jobcounts and result == "":
					jobcounts = False
				line=""
	elif "--cluster" in cmd:
		for line in io.TextIOWrapper(popen.stdout, encoding="utf-8"):
			nline += 1
			if verbose:
				yield line
			else:
				result = line
				if check_for_errors(result):
					yield check_for_errors(result)
				if result.startswith("Singularity"):
					yield result
				if result.startswith("Job counts"):
					jobcounts = True	
				elif jobcounts and result != "" and nline <= 30: #only keep jobcount info from the beginning of the output
					if len(result.split("\t")) == 2: # the line with the total number of jobs has two elements
						yield now()+" Total number of jobs to submit: %s\n" % result.split("\t")[1].strip()
						njobs = int(result.split("\t")[1])
				elif jobcounts and result == "":
					jobcounts = False
				if line.startswith("rule"):
					curr_task += 1
					yield progressbar(range(njobs),curr_task, "Submitting: ", 100)
				line=""
	else:
		for line in io.TextIOWrapper(popen.stdout, encoding="utf-8"):
			nline += 1
			if verbose:
				yield line
			else:
				result = line
				if check_for_errors(result):
					yield check_for_errors(result)
				if result.startswith("Nothing"):
					yield now()+" There is nothing to do. If you want to force a rerun of this step try: -f or --force\n"
				if result.startswith("Job counts"):
					jobcounts = True	
				elif jobcounts and result != "" and nline <= 30: #only keep jobcount info from the beginning of the output
					if len(result.split("\t")) == 2: # the line with the total number of jobs has two elements
						yield now()+" Total number of tasks to run: %s\n" % result.split("\t")[1].strip()
						njobs = int(result.split("\t")[1])
				elif jobcounts and result == "":
					jobcounts = False
				if line.startswith("rule"):
					curr_task += 1
					yield progressbar(range(njobs),curr_task, "Runnning: ", 100)
				line=""

def get_additional_snakemake_flags(flags, rerun):
	if flags:
		flags= flags.strip() #need to remove trailing charcters such as spaces first otherwise the list will be messed up
	if rerun: # add --rerun-incomplete in case it is set
		if flags: #add depending on if flags already contains values or not
			flags += " --rerun-incomplete"
		else:
			flags += "--rerun-incomplete"
	if flags:
		print(now(), "INFO: Additional flags will be passed on to snakemake: ", flags)
		return flags.split(" ")
	else:
		return []

def get_additional_singularity_flags(flags):
	if flags:
		print(now(), "INFO: Additional flags will be passed on to singularity: ", flags)
		return ["--singularity-args"]+[singularity_bindpoints +" " + flags]
	else:
		return ["--singularity-args"]+ [singularity_bindpoints]
	
pars = argparse.ArgumentParser(usage=help_message(default_help))
pars.add_argument('--debug', action='store_true', dest="debug", required=False)
pars.add_argument('-v', '--version', action='store_true', dest='version', required=False)
pars.add_argument('command', action='store', nargs="?")
pars.add_argument('arguments', action='store', nargs=argparse.REMAINDER)

args = pars.parse_args()

def get_commit():
		return str(subprocess.check_output(['git', 'rev-parse', '--short', 'HEAD']).decode('ascii').strip())
if args.version == True:
	commit = get_commit()
	print("Version:", version)
	print("Git commit:", commit)
	sys.exit(0)

if args.debug == True:
	print(now(), "DEBUG: Addditional debugging output enabled.")
	debug=True

if not args.command:
	print(default_help)
	sys.exit(0)

# create own argparse class to include some standard options:
class PhyloParser(argparse.ArgumentParser):
	def __init__(self, **kwargs):
		super().__init__(**kwargs)
		self.add_argument("-f", "--force", action="store_true" )
		self.add_argument("-F", "--FORCE", action="store_true" )
		self.add_argument("-t", "--cluster",  action="store", default="serial")
		self.add_argument("-c", "--cluster-config", action="store")
		self.add_argument("--dry", action="store_true")
		self.add_argument("-h", "--help", action="store_true")
		self.add_argument("--verbose", action="store_true", default=False)	
		self.add_argument("--singularity", action="store",dest="si_args", default="")
		self.add_argument("--snakemake", action="store",dest="sm_args", default="")
		self.add_argument("--rerun-incomplete", action="store_true", dest="rerun", default=False)

class UtilParser(argparse.ArgumentParser):
	def __init__(self, **kwargs):
		super().__init__(**kwargs)
		self.add_argument("-h", "--help", action="store_true")
		self.add_argument("--verbose", action="store_true", default=False)	
	

if args.command == "setup":
	print(now(), "Welcome to phylociraptor setup v%s" % version)
	setup_parser = PhyloParser(usage=help_message(setup_help), add_help=False)
	setup_parser.add_argument("--config-file", action="store", dest="config_file", default="data/config.yaml")
	setup_parser.add_argument("--samples-csv", action="store", dest="samples_file")
	setup_parser.add_argument("--busco-set", action="store", dest="busco_set")
	setup_parser.add_argument("--add_genomes", action="store_true",default=False, dest="add_genomes")

	setup_args = setup_parser.parse_args(args.arguments)
	if setup_args.help or len(sys.argv) <= 2:
		print(help_message(setup_help))
		sys.exit(0)
	
	check_required_software()
	
	if (setup_args.add_genomes):
		print(now(),"Will only add genomes")
		cmd = ["snakemake","-p", "-s" , "rules/setup.smk", "--use-singularity", "-r", "add_genomes", "--latency-wait", latency_wait]
	else:
		cmd = ["snakemake","-p", "-s" , "rules/setup.smk", "--use-singularity", "-r", "setup", "--latency-wait", latency_wait]
	cmd += get_flags(vars(setup_args))
	cmd += determine_submission_mode(setup_args.cluster)

	# check for overrides of config file parameters	
	if setup_args.config_file:
		cmd += check_config_file(setup_args.config_file)
	if setup_args.busco_set or setup_args.samples_file:
		cmd += ["--config"]
		if not os.path.isfile(setup_args.samples_file):
			print(now(), "ERROR: File specified with --samples-csv not found:", setup_args.samples_file)
			sys.exit(0)
		cmd += ["busco={set: %s}" % setup_args.busco_set]
		cmd += ["species='%s'" % setup_args.samples_file]
	
	cmd += get_additional_snakemake_flags(setup_args.sm_args, setup_args.rerun)
	cmd += get_additional_singularity_flags(setup_args.si_args)

	for line in execute_command(cmd, setup_args.verbose):
		print(line, end="\r")
	if debug:
		print(now(),"DEBUG:", cmd)

elif args.command=="orthology":
	print(now(), "Welcome to phylociraptor orthology v%s" % version)
	orthology_parser = PhyloParser(usage=help_message(orthology_help),add_help=False)	
	orthology_parser.add_argument("--busco_threads", action="store")
	orthology_parser.add_argument("--config-file", action="store", dest="config_file", default="data/config.yaml")
	orthology_parser.add_argument("--augustus_species", action="store")
	orthology_parser.add_argument("--additional_params", action="store")
	orthology_args = orthology_parser.parse_args(args.arguments)
	if orthology_args.help or len(sys.argv) <= 2: # help is specified
		print(help_message(orthology_help))
		sys.exit(0)
	check_required_software()
	if check_required_files(args.command):
		print(now(), "ERROR: Some files are missing preventing this part to run.\nDid your run phylociraptor setup already? Missing file:", check_required_files(args.command))
		sys.exit(0)
	cmd = ["snakemake","-p", "-s", "rules/orthology.smk", "--use-singularity", "-r", "orthology", "--latency-wait", latency_wait]
	cmd += determine_submission_mode(orthology_args.cluster)
	cmd += get_flags(vars(orthology_args))

	# check for overrides of config file parameters	
	if orthology_args.config_file:
		cmd += check_config_file(orthology_args.config_file)
	if orthology_args.busco_threads or orthology_args.augustus_species or orthology_args.additional_params:
		cmd += ["--config"]
		arg = "busco={"
		if orthology_args.busco_threads: #evaluation of this string by snakemake is lazy, it does not seem to matter if there are extra , characters...
			arg += "threads: %s, " % orthology_args.busco_threads
		if orthology_args.augustus_species:
			arg += "augustus_species: %s, " % orthology_args.augustus_species
		if orthology_args.additional_params:
			arg += "additional_parameters: %s " % orthology_args.additional_params
		arg+="}"
		cmd += [arg]
	
	cmd += get_additional_snakemake_flags(orthology_args.sm_args, orthology_args.rerun)
	cmd += get_additional_singularity_flags(orthology_args.si_args)

	
	for line in execute_command(cmd, orthology_args.verbose):
                print(line, end="\r")
	if debug:
		print(now(),"DEBUG:", " ".join(cmd))
elif args.command=="filter-orthology":
	print(now(), "Welcome to phylociraptor filter-orthology v%s" % version)
	forthology_parser = PhyloParser(add_help=False)
	forthology_parser.add_argument("--config-file", action="store", dest="config_file", default="data/config.yaml")
	forthology_parser.add_argument("--dupseq", action="store")
	forthology_parser.add_argument("--cutoff", action="store")
	forthology_parser.add_argument("--minsp", action="store")
	forthology_parser.add_argument("--seqtype", action="store")
	
	forthology_args = forthology_parser.parse_args(args.arguments)	
	if forthology_args.help or len(sys.argv) <= 2: # help is specified
		print(help_message(forthology_help))
		sys.exit(0)
	check_required_software()
	if check_required_files(args.command):
		print(now(), "ERROR: Some files are missing preventing this part to run.\nDid your run phylociraptor setup and orthology already? Missing file:", check_required_files(args.command))
		sys.exit(0)
	cmd = ["snakemake", "-p", "-s", "rules/filter-orthology.smk", "--use-singularity", "-r", "filter_orthology", "--latency-wait", latency_wait]
	cmd += determine_submission_mode(forthology_args.cluster)
	cmd += get_flags(vars(forthology_args))
	
	# check for overrides of config file parameters	
	if forthology_args.config_file:
		cmd += check_config_file(forthology_args.config_file)
	if forthology_args.dupseq or forthology_args.cutoff or forthology_args.minsp or forthology_args.seqtype:
		cmd += ["--config"]
		arg = "filtering={"
		if forthology_args.dupseq:
			if forthology_args.dupseq in ["persample", "perfile"]:
				arg += "dupseq: %s, "% forthology_args.dupseq
			else:
				print(now(), "ERROR: Wrong parameter specified for --dupseq: ", forthology_args.dupseq)
				sys.exit(0)
		if forthology_args.cutoff:
			if float(forthology_args.cutoff) >=0 and float(forthology_args.cutoff) <=1:
				arg += "cutoff: %s, " % forthology_args.cutoff
			else:
				print(now(), "ERROR: Wrong parameter specified for --cutoff: ", forthology_args.dubseq)	
				sys.exit(0)
		if forthology_args.minsp:
			arg += "minsp: %s, " %forthology_args.minsp
		if forthology_args.seqtype:
			if forthology_args.seqtype in ["aa", "nu"]:
				arg += "seq_type: %s " % forthology_args.seqtype
			else:
				print(now(), "ERROR: Wrong parameter specified for --seqtype: ", forthology_args.seqtype)
				sys.exit(0)
		arg += "}"
		cmd += [arg]
	
	cmd += get_additional_snakemake_flags(forthology_args.sm_args, forthology_args.rerun)
	cmd += get_additional_singularity_flags(forthology_args.si_args)
	
	for line in execute_command(cmd, forthology_args.verbose):
                print(line, end="\r")
	if debug:
		print(now(),"DEBUG:", cmd)
elif args.command=="align":
	print(now(), "Welcome to phylociraptor align v%s" % version)
	align_parser = PhyloParser(usage=help_message(align_help), add_help=False)
	align_parser.add_argument("--config-file", action="store", dest="config_file", default="data/config.yaml")
	align_parser.add_argument("--method", action="store")
	align_parser.add_argument("--threads", action="store")
	align_parser.add_argument("--parameters", action="store")
	align_args = align_parser.parse_args(args.arguments)	
	if align_args.help or len(sys.argv) <= 2: # help is specified
		print(help_message(align_help))
		sys.exit(0)
	check_required_software()
	if check_required_files(args.command):
		print(now(), "ERROR: Some files are missing preventing this part to run.\nDid your run the required steps before? Missing file:", check_required_files(args.command))
		sys.exit(0)

	cmd = ["snakemake", "-p", "-s", "rules/align.smk", "--use-singularity", "-r", "align", "--latency-wait", latency_wait]
	cmd += determine_submission_mode(align_args.cluster)
	cmd += get_flags(vars(align_args))

	# check for overrides of config file parameters	
	if align_args.config_file:
		cmd += check_config_file(align_args.config_file)
	if align_args.method or align_args.threads or align_args.parameters:
		cmd += ["--config"]
		arg = "alignment={"
		if align_args.method:
			if align_args.method in ["mafft"]:
				arg += "method: %s, " % align_args.method
			else:
				print(now(), "ERROR: Wrong parameter specified with --method:", align_args.method)
				sys.exit(0)	
		if align_args.threads:
			arg += "threads: %s, " % align_args.threads
		if align_args.parameters:
			arg += "parameters: '%s'" % align_args.parameters
		arg += " }"
		cmd += [arg]	
	
	cmd += get_additional_snakemake_flags(align_args.sm_args, align_args.rerun)
	cmd += get_additional_singularity_flags(align_args.si_args)
	
	for line in execute_command(cmd, align_args.verbose):
                print(line, end="\r")
	if debug:
		print(now(),"DEBUG:", " ".join(cmd))
elif args.command=="filter-align":
	print(now(), "Welcome to phylociraptor filter-align v%s" % version )
	falign_parser = PhyloParser(usage=help_message(falign_help), add_help=False)
	falign_parser.add_argument("--config-file", action="store", dest="config_file", default="data/config.yaml")
	falign_parser.add_argument("--method", action="store")
	falign_parser.add_argument("--parameters", action="store")
	falign_parser.add_argument("--min_parsimony_sites", action="store")
	
	falign_args = falign_parser.parse_args(args.arguments)	
	if falign_args.help or len(sys.argv) <= 2: # help is specified
		print(help_message(falign_help))
		sys.exit(0)
	check_required_software()
	if check_required_files(args.command):
		print(now(), "ERROR: Some files are missing preventing this part to run.\nDid your run the required steps before? Missing file:", check_required_files(args.command))
		sys.exit(0)

	cmd = ["snakemake","-p", "-s", "rules/filter-align.smk", "--use-singularity", "-r", "filter_align", "--latency-wait", latency_wait]
	cmd += determine_submission_mode(falign_args.cluster)
	cmd += get_flags(vars(falign_args))
	
	# check for overrides of config file parameters	
	if falign_args.config_file:
		cmd += check_config_file(falign_args.config_file)
	if falign_args.method or falign_args.min_parsimony_sites or falign_args.parameters:
		cmd += ["--config"]
		arg = "trimming={" 
		if falign_args.min_parsimony_sites:
			arg += "min_parsimony_sites: %s, " % falign_args.min_parsimony_sites
		if falign_args.method:
			arg += "method: %s, " % falign_args.method
		if falign_args.parameters:
			arg += "parameters: %s " % falign_args.parameters
		arg += "}"
		cmd += [arg]

	cmd += get_additional_snakemake_flags(falign_args.sm_args, falign_args.rerun)
	cmd += get_additional_singularity_flags(falign_args.si_args)

	for line in execute_command(cmd, falign_args.verbose):
                print(line, end="\r")
	if debug:
		print(now(),"DEBUG:", cmd)
elif args.command=="speciestree":
	print(now(), "Welcome to phylociraptor speciestree v%s" % version)
	sptree_parser = PhyloParser(add_help=False)
	sptree_parser.add_argument("--config-file", action="store", dest="config_file", default="data/config.yaml")
	sptree_args = sptree_parser.parse_args(args.arguments)	
	if sptree_args.help or len(sys.argv) <= 2: # help is specified
		print(help_message(sptree_help))
		sys.exit(0)
	check_required_software()
	if check_required_files(args.command):
		print(now(), "ERROR: Some files are missing preventing this part to run.\nDid your run the required steps before? Missing file:", check_required_files(args.command))
		sys.exit(0)

	cmd = ["snakemake", "-p","-s", "rules/speciestree.smk", "--use-singularity", "-r", "speciestree", "--latency-wait", latency_wait]
	cmd += determine_submission_mode(sptree_args.cluster)
	cmd += get_flags(vars(sptree_args))
	
	# check for overrides of config file parameters	
	if sptree_args.config_file:
		cmd += check_config_file(sptree_args.config_file)
	
	cmd += get_additional_snakemake_flags(sptree_args.sm_args, sptree_args.rerun)
	cmd += get_additional_singularity_flags(sptree_args.si_args)
	
	for line in execute_command(cmd, sptree_args.verbose):
                print(line, end="\r")
	if debug:
		print(now(),"DEBUG:", cmd)
elif args.command=="njtree":
	print(now(), "Welcome to phylociraptor njtree %s" % version)
	njtree_parser = PhyloParser(add_help=False)
	njtree_parser.add_argument("--config-file", action="store", dest="config_file", default="data/config.yaml")
	njtree_args = njtree_parser.parse_args(args.arguments)	
	if njtree_args.help or len(sys.argv) <= 2: # help is specified
		print(help_message(njtree_help))
		sys.exit(0)
	check_required_software()
	if check_required_files(args.command):
		print(now(), "ERROR: Some files are missing preventing this part to run.\nDid your run the required steps before? Missing file:", check_required_files(args.command))
		sys.exit(0)

	cmd = ["snakemake", "-s", "rules/quicktree.smk", "--use-singularity", "-r", "njtree", "--latency-wait", latency_wait]
	cmd += determine_submission_mode(njtree_args.cluster)
	cmd += get_flags(vars(njtree_args))
	
	# check for overrides of config file parameters	
	if njtree_args.config_file:
		cmd += check_config_file(njtree_args.config_file)
	
	cmd += get_additional_snakemake_flags(njtree_args.sm_args, njtree_args.rerun)
	cmd += get_additional_singularity_flags(njtree_args.si_args)
	
	for line in execute_command(cmd, njtree_args.verbose):
                print(line, end="\r")
	if debug:
		print(now(),"DEBUG:", cmd)
elif args.command=="mltree":
	print(now(), "Welcome to phylociraptor mltree v%s" % version)
	tree_parser = PhyloParser(add_help=False)
	tree_parser.add_argument("--config-file", action="store", dest="config_file", default="data/config.yaml")
	tree_args = tree_parser.parse_args(args.arguments)	
	if tree_args.help or len(sys.argv) <= 2: # help is specified
		print(help_message(tree_help))
		sys.exit(0)
	check_required_software()
	if check_required_files(args.command):
		print(now(), "ERROR: Some files are missing preventing this part to run.\nDid your run the required steps before? Missing file:", check_required_files(args.command))
		sys.exit(0)

	cmd = ["snakemake", "-s", "rules/tree.smk", "--use-singularity", "-r", "mltree", "--latency-wait", latency_wait]
	cmd += determine_submission_mode(tree_args.cluster)
	cmd += get_flags(vars(tree_args))
	
	# check for overrides of config file parameters	
	if tree_args.config_file:
		cmd += check_config_file(tree_args.config_file)
	
	cmd += get_additional_snakemake_flags(tree_args.sm_args, tree_args.rerun)
	cmd += get_additional_singularity_flags(tree_args.si_args)
	
	
	for line in execute_command(cmd, tree_args.verbose):
                print(line, end="\r")
	if debug:
		print(now(),"DEBUG:", cmd)
elif args.command=="modeltest":
	print(now(), "Welcome to phylociraptor modeltest v%s" % version)
	model_parser = PhyloParser(add_help=False)
	model_parser.add_argument("--config-file", action="store", dest="config_file", default="data/config.yaml")
	model_args = model_parser.parse_args(args.arguments)	
	if model_args.help or len(sys.argv) <= 2: # help is specified
		print(help_message(model_help))
		sys.exit(0)
	check_required_software()
	if check_required_files(args.command):
		print(now(), "ERROR: Some files are missing preventing this part to run.\nDid your run the required steps before? Missing file:", check_required_files(args.command))
		sys.exit(0)

	cmd = ["snakemake", "-p", "-s", "rules/model.smk", "--use-singularity", "-r", "modeltest", "--latency-wait", latency_wait]
	cmd += determine_submission_mode(model_args.cluster)
	cmd += get_flags(vars(model_args))
	
	# check for overrides of config file parameters	
	if model_args.config_file:
		cmd += check_config_file(model_args.config_file)
	
	cmd += get_additional_snakemake_flags(model_args.sm_args, model_args.rerun)
	cmd += get_additional_singularity_flags(model_args.si_args)

	for line in execute_command(cmd, model_args.verbose):
                print(line, end="\r")
	if debug:
		print(now(),"DEBUG:", cmd)	
elif args.command=="report":
	print(now(), "Welcome to phylociraptor report v%s" % version)
	report_parser = argparse.ArgumentParser(add_help=False)
	report_parser.add_argument("-h", "--help", action="store_true")
	report_parser.add_argument("--verbose", action="store_true", default=False)
	report_parser.add_argument("--figure", action="store_true", default=False)
	report_parser.add_argument("--config-file", action="store", dest="config_file", default="data/config.yaml")
	report_args = report_parser.parse_args(args.arguments)
	print(now(), "Will create a phylociraptor report...")
	
	if report_args.help: # help is specified
		print(help_message(report_help))
		sys.exit(0)
	check_required_software()

	if check_required_files(args.command):
		print(now(), "ERROR: Some files are missing preventing this part to run.\nDid your run the required steps before? Missing file:", check_required_files(args.command))
		sys.exit(0)

	# combine files to produce final statistics file, this is necessary due to the new concurrency setting. 
	#for line in execute_command(cmd, report_args.verbose):
        #        print(line, end="\r")
	#if debug:
	#	print(cmd)
	# check for overrides of config file parameters	
	if report_args.figure:
		cmd = ["singularity", "exec", "-B", os.getcwd(), "docker://reslp/rmarkdown:4.0.3", "Rscript", 'bin/report_figure.R']
		if report_args.config_file:
			cmd += [report_args.config_file]
	else:
		rcommand = """ rmarkdown::render('./bin/report.Rmd') """
		cmd = ["singularity", "exec", "docker://reslp/rmarkdown:4.0.3", "Rscript", "-e", rcommand]
		if report_args.config_file:
			cmd += [report_args.config_file]
	#proc = subprocess.run(cmd)
	for line in execute_command(cmd, report_args.verbose):
                print(line, end="\r")
	if debug:
		print(now(), "DEBUG:", cmd)
	if os.path.isfile("bin/report.html"):
		subprocess.call(["mv", "bin/report.html", "results/report.html"])
		print(now(), "Your phylociraptor report has been created in the results folder.")
	if os.path.isfile("bin/report-figure.pdf"):
		subprocess.call(["mv", "bin/report-figure.pdf", "results/report-figure.pdf"])
		print(now(), "Your phylociraptor report figure has been created in the results folder.")
elif args.command == "check":
	print(now(), "Welcome to phylociraptor check v%s" % version)
	check_parser = argparse.ArgumentParser(add_help=False)
	check_parser.add_argument("-h", "--help", action="store_true")
	check_parser.add_argument("--verbose", action="store_true", default=False)
	check_args = check_parser.parse_args(args.arguments)
	
	if check_args.help: 
		print(help_message(check_help))
		sys.exit(0)

	print(now(), "Checking the status of the pipeline:")

	check_rest = True
	filter_align_done = False
	for mode in steps_to_check:
		if mode == "report":
			continue
		failed = 0
		if filter_align_done and not check_rest: # in case filter_align is already finished, check the rest of the parts.
			check_rest = True 
		if debug:
			print("\n(Debug) Checking runmode:", mode)
		if check_hasrun(mode) and check_rest == True:
			failed = 1
			check_rest = False
			if debug:
				print("After check_hasrun: Failed=", failed, "check_rest=", check_rest, "filter_align_done=", filter_align_done)
		if check_required_files(mode) and check_rest == True:
			failed = 1
			check_rest = False
			if debug:
				print("After check_required files: Failed=", failed, "check_rest=", check_rest, "filter_align_done=", filter_align_done)
		if check_directories(mode) and check_rest == True:
			failed = 1
			check_rest  = False
			if debug:
				print("After check_directories: Failed=", failed, "check_rest=", check_rest, "filter_align_done=", filter_align_done)
		# if non of the checkpoint files are missing, lets check if the dataset has changed:
		if failed == 0 and check_rest == True:
			last_files = []
			for n in outfile_dict[mode]:
				last_files.append(os.path.getmtime(n))
			last = max(last_files)
			current = os.path.getmtime(checkpoint_file_dict[mode])
			if debug: 
				print("Current file timestamp:", current)
				print("Last file timestamp:", last)
			if last > current and mode != "setup":	
				failed = 1
				check_rest = False
				if debug:
					print("Due to timestamp. Failed=", failed, "check_rest=", check_rest, "filter_align_done=", filter_align_done)
		if failed == 0 and check_rest == False and filter_align_done == False:
			failed = 2
		if failed == 1: 
			print(" ", mode, " ...", '\033[91m' , "INCOMPLETE", '\033[0m')
			if mode == "orthology":
				missing=""
				if not os.path.isdir("results/assemblies"):
					print("\tDirectory not found: results/assemblies")
					continue
				nassemblies = 0
				for assembly in os.listdir("results/assemblies"):
					assembly = os.path.splitext(assembly)[0]
					if ".fna" in assembly:
						assembly = assembly.split(".fna")[0]
					if not os.path.isfile("results/checkpoints/busco/busco_"+ assembly + ".done"):
						missing += assembly +","
					nassemblies += 1
				missing = missing.strip(",")
				if missing:
					print("\t Missing orthology results for", len(missing.split(",")), "of", nassemblies, "genomes.")
					if check_args.verbose:
						print("\tOrthology step was not successful for these genomes:")
						print("\t", missing)		
			if mode == "align":
				missing = ""
				genes = []
				# first get all sequences files after filter-orthology
				if not os.path.isdir("results/alignments"):
					continue
				for gene in os.listdir("results/orthology/busco/busco_sequences_deduplicated"):
					gene = os.path.splitext(gene)[0]
					genes.append(gene.split("_")[0])
				# now check if the same files are in align
				total_genes = len(genes)
				for directory in os.listdir("results/alignments/full"):
					for gene in os.listdir("results/alignments/full/"+directory):
						gene = os.path.splitext(gene)[0]
						gene = gene.split("_")[0]
						if gene in genes:
							genes.remove(gene)	
					missing = ",".join(genes)
					missing = missing.strip(",")
					if missing:
						print("\tMissing", directory, "alignments for", len(missing.split(",")), "of", total_genes, "genes.")
						if check_args.verbose:
							print("\t These alignments are missing:")
							print("\t", missing)
			if mode == "filter-align":
				# first look for missing trimmed files
				if not os.path.isdir("results/alignments/trimmed"):
					continue
				aligners = os.listdir("results/alignments/full")
				for aligner in aligners:
					for directory in glob.glob("results/alignments/trimmed/"+aligner+"-*"):
						missing = ""
						genes = []
						for gene in os.listdir("results/alignments/full/"+aligner):
							gene = os.path.splitext(gene)[0]
							genes.append(gene.split("_")[0])
						total_genes = len(genes)
						genes_all = genes
						for gene in glob.glob(directory+"/*.fas"):
							gene = os.path.basename(gene)
							gene = os.path.splitext(gene)[0]
							gene = gene.split("_")[0]
							if gene in genes_all:
								genes_all.remove(gene)	
						missing = ",".join(genes_all)
						missing = missing.strip(",")
						if missing:
							print("\tMissing",aligner,directory.split("-")[-1],"trimmed alignments for", len(missing.split(",")), "of", total_genes, "genes.")
							if check_args.verbose:
								print("\t These trimmed alignments are missing:")
								print("\t", missing)
				# second check missing filtered files (after trimming)
				if not os.path.isdir("results/alignments/filtered"):
					continue
				aligners = os.listdir("results/alignments/trimmed")
				for aligner in aligners:
					for directory in glob.glob("results/alignments/filtered/"+aligner+"-*"):
						missing = ""
						genes = []
						for gene in os.listdir("results/alignments/trimmed/"+aligner):
							gene = os.path.splitext(gene)[0]
							genes.append(gene.split("_")[0])
						total_genes = len(genes)
						genes_all = genes
						for gene in glob.glob(directory+"/*.fas"):
							gene = os.path.basename(gene)
							gene = os.path.splitext(gene)[0]
							gene = gene.split("_")[0]
							if gene in genes_all:
								genes_all.remove(gene)	
						missing = ",".join(genes_all)
						missing = missing.strip(",")
						if missing:
							print("\tMissing",aligner,directory.split("-")[-1],"filtered alignments for", len(missing.split(",")), "of", total_genes, "genes.")
							print("(Automatically filtered alignments will also show up here)")
							if check_args.verbose:
								print("\t These filtered alignments are missing:")
								print("\t", missing)
			if mode == "modeltest":
				if not os.path.isdir("results/checkpoints/modeltest"):
					continue
				units = os.listdir("results/alignments/filtered")
				for unit in units:
					if not os.path.isdir("results/checkpoints/modeltest/"+unit):
						continue
					genes = []
					missing = ""
					for gene in os.listdir("results/alignments/filtered/"+unit):
						gene = os.path.splitext(gene)[0]
						genes.append(gene.split("_")[0])
					total_genes = len(genes)
					for gene in os.listdir("results/checkpoints/modeltest/"+unit):
						gene = os.path.splitext(gene)[0]
						gene = gene.split("_")[0]
						if gene in genes:
							genes.remove(gene)	
					missing = ",".join(genes)
					missing = missing.strip(",")
					if missing:
						print("\tMissing",unit,"modeltest results for", len(missing.split(",")), "of", total_genes, "genes.")
						if check_args.verbose:
							print("\t These modeltest results are missing:")
							print("\t", missing)
			if mode == "speciestree":
				if not os.path.isdir("results/checkpoints/gene_trees"):
					continue
				units = os.listdir("results/alignments/filtered")
				for unit in units:
					missing = ""
					genes = []
					if not os.path.isdir("results/checkpoints/gene_trees/"+unit):
						continue
					for gene in os.listdir("results/alignments/filtered/"+unit):
						gene = os.path.splitext(gene)[0]
						genes.append(gene.split("_")[0])
					total_genes = len(genes)
					for gene in os.listdir("results/checkpoints/gene_trees/"+unit):
						gene = os.path.splitext(gene)[0]
						gene = gene.split("_")[0]
						if gene in genes:
							genes.remove(gene)	
					missing = ",".join(genes)
					missing = missing.strip(",")
					if missing:
						print("\tMissing",unit,"gene tree for", len(missing.split(",")), "of", total_genes, "genes.")
						if check_args.verbose:
							print("\t These gene trees are missing:")
							print("\t", missing)
						
				
		elif failed == 2:	
			print(" ", mode, " ...", '\033[93m', "NOT EVALUATED", '\033[0m', "(preceeding step not finished)")
		else:
			if mode == "filter-align":
				filter_align_done = True
			print(" ", mode," ...", '\033[92m', "DONE", '\033[0m')
	print("\nWARNING: phylociraptor check is just a quick and shallow verification of the run. In case you run into problems, please also check logfiles in the log directory for more in-depth diagnostics.")  

elif args.command == "util":
	print(now(), "Welcome to phylociraptor util v%s" % version)
	#util_parser = argparse.ArgumentParser(add_help=False)
	#util_parser.add_argument("-h", "--help", action="store_true")
	#util_parser.add_argument("command", action="store", nargs="?")
	#util_parser.add_argument("--verbose", action="store_true", default=False)
	#util_args = util_parser.parse_args(args.arguments)

	if  len(args.arguments) == 0: 
		print(help_message(util_help))
		sys.exit(0)
	elif args.arguments[0] == "-h" or args.arguments[0] == "--help":
		print(help_message(util_help))
		sys.exit(0)
	else:
		which_util = args.arguments.pop(0)
	
	if which_util == "get-lineage":
		gl_parser = UtilParser(add_help=False)
		gl_parser.add_argument("-g","--genomefile", action="store", default="results/statistics/downloaded_genomes_statistics.txt")
		gl_parser.add_argument("-o", "--outfile", action="store") 
		gl_parser.add_argument("--quiet", action="store_true", default=False) 
		gl_args = gl_parser.parse_args(args.arguments) 
		if gl_args.help: 
			print(help_message(util_lineage_help))
			sys.exit(0)
		if not gl_args.outfile:
			print(now(), "You need to specify an output file with -o/--outfile")
			sys.exit(1)
		if not os.path.isfile(gl_args.genomefile):
			print(now(), "Genome file not found:",gl_args.genomefile)
			sys.exit(1)
		cmd = ["singularity", "exec", "-B", os.getcwd(), "docker://reslp/biopython_plus:1.77", "python", "bin/get_lineage.py", gl_args.genomefile,  gl_args.outfile]
		if debug:
			print(cmd)
		for line in execute_command(cmd, not gl_args.quiet):
			print(line, end="\r")
	if which_util == "estimate-conflict":
		qs_parser = UtilParser(add_help=False)
		qs_parser.add_argument("-i","--intrees", action="store", default="all")
		qs_parser.add_argument("-o", "--outfile", action="store") 
		qs_parser.add_argument("-s", "--seed", action="store", default="random") 
		qs_parser.add_argument("-n", "--nquartets", action="store", default=None)
		qs_parser.add_argument("-t", "--threads", action="store", default="1")
		qs_parser.add_argument("-l", "--lineagefile", action="store")
		qs_parser.add_argument("-b", "--stopby", action="store", default=None)
		qs_parser.add_argument("-a", "--selecttaxa", action="store")
		qs_parser.add_argument("--quiet", action="store_true", default=False) 
		qs_args = qs_parser.parse_args(args.arguments) 
		if qs_args.help: 
			print(help_message(util_estimate_conflict_help))
			sys.exit(0)
		if not qs_args.outfile:
			print(now(),"You need to specify an output file with -o/--outfile")
			sys.exit(1)
		if not qs_args.intrees:
			print(now(),"No input trees specified")
			sys.exit(1)
		elif qs_args.intrees == "all":
			print(now(), "Will compare all trees, this expects phylociraptor mltree and/or speciestree to be finished.")
			iqtree_trees = glob.glob("results/phylogeny-*/iqtree/*/concat.contree")
			raxml_trees =  glob.glob("results/phylogeny-*/raxml/*/raxmlng.raxml.support")
			astral_trees = glob.glob("results/phylogeny-*/astral/*/species_tree.tre")
			all_trees = ",".join(iqtree_trees + raxml_trees + astral_trees)
		else:
			all_trees = qs_args.intrees	
		cmd = ["singularity", "exec", "-B", os.getcwd(), "docker://reslp/phylopy:2", "python", "bin/estimate_conflict.py", "-i", all_trees,  "-o", qs_args.outfile, "-s", qs_args.seed, "-t", qs_args.threads]

		if qs_args.lineagefile and qs_args.selecttaxa:
			cmd += ["-l", qs_args.lineagefile, "--selecttaxa", qs_args.selecttaxa]
		if qs_args.nquartets:
			cmd += ["-q", qs_args.nquartets]
		elif qs_args.stopby:
			cmd += ["-b", qs_args.stopby]
		else:
			print(now(), "You need to either specify a stopping criterion (-b) or the total number of quartets (-n).")
			sys.exit(1)
			

		if debug:
			print(cmd)
		for line in execute_command(cmd, not qs_args.quiet):
			print(line, end="\r")
	if which_util == "plot-tree":
		qs_parser = UtilParser(add_help=False)
		qs_parser.add_argument("-i","--intrees", action="store", default="all")
		qs_parser.add_argument("-o", "--outprefix", action="store", default="tree-plot") 
		qs_parser.add_argument("-g", "--outgroup", action="store", default="none")
		qs_parser.add_argument("-l", "--lineagefile", action="store", default="none")
		qs_parser.add_argument("-e", "--level", action="store", default="none")
		qs_parser.add_argument("-s", "--seed", action="store", default="random") 
		qs_parser.add_argument("--quiet", action="store_true", default=False) 
		qs_args = qs_parser.parse_args(args.arguments) 
		if qs_args.help: 
			print(help_message(util_plot_tree_help))
			sys.exit(0)
		if not qs_args.intrees:
			print(now(),"No input trees specified")
			sys.exit(1)
		elif qs_args.intrees == "all":
			print(now(), "Will compare all trees, this expects phylociraptor mltree and/or speciestree to be finished.")
			iqtree_trees = glob.glob("results/phylogeny-*/iqtree/*/concat.contree")
			raxml_trees =  glob.glob("results/phylogeny-*/raxml/*/raxmlng.raxml.support")
			astral_trees = glob.glob("results/phylogeny-*/astral/*/species_tree.tre")
			all_trees = ",".join(iqtree_trees + raxml_trees + astral_trees)
		else:
			all_trees = qs_args.intrees	
		
		cmd = ["singularity", "exec", "-B", os.getcwd(), "docker://reslp/rphylogenetics:4.0.3", "Rscript", "bin/plot-tree.R"]
		cmd += ["plot", os.getcwd(), qs_args.seed, all_trees, qs_args.outgroup, qs_args.lineagefile, qs_args.level, qs_args.outprefix]

		if debug:
			print(cmd)
		for line in execute_command(cmd, not qs_args.quiet):
			print(line, end="\r")
	if which_util == "plot-conflict":
		qs_parser = UtilParser(add_help=False)
		qs_parser.add_argument("-i","--intrees", action="store", default=None)
		qs_parser.add_argument("-g", "--outgroup", action="store", default="none")
		qs_parser.add_argument("-l", "--lineagefile", action="store", default="none")
		qs_parser.add_argument("-q", "--quartetfile", action="store", default=None)
		qs_parser.add_argument("-e", "--level", action="store", default="none")
		qs_parser.add_argument("-s", "--seed", action="store", default="random") 
		qs_parser.add_argument("-r", "--treelist", action="store", default=None)
		qs_parser.add_argument("--quiet", action="store_true", default=False) 
		qs_args = qs_parser.parse_args(args.arguments) 
		if qs_args.help: 
			print(help_message(util_plot_conflict_help))
			sys.exit(0)
		if not qs_args.intrees:
			print(now(),"No input trees specified")
			sys.exit(1)
		elif len(qs_args.intrees.split(",")) != 2:
			print(now(), "Need exactly two trees for comparison, seperated by a comma (,).")
			sys.exit(1)
		else:
			all_trees = qs_args.intrees	
		if not qs_args.quartetfile:
			print(now(), "You need to provide a quartet conflict file.")
			sys.exit(1)
		if not qs_args.treelist:
			print(now(), "You need to provide the treelist file.")
			sys.exit(1)
		
		cmd = ["singularity", "exec", "-B", os.getcwd(), "docker://reslp/rphylogenetics:4.0.3", "Rscript", "bin/plot-tree.R"]
		cmd += ["conflicts", os.getcwd(), qs_args.seed, all_trees, qs_args.outgroup, qs_args.lineagefile, qs_args.level, qs_args.quartetfile, qs_args.treelist]

		if debug:
			print(cmd)
		for line in execute_command(cmd, not qs_args.quiet):
			print(line, end="\r")
	if which_util == "plot-similarity":
		qs_parser = UtilParser(add_help=False)
		qs_parser.add_argument("-m","--simmatrix", action="store", default=None)
		qs_parser.add_argument("-r", "--treelist", action="store", default="none")
		qs_parser.add_argument("--quiet", action="store_true", default=False)
		qs_parser.add_argument("-t", "--threads", action="store", default="1") 
		qs_parser.add_argument("-n", "--ndistances", action="store", default="all")
		qs_parser.add_argument("-s", "--seed", action="store", default="random")
		qs_args = qs_parser.parse_args(args.arguments) 
		print(qs_args)
		if qs_args.help: 
			print(help_message(util_plot_similarity_help))
			sys.exit(0)
		if not qs_args.simmatrix:
			print(now(), "You have to the similarity matrix file created with phylociraptor util estimate-conflict")
			sys.exit(1)	
		cmd = ["singularity", "exec", "-B", os.getcwd(), "docker://reslp/rphylogenetics:4.0.3", "Rscript", "bin/plot-similarity-matrix.R", os.getcwd(), qs_args.simmatrix, qs_args.treelist, qs_args.threads, qs_args.ndistances, qs_args.seed]
		if debug:
			print(cmd)
		for line in execute_command(cmd, not qs_args.quiet):
			print(line, end="\r")
else:
	print("Runmode not recognized: %s" % args.command)
	print("Please run phylociraptor -h to see avilable options")

print(now(), "DONE")