<repo_name>jasonsahl/wgfast<file_sep>/README.md
#### *WG-FAST*
The whole genome focused array SNP typing (*WG-FAST*) pipeline
written by: <NAME>
Email: jasonsahl at gmail dot com

for a more comprehensive overview, look at the Manual

#### overview
The goal of *WG-FAST* is to phylogenetically genotype an unknown
sample in the context of a well studied pathogen.  This sample
can either be a metagenomics dataset, a metatranscriptomics dataset,
or a single isolate sequencing dataset

#### Installation
1.  The easiest way to install is through conda:  

```conda create -n wgfast python=3.6```  
```source activate wgfast```  
```conda install -c bioconda gatk4=4.2.6.1 picard=2.27.4 raxml samtools bbmap dendropy minimap2 biopython```  
#You might need to install Bioconda with: pip install Biopython  

2. Download the wgfast github repository, install:  
```git clone https://github.com/jasonsahl/wgfast.git```  
```python setup.py build```  
```python setup.py install```

3. Open the script (wgfast.py) with a text editor and change the path to your *WG-FAST* installation directory.  
For example:  
WGFAST_PATH="/Users/jsahl/wgfast"  

4. To verify your installation, enter the wgfast directory and type the command below.  If everything
is working correctly, all tests should pass:  
```python tests/test_all_functions.py```  
<file_sep>/wg_fast/util.py
from __future__ import division
import os,os.path
import re
import logging
import subprocess
import sys
from subprocess import Popen
import time
try:
    from Bio import SeqIO
    from Bio import Phylo
except:
    print("BioPython is not in your PYTHONPATH, but needs to be")
    sys.exit()
try:
    import dendropy
    from dendropy import treecalc
except:
    print("dendropy is not installed, but needs to be")
    sys.exit()
import glob
from operator import itemgetter
import threading
import collections
import random

#logPrint stuff
OUTSTREAM = sys.stdout
ERRSTREAM = sys.stderr
DEBUG = False

def logPrint(msg, stream=None):
    if stream is None:
        stream = OUTSTREAM
    stream.write('LOG: %s - %s\n' % (timestamp(), removeRecursiveMsg(msg)))
    stream.flush()

def errorPrint(msg, stream=None):
    if stream is None:
        stream = ERRSTREAM

    stream.write('ERROR: %s - %s\n' % (timestamp(), removeRecursiveMsg(msg)))
    stream.flush()

def debugPrint(fmsg, stream=None):
    """In this case msg is a function, so the work is only done if debugging is one"""
    if DEBUG:
        if stream is None:
            stream = ERRSTREAM

        stream.write('DEBUG: %s - %s\n' % (timestamp(), removeRecursiveMsg(fmsg())))
        stream.flush()

def timestamp():
    return time.strftime('%Y/%m/%d %H:%M:%S')

def removeRecursiveMsg(msg):
    """
    This takes a message and if it starts with something that looks like
    a message generated with these tools it chops it off.  Useful if using
    one of these logging functions to print output from a program using
    the same logging functions
    """
    if msg.startswith('ERROR: ') or msg.startswith('DEBUG: ') or msg.startswith('LOG: '):
        return msg.split(' - ', 1)[1]
    else:
        return msg

def mp_shell(func, params, numProc):
    from multiprocessing import Pool
    p = Pool(numProc)
    out = p.map(func, params)
    p.terminate()
    return out

def report_stats(results,name,output):
    outfile = open(output, "w")
    total_size = []
    mapped_size = []
    with open(results) as infile:
        for line in infile:
            newline = line.strip()
            fields = newline.split()
            """this just makes sure that the file looks correct"""
            if len(fields)==3:
                total_size.append(float(fields[1]))
                mapped_size.append(float(fields[2]))
            else:
                print("coverage file for %s is malformed" % name)
                print("-------------------------")
    try:
        total_summed = sum(total_size)
        total_mapped = sum(mapped_size)
        mapped_value = (total_mapped/total_summed)*100
        outfile.write(str(name)+"\t"+str.format('{0:.4f}',mapped_value)+"\n")
        outfile.close()
    except:
        pass

def merge_files_by_column(column,file_1,file_2,out_file):
    """Takes 2 file and merge their columns based on the column. It is assumed
    that line ordering in the files do not match, so we read both files into memory
    and join them"""
    join_map = {}
    with open(file_1) as my_file_1:
        for line in my_file_1:
            line.strip()
            row = line.split()
            column_value = row.pop(column)
            join_map[column_value] = row
    with open(file_2) as my_file_2:
        for line in my_file_2:
            line.strip()
            row = line.split()
            column_value = row.pop(column)
            if column_value in join_map:
                join_map[column_value].extend(row)
    fout = open(out_file, 'w')
    for k, v in join_map.items():
        fout.write('\t'.join([k] + v) + '\n')
    fout.close()

def sum_coverage(coverage,cov,name):
    outfile = open("%s.amount_covered.txt" % name, "w")
    another_outfile = open("%s.sum_covered.txt" % name, "w")
    all = []
    my_dict = {}
    cov_dict = {}
    with open(coverage) as my_file:
        for line in my_file:
            fields=line.split()
            fields = map(lambda s: s.strip(), fields)
            all.append(fields)
    for x, y in all:
        """Here we're only counting it if it is above the given coverage threshold"""
        try:
            cov_dict[x].append(int(y))
        except KeyError:
            cov_dict[x] = [int(y)]
        if int(y)>=int(cov):
           try:
               my_dict[x].append(y)
           except KeyError:
               my_dict[x] = [y]
        else:
               pass
    for k,v in my_dict.items():
        outfile.write(str(k)+"\t"+str(len(v))+"\n")
    for k,v in cov_dict.items():
        another_outfile.write(str(k)+"\t"+str(sum(v))+"\n")
    outfile.close()
    another_outfile.close()

def test_file(option, opt_str, value, parser):
    try:
        with open(value): setattr(parser.values, option.dest, value)
    except IOError:
        print('%s file cannot be opened' % option)
        sys.exit()

def test_dir(option, opt_str, value, parser):
    if os.path.exists(value):
        setattr(parser.values, option.dest, value)
    else:
        print("directory of fastqs cannot be found")
        sys.exit()

def test_filter(option, opt_str, value, parser):
    if "F" in value:
        setattr(parser.values, option.dest, value)
    elif "T" in value:
        setattr(parser.values, option.dest, value)
    else:
        print("option not supported.  Only select from T and F")
        sys.exit()

def test_methods(option, opt_str, value, parser):
    if "ML" in value:
        setattr(parser.values, option.dest, value)
    elif "MP" in value:
        setattr(parser.values, option.dest, value)
    else:
        print("option not supported.  Only select from MP or ML")
        sys.exit()

def test_models(option, opt_str, value, parser):
    if "GTRGAMMA" in value:
        setattr(parser.values, option.dest, value)
    elif "ASC_GTRGAMMA" in value:
        setattr(parser.values, option.dest, value)
    else:
        print("substitution model is not supported")
        sys.exit()

def get_seq_length(ref, name):
    """uses BioPython in order to calculated the length of
    each fasta entry in the reference fasta"""
    outfile = open("%s.tmp.txt" % name, "w")
    for record in SeqIO.parse(open(ref), "fasta"):
        outfile.write(str(record.id)+"\t"+str(len(record.seq))+"\n")
    outfile.close()

def remove_column(temp_file, name):
    outfile = open("%s.coverage.out" % name, "w")
    my_fields = []
    with open(temp_file) as my_file:
        for line in my_file:
            fields=line.split()
            del fields[1]
            my_fields.append(fields)
    for x in my_fields:
        outfile.write("\t".join(x))
        outfile.write("\n")
    outfile.close()

def get_seq_name(in_fasta):
    """used for renaming the sequences - tested"""
    return os.path.basename(in_fasta)

def get_readFile_components(full_file_path):
    """function adapted from:
    https://github.com/katholt/srst2 - tested"""
    (file_path,file_name) = os.path.split(full_file_path)
    m1 = re.match("(.*).gz",file_name)
    ext = ""
    if m1 is not None:
        ext = ".gz"
        file_name = m1.groups()[0]
    (file_name_before_ext,ext2) = os.path.splitext(file_name)
    full_ext = ext2+ext
    return file_path,file_name_before_ext,full_ext

def read_file_sets(dir_path):
    """match up pairs of sequence data, adapted from
    https://github.com/katholt/srst2 will be tough to test
    with variable names and read paths"""
    fileSets = {}
    forward_reads = {}
    reverse_reads = {}
    num_paired_readsets = 0
    num_single_readsets = 0
    for infile in glob.glob(os.path.join(dir_path, "*.fastq.gz")):
        (file_path,file_name_before_ext,full_ext) = get_readFile_components(infile)
        m=re.match("(.*)(_S.*)(_L.*)(_R.*)(_.*)", file_name_before_ext)
        if m is None:
            #m=re.match("(.*)("+"_R1"+")(_.*)$",file_name_before_ext)
            m=re.match("(.*)("+"_R1"+")(.*)$",file_name_before_ext)
            if m is not None:
                (baseName,read) = m.groups()[0], m.groups()[1]
                forward_reads[baseName] = infile
            else:
                #m=re.match("(.*)("+"_R2"+")(_.*)$",file_name_before_ext)
                m=re.match("(.*)("+"_R2"+")(.*)$",file_name_before_ext)
                if m is not None:
                    (baseName,read) = m.groups()[0], m.groups()[1]
                    reverse_reads[baseName] = infile
                else:
                    print("Could not determine forward/reverse read status for input file %s" % infile)
        else:
            baseName, read  = m.groups()[0], m.groups()[3]
            if read == "_R1":
                forward_reads[baseName] = infile
            elif read == "_R2":
                reverse_reads[baseName] = infile
            else:
                print("Could not determine forward/reverse read status for input file")
                fileSets[file_name_before_ext] = infile
                num_single_readsets += 1
    for sample in forward_reads:
        if sample in reverse_reads:
            fileSets[sample] = [forward_reads[sample],reverse_reads[sample]] # store pair
            num_paired_readsets += 1
        else:
            fileSets[sample] = [forward_reads[sample]] # no reverse found
            num_single_readsets += 1
            logging.info('Warning, could not find pair for read:' + forward_reads[sample])
    for sample in reverse_reads:
        if sample not in fileSets:
            fileSets[sample] = reverse_reads[sample] # no forward found
            num_single_readsets += 1
            logging.info('Warning, could not find pair for read:' + reverse_reads[sample])
    if num_paired_readsets > 0:
        logging.info('Total paired readsets found:' + str(num_paired_readsets))
    if num_single_readsets > 0:
        logging.info('Total single reads found:' + str(num_single_readsets))
    return fileSets

def get_sequence_length(fastq_in):
    from itertools import islice
    from gzip import GzipFile
    with GzipFile("%s" % fastq_in) as file:
        head = list(islice(file, 2))
    return len(head[1])

def _perform_workflow_run_loop_dev(data):
    """sample ID"""
    idx = data[0]
    """path to data"""
    f = data[1]
    dir_path = data[2]
    reference = data[3]
    ref_coords = data[4]
    coverage = data[5]
    proportion = data[6]
    matrix = data[7]
    scratch_dir = data[8]
    doc = data[9]
    tmp_dir = data[10]
    """This is now the PATH to bbduk.sh"""
    wgfast_path = data[11]
    processors = data[12]
    ploidy = data[13]
    if os.path.isfile("%s.tmp.xyx.matrix" % idx):
        pass
    else:
        """This means that the data is paired end"""
        if len(f)>1:
            if os.path.isfile("%s.F.paired.fastq.gz" % idx):
                pass
            else:
                length = int(get_sequence_length(f[0])/2)
                try:
                    subprocess.check_call("bbduk.sh in=%s in2=%s ref=%s/bin/illumina_adapters_all.fasta out=%s.F.paired.fastq.gz out2=%s.R.paired.fastq.gz minlen=%s overwrite=true > /dev/null 2>&1" % (f[0],f[1],wgfast_path,idx,idx,length), shell=True)
                except:
                    print("Read trimmer did not finish correctly")
                    sys.exit()
            subprocess.check_call("minimap2 -ax sr %s/reference.fasta %s.F.paired.fastq.gz %s.R.paired.fastq.gz | samtools sort -l 0 -@ %s - | samtools view -F 4 -Su -o %s_renamed.bam -" % (scratch_dir,idx,idx,processors,idx),stdout=open(os.devnull, 'wb'),stderr=open(os.devnull, 'wb'),shell=True)
        else:
            """Single end read support"""
            length = int(get_sequence_length(f[0])/2)
            try:
                subprocess.check_call("bbduk.sh -Xmx2g in=%s ref=%s/bin/illumina_adapters_all.fasta out=%s.single.fastq.gz minlen=%s overwrite=true ignorebadquality" % (f[0],wgfast_path,idx,length), stdout=open(os.devnull, 'wb'),stderr=open(os.devnull, 'wb'),shell=True)
            except:
                print("Read trimmer did not finish correctly")
                sys.exit()
            subprocess.check_call("minimap2 -ax sr %s/reference.fasta %s.single.fastq.gz | samtools sort -l 0 -@ %s - | samtools view -F 4 -Su -o %s_renamed.bam -" % (scratch_dir,idx,processors,idx),stdout=open(os.devnull, 'wb'),stderr=open(os.devnull, 'wb'),shell=True)
        """inserts read group information, required by new versions of GATK"""
        try:
            subprocess.check_call("picard AddOrReplaceReadGroups I=%s_renamed.bam O=%s_renamed_header.bam SORT_ORDER=coordinate RGID=%s RGLB=%s RGPL=illumina RGSM=%s RGPU=name CREATE_INDEX=true VALIDATION_STRINGENCY=SILENT" % (idx,idx,idx,idx,idx),stdout=open(os.devnull, 'wb'),stderr=open(os.devnull, 'wb'),shell=True)
        except:
            time.sleep(5)
            subprocess.check_call("picard AddOrReplaceReadGroups I=%s_renamed.bam O=%s_renamed_header.bam SORT_ORDER=coordinate RGID=%s RGLB=%s RGPL=illumina RGSM=%s RGPU=name CREATE_INDEX=true VALIDATION_STRINGENCY=SILENT" % (idx,idx,idx,idx,idx),stdout=open(os.devnull, 'wb'),stderr=open(os.devnull, 'wb'),shell=True)
        try:
            subprocess.check_call("samtools index %s_renamed_header.bam > /dev/null 2>&1" % idx,stdout=open(os.devnull, 'wb'),stderr=open(os.devnull, 'wb'),shell=True)
        except:
            print("problem indexing file with samtools..exiting")
            sys.exit()
        #TODO: be able to change the ploidy, which will affect the number of calls
        subprocess.check_call("gatk HaplotypeCaller -R %s/reference.fasta -I %s_renamed_header.bam -O %s.test.vcf -ERC BP_RESOLUTION -ploidy %s" % (scratch_dir,idx,idx,ploidy), stdout=open(os.devnull, 'wb'),stderr=open(os.devnull, 'wb'),shell=True)
        subprocess.check_call("gatk GenotypeGVCFs -O %s.vcf.out -R %s/reference.fasta -V %s.test.vcf --include-non-variant-sites" % (idx,scratch_dir,idx), stdout=open(os.devnull, 'wb'),stderr=open(os.devnull, 'wb'),shell=True)
        #subprocess.check_call("bcftools filter -G 5 -g 5 %s.test.vcf2 > %s.vcf.out" % (idx,idx),stdout=open(os.devnull, 'wb'),stderr=open(os.devnull, 'wb'),shell=True)
        if "T" == doc:
            subprocess.check_call("samtools depth -aa %s_renamed_header.bam > %s.coverage" % (idx,idx), shell=True)
            remove_column("%s.coverage" % idx, idx)
            sum_coverage("%s.coverage.out" % idx, coverage, idx)
            merge_files_by_column(0,"ref.genome_size.txt", "%s.amount_covered.txt" % idx, "%s.results.txt" % idx)
            merge_files_by_column(0,"ref.genome_size.txt", "%s.sum_covered.txt" % idx, "%s.cov.results.txt" % idx)
            report_stats("%s.results.txt" % idx, idx, "%s_breadth.txt" % idx)
            report_stats("%s.cov.results.txt" % idx, idx, "%s_sum_cov.txt" % idx)
        else:
            pass
        #filtered.vcf will be created in this function
        good_calls = process_vcf("%s.vcf.out" % idx, ref_coords, coverage, proportion, idx)
        if int(good_calls) > 0:
            make_temp_matrix("%s.filtered.vcf" % idx, matrix, idx)
        else:
            print("sample %s had no SNP calls and will not be inserted into the tree" % idx)
            print("-------------------------")

def run_loop_dev(fileSets,dir_path,reference,processors,ref_coords,coverage,proportion,
    matrix,scratch_dir,doc,tmp_dir,wgfast_path,ploidy):
    files_and_temp_names = []
    for idx, f in fileSets.items():
        files_and_temp_names.append([idx,f,dir_path,reference,ref_coords,coverage,proportion,
                                        matrix,scratch_dir,doc,tmp_dir,wgfast_path,processors,ploidy])
    mp_shell(_perform_workflow_run_loop_dev,files_and_temp_names,processors)

