reads_first.py

#!/usr/bin/env python

import argparse,os,sys,importlib,shutil,subprocess

helptext="""
This script is a wrapper around several scripts in the HybSeqPipeline.
It can check whether you have the appropriate dependencies available (see --check-depend).
It makes sure that the other scripts needed are in the same directory as this one.
Command line options are passed to the other executables.
Unless --prefix is set, output will be put within a directory named after your read files."""

velvet_genefilename = "velvet_genelist.txt"
cap3_genefilename = "cap3_genelist.txt"
exonerate_genefilename = "exonerate_genelist.txt"
spades_genefilename = "spades_genelist.txt"


def py_which(cmd, mode=os.F_OK | os.X_OK, path=None):
    """Given a command, mode, and a PATH string, return the path which
    conforms to the given mode on the PATH, or None if there is no such
    file.

    `mode` defaults to os.F_OK | os.X_OK. `path` defaults to the result
    of os.environ.get("PATH"), or can be overridden with a custom search
    path.

    """
    # Check that a given file can be accessed with the correct mode.
    # Additionally check that `file` is not a directory, as on Windows
    # directories pass the os.access check.
    def _access_check(fn, mode):
        return (os.path.exists(fn) and os.access(fn, mode)
                and not os.path.isdir(fn))

    # If we're given a path with a directory part, look it up directly rather
    # than referring to PATH directories. This includes checking relative to the
    # current directory, e.g. ./script
    if os.path.dirname(cmd):
        if _access_check(cmd, mode):
            return cmd
        return None

    if path is None:
        path = os.environ.get("PATH", os.defpath)
    if not path:
        return None
    path = path.split(os.pathsep)

    if sys.platform == "win32":
        # The current directory takes precedence on Windows.
        if not os.curdir in path:
            path.insert(0, os.curdir)

        # PATHEXT is necessary to check on Windows.
        pathext = os.environ.get("PATHEXT", "").split(os.pathsep)
        # See if the given file matches any of the expected path extensions.
        # This will allow us to short circuit when given "python.exe".
        # If it does match, only test that one, otherwise we have to try
        # others.
        if any([cmd.lower().endswith(ext.lower()) for ext in pathext]):
            files = [cmd]
        else:
            files = [cmd + ext for ext in pathext]
    else:
        # On other platforms you don't have things like PATHEXT to tell you
        # what file suffixes are executable, so just pass on cmd as-is.
        files = [cmd]

    seen = set()
    for dir in path:
        normdir = os.path.normcase(dir)
        if not normdir in seen:
            seen.add(normdir)
            for thefile in files:
                name = os.path.join(dir, thefile)
                if _access_check(name, mode):
                    return name
    return None



def check_dependencies():
	"""Checks for the presence of executables and Python packages"""
	executables = 		["blastx",
						"exonerate",
						#"velvetg",
						#"velveth",
						#"cap3",
						"parallel",
						"makeblastdb",
						"spades.py",
						"bwa",
						"samtools"]

	python_packages = 	["Bio"]
	
	everything_is_awesome = True
	for e in executables:
		e_loc = py_which(e)
		if e_loc:
			print "{} found at {}".format(e,e_loc)
		else:
			print "{} not found in your $PATH!".format(e)
			everything_is_awesome = False
			
	for p in python_packages:
		try:
			i = importlib.import_module(p)
			print "Package {} successfully loaded!".format(p)
		except ImportError:
			print "Package {} not found!".format(p)
			everything_is_awesome=False
	return everything_is_awesome


