reads_first.py

#!/usr/bin/env python

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

helptext="""

HybPiper Version 1.2 (March 2017)

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,exclude=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 exclude:
        target_cmds.append("-- exclude {}".format(exclude))
    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,exclude=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:
        up_bamfile = bamfile.replace(".bam","_unpaired.bam")
        unpaired_cmd = "time python {} {} {}".format(os.path.join(run_dir,"distribute_reads_to_targets_bwa.py"),up_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))
    if unpaired:
        target_cmds.append("--unpaired")
    if exclude:
        target_cmds.append("--exclude {}".format(exclude))
    target_cmd = " ".join(target_cmds)
    print("[DISTRIBUTE]: {}".format(target_cmd))
    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))
    if kvals:
        spades_runner_list.append("--kvals")
        spades_runner_list.append("{}".format(",".join(kvals)))
    
    
    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. Can be a tab-delimited file (one gene per line) or a single sequence name",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.add_argument("--exclude",help="Do not use any sequence with the specified string as a target sequence for exonerate. The sequence will be used for read sorting.",default=None)
    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])
    print(("HybPiper was called with these arguments:\n{}\n".format(" ".join(sys.argv))))
    
    #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:
        pre_existing_fastas = glob.glob("./*/*_interleaved.fasta") + glob.glob("./*/*_unpaired.fasta")
        for fn in pre_existing_fastas:
            os.remove(fn)
        if args.bwa:
            exitcode = distribute_bwa(bamfile,readfiles,baitfile,run_dir,args.target,unpaired_readfile,args.exclude)
        else:
            exitcode=    distribute(blastx_outputfile,readfiles,baitfile,run_dir,args.target,unpaired_readfile,args.exclude)
        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"]
    print("READFILES:\n")
    print(readfiles)
    print(genes)
    print(os.listdir('.'))
    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()