def process_vcf(vcf, ref_coords, coverage, proportion, name):
    """finds SNPs that pass user-defined thresholds
    for coverage and proportion - needs to look at tests"""
    vcf_out = open("%s.filtered.vcf" % name, "w")
    good_snps = []
    mixed_snps = []
    mixed_refs = []
    ref_set = set(ref_coords)
    with open(vcf) as vcf_in:
        for line in vcf_in:
            newline=line.strip("\n")
            if newline.startswith('#'):
               pass
            elif newline.startswith('Java'):
                pass
            elif newline.startswith('/tmp'):
                pass
            elif newline.startswith('Try'):
                pass
            elif newline.startswith('INFO'):
                pass
            elif newline in ['\n', '\r\n']:
                pass
            else:
                fields=newline.split()
                """for GATK, a period signifies a reference call.
                First we want to look at the situation where this is
                not the case"""
                merged_fields=fields[0]+"::"+fields[1]
                if merged_fields in ref_set:
                    #This indicates that the position is a SNP
                    if "." != fields[4] and len(fields[4]) == 1 and len(fields[3])==1 and fields[8]!="GT:AD:PGT:PID:PS":
                        if fields[6] == "LowQual":
                            pass
                        else:
                            snp_fields=fields[9].split(':')
                            if int(len(snp_fields))>2:
                                prop_fields=snp_fields[1].split(',')
                                fixed_coverage = int(float(snp_fields[2]))
                                if fixed_coverage>=coverage:
                                    if int(prop_fields[1])/int(snp_fields[2])>=float(proportion):
                                        vcf_out.write(fields[0]+"::"+fields[1]+"\t"+fields[4]+"\n")
                                        good_snps.append("1")
                                    else:
                                        #Changed out a gap character with an N
                                        vcf_out.write(fields[0]+"::"+fields[1]+"\t"+"N"+"\n")
                                        mixed_snps.append("1")
                                    """if problems are encountered, throw in a gap.  Could be too conservative"""
                                else:
                                    vcf_out.write(fields[0]+"::"+fields[1]+"\t"+"N"+"\n")
                            else:
                                pass
                    #This is if the position is reference
                    elif "." in fields[4] and len(fields[4])==1 and len(fields[3])==1:
                        if fields[6] == "LowQual":
                            pass
                        else:
                            if len(fields) == 10:
                                if "DP" in fields[7]:
                                    nosnp_fields=fields[7].split(';')
                                    #This will provide the coverage
                                    if "DP" in nosnp_fields[0]:
                                        cov_fields=nosnp_fields[0].replace("DP=","")
                                    elif "DP" in nosnp_fields[1]:
                                        cov_fields=nosnp_fields[1].replace("DP=","")
                                    fixed_coverage = int(float(cov_fields))
                                    if fixed_coverage>=coverage:
                                        vcf_out.write(fields[0]+"::"+fields[1]+"\t"+fields[3]+"\n")
                                    else:
                                        vcf_out.write(fields[0]+"::"+fields[1]+"\t"+"N"+"\n")
                                        mixed_refs.append("1")
                                else:
                                    vcf_out.write(fields[0]+"::"+fields[1]+"\t"+"N"+"\n")
                                    mixed_refs.append("1")
                    #If it can't determine the status of the position, add an N
                    else:
                        vcf_out.write(fields[0]+"::"+fields[1]+"\t"+"N"+"\n")
    vcf_out.close()
    print("number of SNPs in genome %s = %s" % (name, str(len(good_snps))))
    print("number of discarded SNPs in genome %s = %s" % (name, str(len(mixed_snps))))
    print("number of discarded Reference positions in genome due to no coverage %s = %s" % (name, str(len(mixed_refs))))
    print("-------------------------")
    return len(good_snps)

def sort_information(x):
    """simple sort - tested"""
    try:
        fields = x.split("::")
        return int(fields[1])
    except:
        raise TypeError("problem encountered parsing fields")

def matrix_to_fasta(matrix_in, outfile):
    """function to convert a SNP matrix to a multi-fasta file - tested"""
    reduced = []
    out_fasta = open(outfile, "w")
    redux = []
    with open(matrix_in) as my_file:
        for line in my_file:
            newline = line.strip()
            fields = newline.split()
            reduced.append(fields[1:])
    test=map(list, zip(*reduced))
    for x in test:
        out_fasta.write(">"+str(x[0])+"\n")
        out_fasta.write("".join(x[1:]))
        out_fasta.write("\n")
        redux.append(">"+str(x[0])+"".join(x[1:3]))
    out_fasta.close()
    return redux

def run_raxml(fasta_in, tree, out_class_file, insertion_method, parameters, model, suffix):
    """untested function, system calls"""
    if "NULL" == parameters:
        if "ASC_GTRGAMMA" == model:
            args = ['raxmlHPC-SSE3', '-f', '%s' % insertion_method,
	         '-s', '%s' % fasta_in, '-m', '%s' % model, '-n', '%s' % suffix, '-t',
	         '%s' % tree, '--asc-corr=lewis', '--no-bfgs', '>', '/dev/null 2>&1']
        else:
            args = ['raxmlHPC-SSE3', '-f', '%s' % insertion_method,
	         '-s', '%s' % fasta_in, '-m', '%s' % model, '-n', '%s' % suffix, '-t',
	         '%s' % tree, '--no-bfgs', '>', '/dev/null 2>&1']
    else:
        if "ASC_GTRGAMMA" == model:
            args = ['raxmlHPC-SSE3', '-f', '%s' % insertion_method,
	         '-s', '%s' % fasta_in, '-m', '%s' % model, '-n', '%s' % suffix, '-R', parameters, '-t',
	         '%s' % tree, '--asc-corr=lewis', '--no-bfgs', '>', '/dev/null 2>&1']
        else:
            args = ['raxmlHPC-SSE3', '-f', '%s' % insertion_method,
	         '-s', '%s' % fasta_in, '-m', '%s' % model, '-n', '%s' % suffix, '-R', parameters, '-t',
	         '%s' % tree, '--no-bfgs', '>', '/dev/null 2>&1']
    try:
        vcf_fh = open('%s.raxml.out' % suffix, 'w')
    except:
        logPrint('could not open raxml file')
    try:
        log_fh = open('%s.raxml.log' % suffix, 'w')
    except:
        logPrint('could not open log file')
    #TODO: See why this doesn't catch many RAxML errors
    try:
        raxml_run = Popen(args, stderr=log_fh, stdout=vcf_fh)
        raxml_run.wait()
    except:
        logPrint("sequence(s) were not inserted into tree!!!!!")
    try:
        os.system("sed 's/\[[^]]*\]//g' RAxML_labelledTree.%s > %s.tree_including_unknowns_noedges.tree" % (suffix, suffix))
        subprocess.check_call("mv RAxML_labelledTree.%s %s_tree_including_unknowns_edges.tree" % (suffix, suffix) , shell=True)
    except:
        logPrint("error encountered running RAxML...check log file")
        sys.exit()
    try:
        subprocess.check_call("cat RAxML_classificationLikelihoodWeights.%s >> %s" % (suffix, out_class_file), shell=True)
    except:
        pass
    os.system("rm RAxML_*.%s" % suffix)
    return suffix

def subsample_snps(matrix, dist_sets, used_snps, subnums):
    """get a list of all possible positions, depending
    on those positions in the original matrix - tested"""
    allSNPs = [ ]
    with open(matrix) as my_file:
        for line in my_file:
            if line.startswith("LocusID"):
                pass
            else:
                fields=line.split()
                allSNPs.append(fields[0])
    for k,v in used_snps.items():
        for z in dist_sets:
            if len(z) == 0:
                pass
            else:
                if z[0]==k:
                    for x in range(1,int(subnums)+1):
                        kept_snps=random.sample(set(allSNPs), int(v))
                        outfile = open("%s.%s.%s.tmp.matrix" % (k,x,z[1]), "w")
                        in_matrix=open(matrix)
                        firstLine = in_matrix.readline()
                        outfile.write(firstLine)
                        outfile.write("\n")
                        first_fields = firstLine.split()
                        fixed_fields = []
                        for y in first_fields:
                            fixed_fields.append(re.sub('[:,]', '', y))
                        gindex=fixed_fields.index(z[1])
                        for line in in_matrix:
                            matrix_fields=line.split()
                            if matrix_fields[0] in kept_snps:
                                outfile.write(line)
                                outfile.write("\n")
                            else:
                                outfile.write("\t".join(matrix_fields[:gindex])+"\t"+"-"+"\t"+"\t".join(matrix_fields[gindex+1:])+"\n")
                        in_matrix.close()
                        outfile.close()
    return allSNPs

def find_used_snps():
    """report how many SNPs were used in a given sample.  This is
    then used for the sub-sampling routine - tested"""
    curr_dir= os.getcwd()
    used_SNPs = {}
    for infile in glob.glob(os.path.join(curr_dir, "*.filtered.vcf")):
        name=get_seq_name(infile)
        reduced=name.replace('.filtered.vcf', '')
        good_snps=[]
        with open(infile) as my_file:
            for line in my_file:
                fields=line.split()
                try:
                    if fields[1] != "N":
                        good_snps.append("1")
                    else:
                        pass
                except:
                    raise TypeError("abnormal number of fields observed")
        used_SNPs[str(reduced)] = int(len(good_snps))
    return used_SNPs

def branch_lengths_2_decimals(str_newick_tree):
    """replaces branch lengths in scientific notation with decimals = tested"""
    colon_s = 0
    comma_back_paren_s = 0
    num = ''
    new_tree = ''
    for count, char in enumerate(str_newick_tree):
        if char == ':':
            colon_s = count
            continue
        if char in (')', ','):
            comma_back_paren_s = 1
            num = '%f' % float(num)
            new_tree += ":" + num
            colon_s = 0
            num = ''
        if colon_s != 0:
            num = num + char
        if colon_s == 0:
            new_tree += char
    new_tree = new_tree.strip('\'').strip('\"').strip('\'') + ";"
    return new_tree

def fasta_to_tab(fasta):
    """tested"""
    outfile = open("out.tab", "w")
    with open(fasta) as my_fasta:
        for record in SeqIO.parse(my_fasta,"fasta"):
            """this list is just for testing,
            and is ok if it's overwritten for each
            fasta"""
            for_test = []
            outfile.write(str(record.id)+"\t"+str(record.seq)+"\n")
            for_test.append(record.id)
            for_test.append(str(record.seq))
    outfile.close()
    return for_test

def tab_to_fasta(new_tab):
    """tested"""
    outfile = open("out.fasta", "w")
    with open(new_tab) as infile:
        for line in infile:
            to_test = []
            fields = line.split()
            outfile.write(">"+fields[0]+"\n")
            outfile.write(fields[1].upper())
            to_test.append(fields[0])
            to_test.append(fields[1].upper())
    outfile.close()
    return to_test

def tab_to_matrix(tab):
    """tested"""
    reduced = []
    out_matrix = open("tab_matrix", "w")
    with open(tab) as my_file:
        for line in my_file:
            tmp_list = []
            fields = line.split()
            tmp_list.append(fields[0])
            for nucs in fields[1]:
                tmp_list.append(nucs.upper())
            reduced.append(tmp_list)
    test=list(map(list,zip(*reduced)))
    to_return = []
    for x in test:
        out_matrix.write("\t".join(x))
        out_matrix.write("\n")
        to_return.append(list(x))
    out_matrix.close()
    return to_return

def filter_alignment(tab):
    """tested"""
    outfile = open("tab.filtered", "w")
    with open(tab) as infile:
        firstLine = infile.readline()
        test_fields = []
        outfile.write(firstLine)
        outfile.write("\n")
        for line in infile:
            valid_fields = []
            fields = line.split()
            for field in fields:
                """skip fields that might be present when missing data are included"""
                if field != "-" and field != "N" and field != "X":
                    valid_fields.append(field)
                else:
                    pass
            counter=collections.Counter(valid_fields)
            values=counter.values()
            #sorted_values = values.sort(key=int)
            sorted_values = sorted(values)
            if len(sorted_values)>int(1):
                outfile.write(line)
                outfile.write("\n")
            else:
                pass
            test_fields.append(valid_fields)
    outfile.close()
    return test_fields

def raxml_calculate_base_tree(in_fasta, model, name):
    """not tested, all system calls"""
    args = ['raxmlHPC-SSE3', '-f', 'd', '-p', '12345',
	     '-s', '%s' % in_fasta, '-m', '%s' % model, '-n', '%s' % name, "--no-bfgs",
	     '>', '/dev/null 2>&1']
    try:
        vcf_fh = open('raxml.out', 'w')
    except:
        logPrint('could not open raxml file')
    try:
        log_fh = open('raxml.log', 'w')
    except:
        logPrint('could not open log file')
    try:
        raxml_run = Popen(args, stderr=log_fh, stdout=vcf_fh)
        raxml_run.wait()
    except:
        print("could not infer base pruned tree")
        sys.exit()

def file_to_fasta(matrix, out_fasta):
    """almost identical to matrix_to_fasta. Not tested"""
    reduced = [ ]
    out_matrix = open(out_fasta, "w")
    with open(matrix) as my_file:
        for line in my_file:
            fields = line.strip().split()
            reduced.append(fields)
    test=map(list, zip(*reduced))
    for x in test:
        out_matrix.write(">"+str(x[0])+"\n")
        out_matrix.write("".join(x[1:]))
        out_matrix.write("\n")
    out_matrix.close()

def prune_fasta(to_prune, infile, outfile):
    """tested"""
    my_out = open(outfile, "w")
    seqrecords = [ ]
    ids = [ ]
    with open(infile) as my_fasta:
        for record in SeqIO.parse(my_fasta, "fasta"):
            if record.id not in to_prune:
                seqrecords.append(record)
                ids.append(record.id)
    SeqIO.write(seqrecords, my_out, "fasta")
    my_out.close()
    return ids

def remove_invariant_sites(in_fasta, out_fasta):
    """only keep invarint sites, all functions are tested"""
    fasta_to_tab(in_fasta)
    tab_to_matrix("out.tab")
    filter_alignment("tab_matrix")
    file_to_fasta("tab.filtered", out_fasta)

def compare_subsample_results(outnames,distances,fudge):
    """needs testing"""
    curr_dir= os.getcwd()
    for infile in glob.glob(os.path.join(curr_dir, "*.subsample.distances.txt")):
        name=get_seq_name(infile)
        split_fields=name.split(".")
        genomes_used = []
        all_dists=[]
        dists_greater_than_true=[]
        dists_equal_to_true=[]
        dists_less_than_true=[]
        """Here, I get the list of the genomes that are being analyzed"""
        try:
            with open(infile) as my_file:
                for line in my_file:
                    fields = line.split()
                    all_dists.append(float(fields[7]))
        except:
           print("problem parsing input file: %s" % infile)
        if len(all_dists)>=1:
            max_dist=max(all_dists)
            print("")
            print("maximum subsample distance between %s and %s = %.2f" % (fields[3],fields[5],float(max_dist)))
        else:
            print("problem grabbing distances- make sure that subsample.distances.txt files aren't empty")
        true_dists = [ ]
        for distance in distances:
            if distance[1] == split_fields[1]:
                true_dists.append(distance[2])
        for all_dist in all_dists:
            if "%.2f" % float(all_dist)> "%.2f" % (float(true_dists[0])+float(fudge)):
                dists_greater_than_true.append("1")
            elif "%.2f" % float(all_dist)<"%.2f" % (float(true_dists[0])-float(fudge)):
                dists_less_than_true.append("1")
            else:
                dists_equal_to_true.append("1")
        try:
            greaters = int(len(dists_greater_than_true))
            equals = int(len(dists_equal_to_true))
            lessers = int(len(dists_less_than_true))
            print("True distance between Reference and %s = %.2f" % (split_fields[1],float(true_dists[0])))
            print("Sample: %s" % split_fields[0])
            print("Subsample distances between Reference and %s greater than true value = %s" % (split_fields[2],greaters))
            print("Subsample distances between Reference and %s equal to true value = %s" % (split_fields[2],equals))
            print("Subsample distances between Reference and %s less than true value = %s" % (split_fields[2],lessers))
            p = (greaters+lessers)/(greaters+lessers+equals)
            print("Placement p value = %.3f" % float(p))
            print("-------------------------")
        except:
            pass

def transform_tree(tree):
    """converts a Newick tree into a Nexus-formatted
    tree that can be visualized with FigTree - needs testing"""
    #infile = open(tree, "U")
    tree_string = []
    with open(tree) as infile:
        for line in infile:
            tree_string.append(line)
    #infile.close()
    mytree = Phylo.read(tree, 'newick')
    tree_names = [ ]
    for clade in mytree.find_clades():
        if clade.name:
            tree_names.append(clade.name)
    outfile = open("transformed.tree", "w")
    outfile.write("#NEXUS"+"\n")
    outfile.write("begin taxa;"+"\n")
    outfile.write("\t"+"dimensions ntax="+str(len(tree_names))+";"+"\n")
    outfile.write("\t"+"taxlabels"+"\n")
    for clade in mytree.find_clades():
        if clade.name:
            tree_names.append(clade.name)
    nr=[x for i, x in enumerate(tree_names) if x not in tree_names[i+1:]]
    for tree_name in nr:
        if "QUERY__" in tree_name:
            outfile.write("\t"+"'%s'" % tree_name+"[&!color=#-3407821]"+"\n")
        else:
            outfile.write("\t"+tree_name+"\n")
    outfile.write(";"+"\n")
    outfile.write("end;"+"\n")
    outfile.write(""+"\n")
    outfile.write("begin trees;"+"\n")
    for x in tree_string:
        outfile.write("\t"+"tree tree_1 = [&R] "+str(x)+"\n")
    outfile.write("end;"+"\n")
    infile.close()
    outfile.close()