def blastx(readfiles,baitfile,evalue,basename,cpu=None,max_target_seqs=10,unpaired=False):
	dna = set("ATCGN")
	if os.path.isfile(baitfile):
		#Quick detection of whether baitfile is DNA. 
		with open(baitfile) as bf:
			header = bf.readline()
			seqline = bf.readline().rstrip().upper()
			if not set(seqline) - dna:
				print "ERROR: only ATCGN characters found in first line. You need a protein bait file for BLASTx!"
				return None

		if os.path.isfile(os.path.split(baitfile)[0]+'.psq'):
			db_file = baitfile
		else:
			print "Making protein blastdb in current directory."
			if os.path.split(baitfile)[0]:
				shutil.copy(baitfile,'.')
			db_file = os.path.split(baitfile)[1]
			makeblastdb_cmd = "makeblastdb -dbtype prot -in {}".format(db_file)
			print makeblastdb_cmd
			exitcode = subprocess.call(makeblastdb_cmd,shell=True)
			if exitcode:
				return None
	else:
		print "Cannot find baitfile at: {}".format(baitfile)
		return None

	#Remove previous blast results if they exist (because we will be appending)
	if os.path.isfile(basename+".blastx"):
		os.remove(basename+".blastx")
	
	if unpaired:
		read_file = readfiles
		pipe_cmd = "cat {} |  awk '{{if(NR % 4 == 1 || NR % 4 == 2) {{sub(/@/, \">\"); print; }} }}'".format(read_file)
		blastx_command = "blastx -db {} -query - -evalue {} -outfmt 6 -max_target_seqs {}".format(db_file,evalue,max_target_seqs)
		if cpu:
			full_command = "time {} | parallel -j {} -k --block 200K --recstart '>' --pipe '{}' >> {}_unpaired.blastx ".format(pipe_cmd,cpu,blastx_command,basename)
		else:
			full_command = "time {} | parallel -k --block 200K --recstart '>' --pipe '{}' >> {}_unpaired.blastx ".format(pipe_cmd,blastx_command,basename)
		print full_command
		exitcode = subprocess.call(full_command,shell=True)
		if exitcode:
			#Concatenate the two blastfiles.
		
			return None
		return basename + "_unpaired.blastx"
		
	else:
		for read_file in readfiles:

		#Piping commands for Fastq -> FASTA	
		# Curly braces must be doubled within a formatted string.
			pipe_cmd = "cat {} |  awk '{{if(NR % 4 == 1 || NR % 4 == 2) {{sub(/@/, \">\"); print; }} }}'".format(read_file)
	
			blastx_command = "blastx -db {} -query - -evalue {} -outfmt 6 -max_target_seqs {}".format(db_file,evalue,max_target_seqs)
			if cpu:
				full_command = "time {} | parallel -j {} -k --block 200K --recstart '>' --pipe '{}' >> {}.blastx ".format(pipe_cmd,cpu,blastx_command,basename)
			else:
				full_command = "time {} | parallel -k --block 200K --recstart '>' --pipe '{}' >> {}.blastx ".format(pipe_cmd,blastx_command,basename)
			print full_command
			exitcode = subprocess.call(full_command,shell=True)
			if exitcode:
				#Concatenate the two blastfiles.
		
				return None
			
	return basename + '.blastx'

def distribute(blastx_outputfile,readfiles,baitfile,run_dir,target=None,unpaired_readfile=None):
	#NEED TO ADD SOMETHING ABOUT DIRECTORIES HERE.
	#print run_dir
	read_cmd = "time python {} {} {}".format(os.path.join(run_dir,"distribute_reads_to_targets.py"),blastx_outputfile," ".join(readfiles))
	exitcode = subprocess.call(read_cmd,shell=True)
	if exitcode:
		print "ERROR: Something went wrong with distributing reads to gene directories."
		return exitcode
	target_cmds = ["time python", os.path.join(run_dir,"distribute_targets.py"),baitfile,"--blastx",blastx_outputfile]
	if target:
		target_cmds.append("--target {}".format(target))
	if unpaired_readfile:
		blastx_outputfile = blastx_outputfile.replace(".blastx","_unpaired.blastx")
		unpaired_cmd = "time python {} {} {}".format(os.path.join(run_dir,"distribute_reads_to_targets.py"),blastx_outputfile,unpaired_readfile)
		print("[CMD] {}\n".format(unpaired_cmd))
		exitcode = subprocess.call(unpaired_cmd,shell=True)
	target_cmd = " ".join(target_cmds)
	exitcode = subprocess.call(target_cmd,shell=True)
	if exitcode:
		print "ERROR: Something went wrong distributing targets to gene directories."
		return exitcode
	return None

def distribute_bwa(bamfile,readfiles,baitfile,run_dir,target=None,unpaired=None):
	#NEED TO ADD SOMETHING ABOUT DIRECTORIES HERE.
	#print run_dir
	read_cmd = "time python {} {} {}".format(os.path.join(run_dir,"distribute_reads_to_targets_bwa.py"),bamfile," ".join(readfiles))
	print("[CMD] {}\n".format(read_cmd))
	exitcode = subprocess.call(read_cmd,shell=True)
	
	if unpaired:
		bamfile = bamfile.replace(".bam","_unpaired.bam")
		unpaired_cmd = "time python {} {} {}".format(os.path.join(run_dir,"distribute_reads_to_targets_bwa.py"),bamfile,unpaired)
		print("[CMD] {}\n".format(unpaired_cmd))
		exitcode = subprocess.call(unpaired_cmd,shell=True)
	
	if exitcode:
		print "ERROR: Something went wrong with distributing reads to gene directories."
		return exitcode
	target_cmds = ["time python", os.path.join(run_dir,"distribute_targets.py"),baitfile,"--bam",bamfile]
	if target:
		target_cmds.append("--target {}".format(target))
	target_cmd = " ".join(target_cmds)
	exitcode = subprocess.call(target_cmd,shell=True)
	if exitcode:
		print "ERROR: Something went wrong distributing targets to gene directories."
		return exitcode
	return None



def make_basename(readfiles,prefix=None):
	"""Unless prefix is set, generate a directory based off the readfiles, using everything up to the first underscore.
        If prefix is set, generate the directory "prefix" and set basename to be the last component of the path.
        
        """
	if prefix:
            if not os.path.exists(prefix):
                os.makedirs(prefix)
#             else:
#                 if os.listdir(prefix):
#                 ## When running under HTCondor, the execution might be restarted
#                 ## However, this code in not tolerant of restarts
#                     print "ERROR: Directory {} not empty, exiting!".format(prefix)
#                     return

            prefixParentDir, prefix = os.path.split(prefix)
            if not prefix:
                # if prefix has a trailing /, prefixParentDir will have the / stripped and prefix will be empty.
                # so try again
                prefix = os.path.split(prefixParentDir)[1]


            return prefixParentDir,prefix

        ## --prefix is not set on cmd line;  Write output to subdir in . 
	basename = os.path.split(readfiles[0])[1].split('_')[0]

	if not os.path.exists(basename):
		os.makedirs(basename)
	return '.', basename

def spades(genes,run_dir,cov_cutoff=8,cpu=None,paired=True,kvals=None,timeout=None,unpaired=False):
	"Run SPAdes on each gene separately using GNU paralell."""
	
#	import spades_runner

 	with open(spades_genefilename,'w') as spadesfile:
 		spadesfile.write("\n".join(genes)+"\n")
	
	if os.path.isfile("spades.log"):
		os.remove("spades.log")
	if os.path.isfile("spades_redo.log"):
		os.remove("spades_redo.log")

	spades_runner_list = ["python","{}/spades_runner.py".format(run_dir),spades_genefilename,"--cov_cutoff",str(cov_cutoff)]
	if cpu:
		spades_runner_list.append("--cpu")
		spades_runner_list.append(str(cpu))
	if not paired:
		spades_runner_list.append("--single")
	if unpaired:
		spades_runner_list.append("--unpaired")
	if timeout:
		spades_runner_list.append("--timeout")
		spades_runner_list.append("{}%".format(timeout))
	
	
	spades_runner_cmd = " ".join(spades_runner_list)	

# 	if cpu:
# 		if paired:
# 			spades_runner_cmd = "python {} {} --cpu {} --cov_cutoff {}".format(os.path.join(run_dir,"spades_runner.py"),spades_genefilename, cpu, cov_cutoff)
# 		else:
# 			spades_runner_cmd = "python {} {} --cpu {} --cov_cutoff {} --single".format(os.path.join(run_dir,"spades_runner.py"),spades_genefilename, cpu, cov_cutoff)
# 	else:
# 		if paired:
# 			spades_runner_cmd = "python {} {} --cov_cutoff {}".format(os.path.join(run_dir,"spades_runner.py"),spades_genefilename, cov_cutoff)
# 		else:
# 			spades_runner_cmd = "python {} {} --cov_cutoff {} --single".format(os.path.join(run_dir,"spades_runner.py"),spades_genefilename, cov_cutoff)
	exitcode = subprocess.call(spades_runner_cmd,shell=True)
	if exitcode:
		sys.stderr.write("WARNING: Something went wrong with the assemblies! Check for failed assemblies and re-run! \n")
		return None
	else:
		if os.path.isfile("spades_duds.txt"):
			spades_duds = [x.rstrip() for x in open("spades_duds.txt")]
		else:
			spades_duds = []



# 
# 	spades_failed = spades_runner.spades_initial(spades_genefilename,cov_cutoff,cpu,kvals)
# 	if len(spades_failed) > 0:
# 		with open("failed_spades.txt",'w') as failed_spadefile:
# 			failed_spadefile.write("\n".join(spades_failed))
# 	
# 		spades_failed_redos,spades_duds = spades_runner.rerun_spades("failed_spades.txt")
# 		if len(spades_failed) == 0:
# 			sys.stderr.write("All redos completed successfully!\n")
# 	
 	spades_genelist = []
	for gene in genes:
#		if gene not in set(spades_failed):
		if gene not in set(spades_duds):
#				if gene not in set(spades_failed_redos):
			spades_genelist.append(gene)
	
	with open(exonerate_genefilename,'w') as genefile:
		genefile.write("\n".join(spades_genelist)+"\n")

	return spades_genelist	