def write_reduced_matrix(matrix):
    """This function takes a NASP formatted
    SNP matrix and writes a temporary matrix
    that can be easily combined with temporary files - tested"""
    in_matrix = open(matrix)
    outfile = open("temp.matrix", "w")
    outdata = [ ]
    firstLine = in_matrix.readline()
    first_fields=firstLine.split()
    last=first_fields.index("#SNPcall")
    outfile.write("\t".join(first_fields[:last]))
    outfile.write("\n")
    for line in in_matrix:
        fields = line.split()
        outfile.write("\t".join(fields[:last]))
        outfile.write("\n")
        outdata.append(len(fields[:last]))
    outfile.close()
    in_matrix.close()
    return outdata

def make_temp_matrix(vcf, matrix, name):
    """these are all of the screened SNPs - tested"""
    matrix_ids=[]
    with open(matrix) as in_matrix:
        firstLine = in_matrix.readline()
        for line in in_matrix:
            newline = line.strip("\n")
            mfields=newline.split()
            matrix_ids.append(mfields[0])
    value_dict={}
    new_dicts={}
    with open(vcf) as my_file:
        for line in my_file:
            newline = line.strip("\n")
            if newline.startswith("#"):
                pass
            else:
                fields=newline.split()
                value_dict.update({fields[0]:fields[1]})
                #value_dict could only contain Ns
    if len(value_dict)>=1:
        """here is where value_dict_set is populated"""
        keys = []
        for k,v in value_dict.items():
            keys.append(k)
        value_dict_set = set(keys)
        new_dicts = value_dict
        for x in matrix_ids:
            if x not in value_dict_set:new_dicts.update({x:"N"})
    else:
        print("no usable information in vcf file, did you use the correct reference?")
    #new_dicts contains all calls, good and bad
    value_dict = {}
    """variety will contain a complete set of SNPs"""
    variety = []
    #This ensures that the values are in order
    for x in matrix_ids:
        for k,v in new_dicts.items():
            if x == k:
                variety.append(v)
    variety_set = set(variety)
    #Changed this from a to w on 12/5/19
    out_matrix = open("%s.tmp.xyx.matrix" % name,"w")
    if len(variety)>=1:
        if "A" or "T" or "G" or "C" in value_dict_set:
            out_matrix.write('%s\n' % name)
            for x in matrix_ids:
                if x in new_dicts:
                    out_matrix.write("%s\n" % new_dicts.get('%s' % x))
        else:
            print("sample %s had no usable positions!!!" % name)
    else:
        print("sample %s has problems" % name)
        print("-------------------------")
    out_matrix.close()
    value_dict_set = []
    return new_dicts

def grab_names():
    """These names will be used for future iterations - tested"""
    curr_dir= os.getcwd()
    outnames = [ ]
    for infile in glob.glob(os.path.join(curr_dir, "*.filtered.vcf")):
        name=get_seq_name(infile)
        reduced=name.replace(".filtered.vcf","")
        outnames.append(reduced)
    return outnames

def parse_likelihoods(infile):
    """This function parses the likelihoods output from RAxML - tested"""
    like_dict = {}
    with open(infile) as my_file:
        for line in my_file:
            fields = line.split()
            try:
                like_dict[fields[0]].append(fields[2])
            except KeyError:
                like_dict[fields[0]] = [fields[2]]
    print("sample_name"+"\t"+"insertion_likelihood"+"\t"+"number of potential insertion nodes")
    for k,v in like_dict.items():
        print(k+"\t"+v[0]+"\t"+str(len(v)))
    return like_dict

def calculate_pairwise_tree_dists(intree, output):
    """uses dendropy function to calculate all pairwise distances between tree - tested"""
    tree = dendropy.Tree.get_from_path(intree, "newick", preserve_underscores=True)
    outfile = open("%s" % output, "w")
    distances = tree.phylogenetic_distance_matrix()
    distances_sets = [ ]
    for i,t1 in enumerate(tree.taxon_namespace):
        for t2 in tree.taxon_namespace[i+1:]:
            distances_sets.append(distances(t1, t2))
    try:
        for i, t1 in enumerate(tree.taxon_namespace):
            for t2 in tree.taxon_namespace[i+1:]:
                outfile.write("Distance between '%s' and '%s': %s\n" % (t1.label, t2.label, distances(t1, t2)))
    except:
        print("problem iterating through tree. Tree is empty or not Newick format")
    outfile.close()
    return distances_sets

def get_closest_dists_new(final_sets, outnames):
    """tested"""
    results = []
    for final_set in final_sets:
        if len(final_set) == 0:
            pass
        results.append(final_set[1]+final_set[2])
    return results

def find_two_new(infile,outnames):
    """find two closest genomes to each query genome,
    return the names and distances (sorted), for the
    two genomes - tested"""
    distances = ()
    output_tuples = ()
    for outname in outnames:
        outname_tuple = ()
        with open(infile) as my_file:
            for line in my_file:
                fields=line.split()
                new_fields=[ ]
                for x in fields:
                    new_fields.append(re.sub('[:,]', '', x))
                final_fields=[ ]
                for y in new_fields:
                    final_fields.append(y.replace("QUERY___",""))
                """We need the distance of all samples, compared to the reference"""
                if "Reference" in final_fields[2].replace("'",""):
                    distances=((final_fields[2].replace("'",""),final_fields[4].replace("'",""),final_fields[5].replace("'","")),)+distances
                elif "Reference" in final_fields[4].replace("'",""):
                    distances=((final_fields[4].replace("'",""),final_fields[2].replace("'",""),final_fields[5].replace("'","")),)+distances
                if final_fields[4].replace("'","") == outname:
                    if "Reference" not in final_fields[2].replace("'","") and final_fields[2].replace("'","") not in outnames:
                        outname_tuple=((final_fields[4].replace("'",""),final_fields[2].replace("'",""),final_fields[5].replace("'","")),)+outname_tuple
                elif final_fields[2].replace("'","") == outname:
                    if "Reference" not in final_fields[4].replace("'","") and final_fields[4].replace("'","") not in outnames:
                        outname_tuple=((final_fields[2].replace("'",""),final_fields[4].replace("'",""),final_fields[5].replace("'","")),)+outname_tuple
        for result in sorted(outname_tuple,key=lambda x: float(x[2]))[:2]:
            output_tuples=((result),)+output_tuples
    return output_tuples,distances

def _perform_workflow_subsample_snps_dev(data):
    final_set = data[0]
    """final_set is a tuple of sample, NN, value"""
    used_snps = data[1]
    """used_snps is a dictionary"""
    subnums = data[2]
    """subnums is an integer"""
    allsnps = data[3]
    """allsnps is a list of all SNPs used"""
    matrix = data[4]
    for k,v in used_snps.items():
        if final_set[0]==k:
            for x in range(1,int(subnums)+1):
                kept_snps=random.sample(set(allsnps), int(v))
                solids = set(kept_snps)
                if os.path.isfile("%s.%s.%s.tmp.matrix" % (k,x,final_set[1])):
                    pass
                else:
                    outfile = open("%s.%s.%s.tmp.matrix" % (k,x,final_set[1]), "w")
                    with open(matrix) as my_matrix:
                        firstLine = my_matrix.readline()
                        outfile.write(firstLine)
                        first_fields = firstLine.split()
                        fixed_fields = []
                        for y in first_fields:
                            fixed_fields.append(re.sub('[:,]', '', y))
                        gindex=fixed_fields.index(final_set[1])
                        for line in my_matrix:
                            matrix_fields=line.split()
                            if matrix_fields[0] in solids:
                                outfile.write(line)
                            else:
                                outfile.write("\t".join(matrix_fields[:gindex])+"\t"+"-"+"\t"+"\t".join(matrix_fields[gindex+1:])+"\n")
                    outfile.close()
        else:
            pass

def subsample_snps_2(final_sets,used_snps,subnums,allsnps,processors,matrix):
    files_and_temp_names = []
    for f in final_sets:
        files_and_temp_names.append([f,used_snps,subnums,allsnps,matrix])
    mp_shell(_perform_workflow_subsample_snps_dev, files_and_temp_names, processors)

def subsample_snps_dev(matrix, final_set, used_snps, subnums, allsnps):
    """needs testing"""
    for k,v in used_snps.items():
        if final_set[0]==k:
            for x in range(1,int(subnums)+1):
                kept_snps=random.sample(set(allsnps), int(v))
                solids = set(kept_snps)
                if os.path.isfile("%s.%s.%s.tmp.matrix" % (k,x,final_set[1])):
                    pass
                else:
                    outfile = open("%s.%s.%s.tmp.matrix" % (k,x,final_set[1]), "w")
                    in_matrix=open(matrix)
                    firstLine = in_matrix.readline()
                    outfile.write(firstLine)
                    first_fields = firstLine.split()
                    fixed_fields = []
                    for y in first_fields:
                        fixed_fields.append(re.sub('[:,]', '', y))
                    gindex=fixed_fields.index(final_set[1])
                    for line in in_matrix:
                        matrix_fields=line.split()
                        if matrix_fields[0] in solids:
                            outfile.write(line)
                        else:
                            outfile.write("\t".join(matrix_fields[:gindex])+"\t"+"-"+"\t"+"\t".join(matrix_fields[gindex+1:])+"\n")
                    in_matrix.close()
                    outfile.close()
        else:
            pass

def get_all_snps(matrix):
    """tested"""
    allSNPs = [ ]
    with open(matrix) as my_matrix:
        for line in my_matrix:
            if line.startswith("LocusID"):
                pass
            else:
                fields=line.split()
                allSNPs.append(fields[0])
    return allSNPs

def create_params_files(id, to_prune_set, full_tree, full_matrix, dist_sets, processors):
    """not currently tested, but needs to be"""
    if int(processors)<=2:
        my_processors = 2
    else:
        my_processors = int(int(processors)/2)
    for item in to_prune_set:
        new_name = str(id)+str(item)
        if os.path.isfile("%s-PARAMS" % new_name):
            continue
        else:
            tmptree = open("%s.tmp.tree" % new_name, "w")
            to_prune = []
            for x in dist_sets:
                if x[0] == id:
                    if x[1] == item or x[2] == item:
                        to_prune.append(x[1])
            to_prune_fixed=[]
            for x in to_prune:
                to_prune_fixed.append(re.sub('[:,]', '', x))
            tree_full = dendropy.Tree.get_from_path(full_tree,schema="newick",preserve_underscores=True)
            tree_full.prune_taxa_with_labels(to_prune_fixed)
            final_tree = branch_lengths_2_decimals(tree_full.as_string("newick"))
            tmptree.write(final_tree)
            tmptree.close()
            tmptree2 = open("%s.tree" % new_name, "w")
            with open("%s.tmp.tree" % new_name) as my_file:
                for line in my_file:
                    if line.startswith("[&U]"):
                        fields = line.split()
                        fixed_fields = [ ]
                        for x in fields:
                            fixed_fields.append(x.replace("'",""))
                        tmptree2.write(fixed_fields[1])
                    else:
                        pass
            tmptree2.close()
            """result is a pruned tree that is ready for RAxML"""
            matrix_to_fasta(full_matrix, "%s.fasta" % new_name)
            os.system("sed 's/://g' %s.fasta | sed 's/,//g' > %s_in.fasta" % (new_name, new_name))
            prune_fasta(to_prune, "%s_in.fasta" % new_name, "%s_pruned.fasta" % new_name)
            try:
                subprocess.check_call("rm RAxML*%s-PARAMS" % new_name, shell=True, stderr=open(os.devnull, 'w'))
            except:
                pass
            #forcing GTRGAMMA
            try:
                subprocess.check_call("raxmlHPC-PTHREADS-SSE3 -T %s -f e -m GTRGAMMA -s %s_pruned.fasta -t %s.tree -n %s-PARAMS > /dev/null 2>&1" % (my_processors, new_name, new_name, new_name), shell=True)
                os.system("mv RAxML_binaryModelParameters.%s-PARAMS %s-PARAMS" % (new_name, new_name))
            except:
                continue

def process_temp_matrices_dev(dist_sets, sample, tree, processors, patristics, insertion_method, parameters, model):
    """not currently tested, but needs to be"""
    name=get_seq_name(sample)
    split_fields=name.split(".")
    outfile=open("%s.%s.subsample.distances.txt" % (split_fields[0],split_fields[2]), "a")
    name_fixed = []
    name_fixed.append(re.sub('[:,]', '', split_fields[2]))
    to_prune = []
    for x in dist_sets:
        if x[0] == split_fields[0]:
            if x[1] == split_fields[2] or x[2] == split_fields[2]:
                to_prune.append(x[1])
    to_prune_fixed=[]
    for x in to_prune:
        to_prune_fixed.append(re.sub('[:,]', '', x))
    #full_context includes sample, near neighbor, and replicate
    full_context = split_fields[0]+split_fields[1]+split_fields[2]
    new_name = split_fields[0]+split_fields[2]
    if os.path.isfile("%s.tree_including_unknowns_noedges.tree" % full_context):
        print("tree already present, skipping")
        pass
    else:
        tree_full = dendropy.Tree.get_from_path(tree,schema="newick",preserve_underscores=True)
        tree_full.prune_taxa_with_labels(to_prune_fixed)
        tmptree = open("%s.tmp.tree" % full_context, "w")
        final_tree = branch_lengths_2_decimals(tree_full.as_string("newick"))
        tmptree.write(final_tree)
        tmptree.close()
        tmptree2 = open("%s.tree" % full_context, "w")
        #for line in open("%s.tmp.tree" % full_context, "U"):
        with open("%s.tmp.tree" % full_context) as my_file:
            for ine in my_file:
                if line.startswith("[&U]"):
                    fields = line.split()
                    fixed_fields = [ ]
                    for x in fields:
                        fixed_fields.append(x.replace("'",""))
                    tmptree2.write(fixed_fields[1])
                else:
                    pass
        tmptree2.close()
        matrix_to_fasta(sample, "%s.fasta" % full_context)
        os.system("sed 's/://g' %s.fasta | sed 's/,//g' > %s_in.fasta" % (full_context, full_context))
        if os.path.isfile("%s-PARAMS" % new_name):
            try:
                run_raxml("%s_in.fasta" % full_context, "%s.tree" % full_context, "%s.subsampling_classifications.txt" % full_context, insertion_method, "%s-PARAMS" % new_name, "GTRGAMMA", "%s" % full_context)
            except:
                pass
        else:
            try:
                subprocess.check_call("raxmlHPC-PTHREADS-SSE3 -T %s -f e -m GTRGAMMA -s %s_pruned.fasta -t %s.tree -n %s-PARAMS --no-bfgs > /dev/null 2>&1" % (my_processors, new_name, new_name, new_name), shell=True)
                os.system("mv RAxML_binaryModelParameters.%s-PARAMS %s-PARAMS" % (new_name, new_name))
                run_raxml("%s_in.fasta" % full_context, "%s.tree" % full_context, "%s.subsampling_classifications.txt" % full_context, insertion_method, "%s-PARAMS" % (split_fields[0]+split_fields[2]), "GTRGAMMA", "%s" % full_context)
            except:
                pass
    if os.path.isfile("%s.resampling_distances.txt" % full_context):
        pass
    else:
        try:
            calculate_pairwise_tree_dists("%s.tree_including_unknowns_noedges.tree" % full_context, "%s.resampling_distances.txt" % full_context)
            for line in open("%s.resampling_distances.txt" % full_context,"U"):
                resample_fields = line.split()
                myid = re.sub("[:']", "",resample_fields[4])
                fixedid = myid.replace("QUERY___","")
                newid = re.sub("[:']","",resample_fields[2])
                fixedid2 = newid.replace("QUERY___","")
                if resample_fields[2] == "'Reference'" and fixedid in name_fixed:
                    outfile.write("resampled distance between Reference and %s = %s\n" % (fixedid, resample_fields[5]))
                elif resample_fields[4] == "'Reference':" and fixedid2 in name_fixed:
                    outfile.write("resampled distance between Reference and %s = %s\n" % (fixedid2, resample_fields[5]))
                else:
                    pass
        except:
            pass

def check_input_files(matrix,reference):
    ref_names = []
    with open(reference) as my_ref:
        for record in SeqIO.parse(my_ref,"fasta"):
            ref_names.append(record.id)
    with open(matrix) as f:
        for line in f.readlines()[:2]:
            if line.startswith("LocusID"):
                pass
            else:
                fields = line.split()
                name_fields=fields[0].split("::")
                if name_fields[0] in ref_names:
                    pass
                else:
                    print("The IDs in your Reference don't match the names in your SNP matrix! Please fix and re-start...exiting...")
                    sys.exit()

def create_merged_vcf():
    out_file = open("merged.vcf", "w")
    start_dir = os.getcwd()
    lists = []
    for infile in glob.glob(os.path.join(start_dir, "*.tmp.xyx.matrix")):
        data = open(infile).read().splitlines()
        lists.append(data)
    test=map(list, zip(*lists))
    for x in test:
        out_file.write("\t".join(x))
        out_file.write("\n")
    out_file.close()

def _perform_workflow_create_params(data):
    id = data[0]
    to_prune_set = data[1]
    full_tree = data[2]
    full_matrix = data[3]
    dist_sets = data[4]
    processors = data[5]
    if int(processors)<=2:
        my_processors = 2
    else:
        my_processors = int(int(processors)/2)
    for item in to_prune_set:
        new_name = str(id)+str(item)
        if os.path.isfile("%s-PARAMS" % new_name):
            pass
        else:
            tmptree = open("%s.tmp.tree" % new_name, "w")
            to_prune = []
            for x in dist_sets:
                if x[0] == id:
                    if x[1] == item or x[2] == item:
                        to_prune.append(x[1])
            to_prune_fixed=[]
            for x in to_prune:
                to_prune_fixed.append(re.sub('[:,]', '', x))
            tree_full = dendropy.Tree.get_from_path(full_tree,schema="newick",preserve_underscores=True)
            tree_full.prune_taxa_with_labels(to_prune_fixed)
            final_tree = branch_lengths_2_decimals(tree_full.as_string("newick"))
            tmptree.write(final_tree)
            tmptree.close()
            tmptree2 = open("%s.tree" % new_name, "w")
            with open("%s.tmp.tree" % new_name) as my_file:
                for line in my_file:
                    line = line.replace("'","")
                    tmptree2.write(line)
            tmptree2.close()
            """result is a pruned tree that is ready for RAxML"""
            matrix_to_fasta(full_matrix, "%s.fasta" % new_name)
            os.system("sed 's/://g' %s.fasta | sed 's/,//g' > %s_in.fasta" % (new_name, new_name))
            prune_fasta(to_prune,"%s_in.fasta" % new_name,"%s_pruned.fasta" % new_name)
            try:
                subprocess.check_call("rm RAxML*%s-PARAMS" % new_name, shell=True, stderr=open(os.devnull, 'w'))
            except:
                pass
            #forcing GTRGAMMA
            try:
                subprocess.check_call("raxmlHPC-PTHREADS-SSE3 -T %s -f e -m GTRGAMMA -s %s_pruned.fasta -t %s.tree -n %s-PARAMS > /dev/null 2>&1" % (my_processors, new_name, new_name, new_name), shell=True)
                os.system("mv RAxML_binaryModelParameters.%s-PARAMS %s-PARAMS" % (new_name, new_name))
            except:
                print("problem creating PARAMS file for %s.tree, exiting" % new_name)
                sys.exit()

def create_params_files_dev(new_sample_dicts,tree,matrix,final_sets,processors):
    to_run = []
    for k,v in new_sample_dicts.items():
        to_run.append([k,v,tree,matrix,final_sets,processors])
    mp_shell(_perform_workflow_create_params,to_run,processors)

def process_temp_matrices_2(final_sets,final_matrices,tree,processors,patristic_distances, V, parameters, model):
    to_run = []
    for matrix in final_matrices:
        to_run.append([final_sets,matrix,tree,processors,patristic_distances,V,parameters,model])
    mp_shell(_perform_workflow_temp_matrices, to_run, processors)

def _perform_workflow_temp_matrices(data):
    dist_sets = data[0]
    sample = data[1]
    tree = data[2]
    processors = data[3]
    patristics = data[4]
    insertion_method = data[5]
    parameters = data[6]
    model = data[7]
    name=get_seq_name(sample)
    split_fields=name.split(".")
    outfile=open("%s.%s.subsample.distances.txt" % (split_fields[0],split_fields[2]), "a")
    name_fixed = []
    name_fixed.append(re.sub('[:,]', '', split_fields[2]))
    to_prune = []
    for x in dist_sets:
        if x[0] == split_fields[0]:
            if x[1] == split_fields[2] or x[2] == split_fields[2]:
                to_prune.append(x[1])
    to_prune_fixed=[]
    for x in to_prune:
        to_prune_fixed.append(re.sub('[:,]', '', x))
    #full_context includes sample, near neighbor, and replicate
    full_context = split_fields[0]+split_fields[1]+split_fields[2]
    new_name = split_fields[0]+split_fields[2]
    if os.path.isfile("%s.tree_including_unknowns_noedges.tree" % full_context):
        print("tree already present, skipping")
        pass
    else:
        tree_full = dendropy.Tree.get_from_path(tree,schema="newick",preserve_underscores=True)
        tree_full.prune_taxa_with_labels(to_prune_fixed)
        tmptree = open("%s.tmp.tree" % full_context, "w")
        final_tree = branch_lengths_2_decimals(tree_full.as_string("newick"))
        tmptree.write(final_tree)
        tmptree.close()
        tmptree2 = open("%s.tree" % full_context, "w")
        with open("%s.tmp.tree" % full_context) as my_file:
            for line in my_file:
                line = line.replace("'","")
                tmptree2.write(line)
        tmptree2.close()
        matrix_to_fasta(sample, "%s.fasta" % full_context)
        os.system("sed 's/://g' %s.fasta | sed 's/,//g' > %s_in.fasta" % (full_context, full_context))
        if os.path.isfile("%s-PARAMS" % new_name):
            try:
                run_raxml("%s_in.fasta" % full_context, "%s.tree" % full_context, "%s.subsampling_classifications.txt" % full_context, insertion_method, "%s-PARAMS" % new_name, "GTRGAMMA", "%s" % full_context)
            except:
                pass
        else:
            try:
                subprocess.check_call("raxmlHPC-PTHREADS-SSE3 -T %s -f e -m GTRGAMMA -s %s_pruned.fasta -t %s.tree -n %s-PARAMS --no-bfgs > /dev/null 2>&1" % (my_processors, new_name, new_name, new_name), shell=True)
                os.system("mv RAxML_binaryModelParameters.%s-PARAMS %s-PARAMS" % (new_name, new_name))
                run_raxml("%s_in.fasta" % full_context, "%s.tree" % full_context, "%s.subsampling_classifications.txt" % full_context, insertion_method, "%s-PARAMS" % (split_fields[0]+split_fields[2]), "GTRGAMMA", "%s" % full_context)
            except:
                pass
    if os.path.isfile("%s.resampling_distances.txt" % full_context):
        pass
    else:
        try:
            calculate_pairwise_tree_dists("%s.tree_including_unknowns_noedges.tree" % full_context, "%s.resampling_distances.txt" % full_context)
            for line in open("%s.resampling_distances.txt" % full_context,"U"):
                resample_fields = line.split()
                myid = re.sub("[:']", "",resample_fields[4])
                fixedid = myid.replace("QUERY___","")
                newid = re.sub("[:']","",resample_fields[2])
                fixedid2 = newid.replace("QUERY___","")
                if resample_fields[2] == "'Reference'" and fixedid in name_fixed:
                    outfile.write("resampled distance between Reference and %s = %s\n" % (fixedid, resample_fields[5]))
                elif resample_fields[4] == "'Reference':" and fixedid2 in name_fixed:
                    outfile.write("resampled distance between Reference and %s = %s\n" % (fixedid2, resample_fields[5]))
                else:
                    pass
        except:
            pass
    outfile.close()

def qc_files(fasta,tree):
    #first grab the IDs from the FASTA file:
    fasta_ids = []
    with open(fasta) as my_fasta:
        for record in SeqIO.parse(fasta,"fasta"):
            fasta_ids.append(record.id)
    #Now parse the names from the tree:
    tree_ids = []
    mytree = Phylo.read(tree,'newick')
    for clade in mytree.find_clades():
        if clade.name:
            tree_ids.append(clade.name)
    tree_set = set(tree_ids)
    fasta_set = set(fasta_ids)
    if len(fasta_set) == len(tree_set):
        print("new samples didn't get added correctly...exiting")
        sys.exit()
    else:
        pass
<file_sep>/tools/subsample_reads_and_place.py
#!/usr/bin/env python

"""From a SNP matrix, sub-sample
SNPs from a given genome, then find
the SNP level where the unknown genome
is placed correctly >90% of the time"""

from optparse import OptionParser
import sys
import random
import re
import os
import subprocess
import dendropy
import collections

try:
    from dendropy import treecalc
except:
    print("Dendropy needs to be installed for this script to run")
    sys.exit()
from subprocess import Popen
try:
    from Bio import SeqIO
except:
    print("BioPython needs to be installed for this script to run")
    sys.exit()

def test_file(option, opt_str, value, parser):
    try:
        with open(value): setattr(parser.values, option.dest, value)
    except IOError:
        print('%s file cannot be opened' % option)
        sys.exit()

def subsample_snps(matrix, name, start):
    """get a list of all possible positions, depending
    on those positions in the original matrix.  Similar
    to method in the main script"""
    allSNPs = []
    with open(matrix) as my_matrix:
        for line in my_matrix:
            if line.startswith("LocusID"):
                pass
            else:
                fields=line.split()
                allSNPs.append(fields[0])
    kept_snps=random.sample(set(allSNPs), int(start))
    outfile = open("%s.%s.tmp.matrix" % (name,start), "w")
    with open(matrix) as my_matrix:
        firstLine = my_matrix.readline()
        outfile.write(firstLine)
        first_fields = firstLine.split()
        last=first_fields.index("#SNPcall")
        mygenome=first_fields.index(name)
        fixed_fields = []
        for x in first_fields[:last]:
            fixed_fields.append(re.sub('[:,]', '', x))
        gindex = []
        for x in first_fields[:last]:
            gindex.append(first_fields.index(x))
        for line in my_matrix:
            matrix_fields=line.split()
            if matrix_fields[0] in kept_snps:
                outfile.write(line)
            else:
                outfile.write("\t".join(matrix_fields[:mygenome])+"\t"+"-"+"\t"+"\t".join(matrix_fields[mygenome+1:])+"\n")
    outfile.close()
    return last

def subsample_snps_keep(matrix, name, start):
    """get a list of all possible positions, depending
    on those positions in the original matrix.  Similar
    to method in the main script"""
    allSNPs = []
    with open(matrix) as my_matrix:
        for line in my_matrix:
            if line.startswith("LocusID"):
                pass
            else:
                fields=line.split()
                allSNPs.append(fields[0])
    kept_snps=random.sample(set(allSNPs), int(start))
    outfile = open("%s.%s.tmp.matrix" % (name,start), "w")
    with open(matrix) as in_matrix:
        firstLine = in_matrix.readline()
        first_fields = firstLine.split()
        last=first_fields.index("#SNPcall")
        first_fields.insert(last,"QUERY_%s" % name)
        outfile.write("\t".join(first_fields)+"\n")
        """mygenome is the index of the genome that we want to subsample"""
        mygenome=first_fields.index(name)
        fixed_fields = []
        for x in first_fields[:last]:
            fixed_fields.append(re.sub('[:,]', '', x))
        gindex = [ ]
        for x in first_fields[:last]:
            gindex.append(first_fields.index(x))
        for line in in_matrix:
            matrix_fields=line.split()
            if matrix_fields[0] in kept_snps:
                outfile.write("\t".join(matrix_fields[:last])+"\t"+"\t".join(matrix_fields[mygenome])+"\t"+"".join(matrix_fields[last:])+"\n")
            else:
                outfile.write("\t".join(matrix_fields[:last])+"\t"+"-"+"\t"+"\t".join(matrix_fields[last:])+"\n")
    outfile.close()
    return last

def matrix_to_fasta(matrix_in,name,last):
    """converts a SNP matrix to fasta format.
    Again, slightly different output compared to tested
    function"""
    reduced = []
    out_fasta = open("%s.fasta" % name, "w")
    with open(matrix_in) as my_matrix:
        for line in my_matrix:
            fields = line.split("\t")
            reduced.append(fields[1:last])
    test=map(list, zip(*reduced))
    for x in test:
        out_fasta.write(">"+str(x[0])+"\n"+"".join(x[1:])+"\n")
    out_fasta.close()

def branch_lengths_2_decimals(str_newick_tree):
    """replaces branch lengths in scientific notation with decimals
    Identical to tested function in main script"""
    colon_s = 0
    comma_back_paren_s = 0
    num = ''
    new_tree = ''
    for count, char in enumerate(str_newick_tree):
        if char == ':':
            colon_s = count
            continue
        if char in (')', ','):
            comma_back_paren_s = 1
            num = '%f' % float(num)
            new_tree += ":" + num
            colon_s = 0
            num = ''
        if colon_s != 0:
            num = num + char
        if colon_s == 0:
            new_tree += char
    new_tree = new_tree.strip('\'').strip('\"').strip('\'') + ";"
    return new_tree

def insert_sequence(in_fasta, tree, fixed_name, parameters, processors):
    args = ['raxmlHPC-PTHREADS-SSE3', '-f', 'V',
            '-s', '%s' % in_fasta, '-m', 'GTRGAMMA', '-n', '%s' % fixed_name, '-t',
            '%s' % tree, '-T', '%s' % processors, '-R', '%s' % parameters, '--no-bfgs', '>', '/dev/null 2>&1']
    try:
        vcf_fh = open('raxml.out', 'w')
    except:
        print('could not open raxml file')
    try:
        log_fh = open('raxml.log', 'w')
    except:
        print('could not open log file')
    try:
        raxml_run = Popen(args, stderr=log_fh, stdout=vcf_fh)
        raxml_run.wait()
        os.system("sed 's/\[[^]]*\]//g' RAxML_labelledTree.%s > %s.tree_including_unknowns_noedges.tree" % (fixed_name,fixed_name))
        subprocess.check_call("rm RAxML*.%s" % fixed_name, shell=True)
    except:
        print("sequence(s) were not inserted into tree!!!!!")

def prune_tree(fixed_name,tree):
    tree_full = dendropy.Tree.get_from_path(tree,schema="newick",preserve_underscores=True)
    tree_full.prune_taxa_with_labels(["%s" % fixed_name])
    final_tree = branch_lengths_2_decimals(tree_full.as_string("newick"))
    tmptree = open("%s.tmpx.tree" % fixed_name, "w")
    tmptree.write(final_tree)
    tmptree.close()
    tmptree2 = open("%s.tmpxz.tree" % fixed_name, "w")
    with open("%s.tmpx.tree" % fixed_name) as my_tree:
        for line in my_tree:
            fields = line.split()
            fixed_fields = []
            for x in fields:
                fixed_fields.append(x.replace("'",""))
            else:
                pass
            tmptree2.write("".join(fixed_fields))
    tmptree2.close()

def calculate_pairwise_tree_dists(intree, output):
    tree = dendropy.Tree.get_from_path(intree, "newick", preserve_underscores=True)
    outfile = open("%s" % output, "w")
    distances = tree.phylogenetic_distance_matrix()
    for i, t1 in enumerate(tree.taxon_namespace):
        for t2 in tree.taxon_namespace[i+1:]:
            outfile.write("Distance between '%s' and '%s': %s\n" % (t1.label, t2.label, distances(t1, t2)))
    outfile.close()

def parse_distances(distance_file,fixed_name):
    true_value = []
    with open(distance_file) as my_file:
        for line in my_file:
            if len(fixed_name)>1:
                for name in fixed_name:
                    #print(name)
                    fields = line.split()
                    if fields[2] == "'%s'" % name and fields[4] == "'Reference':":
                        print(fields)
                        true_value.append(fields[5])
                    elif fields[2] =="'Reference'" and fields[4] =="'%s':" % name:
                        true_value.append(fields[5])
            else:
                fields = line.split()
                if fields[2] == "'%s'" % ''.join(fixed_name) and fields[4] == "'Reference':":
                    true_value.append(fields[5])
                elif fields[2] =="'Reference'" and fields[4] =="'%s':" % ''.join(fixed_name):
                    true_value.append(fields[5])
                else:
                    pass
    return true_value

def get_name_by_ID(in_fasta, ID, out_fasta):
    output_handle = open(out_fasta, "w")
    seqrecords=[]
    with open(in_fasta) as my_fasta:
        for record in SeqIO.parse(my_fasta, "fasta"):
            if record.id == ID:
                seqrecords.append(record)
    SeqIO.write(seqrecords, output_handle, "fasta")
    output_handle.close()

def rename_fasta(in_fasta, name, out_fasta):
    output_handle = open(out_fasta, "w")
    seqrecords=[]
    with open(in_fasta) as my_fasta:
        for record in SeqIO.parse(my_fasta, "fasta"):
            seqrecords.append(record.seq)
    output_handle.write(">"+name+"\n")
    for seqrecord in seqrecords:
        output_handle.write(str(seqrecord)+"\n")
    output_handle.close()

def get_field_index(matrix_in):
    firstLine = open(matrix_in).readline()
    first_fields = firstLine.split("\t")
    last=first_fields.index("#SNPcall")
    return last

def remove_sequence(in_fasta,name,out_fasta):
    seqrecords = []
    output_handle = open(out_fasta, "w")
    with open(in_fasta) as my_fasta:
        for record in SeqIO.parse(my_fasta, "fasta"):
            if record.id != name:
                seqrecords.append(record)
    SeqIO.write(seqrecords, output_handle, "fasta")
    output_handle.close()