# def spades(genes,cov_cutoff=8,cpu=None,paired=True,kvals=None):
# 	"Run SPAdes on each gene separately using GNU paralell."""
# 	if os.path.isfile("spades.log"):
# 		os.remove("spades.log")
# 	
# 	with open(spades_genefilename,'w') as spadesfile:
# 		spadesfile.write("\n".join(genes)+"\n")
# 	
# 	if paired:
# 		fileflag = "--12"
# 	else:
# 		fileflag = "-s"
# 	
# 	if kvals:
# 		kvals = ",".join(kvals)
# 	
# 	if cpu:
# 		if kvals:
# 			spades_cmd = "time parallel -j {} --eta spades.py --only-assembler -k {} --threads 1 --cov-cutoff {} {} {{}}/{{}}_interleaved.fasta -o {{}}/{{}}_spades :::: {} > spades.log".format(cpu,kvals,cov_cutoff,fileflag,spades_genefilename)
# 		else:
# 			spades_cmd = "time parallel -j {} --eta spades.py --only-assembler --threads 1 --cov-cutoff {} {} {{}}/{{}}_interleaved.fasta -o {{}}/{{}}_spades :::: {} > spades.log".format(cpu,cov_cutoff,fileflag,spades_genefilename)
# 	else:
# 		if kvals:
# 			spades_cmd = "time parallel --eta spades.py --only-assembler -k {} --threads 1 --cov-cutoff {} {} {{}}/{{}}_interleaved.fasta -o {{}}/{{}}_spades :::: {} > spades.log".format(kvals,cov_cutoff,fileflag,spades_genefilename)
# 		else:
# 			spades_cmd = "time parallel --eta spades.py --only-assembler --threads 1 --cov-cutoff {} {} {{}}/{{}}_interleaved.fasta -o {{}}/{{}}_spades :::: {} > spades.log".format(cov_cutoff,fileflag,spades_genefilename)
# 		
# 	sys.stderr.write("Running SPAdes on {} genes\n".format(len(genes)))
# 	sys.stderr.write(spades_cmd + "\n")
# 	exitcode = subprocess.call(spades_cmd,shell=True)
# 	
# 	#Need to handle an error differently with SPAdes, which can fail if there are simply not enough reads. 
# 	
# 	if exitcode:
# 		sys.stderr.write("ERROR: One or more genes had an error with SPAdes assembly. This may be due to low coverage. No contigs found for the following genes:\n")
# 	
# 	spades_genelist = []
# 	
# 	for gene in genes:
# 		if os.path.isfile("{}/{}_spades/contigs.fasta".format(gene,gene)):
# 			shutil.copy("{}/{}_spades/contigs.fasta".format(gene,gene),"{}/{}_contigs.fasta".format(gene,gene))
# 			spades_genelist.append(gene)
# 		else:
# 			sys.stderr.write("{}\n".format(gene))
# 	
# 	with open(exonerate_genefilename,'w') as genefile:
# 		genefile.write("\n".join(spades_genelist)+"\n")
# 
# 	return spades_genelist	
		