def main(matrix,tree,name,start,step,end,processors,iterations,deviation,remove):
    aa = subprocess.call(['which', 'raxmlHPC-PTHREADS-SSE3'])
    if aa == 0:
        pass
    else:
        print("RAxML must be in your path as raxmlHPC-PTHREADS-SSE3")
        sys.exit()
    """get starting information"""
    start_dir = os.getcwd()
    start_path = os.path.abspath("%s" % start_dir)
    matrix_path = os.path.abspath("%s" % matrix)
    tree_path = os.path.abspath("%s" % tree)
    """done with getting starting information"""
    os.system("mkdir %s/%s.tmp" % (start_path,name))
    fixed_name = []
    #Changing directory into temporary one
    os.chdir("%s/%s.tmp" % (start_path,name))
    os.system("sed 's/://g' %s | sed 's/,//g' > REF.matrix" % matrix_path)
    fixed_name.append(re.sub('[:,]', '', name))
    calculate_pairwise_tree_dists(tree_path, "%s.all_snps_patristic_distances.txt" % "".join(fixed_name))
    last=get_field_index(matrix_path)
    matrix_to_fasta("REF.matrix","REF",last)
    remove_sequence("REF.fasta", "".join(fixed_name), "REF_pruned.fasta")
    true_value = parse_distances("%s.all_snps_patristic_distances.txt" % "".join(fixed_name),fixed_name)
    outfile = open("%s.results.out" % ''.join(fixed_name), "w")
    if remove == "T":
        prune_tree(''.join(fixed_name),tree_path)
    else:
        pass
    print("creating parameters file")
    if remove == "T":
        subprocess.check_call("raxmlHPC-PTHREADS-SSE3 -f e -m GTRGAMMA -s REF_pruned.fasta -t %s.tmpxz.tree -n PARAMS --no-bfgs -T %s > /dev/null 2>&1" % ("".join(fixed_name),processors) , shell=True)
        subprocess.check_call("mv RAxML_binaryModelParameters.PARAMS %s.PARAMS" % "".join(fixed_name), shell=True)
    elif remove == "F":
        subprocess.check_call("raxmlHPC-PTHREADS-SSE3 -f e -m GTRGAMMA -s REF.fasta -t %s -n PARAMS --no-bfgs -T %s > /dev/null 2>&1" % (tree_path,processors) , shell=True)
        subprocess.check_call("mv RAxML_binaryModelParameters.PARAMS %s.PARAMS" % "".join(fixed_name), shell=True)
    else:
        print("you need to choose between 'T' and 'F' for remove value. Exiting...")
        sys.exit()
    print("starting loop")
    for i in range(start, end+1, step):
        hits = []
        for j in range(1,iterations+1):
            if remove == "T":
                last=subsample_snps("REF.matrix", "".join(fixed_name), i)
            else:
                last=subsample_snps_keep("REF.matrix", "".join(fixed_name), i)
            os.system("sed 's/://g' %s.%s.tmp.matrix | sed 's/,//g' > %s.%s.tmp.fixed.matrix" % ("".join(fixed_name),i,"".join(fixed_name),i))
            last=get_field_index("%s.%s.tmp.fixed.matrix" % ("".join(fixed_name),i))
            matrix_to_fasta("%s.%s.tmp.fixed.matrix" % ("".join(fixed_name),i), "%s.%s" % ("QUERY_"+"".join(fixed_name),i), last)
            get_name_by_ID("%s.%s.fasta" % ("QUERY_"+"".join(fixed_name),i), ''.join(fixed_name), "%s.%s.%s.tmp.fasta" % ("QUERY_"+"".join(fixed_name),i,j))
            tmp_name = "QUERY_"+''.join(fixed_name)+str(j)
            rename_fasta("%s.%s.%s.tmp.fasta" % ("QUERY_"+"".join(fixed_name),i,j), tmp_name,"%s.%s.%s.zzyzz.fasta" % ("QUERY_"+"".join(fixed_name),i,j))
        if remove == "T":
            os.system("cat %s.*.zzyzz.fasta REF.fasta > %s.joined.fasta" % ("QUERY_"+"".join(fixed_name),"".join(fixed_name)))
        else:
            os.system("cat %s.*.zzyzz.fasta REF.fasta > %s.joined.fasta" % ("QUERY_"+"".join(fixed_name),"QUERY_"+"".join(fixed_name)))
        os.system("rm %s.*.tmp.fasta %s.*.zzyzz.fasta" % ("QUERY_"+"".join(fixed_name),"QUERY_"+"".join(fixed_name)))
        if remove =="T":
            insert_sequence("%s.joined.fasta" % "".join(fixed_name), "%s.tmpxz.tree" % "".join(fixed_name), ''.join(fixed_name), "%s.PARAMS" % "".join(fixed_name), processors)
        else:
            insert_sequence("%s.joined.fasta" % ("QUERY_"+"".join(fixed_name)), "%s" % tree_path, ''.join(fixed_name), "%s.PARAMS" % "".join(fixed_name), processors)
        calculate_pairwise_tree_dists("%s.tree_including_unknowns_noedges.tree" % "".join(fixed_name),"%s.all_patristic_distances.txt" % "".join(fixed_name))
        os.system("cp %s.tree_including_unknowns_noedges.tree %s.%s.%s.tree" % ("".join(fixed_name),"".join(fixed_name),i,j))
        query_names = []
        for j in range(1,iterations+1):
            #query_names.append("QUERY__"+"QUERY_"+"".join(fixed_name)+str(j))
            query_names.append("QUERY___"+"QUERY_"+"".join(fixed_name)+str(j))
        subsampled_values = parse_distances("%s.all_patristic_distances.txt" % "".join(fixed_name), query_names)
        for value in subsampled_values:
            if (float(value)/float(''.join(true_value)))<(1+float(deviation)):
                hits.append("1")
            else:
                pass
        outfile.write(str(i)+"\t"+str(len(hits))+"\n")
        if int(len(hits))>=int(0.95*iterations):
            print("optimal value is for %s is %s" % ("".join(fixed_name),i))
            break
    if len(hits) == 0:
        print("no results at your level of sampling. Change settings and try again")
    os.system("mv %s.results.out %s" % (''.join(fixed_name), start_path))
    os.chdir("%s" % start_path)
    os.system("rm -rf %s/%s.tmp" % (start_path,name))

if __name__ == "__main__":
    usage="usage: %prog [options]"
    parser = OptionParser(usage=usage)
    parser.add_option("-m", "--snp_matrix", dest="matrix",
                      help="path to NASP snp_matrix [REQUIRED]",
                      action="callback", callback=test_file, type="string")
    parser.add_option("-t", "--tree", dest="tree",
                      help="path to input tree [REQUIRED]",
                      action="callback", callback=test_file, type="string")
    parser.add_option("-n", "--name", dest="name",
                      help="name of genome to test [REQUIRED]",
                      action="store", type="string")
    parser.add_option("-s", "--start", dest="start",
                      help="starting number of SNPs to sample, defaults to 50",
                      action="store", type="int", default="50")
    parser.add_option("-p", "--step", dest="step",
                      help="step through SNPs at this level, defaults to 100",
                      action="store", type="int", default="100")
    parser.add_option("-e", "--end", dest="end",
                      help="ending number of SNPs to sample, defaults to 100000",
                      action="store", type="int", default="100000")
    parser.add_option("-o", "--processors", dest="processors",
                      help="number of processors to use with RAxML, defaults to 4",
                      action="store", type="int", default="4")
    parser.add_option("-i", "--iterations", dest="iterations",
                      help="number of iterations at each level, defaults to 10",
                      action="store", type="int", default="10")
    parser.add_option("-d", "--deviation", dest="deviation",
                      help="deviation from 1, to determine correct placement, defaults to 0.05",
                      action="store", type="float", default="0.05")
    parser.add_option("-r", "--remove", dest="remove",
                      help="remove original from tree and place?  Defaults to T",
                      action="store", type="string", default="T")

    options, args = parser.parse_args()

    mandatories = ["matrix", "tree", "name"]
    for m in mandatories:
        if not options.__dict__[m]:
            print("\nMust provide %s.\n" %m)
            parser.print_help()
            exit(-1)

    main(options.matrix,options.tree,options.name,options.start,options.step,options.end,
         options.processors,options.iterations,options.deviation,options.remove)
<file_sep>/tests/test_all_functions.py
#!/usr/bin/env python

"""test each function in the wg-fast
code"""

import unittest
import os
import tempfile
import shutil
import re
import sys
from wg_fast.util import *

curr_dir=os.getcwd()

class Test1(unittest.TestCase):
    def test_get_seq_name_basic_function(self):
        self.assertEqual(get_seq_name("/path/to/test.fasta"), "test.fasta")
    """tests the condition where you use a tilda instead of full path"""
    def test_get_seq_name_tilda(self):
        self.assertEqual(get_seq_name("~/test.fasta"), "test.fasta")
    """tests the case where no path is passed"""
    def test_get_seq_name_empty(self):
        self.assertEqual(get_seq_name(""), "")
    """tests the case where something weird is passed"""
    def test_get_seq_name_wrong_slash(self):
        self.assertEqual(get_seq_name("\wrong\way"), "\\wrong\\way")

class Test2(unittest.TestCase):
    def test_get_readFile_components_basic_function(self):
        self.assertEqual(get_readFile_components("/path/to/file.gz"), ('/path/to', 'file', '.gz'))
    def test_get_readFile_components_tilda(self):
        self.assertEqual(get_readFile_components("~/path/to/file.gz"), ('~/path/to', 'file', '.gz'))
    def test_get_readFile_components_non_gz(self):
        self.assertEqual(get_readFile_components("~/path/to/file.fasta"), ('~/path/to', 'file', '.fasta'))
    def test_get_readFile_components_wrong_slash(self):
        self.assertEqual(get_readFile_components("~\path\to\file.fasta"), ('', '~\\path\to\x0cile', '.fasta'))
    def test_get_readFile_components_empty(self):
        self.assertEqual(get_readFile_components(""), ('', '', ''))

class Test5(unittest.TestCase):
    def test_sort_information_basic_function(self):
        self.assertEqual(sort_information("ADK1::460"), 460)
    def test_sort_information_cant_parse(self):
        self.assertRaises(TypeError, sort_information, "ADK1__460")
    def test_sort_information_no_input(self):
        self.assertRaises(TypeError, sort_information, None)

class Test6(unittest.TestCase):
    def test_matrix_to_fasta_basic_function(self):
        tdir = tempfile.mkdtemp(prefix="filetest_",)
        fpath = os.path.join(tdir,"testfile.filtered.vcf")
        fp = open(fpath, "w")
        fp.write("LocusID\tReference\tgenome1\tgenome2\n")
        fp.write("ADK::1\tA\tT\tT\n")
        fp.write("ADK::2\tT\tT\tT\n")
        fp.close()
        self.assertEqual(matrix_to_fasta(fpath, "%s/outfile.txt" % tdir), [">ReferenceAT", ">genome1TT", ">genome2TT"])
        shutil.rmtree(tdir)
    def test_matrix_to_fasta_unequal_fields(self):
        tdir = tempfile.mkdtemp(prefix="filetest_",)
        fpath = os.path.join(tdir,"testfile.filtered.vcf")
        fp = open(fpath, "w")
        fp.write("LocusID\tReference\tgenome1\tgenome2\n")
        fp.write("ADK::1\tA\tT\tT\n")
        fp.write("ADK::2\tT\tT\n")
        fp.close()
        self.assertEqual(matrix_to_fasta(fpath, "%s/outfile.txt" % tdir), [">ReferenceAT", ">genome1TT"])
        shutil.rmtree(tdir)
    def test_matrix_to_fasta_multiple_states(self):
        tdir = tempfile.mkdtemp(prefix="filetest_",)
        fpath = os.path.join(tdir,"testfile.filtered.vcf")
        fp = open(fpath, "w")
        fp.write("LocusID\tReference\tgenome1\tgenome2\n")
        fp.write("ADK::1\tA\tT\tT\n")
        fp.write("ADK::2\tT\tT\tTT\n")
        fp.close()
        self.assertEqual(matrix_to_fasta(fpath, "%s/outfile.txt" % tdir), [">ReferenceAT", ">genome1TT", ">genome2TTT"])
        shutil.rmtree(tdir)

class Test7(unittest.TestCase):
    def test_write_reduced_matrix_basic_function(self):
        tdir = tempfile.mkdtemp(prefix="filetest_",)
        fpath = os.path.join(tdir,"testfile.filtered.vcf")
        fp = open(fpath, "w")
        fp.write("LocusID\tReference\tgenome1\tgenome2\t#SNPcall\n")
        fp.write("ADK::1\tA\tT\tT\n")
        fp.write("ADK::2\tT\tT\tT\n")
        fp.close()
        self.assertEqual(write_reduced_matrix(fpath), [4, 4])
        shutil.rmtree(tdir)
    def test_write_reduced_matrix_odd_field_numbers(self):
        tdir = tempfile.mkdtemp(prefix="filetest_",)
        fpath = os.path.join(tdir,"testfile.filtered.vcf")
        fp = open(fpath, "w")
        fp.write("LocusID\tReference\tgenome1\tgenome2\t#SNPcall\n")
        fp.write("ADK::1\tA\tT\tT\n")
        fp.write("ADK::2\tT\tT\n")
        fp.close()
        self.assertEqual(write_reduced_matrix(fpath), [4, 3])
        shutil.rmtree(tdir)

class Test8(unittest.TestCase):
    def test_make_temp_matrix_basic_function(self):
        tdir = tempfile.mkdtemp(prefix="filetest_",)
        fpath = os.path.join(tdir,"testfile.matrix")
        fp = open(fpath, "w")
        fp.write("LocusID\tReference\tgenome1\tgenome2\n")
        fp.write("ADK::1\tA\tT\tT\n")
        fp.write("ADK::2\tT\tT\tT\n")
        fp.write("ADK::3\tG\tG\tT\n")
        fp.close()
        vpath = os.path.join(tdir,"testfile.filtered.vcf")
        vp = open(vpath, "w")
        vp.write("ADK::1\tG\n")
        vp.close()
        self.assertEqual(make_temp_matrix(vpath,fpath,"test"), {'ADK::1': 'G', 'ADK::2': 'N', 'ADK::3': 'N'})
        shutil.rmtree(tdir)
    def test_make_temp_matrix_no_matches(self):
        tdir = tempfile.mkdtemp(prefix="filetest_",)
        fpath = os.path.join(tdir,"testfile.matrix")
        fp = open(fpath, "w")
        fp.write("LocusID\tReference\tgenome1\tgenome2\n")
        fp.write("ADK::1\tA\tT\tT\n")
        fp.write("ADK::2\tT\tT\tT\n")
        fp.write("ADK::3\tG\tG\tT\n")
        fp.close()
        vpath = os.path.join(tdir,"testfile.filtered.vcf")
        vp = open(vpath, "w")
        vp.write("ADK::6\tA\n")
        vp.close()
        self.assertEqual(make_temp_matrix(vpath,fpath,"test"), {'ADK::1': 'N', 'ADK::2': 'N', 'ADK::3': 'N', 'ADK::6': 'A'})
        shutil.rmtree(tdir)

class Test9(unittest.TestCase):
    def test_grab_names_basic_function(self):
        tdir = tempfile.mkdtemp(prefix="filetest_",)
        fpath = os.path.join(tdir,"test.filtered.vcf")
        fpath2 = os.path.join(tdir,"name_1.filtered.vcf")
        os.chdir("%s" % tdir)
        fp = open(fpath, "w")
        fp2 = open(fpath2, "w")
        fp.close()
        fp2.close()
        self.assertEqual(grab_names(), ['name_1', 'test'])
        os.chdir("%s" % curr_dir)
        shutil.rmtree(tdir)

#class Test10(unittest.TestCase):
#    def test_process_coverage_basic_function(self):
#        tdir = tempfile.mkdtemp(prefix="filetest_",)
#        fpath = os.path.join(tdir,"ECOLI_coverage.sample_summary")
#        os.chdir("%s" % tdir)
#        fp = open(fpath,"w")
#        fp.write("sample_id       total   mean    granular_third_quartile granular_median granular_first_quartile %_bases_above_15\n")
#        fp.write("ECOLI    2050    3.82    6       5       4       0.0\n")
#        fp.write("Total    2050    3.82    N/A     N/A     N/A")
#        fp.close()
#        self.assertEqual(process_coverage("ECOLI"), {'ECOLI':'3.82'})
#        os.chdir("%s" % curr_dir)
#        shutil.rmtree(tdir)
#    def test_process_coverage_missing_match(self):
#        tdir = tempfile.mkdtemp(prefix="filetest_",)
#        fpath = os.path.join(tdir,"ECOLI_coverage.sample_summary")
#        os.chdir("%s" % tdir)
#        fp = open(fpath,"w")
#        fp.write("sample_id       total   mean    granular_third_quartile granular_median granular_first_quartile %_bases_above_15\n")
#        fp.write("EOLI    2050    3.82    6       5       4       0.0\n")
#        fp.write("Total    2050    3.82    N/A     N/A     N/A")
#        fp.close()
#        self.assertRaises(TypeError, process_coverage, "ECOLI")
#        os.chdir("%s" % curr_dir)
#        shutil.rmtree(tdir)