def velvet(genes,cov_cutoff=5,ins_length=200,kvals = ["21","31","41","51","61"],cpu=None,paired=True):
	"""Use parallel to run velveth and velvetg on a set of k values on every gene with blastx hits from the previous steps."""
	if os.path.isfile('velveth.log'):
		os.remove('velveth.log')
	if os.path.isfile('velvetg.log'):
		os.remove('velvetg.log')
	
	#Write a file with the list of genes for velvet, to avoid screen spam and hitting the argument cap.
	with open(velvet_genefilename,'w') as genefile:
		genefile.write("\n".join(genes)+"\n")
	
	if paired:
		if cpu:
			velveth_cmd = "time parallel  -j {} --eta velveth {{1}}/velvet{{2}} {{2}} -shortPaired {{1}}/{{1}}_interleaved.fasta '>>' velveth.log :::: {} ::: {}".format(cpu,velvet_genefilename," ".join(kvals))
		else:
			velveth_cmd = "time parallel  --eta velveth {{1}}/velvet{{2}} {{2}} -shortPaired {{1}}/{{1}}_interleaved.fasta '>>' velveth.log :::: {} ::: {}".format(velvet_genefilename," ".join(kvals))
	else:
		if cpu:
			velveth_cmd = "time parallel  -j {} --eta velveth {{1}}/velvet{{2}} {{2}} -short {{1}}/{{1}}_interleaved.fasta '>>' velveth.log :::: {} ::: {}".format(cpu,velvet_genefilename," ".join(kvals))
		else:
			velveth_cmd = "time parallel  --eta velveth {{1}}/velvet{{2}} {{2}} -short {{1}}/{{1}}_interleaved.fasta '>>' velveth.log :::: {} ::: {}".format(velvet_genefilename," ".join(kvals))
	

	#print os.getcwd()
	#print os.listdir(".")
	
	
	print "Running velveth on {} genes".format(len(genes))
	print velveth_cmd
	exitcode = subprocess.call(velveth_cmd,shell=True)
	if exitcode:
		print "ERROR: Something went wrong with velveth!"
		return exitcode
	if cpu:
		velvetg_cmd = "time parallel -j {} --eta velvetg {{1}}/velvet{{2}} -ins_length {} -cov_cutoff {} '>>' velvetg.log :::: {} ::: {}".format(cpu,ins_length,cov_cutoff,velvet_genefilename," ".join(kvals))
	else:
		velvetg_cmd = "time parallel --eta velvetg {{1}}/velvet{{2}} -ins_length {} -cov_cutoff {} '>>' velvetg.log :::: {} ::: {}".format(ins_length,cov_cutoff,velvet_genefilename," ".join(kvals))
	print "Running velvetg on {} genes".format(len(genes))
	print velvetg_cmd
	exitcode = subprocess.call(velvetg_cmd,shell=True)
	if exitcode:
		print "ERROR: Something went wrong with velvetg!"
		return exitcode
	with open(exonerate_genefilename,'w') as genefile:
		genefile.write("\n".join([x for x in genes if os.path.getsize(os.path.join(x,'velvet_contigs.fa')) > 0]))	
	return None
	
def cap3(genes,cpu=None):
	print "Concatenating velvet output"
	
	if cpu:
		cat_cmd = "time parallel -j {} cat {{1}}/*/contigs.fa '>' {{1}}/velvet_contigs.fa :::: {}".format(cpu,velvet_genefilename)
	else:
		cat_cmd = "time parallel cat {{1}}/*/contigs.fa '>' {{1}}/velvet_contigs.fa :::: {}".format(velvet_genefilename)
	print cat_cmd
	exitcode = subprocess.call(cat_cmd,shell=True)
	if exitcode:
		print "ERROR: Something went wrong while concatenating the velvet output files."
		return exitcode
	#Check that the velvet output actually has data in it.
	genes = [x for x in genes if os.path.getsize(os.path.join(x,'velvet_contigs.fa')) > 0]
	
	with open(cap3_genefilename,'w') as genefile:
		genefile.write("\n".join(genes)+"\n")
	
	if len(genes) == 0:
		print "No genes left! Exiting!"
		return 1
		
	print "Running Cap3 on {} genes".format(len(genes))
	if cpu:
		cap3_cmd = "time parallel -j {} --eta cap3 {{1}}/velvet_contigs.fa -o 20 -p 99 '>' {{1}}/{{1}}_cap3.log :::: {}".format(cpu,cap3_genefilename)
	else:
		cap3_cmd = "time parallel --eta cap3 {{1}}/velvet_contigs.fa -o 20 -p 99 '>' {{1}}/{{1}}_cap3.log :::: {}".format(cap3_genefilename)
	print cap3_cmd
	exitcode = subprocess.call(cap3_cmd,shell=True)
	if exitcode:
		print "ERROR: Something went wrong with CAP3!"
		return exitcode

	print "Joining CAP3 contigs and singletons"
	join_cmd = "time parallel cat {{1}}/velvet_contigs.fa.cap.contigs {{1}}/velvet_contigs.fa.cap.singlets '>' {{1}}/{{1}}_cap3ed.fa :::: {}".format(cap3_genefilename)
	print join_cmd
	exitcode = subprocess.call(join_cmd,shell=True)
	if exitcode:
		print "ERROR: Something went wrong joining the CAP3 output files!"
		return exitcode
		
	return None	

def exonerate(genes,basename,run_dir,replace=True,cpu=None,thresh=55,use_velvet=False,depth_multiplier=0,length_pct=100,timeout=None):
	#Check that each gene in genes actually has CAP3 output
	#cap3_sizes = [os.stat(os.path.join(x,x+"_cap3ed.fa")).st_size for x in genes]
	#print cap3_sizes
	
	if replace:
		for g in genes:
			if os.path.isdir(os.path.join(g,basename)):
				shutil.rmtree(os.path.join(g,basename))
	#genes = [x for x in genes if os.stat(os.path.join(x,x+"_cap3ed.fa")).st_size > 0]
	if len(genes) == 0:
		print "ERROR: No genes recovered for {}!".format(basename)
		return 1
	
	if os.path.isfile("genes_with_seqs.txt"):
		os.remove("genes_with_seqs.txt")
	
	if use_velvet:
		file_stem = "cap3ed.fa"
	else:
		file_stem = "contigs.fasta"
	
	print "Running Exonerate to generate sequences for {} genes".format(len(genes))
	
	parallel_cmd_list = ["time parallel","--eta"]
	if cpu:
		parallel_cmd_list.append("-j {}".format(cpu))
	if timeout:
		parallel_cmd_list.append("--timeout {}%".format(timeout))	
	
	exonerate_cmd_list = ["python","{}/exonerate_hits.py".format(run_dir),
				"{}/{}_baits.fasta","{{}}/{{}}_{}".format(file_stem),
				"--prefix {{}}/{}".format(basename), 
				"-t {}".format(thresh),
				"--depth_multiplier {}".format(depth_multiplier),
				"--length_pct {}".format(length_pct),
				"::::",
				exonerate_genefilename,
				"> genes_with_seqs.txt"]
				
	exonerate_cmd = " ".join(parallel_cmd_list) + " " + " ".join(exonerate_cmd_list)			
	
# 	if cpu:
# 		exonerate_cmd = "time parallel -j {} python {} {{}}/{{}}_baits.fasta {{}}/{{}}_{} --prefix {{}}/{} -t {} --depth_multiplier {} --length_pct {} :::: {} > genes_with_seqs.txt".format(cpu,os.path.join(run_dir,"exonerate_hits.py"),file_stem,basename,thresh,depth_multiplier,length_pct,exonerate_genefilename)
# 	else:
# 		exonerate_cmd = "time parallel python {} {{}}/{{}}_baits.fasta {{}}/{{}}_{} --prefix {{}}/{} -t {} --depth_multiplier {} --length_pct {} :::: {} > genes_with_seqs.txt".format(os.path.join(run_dir,"exonerate_hits.py"),file_stem,basename,thresh,depth_multiplier,length_pct,exonerate_genefilename)
	print exonerate_cmd
	exitcode = subprocess.call(exonerate_cmd,shell=True)
	if exitcode:
		print "ERROR: Something went wrong with Exonerate!"
		return exitcode
	return

def bwa(readfiles,baitfile,basename,cpu,unpaired=None):
	"""Conduct BWA search of reads against the baitfile.
	Returns an error if the second line of the baitfile contains characters other than ACTGN"""
	dna = set("ATCGN")
	if os.path.isfile(baitfile):
		#Quick detection of whether baitfile is DNA. 
		with open(baitfile) as bf:
			header = bf.readline()
			seqline = bf.readline().rstrip().upper()
			if set(seqline) - dna:
				print "ERROR: characters other than ACTGN found in first line. You need a nucleotide bait file for BWA!"
				return None
	
		if os.path.isfile(os.path.split(baitfile)[0]+'.amb'):
			db_file = baitfile
		else:
			print "Making nucleotide bwa index in current directory."
                        baitfileDir = os.path.split(baitfile)[0]
                        if baitfileDir:
                                if os.path.realpath(baitfileDir) != os.path.realpath('.'):
				    shutil.copy(baitfile,'.')
			db_file = os.path.split(baitfile)[1]
			make_bwa_index_cmd = "bwa index {}".format(db_file)
			print "[CMD]: {}".format(make_bwa_index_cmd)
			exitcode = subprocess.call(make_bwa_index_cmd,shell=True)
			if exitcode:
				return None
	else:
		print "ERROR: Cannot find baitfile at: {}".format(baitfile)
		return None
	
	if not cpu:
		import multiprocessing
		cpu = multiprocessing.cpu_count()
	
	if len(readfiles) < 3:
		bwa_fastq = " ".join(readfiles)
	else:
		bwa_fastq = readfiles
			
	bwa_commands = ["time bwa mem","-t",str(cpu),db_file,bwa_fastq," | samtools view -h -b -S - > "]
	if unpaired:
		bwa_commands.append(basename+"_unpaired.bam")
	else:
		bwa_commands.append(basename+".bam")
	full_command = " ".join(bwa_commands)
	print "[CMD]: {}".format(full_command)
	exitcode = subprocess.call(full_command,shell=True)
	if exitcode:
		return None
			
	return basename + '.bam'


			