class Test11(unittest.TestCase):
    def test_find_two_report_error(self):
        tdir = tempfile.mkdtemp(prefix="filetest_",)
        fpath = os.path.join(tdir,"all_patristic_distances.txt")
        os.chdir("%s" % tdir)
        fp = open(fpath,"w")
        fp.write("Distance between 'E2348_69_allexternalnucmer' and 'H10407_allexternalnucmer': 1.39030683167\n")
        fp.close()
        self.assertRaises(TypeError, find_two_new, ['E2348_69_allexternalnucmer'], ['H10407_allexternalnucmer','E2348_69_allexternalnucmer'])
        os.chdir("%s" % curr_dir)
        shutil.rmtree(tdir)
    def test_find_two_basic_function(self):
        tdir = tempfile.mkdtemp(prefix="filetest_",)
        fpath = os.path.join(tdir,"all_patristic_distances.txt")
        os.chdir("%s" % tdir)
        fp = open(fpath,"w")
        fp.write("Distance between 'E2348_69_allexternalnucmer' and 'H10407_allexternalnucmer': 1.39030683167\n")
        fp.write("Distance between 'E2348_69_allexternalnucmer' and 'O157_H7_sakai_allexternalnucmer': 4.53192608862\n")
        fp.write("Distance between 'E2348_69_allexternalnucmer' and 'Reference': 1.29611949657")
        fp.close()
        self.assertEqual(find_two_new(fpath, ["E2348_69_allexternalnucmer"]),((('E2348_69_allexternalnucmer', 'O157_H7_sakai_allexternalnucmer', '4.53192608862'), ('E2348_69_allexternalnucmer', 'H10407_allexternalnucmer', '1.39030683167')),(('Reference', 'E2348_69_allexternalnucmer', '1.29611949657'),)))
        os.chdir("%s" % curr_dir)
        shutil.rmtree(tdir)
    def test_find_two_reversed(self):
        tdir = tempfile.mkdtemp(prefix="reversed_",)
        fpath = os.path.join(tdir,"all_patristic_distances_reversed.txt")
        os.chdir("%s" % tdir)
        fp = open(fpath,"w")
        fp.write("Distance between 'E2348_69_allexternalnucmer' and 'H10407_allexternalnucmer': 1.39030683167\n")
        fp.write("Distance between 'E2348_69_allexternalnucmer' and 'O157_H7_sakai_allexternalnucmer': 4.53192608862\n")
        fp.write("Distance between 'Reference' and 'E2348_69_allexternalnucmer': 0.0941892612547")
        fp.close()
        self.assertEqual(find_two_new(fpath, ["E2348_69_allexternalnucmer"]),((('E2348_69_allexternalnucmer', 'O157_H7_sakai_allexternalnucmer', '4.53192608862'), ('E2348_69_allexternalnucmer', 'H10407_allexternalnucmer', '1.39030683167')),(('Reference', 'E2348_69_allexternalnucmer', '0.0941892612547'),)))
        os.chdir("%s" % curr_dir)
        shutil.rmtree(tdir)
class Test12(unittest.TestCase):
    def test_get_closest_dists_basic_function(self):
        self.assertEqual(get_closest_dists_new((('ECOLI_ISO2', 'SSON_046_allexternalnucmer', '0.08198920048'), ('ECOLI_ISO2', 'H10407_allexternalnucmer', '1.3087194675e-06')),['ECOLI', 'ECOLI_IS03_L007', 'ECOLI_ISO2']),(['SSON_046_allexternalnucmer0.08198920048', 'H10407_allexternalnucmer1.3087194675e-06']))

class Test13(unittest.TestCase):
    def test_calculate_pairwise_tree_dists_basic_function(self):
        """distances were taken directly from Dendropy, run outside of the pipeline"""
        tdir = tempfile.mkdtemp(prefix="filetest_",)
        fpath = os.path.join(tdir,"test.tree")
        os.chdir("%s" % tdir)
        fp = open(fpath,"w")
        fp.write("((H10407_all:0.00000095154411648831,SSON_046_all:0.08004748973386473232):0.09609983129754622044,(E2348_69_all:1.63492542157114217893,O157_H7_sakai_all:4.52711175943011490119):0.00000095154411648831,Reference:0.00000095154411648831):0.0;")
        fp.close()
        self.assertEqual(calculate_pairwise_tree_dists(fpath, "tmp.out"),[0.08004844127798122, 1.7310271559569212, 4.623213493815894, 0.0961017343857792, 1.8110736941466694, 4.703260032005642, 0.17614827257552745, 6.162037181001257, 1.634927324659375, 4.527113662518348])
        os.chdir("%s" % curr_dir)
        shutil.rmtree(tdir)

class Test14(unittest.TestCase):
    def test_subsample_snps(self):
        tdir = tempfile.mkdtemp(prefix="filetest_",)
        fpath = os.path.join(tdir,"test.matrix")
        os.chdir("%s" % tdir)
        fp = open(fpath,"w")
        fp.write("LocusID\tReference\tgenome1\tgenome2\n")
        fp.write("ADK::1\tA\tT\tT\n")
        fp.write("ADK::2\tT\tT\tT\n")
        fp.write("ADK::3\tG\tG\tT\n")
        fp.close()
        self.assertEqual(subsample_snps(fpath,{'ECOLI': ['SSON_046_allexternalnucmer', 'H10407_allexternalnucmer', 'E2348_69_allexternalnucmer', 'O157_H7_sakai_allexternalnucmer'], 'ECOLI_ISO2': ['H10407_allexternalnucmer', 'SSON_046_allexternalnucmer', 'E2348_69_allexternalnucmer', 'O157_H7_sakai_allexternalnucmer'], 'ECOLI_IS03_L007': ['SSON_046_allexternalnucmer', 'H10407_allexternalnucmer', 'E2348_69_allexternalnucmer', 'O157_H7_sakai_allexternalnucmer']}, {'ECOLI': 11, 'ECOLI_ISO2': 14, 'ECOLI_IS03_L007': 11},4),['ADK::1','ADK::2','ADK::3'])
        os.chdir("%s" % curr_dir)
        shutil.rmtree(tdir)

class Test15(unittest.TestCase):
    def test_process_temp_matrices_dev(self):
        tdir = tempfile.mkdtemp(prefix="filetest_",)
        fpath = os.path.join(tdir, "tmp.tree")
        fpath2 = os.path.join(tdir, "all.distances")
        os.chdir("%s" % tdir)
        fp = open(fpath,"w")
        fp2 = open(fpath2, "w")
        fp.write("(E2348_69_allexternalnucmer:1.29611826169204280568,O157_H7_sakai_allexternalnucmer:3.23580782692392832089,(Reference:0.00000061743893499046,((QUERY___ECOLI_ISO2:0.000001,H10407_allexternalnucmer:0.00000030871946749523):0.00000030871946749523,((QUERY___ECOLI:0.000001,QUERY___ECOLI_IS03_L007:0.000001):0.0,SSON_046_allexternalnucmer:0.04099394588025907088):0.04099394588025907088):0.09418733509629113876):0.00000061743893499046);")
        fp.close()
        fp2.write("Distance between 'E2348_69_allexternalnucmer' and 'O157_H7_sakai_allexternalnucmer': 4.53192608862\n")
        fp2.write("Distance between 'E2348_69_allexternalnucmer' and 'Reference': 1.29611949657\n")
        fp2.write("Distance between 'E2348_69_allexternalnucmer' and 'QUERY___ECOLI_ISO2': 1.39030752295\n")
        fp2.write("Distance between 'E2348_69_allexternalnucmer' and 'H10407_allexternalnucmer': 1.39030683167\n")
        fp2.write("Distance between 'E2348_69_allexternalnucmer' and 'QUERY___ECOLI': 1.43130116011\n")
        fp2.write("Distance between 'E2348_69_allexternalnucmer' and 'QUERY___ECOLI_IS03_L007': 1.43130116011\n")
        fp2.write("Distance between 'E2348_69_allexternalnucmer' and 'SSON_046_allexternalnucmer': 1.47229410599\n")
        fp2.write("Distance between 'O157_H7_sakai_allexternalnucmer' and 'Reference': 3.2358090618\n")
        fp2.write("Distance between 'O157_H7_sakai_allexternalnucmer' and 'QUERY___ECOLI_ISO2': 3.32999708818\n")
        fp2.write("Distance between 'O157_H7_sakai_allexternalnucmer' and 'H10407_allexternalnucmer': 3.3299963969\n")
        fp2.write("Distance between 'O157_H7_sakai_allexternalnucmer' and 'QUERY___ECOLI': 3.37099072534\n")
        fp2.write("Distance between 'O157_H7_sakai_allexternalnucmer' and 'QUERY___ECOLI_IS03_L007': 3.37099072534\n")
        fp2.write("Distance between 'O157_H7_sakai_allexternalnucmer' and 'SSON_046_allexternalnucmer': 3.41198367122\n")
        fp2.write("Distance between 'Reference' and 'QUERY___ECOLI_ISO2': 0.0941892612547\n")
        fp2.write("Distance between 'Reference' and 'H10407_allexternalnucmer': 0.0941885699742\n")
        fp2.write("Distance between 'Reference' and 'QUERY___ECOLI': 0.135182898415\n")
        fp2.write("Distance between 'Reference' and 'QUERY___ECOLI_IS03_L007': 0.135182898415\n")
        fp2.write("Distance between 'Reference' and 'SSON_046_allexternalnucmer': 0.176175844296\n")
        fp2.write("Distance between 'QUERY___ECOLI_ISO2' and 'H10407_allexternalnucmer': 1.3087194675e-06\n")
        fp2.write("Distance between 'QUERY___ECOLI_ISO2' and 'QUERY___ECOLI': 0.0409962545997\n")
        fp2.write("Distance between 'QUERY___ECOLI_ISO2' and 'QUERY___ECOLI_IS03_L007': 0.0409962545997\n")
        fp2.write("Distance between 'QUERY___ECOLI_ISO2' and 'SSON_046_allexternalnucmer': 0.08198920048\n")
        fp2.write("Distance between 'H10407_allexternalnucmer' and 'QUERY___ECOLI': 0.0409955633192\n")
        fp2.write("Distance between 'H10407_allexternalnucmer' and 'QUERY___ECOLI_IS03_L007': 0.0409955633192\n")
        fp2.write("Distance between 'H10407_allexternalnucmer' and 'SSON_046_allexternalnucmer': 0.0819885091995\n")
        fp2.write("Distance between 'QUERY___ECOLI' and 'QUERY___ECOLI_IS03_L007': 2e-06\n")
        fp2.write("Distance between 'QUERY___ECOLI' and 'SSON_046_allexternalnucmer': 0.0409949458803\n")
        fp2.write("Distance between 'QUERY___ECOLI_IS03_L007' and 'SSON_046_allexternalnucmer': 0.0409949458803")
        fp2.close()
        #self.assertEqual(process_temp_matrices(({'ECOLI': ['SSON_046_allexternalnucmer', 'H10407_allexternalnucmer', 'E2348_69_allexternalnucmer', 'O157_H7_sakai_allexternalnucmer'], 'ECOLI_ISO2': ['H10407_allexternalnucmer', 'SSON_046_allexternalnucmer', 'E2348_69_allexternalnucmer', 'O157_H7_sakai_allexternalnucmer'], 'ECOLI_IS03_L007': ['SSON_046_allexternalnucmer', 'H10407_allexternalnucmer', 'E2348_69_allexternalnucmer', 'O157_H7_sakai_allexternalnucmer']},fpath,2,fpath2),(

class Test16(unittest.TestCase):
    def test_get_all_snps(self):
        tdir = tempfile.mkdtemp(prefix="filetest_",)
        fpath = os.path.join(tdir,"testfile.filtered.vcf")
        fp = open(fpath, "w")
        fp.write("LocusID\tReference\tgenome1\tgenome2\n")
        fp.write("ADK::1\tA\tT\tT\n")
        fp.write("ADK::2\tT\tT\n")
        fp.close()
        self.assertEqual(get_all_snps(fpath), ['ADK::1','ADK::2'])
        shutil.rmtree(tdir)

class Test17(unittest.TestCase):
    def test_find_used_snps_basic(self):
        tdir = tempfile.mkdtemp(prefix="filetest_",)
        fpath = os.path.join(tdir,"testfile.filtered.vcf")
        fp = open(fpath, "w")
        fp.write("gi|16120353|ref|NC_003143.1|::35501\tG\n")
        fp.write("gi|16120353|ref|NC_003143.1|::52924\tC\n")
        fp.write("gi|16120353|ref|NC_003143.1|::55551\tN\n")
        fp.close()
        os.chdir("%s" % tdir)
        self.assertEqual(find_used_snps(),{'testfile':2})
        os.chdir("%s" % curr_dir)
        shutil.rmtree(tdir)
    def test_find_used_snps_all_missing(self):
        tdir = tempfile.mkdtemp(prefix="filetest_",)
        fpath = os.path.join(tdir,"testfile.filtered.vcf")
        fp = open(fpath, "w")
        fp.write("gi|16120353|ref|NC_003143.1|::35501\tN\n")
        fp.write("gi|16120353|ref|NC_003143.1|::52924\tN\n")
        fp.write("gi|16120353|ref|NC_003143.1|::55551\tN\n")
        fp.close()
        os.chdir("%s" % tdir)
        self.assertEqual(find_used_snps(),{'testfile':0})
        os.chdir("%s" % curr_dir)
        shutil.rmtree(tdir)
    def test_find_used_snps_blank_line(self):
        tdir = tempfile.mkdtemp(prefix="filetest_",)
        fpath = os.path.join(tdir,"testfile.filtered.vcf")
        fp = open(fpath, "w")
        fp.write("gi|16120353|ref|NC_003143.1|::35501\tG\n")
        fp.write("gi|16120353|ref|NC_003143.1|::52924\tC\n")
        fp.write("gi|16120353|ref|NC_003143.1|::55551\t-\n")
        fp.write(" n")
        fp.close()
        os.chdir("%s" % tdir)
        self.assertRaises(TypeError, find_used_snps, ())
        os.chdir("%s" % curr_dir)
        shutil.rmtree(tdir)

class Test18(unittest.TestCase):
    def test_branch_lengths_to_decimals(self):
        self.assertEqual(branch_lengths_2_decimals("('E2348_69_allexternalnucmer':1.63492542157,'O157_H7_sakai_allexternalnucmer':4.52711175943):9.51544116488e-07,Reference:9.51544116488e-07)"),("('E2348_69_allexternalnucmer':1.634925,'O157_H7_sakai_allexternalnucmer':4.527112):0.000001,Reference:0.000001);"))

class Test19(unittest.TestCase):
    def test_parse_likelihoods(self):
        tdir = tempfile.mkdtemp(prefix="filetest_",)
        fpath = os.path.join(tdir,"testfile.filtered.vcf")
        fp = open(fpath, "w")
        fp.write("ECOLI\tI6\t0.298716\t0.298716\n")
        fp.write("ECOLI\tI5\t0.298714\t0.597430\n")
        fp.write("ECOLI\tI4\t0.298714\t0.896143\n")
        fp.close()
        self.assertEqual(parse_likelihoods(fpath), {'ECOLI':['0.298716','0.298714','0.298714']})
        shutil.rmtree(tdir)

class Test20(unittest.TestCase):
    def test_fasta_to_tab(self):
        tdir = tempfile.mkdtemp(prefix="filetest_",)
        fpath = os.path.join(tdir,"testfile.fasta")
        fp = open(fpath, "w")
        fp.write(">id\n")
        fp.write("ATGC")
        fp.close()
        self.assertEqual(fasta_to_tab(fpath), ["id","ATGC"])
        shutil.rmtree(tdir)

class Test21(unittest.TestCase):
    def test_tab_to_fasta(self):
        tdir = tempfile.mkdtemp(prefix="filetest_",)
        fpath = os.path.join(tdir,"testfile.tab")
        fp = open(fpath, "w")
        fp.write(">id\tATCG\n")
        fp.close()
        self.assertEqual(tab_to_fasta(fpath), [">id","ATCG"])
        shutil.rmtree(tdir)

class Test22(unittest.TestCase):
    def test_tab_to_matrix(self):
        tdir = tempfile.mkdtemp(prefix="filetest_",)
        fpath = os.path.join(tdir,"testfile.tab")
        fp = open(fpath, "w")
        fp.write("id\tATCG\n")
        fp.write("id2\tTTCG\n")
        fp.close()
        self.assertEqual(tab_to_matrix(fpath), [['id', 'id2'], ['A', 'T'], ['T', 'T'], ['C', 'C'], ['G', 'G']])
        shutil.rmtree(tdir)

class Test23(unittest.TestCase):
    def test_prune_fasta(self):
        tdir = tempfile.mkdtemp(prefix="filetest_",)
        fpath = os.path.join(tdir,"testfile.tab")
        fp = open(fpath, "w")
        fp.write(">id\nATCG\n")
        fp.write(">id2\nTTCG\n")
        fp.close()
        self.assertEqual(prune_fasta("id",fpath,"%s/outfile.txt" % tdir), ['id2'])
        shutil.rmtree(tdir)

class Test24(unittest.TestCase):
    def test_filter_alignment(self):
        """tests to make sure that this function
        throws out all characters that it is supposed to"""
        tdir = tempfile.mkdtemp(prefix="filetest_",)
        fpath = os.path.join(tdir,"testfile.tab")
        fp = open(fpath, "w")
        fp.write("throw away line\n")
        fp.write("X\tA\tN\tX\t-\n")
        fp.close()
        self.assertEqual(filter_alignment(fpath), [['A']])
        shutil.rmtree(tdir)

class Test25(unittest.TestCase):
    def test_compare_subsample_results(self):
        tdir = tempfile.mkdtemp(prefix="filetest_",)
        os.chdir("%s" % tdir)
        fpath = os.path.join(tdir,"sample..testfile.tab")
        fp = open(fpath, "w")
        fp.write("Cluster0	Cluster0	100.00	15	0	0	1	15	1	15	1e-07	30.2")
        fp.close()
        os.chdir("%s" % curr_dir)
        shutil.rmtree(tdir)
        os.system("rm tab_matrix tab.filtered out.tab out.fasta")