def main():
	parser = argparse.ArgumentParser(description=helptext,formatter_class=argparse.RawTextHelpFormatter)
	parser.add_argument("--check-depend",dest='check_depend',help="Check for dependencies (executables and Python packages) and exit. May not work at all on Windows.",action='store_true')
	parser.add_argument("--bwa",dest="bwa",action='store_true',help="Use BWA to search reads for hits to target. Requires BWA and a bait file that is nucleotides!",default=False)
	parser.add_argument("--no-blast",dest="blast",action="store_false",help="Do not run the blast step. Downstream steps will still depend on the *_all.blastx file. \nUseful for re-runnning assembly/exonerate steps with different options.")
	parser.add_argument("--no-distribute",dest="distribute",action="store_false",help="Do not distribute the reads and bait sequences to sub-directories.")
	parser.add_argument("--no-velvet",dest="velvet",action="store_false",help="Do not run the velvet stages (velveth and velvetg)")
	parser.add_argument("--no-cap3",dest="cap3",action="store_false",help="Do not run CAP3, which joins the output of the different velvet runs")
	parser.add_argument("--no-exonerate",dest="exonerate",action="store_false",help="Do not run the Exonerate step, which assembles full length CDS regions and proteins from each gene")
	parser.add_argument("--velvet-mode",dest="use_velvet",action="store_true",help="Backwards compatability for velvet mode. NOT RECOMMENDED, VELVET MAKES ERROR PRONE ASSEMBLIES!")
	parser.add_argument("--no-assemble",dest="assemble",action="store_false",help="Skip the SPAdes assembly stage.")
	
	parser.add_argument('-r',"--readfiles",nargs='+',help="One or more read files to start the pipeline. If exactly two are specified, will assume it is paired Illumina reads.",default=[])
	parser.add_argument('-b','--baitfile',help="FASTA file containing bait sequences for each gene. If there are multiple baits for a gene, the id must be of the form: >Taxon-geneName",default=None)
	
	parser.add_argument('--cpu',type=int,default=0,help="Limit the number of CPUs. Default is to use all cores available.")
	parser.add_argument('--evalue',type=float,default=1e-10,help="e-value threshold for blastx hits, default: %(default)s")
	parser.add_argument('--max_target_seqs',type=int,default=10,help='Max target seqs to save in blast search, default: %(default)s')
	parser.add_argument('--cov_cutoff',type=int,default=8,help="Coverage cutoff for velvetg. default: %(default)s")
	parser.add_argument('--ins_length',type=int,default=200,help="Insert length for velvetg. default: %(default)s")
	parser.add_argument("--kvals",nargs='+',help="Values of k for velvet assemblies. Velvet needs to be compiled to handle larger k-values! Default auto-dectection by SPAdes.",default=None)
	parser.add_argument("--thresh",type=int,help="Percent Identity Threshold for stitching together exonerate results. Default is 55, but increase this if you are worried about contaminant sequences.",default=65)
	parser.add_argument("--length_pct",help="Include an exonerate hit if it is at least as long as X percentage of the reference protein length. Default = 90%%",default=90,type=int)
	parser.add_argument("--depth_multiplier",help="Accept any full-length exonerate hit if it has a coverage depth X times the next best hit. Set to zero to not use depth. Default = 10",default=10,type=int)

	parser.add_argument('--prefix',help="Directory name for pipeline output, default is to use the FASTQ file name.",default=None)
	parser.add_argument("--timeout",help="Use GNU Parallel to kill long-running processes if they take longer than X percent of average.",default=0)
	
	parser.add_argument("--target",help="Use this target to align sequences for each gene. Other targets for that gene will be used only for read sorting.",default=None)
	parser.add_argument("--unpaired",help="Include a single FASTQ file with unpaired reads along with the two paired read files",default=False)
	
	parser.set_defaults(check_depend=False,blast=True,distribute=True,velvet=False,cap3=False,assemble=True,use_velvet=False,exonerate=True)
	if len(sys.argv) == 1:
		parser.print_help()
		sys.exit(1)
	args = parser.parse_args()
	
	run_dir = os.path.realpath(os.path.split(sys.argv[0])[0])

	#Check dependencies
	if args.check_depend:
		if check_dependencies():
			other_scripts = ["distribute_reads_to_targets.py","distribute_targets.py","exonerate_hits.py"]
			for script in other_scripts:
				if os.path.isfile(os.path.join(run_dir,script)):
					pass
				else:
					print "ERROR: Script {} not found! Please make sure it is in the same directory as this one!".format(script)
					return
			print "Everything looks good!"
			return
		else:
			print "ERROR: One or more dependencies not found!"
			return

	if args.baitfile:
		baitfile = os.path.abspath(args.baitfile)
	else:
		parser.print_help()
		return
	readfiles = [os.path.abspath(x) for x in args.readfiles]	
	if args.unpaired:
		unpaired_readfile = os.path.abspath(args.unpaired)
	else:
		unpaired_readfile = None
	if len(args.readfiles) < 1:
		print "ERROR: Please specify readfiles with -r"
		return
	if not args.baitfile:
		print "ERROR: Please specify a FASTA file containing target sequences."
		return
	
	#Generate directory
	basedir, basename = make_basename(args.readfiles,prefix=args.prefix)
	os.chdir(os.path.join(basedir,basename))
	
	#BWA
	if args.bwa:
		if args.blast:
			args.blast=False
			bamfile = bwa(readfiles,baitfile,basename,cpu=args.cpu)
			if args.unpaired:
				unpaired_bamfile = bwa(unpaired_readfile,baitfile,basename,cpu=args.cpu,unpaired=True)
			if not bamfile:
				print "ERROR: Something went wrong with the BWA step, exiting!"
				return
		else:
			bamfile = basename + ".bam"
	
	#BLAST
	if args.blast:
		if args.unpaired:
			unpaired_blastxfile = blastx(unpaired_readfile,baitfile,args.evalue,basename,cpu=args.cpu,max_target_seqs=args.max_target_seqs,unpaired=True)
		blastx_outputfile = blastx(readfiles,baitfile,args.evalue,basename,cpu=args.cpu,max_target_seqs=args.max_target_seqs)
		if not blastx_outputfile:
			print "ERROR: Something is wrong with the Blastx step, exiting!"
			return
	else:
		blastx_outputfile = basename+".blastx"
	#Distribute
	
	if args.distribute:
		if args.bwa:
			exitcode = distribute_bwa(bamfile,readfiles,baitfile,run_dir,args.target,unpaired_readfile)
		else:
			exitcode=	distribute(blastx_outputfile,readfiles,baitfile,run_dir,args.target,unpaired_readfile)
		if exitcode:
			sys.exit(1)
	if len(readfiles) == 2:
		genes = [x for x in os.listdir(".") if os.path.isfile(os.path.join(x,x+"_interleaved.fasta"))]
	else:
		genes = [x for x in os.listdir(".") if os.path.isfile(os.path.join(x,x+"_unpaired.fasta"))]
	#genes = ["gene008"]
	if len(genes) == 0:
		print "ERROR: No genes with BLAST hits! Exiting!"
		return
	
	if args.use_velvet:
		#Velvet
		if args.velvet:
			exitcode = velvet(genes,cov_cutoff=args.cov_cutoff,ins_length=args.ins_length,kvals=args.kvals,cpu=args.cpu)
			if exitcode:
				return
			
		#CAP3
		if args.cap3:
			exitcode = cap3(genes,cpu=args.cpu)
			if exitcode:
				return

	if args.assemble:
		if len(readfiles) == 1:	
			spades_genelist = spades(genes,run_dir,cov_cutoff=args.cov_cutoff,cpu=args.cpu,kvals=args.kvals,paired=False,timeout=args.timeout)
		elif len(readfiles) == 2:
			if unpaired_readfile:
				spades_genelist = spades(genes,run_dir,cov_cutoff=args.cov_cutoff,cpu=args.cpu,kvals=args.kvals,timeout=args.timeout,unpaired=True)
			else:
				spades_genelist = spades(genes,run_dir,cov_cutoff=args.cov_cutoff,cpu=args.cpu,kvals=args.kvals,timeout=args.timeout)

		else:
			print "ERROR: Please specify either one (unpaired) or two (paired) read files! Exiting!"
			return
		if not spades_genelist:
			print "ERROR: No genes had assembled contigs! Exiting!"
			return
	
	#Exonerate hits
	if args.exonerate:
		genes = [x.rstrip() for x in open(exonerate_genefilename).readlines()]
		exitcode = exonerate(genes,basename,run_dir,cpu=args.cpu,thresh=args.thresh,length_pct = args.length_pct,depth_multiplier=args.depth_multiplier,timeout=args.timeout)
		if exitcode:
			return
	
	sys.stderr.write("Generated sequences from {} genes!\n".format(len(open("genes_with_seqs.txt").readlines())))
	
	paralog_warnings = [x for x in os.listdir(".") if os.path.isfile(os.path.join(x,basename,"paralog_warning.txt"))]
	with open("genes_with_paralog_warnings.txt",'w') as pw:
		pw.write("\n".join(paralog_warnings))
	sys.stderr.write("WARNING: Potential paralogs detected for {} genes!".format(len(paralog_warnings)))

if __name__ == "__main__":main()