if __name__ == "__main__":
    unittest.main()
    main()
<file_sep>/setup.py
#!/usr/bin/env python
from distutils.core import setup
import re

__author__ = "<NAME>"
__credits__ = ["<NAME>"]
__license__ = "GPL v3"
__version__ = "1.0"
__maintainer__ = "<NAME>"
__email__ = "<EMAIL>"
__status__ = "Development"
 
long_description = """WG-FAST, genotype
phylogenetically from a NASP formatted SNP
matrix
"""

setup(name='wg_fast',
      version=__version__,
      description='whole genome focused array SNP typing',
      author=__maintainer__,
      author_email=__email__,
      maintainer=__maintainer__,
      maintainer_email=__email__,
      packages=['wg_fast'],
      long_description=long_description
)
<file_sep>/wgfast_manual.md
## The whole genome focused array SNP typing (*WG-FAST*) pipeline  

### Updated 7/16/2019  

#### Citation:  
<NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME> (2015).
Phylogenetically typing bacterial strains from partial SNP genotypes observed from direct
sequencing of clinical specimen metagenomic data  

#### Contact:  
Please address queries, conserns, improvements to jasonsahl at gmail dot com  

#### What does *WG-FAST* do?  

WG-FAST was designed as a tool to phylogenetically genotype unknown samples, even those
with extremely low read coverage, in the context of a well-studied dataset.

#### What does *WG-FAST* not do?  

*WG-FAST* is not intended to identify new SNPs in a dataset (i.e replace a de novo analysis).
If too many samples are processed with WG-FAST, a phylogenetic discovery bias can exist.

#### Installation  
See Readme

#### Unit tests  
-To test that all functions are working correctly in *WG-FAST*, type:  
```python tests/test_all_functions.py```

#### Required input files  
1. **Directory of sequence reads ("-d")**. The reads must be named according to Illumina HiSeq
or MiSeq conventions. Reads must be in the Illumina 1.9+ FastQ format. If you have old
Illumina FASTQ encodings, they must be converted before running *WG-FAST*. **Important: names
must not have periods ".", brackets "[]", or other weird characters "=:" in the header.**  
2. **Directory of reference files ("-r")**. This directory should only contain the following
files:  
a. **SNP matrix (must end in ".tsv")**. The easiest way to generate this is by using NASP
(https://github.com/TGenNorth/NASP). If other SNP matrix formats are used, they must conform
to having hte first column including (contig::coordinate) and the column following the SNP calls
must be (#SNPcall). For the sub-sampling routine to complete, a genome must be present in your matrix
that is called 'Reference'. **Important: sample names must not have periods in the header.**  
b. **Phylogeny (must end in .tree)**. A script is included with *WG-FAST* that can generate an appropriate
phylogeny for a NASP matrix (see below). This script also generates a 'Parameters' file, which can be used
with *WG-FAST* and cuts down on the computational time required for each subsequent run. The best way to guarantee downstream
compatibility is to generate the phylogeny with RAxML, which is used by the wgfast_prep.py script described below.  
c. **Reference genome in FASTA format (must end in .fasta)**. This should be the same FASTA that was used to call SNPs with NASP.  
d. **RAxML paramters file (must end in .PARAMS) (optional)**. If provided, RAxML will run faster because the likelihood has already
been calculated. Make sure that you use the conda version of RAxML for full compatibility.  

#### Complete list of arguments:  

-h, --help  show this help message and exist  
-r, REFERENCE_DIR, path to reference directory [REQUIRED]  
-d, READ_DIRECTORY, path to directory of fastq.gz files [REQUIRED]  
-p, PROCESSORS, # of processors to use - defaults to 2  
-c, COVERAGE, minimum SNP coverage required to be called, defaults to 3x  
-o, PROPORTION, proportion of alleles to be corrected, defaults to 0.9  
-k, KEEP, keep temporary files? Defaults to F (T or F)  
-s, SUBSAMPLE, run subsample routine? Defaults to T (T or F)  
-n, SUBNUMS, number of subsamples to process, defaults to 100  
-g, DOC, run depth of coverage on all files? Defaults to T (T or F)  
-e, TMP_DIR, temporary directory path for GATK, default = '/tmp'  
-z, insertion method (MP or ML), defaults to ML, MP not well tested  
-f, How close does a subsample have to be from true placement, Defaults to 0.1 (Float)  
-y, Only run sub-sample routine and exit? Defaults to F (T or F)  
-j, MODEL, which model to run with raxml (GTRGAMMA, ASC_GTRGAMMA)  

#### Test Data:  

1. To give *WG-FAST* a try to make sure everything is installed correctly, check out the test_data
directory.  
2. The following command was run from a “run” directory created within the WG-FAST
installation:  
```python ../wgfast.py -r ../test_data -d ../test_data/reads/ -p 4 -c 1 -s F```  

#### What is *WG-FAST* doing?  

*WG-FAST* is a pipeline that wraps other tools together to place the samples into a phylogenetic
context, based on a well-characterized dataset. The order of functions conducted by *WG-FAST*
include:  
1. Adapters are trimmed using BBduk. A list of adapters is included in the “bin” directory
(“illumina_adapters_all.txt”) within the *WG-FAST* distribution. This file includes all of the
standard Illumina adapters that I could find, plus additional ones that we use in our
laboratory. If you have different adapters than these, add the sequences into this file when
you run *WG-FAST*. The minimum length of sequences to keep is hard coded as 50, but I
could change this value to reflect different lengths if needed.  
2. A dictionary is created from the Reference fasta with Picard Tools and the reference fasta is
indexed with SAMtools.  
3. Reads are paired into a single sample based on their names, although single end reads are
also supported. *WG-FAST* assumes that names be something like: "S1_R1_001.fastq.gz” or
“R1_001.fastq.gz”. If the name pairings aren’t recognized, they will be run as single ended.  
4. Reads are mapped to the reference with minimap2, using default settings. SNPs are then
called with the Haplotype caller method in GATK v4.  
5. If selected, coverage across the reference is calculated with the “DepthOfCoverage” method
in SAMtools.  
6. SNPs are compared against the SNPs in your initial SNP matrix. If the SNP is missing from
your VCF file, a “-“ is inserted. Furthermore, if a position fails either the user-defined
thresholds for depth or proportion, the position is replaced with a “-“. The number of
discarded SNPs, or those that didn’t pass the filters, is reported.  
7. The new matrix is converted into a multi-fasta and processed with the evolutionary placement
algorithm in RAxML v8, placing the unknown into a user-provided phylogeny. The resulting
tree is converted into Nexus, such that the unknowns can be easily visualized (Red) in
FigTree (http://tree.bio.ed.ac.uk/software/figtree/).  
8. For the optional sub-sampling routine, the following functions are performed:  
**a.** For each query, the two closest genomes are identified, based on lowest pairwise
patristic distances calculated by DendroPy  
**b.** Two matrices are created for each query, one for each of the nearest neighbors. Each
of these contains the name of the query and the name of the neighbor, ending in
“tmp.matrix”. If you generate these files (and keep them), then run *WG-FAST* in “-y T”
mode, these files will be skipped and will not need to be created again.  
**c.** Each matrix is then converted into a FASTA file and a parameters file is created with
RAxML. The thought process here is that once the parameters files are created, they
can be used for future comparisons. If you run *WG-FAST* with “-y T”, you can re-use
these parameters files and won’t need to generate again. This method uses the
PTHREADS method to take advantage of multiple processors. Currently, only four
separate parameters files can be generated concurrently; this value is hardcoded into
the script. This is due to the large computational demand required by this method, but
could be made into a tunable parameter, if needed.  
**d.** The initial tree is pruned of the genome to be re-inserted. Each sub-sampled neighbor
is then inserted back into the pruned phylogeny, using the GTRGAMMA method, to
stay consistent with the parameters file. The patristic distance is then calculated
between the neighbor and the “Reference”. This is done because if a NASP-formatted
matrix is generated, the reference FASTA file will always be called “Reference”. The
“true” distance is also identified based on the original SNP matrix using all of the
original source data. The sub-sampled distance is then divided by the true distance in
order to get a ratio of placement replication. If this value falls within a user-defined
threshold, then the placement is considered to be correct. For each unknown, the
number of incorrect placements is divided by the total number of placements in order
to get an idea of how stable the placement is in that region of the tree at the level of
SNP density provided by the unknown. The resulting p-value can then be used to
assess the robustness of a placement, based on significance values chosen by the
user.  

Try running the script with the subsample routine with;  
```python ../wgfast.py -r ../test_data -d ../test_data/reads/ -p 4 –c 1```

#### Output printed to screen:  
1. The parameters that you called. Default values are also printed so the run can be
reproduced.  
2. Number of callable positions, including polymorphic and monomorphic positions. These are
all positions called in each sample, compared to the reference. This is the number prior to
any filtering due to mixed SNP positions.  
3. Number of SNPs. These are the number of observed polymorphisms, based on calls made
by GATK.  
4. Number of discarded SNPs. These are polymorphisms that were called by GATK, but were
thrown out because they failed to meet the depth and/or proportion filters.  
5. Insertion likelihood values. The higher the likelihood value and the fewer the number of
possible insertion nodes, the more trusted the placement, although caveats exist (see
Manuscript).  
6. If sub-sample routine is invoked, information is also available for how often the sub-sample
was placed correctly. A sub-sample is considered to be “correct” by comparing the patristic
distance from the sub-sampled genome to the “Reference”, then comparing that to the
distance between the un-sub-sampled genome to the “Reference”. If this ratio falls within the
“fudge-factor” range, then it is considered to be correct. The number of times that the subsample
falls within this range is divided by the total number of iterations and a p-value is
reported.  

### Files generated by *WG-FAST*:  

1. “transformed.tree”. This file is a nexus file of the starting tree with the unknowns placed. If
this file is opened with figtree (http://tree.bio.ed.ac.uk/software/figtree/), the unknown isolates
are shown in red.  
2. “nasp_matrix_with_unknowns.txt”. This file is a modification of the original NASP matrix with
the unknowns inserted, so the user can identify specific SNPs. Information is lost from this
matrix, so it cannot be used with additional *WG-FAST* runs.  
3. “coverage_out.txt”. If you choose to have genome coverage information, the coverage across
the reference genome is provided in this file.  
4. “all_patristic_distances.txt”. This is the patristic distance, or tree path distance, between all
pairwise comparisons in your analysis.  

### Additional scripts included with *WG-FAST*:  

1. wgfast_prep.py  

-What does it do? Given a NASP matrix, this script will generate a maximum likelihood
phylogeny with RAxML and will also generate a “parameters” file that can be used for future
*WG-FAST* runs. The use of a parameters cuts down on the computational requirements when
running additional *WG-FAST* runs using the same input files.  

-What do you need for the script to run? Requirements include a NASP matrix  

-What does the output look like? Two files are produced:  

a. “nasp_raxml.tree”. Your tree. Names have been fixed to work with WG-FAST. Make sure
that names do not include periods.  
b. “nasp.PARAMS”. Parameters file. Use this with the “-X” flag described above.  
c. Example usage:  
```python wgfast_prep.py -m nasp.matrix```  

-Note: if you run this script with a model separate from ASC_GTRGAMMA, you will also need to
use this model when running the main wgfast script, if you use the parameters file. If you have a
large file, you can run the script using the GTRGAMMA model and multiple processors.  

```python wgfast_prep.py -m nasp_matrix.tsv -o GTRGAMMA -p 4```  

2. subsample_snps_pearson.py  

-What does it do? Given a NASP matrix, the script generates a new matrix over a given number
of iterations at a given level of SNP sampling. The script then conducts a Mantel test using the
Pearson correlation between the original similarity matrix and a matrix sampled at a given SNP
density.  

-What do you need for the script to run?
a. NASP nasp_matrix
b. 'mothur executable in your path'. Mothur can be freely obtained from:  
http://www.mothur.org/wiki/Download_mothur  

-What does the output look like? One file is generated “results.txt”, that is new-line delimited,
with each line containing a Pearson correlation value (0 to 1).  

-Example usage:  
```python subsample_snps_pearson.py -m nasp.matrix -s 100```  

3. subsample_reads_and_place.py  

-What does it do? This script helps sub-sample your SNP matrix and identify how robust specific
regions of the phylogeny are. This script will prune the genome from the phylogeny, randomly
subsample the SNPs at a given level n separate times, then will re-insert into the tree. The
distance between the “correct” placement will be compared to the re-sampled placement. This
is a computationally intensive method and works best with a job management system (see
below).  

-What do you need for the script to run? Requirements include:  
a. NASP matrix  
b. Corresponding phylogeny  
c. Name of genome to test  

-What does the output look like? For each genome, you will get an output file, showing each
iteration (1st column), the number of correct placements (2nd column) and the best SNP level, if
one is identified. For example:  
```100  1```  
```200  10```  
```optimal value is for Salmonella_enterica_subsp_enterica_serovar_Uganda_str_R8-3404nucmer is 200```  

-Example usage:  

```python subsample_reads_and_place.py -m matrix.txt -t nasp_raxml.tree -n name -s 100 -o 6 -e 10000```  

<file_sep>/wgfast.py
#!/usr/bin/env python

"""
WG-FAST
written by <NAME>
correspondence: <EMAIL>
"""
from optparse import OptionParser
import subprocess
import os
import sys
import errno
import glob
import tempfile

"""modify line below to reflect your installation directory"""
WGFAST_PATH="/Users/jasonsahl/tools/wgfast"

if os.path.exists(WGFAST_PATH):
    sys.path.append("%s" % WGFAST_PATH)
else:
    print("your WG-FAST path is not correct.  Edit the path in wgfast.py and try again")
    sys.exit()
try:
    from wg_fast.util import *
except:
    print("wgfast path needs to be modified in the wgfast.py file")
    sys.exit()

def main(reference_dir,read_directory,processors,coverage,proportion,keep,subsample,
    subnums,doc,tmp_dir,fudge,only_subs,model,ploidy):
    ref_path=os.path.abspath("%s" % reference_dir)
    dir_path=os.path.abspath("%s" % read_directory)
    """Test to make sure all required files are present"""
    tree = "".join(glob.glob(os.path.join(ref_path, "*.tree")))
    tree_list = glob.glob(os.path.join(ref_path, "*.tree"))
    if len(tree_list) == 0:
        print("You need to provide a tree in your reference directory ending in '.tree'")
        sys.exit()
    elif len(tree_list) == 1:
        pass
    elif len(tree_list) > 1:
        print("More than one tree ending in '.tree' found in your reference directory...exiting")
        sys.exit()
    matrix = "".join(glob.glob(os.path.join(ref_path, "*.tsv")))
    matrix_list = glob.glob(os.path.join(ref_path, "*.tsv"))
    if len(matrix_list) == 0:
        print("You need to provide a NASP formatted matrix ending in '.tsv'")
        sys.exit()
    elif len(matrix_list) > 1:
        print("More than one file in your reference directory ending in '.tsv'...exiting")
        sys.exit()
    reference = "".join(glob.glob(os.path.join(ref_path, "*.fasta")))
    reference_list = glob.glob(os.path.join(ref_path, "*.fasta"))
    if len(reference_list) == 0:
        print("You must provide a REFERENCE FASTA file")
        sys.exit()
    elif len(reference_list) > 1:
        print("More than one file in your reference directory ending in '.fasta'...exiting")
        sys.exit()
    """Get the reference information. Outfile is named $name.tmp.txt"""
    get_seq_length(reference,"ref")
    subprocess.check_call("tr ' ' '\t' < ref.tmp.txt > ref.genome_size.txt", shell=True)
    try:
        parameters = "".join(glob.glob(os.path.join(ref_path, "*.PARAMS")))
        parameters_list = glob.glob(os.path.join(ref_path, "*.PARAMS"))
        if len(parameters_list)>1:
            print("More than one RAxML parameters file found in your reference directory...exiting")
            sys.exit()
        elif len(parameters_list) == 0:
            parameters = "NULL"
    except:
        parameters = "NULL"
    #check for binary dependencies
    logPrint('testing the paths of all dependencies')
    ap=os.path.abspath("%s" % os.getcwd())
    #This is now part of the conda install
    aa = subprocess.call(['which', 'raxmlHPC-SSE3'])
    if aa == 0:
        pass
    else:
        print("RAxML must be in your path as raxmlHPC-SSE3")
        sys.exit()
    print("*citation: '<NAME>. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics (2014).'")
    print("*citation: '<NAME>, <NAME>, <NAME>. Performance, accuracy, and Web server for evolutionary placement of short sequence reads under maximum likelihood. Syst Biol. 2011;60(3):291-302'")
    ab = subprocess.call(['which', 'samtools'])
    if ab == 0:
        pass
    else:
        print("samtools must be in your path")
        sys.exit()
    print("*citation: '<NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, Genome Project Data Processing S. The Sequence Alignment/Map format and SAMtools. Bioinformatics. 2009;25(16):2078-9'")
    ac = subprocess.call(['which','minimap2'])
    if ac == 0:
        pass
    else:
        print("minimap2 must be in your path")
        sys.exit()
    """This is the new test for bbduk.sh"""
    ac = subprocess.call(['which','bbduk.sh'])
    if ac == 0:
        pass
    else:
        print("bbduk need to be in your path as bbduk.sh")
        sys.exit()
    #test for new dependencies
    variant_tools = ['gatk','picard']
    for tool in variant_tools:
        ab = subprocess.call(['which', tool])
        if ab == 0:
            pass
        else:
            print("%s must be in your path" % tool)
            sys.exit()
    print("*citation: 'Li H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXivorg. 2013(arXiv:1303.3997 [q-bio.GN])'")
    print("Patristic distances calculated with DendroPy")
    print("*citation: 'Sukumaran J, Holder MT. DendroPy: a Python library for phylogenetic computing. Bioinformatics. 2010;26(12):1569-71. Epub 2010/04/28. doi: 10.1093/bioinformatics/btq228. PubMed PMID: 20421198'")
    print("Also uses GATK for variant calling")
    print("*citation: '<NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, De<NAME>. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome research. 2010;20(9):1297-303'")
    print("Uses BioPython for FASTA parsing")
    print("*citation :Cock PJ, <NAME>, <NAME>T, <NAME>, <NAME>J, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>J. Biopython: freely available Python tools for computational molecular biology and bioinformatics. Bioinformatics. 2009;25(11):1422-3")
    print("")
    #done checking for dependencies"""
    logPrint('WG-FAST pipeline starting')
    logPrint("WG-FAST was invoked with the following parameters:")
    print("-m %s \\" % "".join(matrix))
    print("-t %s \\" % "".join(tree))
    print("-r %s \\" % "".join(reference))
    print("-d %s \\" % read_directory)
    print("-x %s \\" % "".join(parameters))
    print("-p %s \\" % processors)
    print("-c %s \\" % coverage)
    print("-o %s \\" % proportion)
    print("-k %s \\" % keep)
    print("-s %s \\" % subsample)
    print("-n %s \\" % subnums)
    print("-g %s \\" % doc)
    print("-e %s \\" % tmp_dir)
    print("-f %s \\" % fudge)
    print("-y %s \\" % only_subs)
    print("-j %s \\" % model)
    print("-l %s" % ploidy)
    print("-------------------------")
    #makes a temporary directory
    scratch_dir = tempfile.mkdtemp()
    check_input_files(matrix,reference)
    ########Real work starts here############
    if only_subs == "T":
        pass
    else:
        subprocess.check_call("cp %s %s/reference.fasta" % (reference,scratch_dir), stdout=open(os.devnull, 'wb'),stderr=open(os.devnull, 'wb'),shell=True)
        #index reference file.  GATK appears to do this incorrectly"""
        subprocess.check_call("samtools faidx %s/reference.fasta" % scratch_dir, stdout=open(os.devnull, 'wb'),stderr=open(os.devnull, 'wb'),shell=True)
        subprocess.check_call("picard CreateSequenceDictionary R=%s/reference.fasta" % scratch_dir, stdout=open(os.devnull, 'wb'),stderr=open(os.devnull, 'wb'),shell=True)
    #First checkpoint, not sure this saves much time
    if os.path.isfile("temp.matrix"):
        pass
    else:
        #write reduced matrix with only the SNP data"""
        write_reduced_matrix(matrix)
    ref_name=get_seq_name(reference)
    if only_subs == "T":
        pass
    else:
        fileSets=read_file_sets(dir_path)
        if len(fileSets) == 0:
            print("No usable file sets found...exiting")
            try:
                os.system("rm temp.matrix")
            except:
                pass
            sys.exit()
        else:
            ref_coords = get_all_snps(matrix)
            logPrint("Loop starting")
            print("-------------------------")
            run_loop_dev(fileSets,dir_path,"%s/reference.fasta" % scratch_dir,processors,
            ref_coords,coverage,proportion,matrix,scratch_dir,doc,tmp_dir,WGFAST_PATH,ploidy)
            logPrint("Loop finished")
            print("-------------------------")
    """will subsample based on the number of SNPs reported by the following function"""
    if "T" in doc:
        os.system("cat *breadth.txt > breadth_over_%sx_out.txt" % coverage)
        os.system("cat *sum_cov.txt> coverage_out.txt")
    else:
        pass
    used_snps=find_used_snps()
    #Outnames is required for the sub-sampling routine, even with -y T
    outnames=grab_names()
    for name in outnames:
        for k,v in used_snps.items():
            if name == k:
                logPrint("number of callable positions in genome %s = %s" % (k,v))
    if only_subs == "T":
        try:
            #Starts with a clean slate, to replace with new EPA algorithm
            os.system("rm RAxML*")
        except:
            pass
        pass
    else:
        create_merged_vcf()
        subprocess.check_call("paste temp.matrix merged.vcf > combined.matrix", shell=True)
        matrix_to_fasta("combined.matrix", "all.fasta")
        os.system("mv combined.matrix %s/nasp_matrix.with_unknowns.txt" % ap)
        """this fixes the SNP output to conform with RAxML"""
        os.system("sed 's/://g' all.fasta | sed 's/,//g' > out.fasta")
        #QC step here
        qc_files("out.fasta",tree)
        suffix = run_raxml("out.fasta",tree,"out.classification_results.txt","V",parameters,model,"out")
        transform_tree("%s.tree_including_unknowns_noedges.tree" % suffix)
        print("")
        logPrint("Insertion likelihood values:")
        parse_likelihoods("out.classification_results.txt")
        print("")
        calculate_pairwise_tree_dists("%s.tree_including_unknowns_noedges.tree" % suffix,"all_patristic_distances.txt")
    if subsample=="T":
        aa = subprocess.call(['which','raxmlHPC-PTHREADS-SSE3'])
        if aa == 0:
            pass
        else:
            print("for sub-sample routine, RAxML must be in your path as raxmlHPC-PTHREADS-SSE3")
            sys.exit()
        print("*citation: '<NAME>. RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models. Bioinformatics. 2006;22(21):2688-90'")
        try:
            if os.path.isfile("tmp_patristic_distances.txt"):
                pass
            else:
                os.system("sort -g -k 6 all_patristic_distances.txt | sed 's/://g' > tmp_patristic_distances.txt")
        except:
            print("all_patrisitic_distances.txt must be in your analysis directory!")
            sys.exit()
        final_sets,distances=find_two_new("tmp_patristic_distances.txt", outnames)
        results = get_closest_dists_new(final_sets,outnames)
        logPrint("running subsample routine, forcing GTRGAMMA model")
        """mpshell on this function"""
        allsnps = get_all_snps(matrix)
        subsample_snps_2(final_sets,used_snps,subnums,allsnps,processors,"temp.matrix")
        temp_matrices = glob.glob(os.path.join(ap, "*tmp.matrix"))
        final_matrices = []
        for matrix in temp_matrices:
            final_matrices.append(matrix.replace("%s/" % ap,""))
        sample_sets = {}
        for matrix in final_matrices:
            entries = matrix.split(".")
            if entries[0] in sample_sets:
                sample_sets[entries[0]].append(entries[2])
            else:
                sample_sets[entries[0]]=[entries[2]]
        new_sample_dicts = {}
        for k,v in sample_sets.items():
            uniques = []
            [uniques.append(item) for item in v if item not in uniques]
            new_sample_dicts.update({k:uniques})
        logPrint('creating PARAMS file')
        create_params_files_dev(new_sample_dicts,tree,"temp.matrix",final_sets,processors)
        try:
            """Must make sure that remove previous RAxML files"""
            subprocess.check_call("rm RAxML*", shell=True, stderr=open(os.devnull, 'w'))
        except:
            pass
        """final_matrices does indeed have all of the temp matrices loaded"""
        logPrint('adding unknowns to tree')
        #TODO: enter progress bar instead of printing out lots of text
        process_temp_matrices_2(final_sets,final_matrices,tree,processors,"all_patristic_distances.txt", "V", parameters, model)
        print("-------------------------")
        compare_subsample_results(outnames,distances,fudge)
    else:
        pass
    #Clean up temporary files
    if keep == "T":
        pass
    else:
        for outname in outnames:
            try:
                subprocess.check_call("rm %s* RAxML* temp.matrix all.fasta out*" % outname, shell=True, stderr=open(os.devnull, 'w'))
            except:
                pass
            os.chdir("%s" % ap)
            subprocess.check_call("rm -rf %s" % scratch_dir, shell=True)
    logPrint("all done")

if __name__ == "__main__":
    parser = OptionParser(usage="usage: %prog [options]",version="%prog 1.3")
    parser.add_option("-r", "--reference_directory", dest="reference_dir",
                      help="path to reference file directory [REQUIRED]",
                      action="callback", callback=test_dir, type="string")
    parser.add_option("-d", "--read_directory", dest="read_directory",
                      help="path to directory of fastq files [REQUIRED]",
                      action="callback", callback=test_dir, type="string")
    parser.add_option("-p", "--processors", dest="processors",
                      help="# of processors to use - defaults to 2",
                      default="2", type="int")
    parser.add_option("-c", "--coverage", dest="coverage",
                      help="minimum SNP coverage required to be called a SNP; defaults to 3",
                      default="3", type="int")
    parser.add_option("-o", "--proportion", dest="proportion",
                      help="proportion of alleles to be called a SNP, defaults to 0.9",
                      default="0.9", type="float")
    parser.add_option("-k", "--keep", dest="keep",
                      help="keep temp files?  Defaults to F",
                      action="callback", callback=test_filter, type="string", default="F")
    parser.add_option("-s", "--subsample", dest="subsample",
                      help="Run subsample routine? Defaults to T",
                      action="callback", callback=test_filter, type="string", default="T")
    parser.add_option("-n", "--subnums", dest="subnums",
                      help="number of subsamples to process, defaults to 100",
                      action="store", type="int", default="100")
    parser.add_option("-g", "--doc", dest="doc",
                      help="run depth of coverage on all files?  Defaults to T",
                      action="callback", callback=test_filter, type="string", default="T")
    parser.add_option("-e", "--temp", dest="tmp_dir",
                      help="temporary directory for GATK analysis, defaults to /tmp",
                      action="store", type="string", default="/tmp")
    parser.add_option("-f", "--fudge_factor", dest="fudge",
                      help="How close does a subsample have to be from true placement?  Defaults to 0.1",
                      action="store", type="float", default="0.1")
    parser.add_option("-y", "--only_subs", dest="only_subs",
                      help="Only run sub-sample routine and exit? Defaults to F",
                      action="callback", callback=test_filter, type="string", default="F")
    parser.add_option("-j", "--model", dest="model",
                      help="which model to run with raxml:GTRGAMMA,ASC_GTRGAMMA",
                      action="callback", callback=test_models, type="string", default="ASC_GTRGAMMA")
    parser.add_option("-l", "--ploidy", dest="ploidy",
                      help="ploidy to use with GATK, choose from 1 or 2 [DEFAULT]",
                      action="store", type="int", default="2")
    options, args = parser.parse_args()

    mandatories = ["reference_dir","read_directory"]
    for m in mandatories:
        if not options.__dict__[m]:
            print("\nMust provide %s.\n" %m)
            parser.print_help()
            exit(-1)

    main(options.reference_dir,options.read_directory,
         options.processors,options.coverage,options.proportion,options.keep,options.subsample,
         options.subnums,options.doc,options.tmp_dir,options.fudge,
         options.only_subs,options.model,options.ploidy)
<file_sep>/tools/wgfast_prep.py
#!/usr/bin/env python

"""from a NASP matrix, creates the tree
necessary for wg-fast to run"""

from optparse import OptionParser
import subprocess
import os
import sys

def test_file(option, opt_str, value, parser):
    try:
        with open(value): setattr(parser.values, option.dest, value)
    except IOError:
        print('%s file cannot be opened' % option)
        sys.exit()

def test_models(option, opt_str, value, parser):
    if "GTRGAMMA" in value:
        setattr(parser.values, option.dest, value)
    elif "ASC_GTRGAMMA" in value:
        setattr(parser.values, option.dest, value)
    else:
        print("substitution model is not supported")
        sys.exit()

def get_field_index(matrix_in):
    """untested function"""
    with open(matrix_in) as my_matrix:
        firstLine = open(matrix_in).readline().rstrip()
        first_fields = firstLine.split("\t")
        last=first_fields.index("#SNPcall")
    return last

def matrix_to_fasta(matrix_in, last):
    """converts a NASP matrix to fasta format.
    Similar to tested function in main script,
    but slightly different output"""
    reduced = [ ]
    out_fasta = open("all.fasta", "w")
    with open(matrix_in) as my_matrix:
        for line in my_matrix:
            fields = line.split("\t")
            reduced.append(fields[1:last])
    test=map(list, zip(*reduced))
    for x in test:
        out_fasta.write(">"+str(x[0])+"\n")
        out_fasta.write("".join(x[1:])+"\n")
    out_fasta.close()

def main(matrix,model,processors,algorithm):
    """determines whether or not raxml is in your path"""
    if algorithm == "raxml-ng":
        ab = subprocess.call(['which', 'raxml-ng'])
        if ab == 0:
            pass
        else:
            print("RAxML must be in your path as raxml-ng")
            sys.exit()
    elif algorithm == "raxml-HPC":
        ab = subprocess.call(['which', 'raxmlHPC-PTHREADS-SSE3'])
        if ab == 0:
            pass
        else:
            print("RAxML must be in your path as raxmlHPC-PTHREADS-SSE3")
            sys.exit()
    last=get_field_index(matrix)
    matrix_to_fasta(matrix, last)
    #Prep the creation of the FASTA file, removing odd characters
    os.system("sed 's/://g' all.fasta | sed 's/,//g' > out.fasta")
    if model == "ASC_GTRGAMMA":
        subprocess.check_call("raxmlHPC-SSE3 -f d -p 12345 -m %s -s out.fasta -n nasp --asc-corr=lewis --no-bfgs > /dev/null 2>&1" % model, stdout=open(os.devnull, 'wb'),stderr=open(os.devnull, 'wb'),shell=True)
        subprocess.check_call("raxmlHPC-SSE3 -f e -m %s -s out.fasta -t RAxML_bestTree.nasp -n PARAMS --asc-corr=lewis --no-bfgs > /dev/null 2>&1" % model, stdout=open(os.devnull, 'wb'),stderr=open(os.devnull, 'wb'),shell=True)
    else:
        if algorithm == "raxml-HPC":
            subprocess.check_call("raxmlHPC-PTHREADS-SSE3 -T %s -f d -p 12345 -m %s -s out.fasta -n nasp --no-bfgs > /dev/null 2>&1" % (processors,model), stdout=open(os.devnull, 'wb'),stderr=open(os.devnull, 'wb'),shell=True)
            subprocess.check_call("raxmlHPC-PTHREADS-SSE3 -T %s -f e -m %s -s out.fasta -t RAxML_bestTree.nasp -n PARAMS --no-bfgs > /dev/null 2>&1" % (processors,model), stdout=open(os.devnull, 'wb'),stderr=open(os.devnull, 'wb'),shell=True)
        elif algorithm == "raxml-ng":
            subprocess.check_call("raxml-ng --msa out.fasta --model GTR+G --threads %s --prefix nasp" % processors,stdout=open(os.devnull, 'wb'),stderr=open(os.devnull, 'wb'),shell=True)
    if algorithm == "raxml-HPC":
        subprocess.check_call("mv RAxML_bestTree.nasp nasp_raxml.tree", shell=True)
        subprocess.check_call("mv RAxML_binaryModelParameters.PARAMS nasp.PARAMS", shell=True)
        subprocess.check_call("rm RAxML_* out.fasta all.fasta", shell=True)
    else:
        subprocess.check_call("mv nasp.raxml.bestTree nasp_raxml.tree", stdout=open(os.devnull, 'wb'),stderr=open(os.devnull, 'wb'),shell=True)
        subprocess.check_call("rm nasp.raxml.startTree out.fasta all.fasta", stdout=open(os.devnull, 'wb'),stderr=open(os.devnull, 'wb'),shell=True)
    print("Model used: %s" % model)

if __name__ == "__main__":
    usage="usage: %prog [options]"
    parser = OptionParser(usage=usage)
    parser.add_option("-m", "--snp_matrix", dest="matrix",
                      help="path to NASP snp_matrix [REQUIRED]",
                      action="callback", callback=test_file, type="string")
    parser.add_option("-o", "--model", dest="model",
                      help="model for RAxML, can be ASC_GTRGAMMA or GTRGAMMA [DEFAULT]",
                      action="callback", callback=test_models, type="string", default="GTRGAMMA")
    parser.add_option("-p", "--processors", dest="processors",
                      help="number of processors to use with GTRGAMMA, defaults to 2",
                      action="store", type="int", default="2")
    parser.add_option("-a", "--algorithm", dest="algorithm",
                      help="algorithm to use, either raxml-ng or raxml-HPC",
                      action="store", type="string", default="raxml-ng")
    options, args = parser.parse_args()

    mandatories = ["matrix"]
    for m in mandatories:
        if not options.__dict__[m]:
            print("\nMust provide %s.\n" %m)
            parser.print_help()
            exit(-1)

    main(options.matrix,options.model,options.processors,options.algorithm)
<file_sep>/wg_fast/__init__.py
#!/usr/bin/env python

__author__ = "<NAME>"
__credits__ = ["<NAME>"]
__license__ = "GPL"
__version__ = "3"
__maintainer__ = "<NAME>"
__email__ = "<EMAIL>"
__status__ = "Development"

__all__ = ['util']
