# HG changeset patch
# User jaredgk
# Date 1539641747 14400
# Node ID 3830d29fca6a936295e21a9b7048c7ab31f7d943

Uploaded

diff -r 000000000000 -r 3830d29fca6a bcftools.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/bcftools.py	Mon Oct 15 18:15:47 2018 -0400
@@ -0,0 +1,131 @@
+import os
+import sys
+import logging
+import subprocess
+
+sys.path.insert(0, os.path.abspath(os.path.join(os.pardir,'jared')))
+
+from vcf_reader_func import checkFormat
+
+def check_bcftools_for_errors (bcftools_stderr):
+    '''
+        Checks the bgzip stderr for errors
+
+        Parameters
+        ----------
+        bcftools_stderr : str
+            bcftools stderr
+
+        Raises
+        ------
+        IOError
+            If bcftools stderr returns an error
+    '''
+
+    # Expand as errors are discovered
+    if bcftools_stderr:
+        logging.error(vcftools_stderr)
+        raise Exception(vcftools_stderr)
+
+def call_bcftools (bcftools_call_args):
+
+    # bcftools subprocess call
+    bcftools_call = subprocess.Popen(['bcftools'] + list(map(str, bcftools_call_args)), stdout=subprocess.PIPE, stderr=subprocess.PIPE)
+
+    # Wait for bcftools to finish
+    bcftools_out, bcftools_err = bcftools_call.communicate()
+
+    check_bcftools_for_errors(bcftools_err)
+
+    logging.info('bcftools call complete')
+
+def check_for_index (filename):
+
+    # Assign the file format
+    file_format = checkFormat(filename)
+
+    # Check if the file to be indexed is a vcf.gz
+    if file_format == 'bgzip':
+        # Check if the index (.tbi) exists
+        if os.path.isfile(filename + '.tbi'):
+            return True
+
+    # Check if the file to be indexed is a bcf
+    elif file_format == 'bcf':
+        # Check if the index (.csi) exists
+        if os.path.isfile(filename + '.csi'):
+            return True
+
+    # Return false if no index is found
+    return False
+
+def create_index (filename):
+
+    # Assign the file format
+    file_format = checkFormat(filename)
+
+    # Check if the file to be indexed is a vcf.gz
+    if file_format == 'bgzip':
+        # Create a index (.tbi)
+        call_bcftools(['index', '-t', filename])
+
+    # Check if the file to be indexed is a bcf
+    elif file_format == 'bcf':
+        # Create a index (.csi)
+        call_bcftools(['index', '-c', filename])
+
+    # Report if file cannot be indexed
+    else:
+        raise Exception('Error creating index for: %s. Only .bcf and .vcf.gz (bgzip) files are supported.' % filename)
+
+def convert_to_bcf (filename, output_prefix):
+
+    # Holds the arguments to convert to BCF format
+    convert_args = ['convert', '-O', 'b']
+
+    # Stores the specified output_prefix to the BCF file
+    bcf_output = '%s.bcf' % output_prefix
+
+    # Assigns the output file to the arguments
+    convert_args.extend(['-o', bcf_output])
+
+    # Assigns the specified input to the arguments
+    convert_args.append(filename)
+
+    # Call bcftools
+    call_bcftools(convert_args)
+
+
+def convert_to_vcf (filename, output_prefix):
+
+    # Holds the arguments to convert to VCF format
+    convert_args = ['view', '-O', 'v']
+
+    # Stores the specified output_prefix to the VCF file
+    vcf_output = '%s.vcf' % output_prefix
+
+    # Assigns the output file to the arguments
+    convert_args.extend(['-o', vcf_output])
+
+    # Assigns the specified input to the arguments
+    convert_args.append(filename)
+
+    # Call bcftools
+    call_bcftools(convert_args)
+
+def convert_to_vcfgz (filename, output_prefix):
+
+    # Holds the arguments to convert to VCFGZ format
+    convert_args = ['view', '-O', 'z']
+
+    # Stores the specified output_prefix to the VCFGZ file
+    vcfgz_output = '%s.vcf.gz' % output_prefix
+
+    # Assigns the output file to the arguments
+    convert_args.extend(['-o', vcfgz_output])
+
+    # Assigns the specified input to the arguments
+    convert_args.append(filename)
+
+    # Call bcftools
+    call_bcftools(convert_args)
diff -r 000000000000 -r 3830d29fca6a beagle.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/beagle.py	Mon Oct 15 18:15:47 2018 -0400
@@ -0,0 +1,163 @@
+import os
+import sys
+import subprocess
+import shutil
+import argparse
+import glob
+import logging
+
+sys.path.insert(0, os.path.abspath(os.path.join(os.pardir, 'jared')))
+
+from vcf_reader_func import checkFormat
+from logging_module import initLogger, logArgs
+from vcftools import bgzip_decompress_vcfgz
+from bcftools import convert_to_bcf, check_for_index, create_index
+
+def delete_beagle_log (output_prefix):
+    '''
+        Delete beagle log file
+
+        This function is used to delete beagle's log file if an error is
+        encountered. A warning is produced if the log file cannot be found.
+
+        Parameters
+        ----------
+        output_prefix : str
+            Output file prefix
+    '''
+
+    # Check that log file exists, if not return warning
+    if not os.path.isfile(output_prefix + '.log'):
+        logging.warning('beagle log file %s.log does not exist' % output_prefix)
+    else:
+        os.remove(output_prefix + '.log')
+
+def check_beagle_for_errors (beagle_stderr, output_prefix):
+    '''
+        Checks the beagle stdout for errors
+
+        Parameters
+        ----------
+        beagle_stderr : str
+            beagle stderr
+        output_prefix : str
+            Output file prefix
+
+        Raises
+        ------
+        Exception
+            If beagle stdout returns an error
+    '''
+
+    # Check if beagle completed without an error
+    if not beagle_stderr.strip():
+        pass
+
+    # Print missing data message if that is likely
+    elif 'ERROR: genotype is missing allele separator:' in str(beagle_stderr):
+        # Delete the beagle log file
+        delete_beagle_log(output_prefix)
+
+        # Store reported error
+        error_reported = 'ERROR: genotype is missing allele separator'
+        # Store message for user about error
+        user_message = 'Please confirm the input has no missing data.'
+        # Report on the error
+        raise Exception(error_reported + '\n' + user_message)
+
+    # Print output for beagle if error is detected
+    elif 'ERROR:' in str(beagle_stderr):
+        # Delete the beagle log file
+        delete_beagle_log(output_prefix)
+
+        # Splits log into list of lines
+        beagle_stderr_lines = beagle_stderr.splitlines()
+        # Prints the error(s)
+        raise Exception('\n'.join((output_line for output_line in beagle_stderr_lines if output_line.startswith('ERROR:'))))
+
+    # Print output if not completed and no error found. Unlikely to be used, but included.
+    else:
+        # Delete the beagle log file
+        delete_beagle_log(output_prefix)
+
+        raise Exception(beagle_stderr)
+
+
+def standard_beagle_call (beagle_path, beagle_call_args, output_prefix):
+    '''
+        Calls beagle using subprocess
+
+        This function is used to call beagle under standard conditions. The
+        functions then passes the stderr to check_beagle_for_errors to check
+        for errors.
+
+        Parameters
+        ----------
+        beagle_path : str
+            Path to beagle.jar
+        beagle_call_args : list
+            Argument list for beagle
+        output_prefix : str
+            Output file prefix
+    '''
+
+    # Assign location of beagle jar file
+    beagle_jar = os.path.join(beagle_path, 'beagle.jar')
+
+    # Check that beagle.jar exists
+    if not os.path.isfile(beagle_jar):
+        raise IOError('beagle.jar not found. Path specified: %s' % beagle_path)
+
+    logging.info('beagle phasing parameters assigned')
+
+    # Phasing subprocess call
+    phase_call = subprocess.Popen(['java', '-jar', beagle_jar] + beagle_call_args, stdout = subprocess.PIPE, stderr = subprocess.PIPE)
+    phase_stdout, phase_stderr = phase_call.communicate()
+
+    # Check if code is running in python 3
+    if sys.version_info[0] == 3:
+        # Convert bytes to string
+        phase_stderr = phase_stderr.decode()
+
+    # Check beagle call for errors
+    check_beagle_for_errors(phase_stderr, output_prefix)
+
+    logging.info('beagle phasing complete')
+
+def call_beagle (beagle_path, beagle_call_args, output_prefix, output_format):
+    '''
+        Automates beagle calls
+
+        This function passes the argument list to standard_beagle_call. Once the
+        beagle call has finished, the function will automatically convert the
+        bgzip compressed output of beagle to BCF and VCF, if either format is
+        specified.
+
+        Parameters
+        ----------
+        beagle_path : str
+            Path to beagle.jar
+        beagle_call_args : list
+            Argument list for beagle
+        output_prefix : str
+            Output file prefix
+        output_format : str
+            Output file format
+    '''
+    print beagle_call_args
+    # Standard call to beagle
+    standard_beagle_call(beagle_path, beagle_call_args, output_prefix)
+
+    # Decompress if a VCF files is requested
+    if output_format == 'vcf':
+        bgzip_decompress_vcfgz(output_prefix + '.vcf.gz')
+
+    # Convert to BCF if requested
+    elif output_format == 'bcf':
+
+        # Check if there is an index file
+        if check_for_index(output_prefix + '.vcf.gz') == False:
+            # Create an index if not found
+            create_index(output_prefix + '.vcf.gz')
+        # Convert vcf.gz to bcf
+        convert_to_bcf(output_prefix + '.vcf.gz', output_prefix)
diff -r 000000000000 -r 3830d29fca6a logging_module.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/logging_module.py	Mon Oct 15 18:15:47 2018 -0400
@@ -0,0 +1,97 @@
+import sys
+import logging
+
+
+
+def initLogger(filename='pipeline.log', filelevel='INFO',
+               streamlevel='WARNING', resetlog=True):
+    """General logger initialization for PPP functions.
+
+    Messages from WARNING level and higher will be logged to file and to
+    stderr. By default, INFO level will also be written to logfile. Both
+    levels are flexible.
+
+    Level groupings:
+        ERROR: Error messages should be generated by an exception call
+        WARNING: Non-terminal behavior that may be unusual (i.e. lists
+            with no values, blank strings)
+        INFO: Every major (user-facing) function should have the following:
+            -Message for function start
+            -List of input arguments and options
+            -Basic sanity checks (dimensions of input data)
+            -Statements before or after major function calls
+            -Message for function end
+        DEBUG: Mainly for developer use/debugging. Generate logs for
+            sub-functions that match INFO level for major functions.
+            Possible care should be used if there are lots of loops called.
+
+    Use: Call with either the individual function (in __name__=="__main__"
+    statement) or in pipeline file.
+
+    Parameters
+    ----------
+    filename : str ("pipeline.log")
+        Name of file that log will be written to
+
+    filelevel : {'INFO','DEBUG','WARNING','ERROR'}
+        Set minimum level of log messages that are written to log file.
+        Note that this acts as a de facto minumum for 'streamlevel' as well.
+
+    streamlevel : {'WARNING','DEBUG','INFO','ERROR'}
+        Set minimum level of log messages that are output to stream.
+
+    resetlog : bool (True)
+        If true, will overwrite logfile when opening. Set to false if log is
+        being initialized multiple times
+
+    Returns
+    -------
+    None
+
+    Exceptions
+    ----------
+    If filelevel or streamlevel are not a valid logging level
+
+    """
+    log_levels = ['DEBUG','INFO','WARNING','ERROR']
+    if filelevel is not None and filelevel.upper() not in log_levels:
+        raise Exception('filelevel value %s is not a valid level' %
+                         filelevel)
+    if streamlevel is not None and streamlevel.upper() not in log_levels:
+        raise Exception('streamlevel value %s is not a valid level' %
+                         streamlevel)
+    fmt_def = "%(asctime)s - %(funcName)s - %(levelname)s: %(message)s"
+    fmt_notime = "%(funcName)s - %(levelname)s: %(message)s"
+    fmtr = logging.Formatter(fmt=fmt_def)
+    fmtr_notime = logging.Formatter(fmt=fmt_notime)
+    filelogger = logging.getLogger()
+    filelogger.setLevel('INFO')
+    if streamlevel is not None:
+        s_handler = logging.StreamHandler()
+        s_handler.setFormatter(fmtr_notime)
+        s_handler.setLevel(streamlevel)
+        filelogger.addHandler(s_handler)
+    logmode = 'a'
+    if resetlog:
+        logmode = 'w'
+    if filelevel is not None:
+        f_handler = logging.FileHandler(filename,mode=logmode)
+        f_handler.setFormatter(fmtr)
+        f_handler.setLevel(filelevel)
+        #filelogger.setLevel(filelevel)
+        filelogger.addHandler(f_handler)
+    #Formats exception messages to be sent to appropriate loggers
+    def exp_handler(etype,val,tb):
+        logging.error("%s" % (val), exc_info=(etype,val,tb))
+
+    sys.excepthook = exp_handler
+
+
+def logArgs(args, func_name=None, print_nones=False):
+    header = "Arguments"
+    if func_name is not None:
+        header+=" for"+func_name
+    for arg in vars(args):
+        val = vars(args)[arg]
+        if val is not None or print_nones:
+            logging.info('Argument %s: %s' % (arg,val))
diff -r 000000000000 -r 3830d29fca6a model.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/model.py	Mon Oct 15 18:15:47 2018 -0400
@@ -0,0 +1,236 @@
+import os
+import sys
+import json
+import subprocess
+import argparse
+import logging
+import itertools
+
+from collections import defaultdict
+
+# Insert Jared's directory path, required for calling Jared's functions. Change when directory structure changes.
+sys.path.insert(0, os.path.abspath(os.path.join(os.pardir, 'jared')))
+
+from logging_module import initLogger
+
+class ModelFile(dict):
+    def __init__(self, *arg, **kw):
+        super(ModelFile, self).__init__(*arg, **kw)
+        self.inds = []
+        self.ind_file = ''
+        self.exclude_file = ''
+
+    def assign_inds (self, inds = []):
+        # Return error if inds is empty
+        if not inds:
+            raise IOError('No individuals found in the model file.')
+        # Store the individuals
+        self.inds = [str(ind) for ind in inds]
+
+    def create_ind_file (self, file_ext = '', file_path = '', overwrite = False):
+        # Assign the filename for the population file
+        ind_filename = 'unique_individuals' + file_ext
+
+        # If a path is assigned, create the file at the specified location
+        if file_path:
+            ind_filename = os.path.join(file_path, ind_filename)
+
+        # Check if previous files should be overwriten
+        if not overwrite:
+            # Check if the file already exists
+            if os.path.isfile(ind_filename):
+                raise IOError('Individuals file exists.')
+
+        # Create the population file
+        ind_file = open(ind_filename, 'w')
+        ind_file.write('%s\n' %'\n'.join(self.inds))
+        ind_file.close()
+
+        # Save the individuals filename
+        self.ind_file = ind_filename
+
+    def delete_ind_file (self):
+        # Check if an individuals file was created
+        if self.ind_file:
+
+            # Delete the individuals file
+            os.remove(self.ind_file)
+
+            # Remove the filename
+            self.ind_file = ''
+
+    def create_exclude_ind_file (self, inds_to_include = [], file_ext = '', file_path = '', overwrite = False):
+        # Assign the filename for the population file
+        ind_filename = 'exclude_individuals' + file_ext
+
+        # If a path is assigned, create the file at the specified location
+        if file_path:
+            ind_filename = os.path.join(file_path, ind_filename)
+
+        # Check if previous files should be overwriten
+        if not overwrite:
+            # Check if the file already exists
+            if os.path.isfile(ind_filename):
+                raise IOError('Individuals file exists.')
+
+        # Create exclude list by removing included individuals
+        exclude_inds = list(set(self.inds) - set(inds_to_include))
+
+        # Create the population file
+        ind_file = open(ind_filename, 'w')
+        ind_file.write('%s\n' %'\n'.join(exclude_inds))
+        ind_file.close()
+
+        # Save the individuals filename
+        self.exclude_file = ind_filename
+
+    def delete_ind_file (self):
+        # Check if an individuals file was created
+        if self.exclude_file:
+
+            # Delete the individuals file
+            os.remove(self.exclude_file)
+
+            # Remove the filename
+            self.exclude_file = ''
+
+class Model:
+    def __init__ (self, name):
+        self.name = name
+        self.tree = ''
+        self.npop = 0
+        self.pop_list = []
+        self.nind = defaultdict(int)
+        self.ind_dict = defaultdict(list)
+        self.pop_files = []
+        self.ind_file = ''
+
+    @property
+    def inds(self):
+        return list(itertools.chain.from_iterable(self.ind_dict.values()))
+
+    def assign_tree (self, tree):
+        self.tree = str(tree)
+
+    def assign_pop (self, pop, inds = []):
+        self.npop += 1
+        self.pop_list.append(str(pop))
+        if inds:
+            self.nind[pop] = len(inds)
+            self.ind_dict[pop] = [str(ind) for ind in inds]
+
+    def create_pop_files (self, file_ext = '', file_path = '', overwrite = False):
+        for pop in self.pop_list:
+            # Assign the filename for the population file
+            pop_filename = pop + file_ext
+
+            # If a path is assigned, create the file at the specified location
+            if file_path:
+                pop_filename = os.path.join(file_path, pop_filename)
+
+            # Check if previous files should be overwriten
+            if not overwrite:
+                # Check if the file already exists
+                if os.path.isfile(pop_filename):
+                    raise IOError('Population file exists.')
+
+            # Create the population file
+            pop_file = open(pop_filename, 'w')
+            pop_file.write('%s\n' %'\n'.join(self.ind_dict[pop]))
+            pop_file.close()
+
+            # Save the population filename
+            self.pop_files.append(pop_filename)
+
+    def delete_pop_files (self):
+        # Check if pop files were created
+        if len(self.pop_files) != 0:
+
+            # Loop the created pop files
+            for pop_file in self.pop_files:
+                # Delete the pop file
+                os.remove(pop_file)
+
+            # Remove the filenames
+            self.pop_files = []
+
+    def create_ind_file (self, file_ext = '', file_path = '', overwrite = False):
+        # Assign the filename for the population file
+        ind_filename = 'individual.keep' + file_ext
+
+        # If a path is assigned, create the file at the specified location
+        if file_path:
+            ind_filename = os.path.join(file_path, ind_filename)
+
+        # Check if previous files should be overwriten
+        if not overwrite:
+            # Check if the file already exists
+            if os.path.isfile(ind_filename):
+                raise IOError('Individuals file exists.')
+
+        # Create the population file
+        ind_file = open(ind_filename, 'w')
+        ind_file.write('%s\n' %'\n'.join(self.inds))
+        ind_file.close()
+
+        # Save the individuals filename
+        self.ind_file = ind_filename
+
+    def delete_ind_file (self):
+        # Check if an individuals file was created
+        if self.ind_file:
+
+            # Delete the individuals file
+            os.remove(self.ind_file)
+
+            # Remove the filename
+            self.ind_file = ''
+
+def read_model_file (model_filename):
+
+    # Check that the file exists
+    if not os.path.isfile(model_filename):
+        raise IOError
+
+    # Create ModelFile object
+    models_to_return = ModelFile()
+
+    # Check if using python 2 or 3
+    if sys.version_info[0] == 2:
+        # Open the model file in python 2
+        model_file = open(model_filename, 'rU')
+    else:
+        # Open the model file in python 3
+        model_file = open(model_filename, 'r', newline=None)
+
+    # Parse the model file using the json reader
+    models_dict = json.load(model_file)
+
+    # List to store all unique individuals (i.e. individuals in all models)
+    individual_list = []
+
+    # Loop the parsed models
+    for model_dict in models_dict:
+
+        # Create the model
+        model = Model(model_dict['name'])
+
+        # Loop the populations in the model
+        for pop, pop_dict in model_dict['pops'].items():
+
+            # Assign the population ans it's individuals to the model
+            model.assign_pop(pop, pop_dict['inds'])
+            # Assign the individuals to the unique individual list
+            individual_list.extend(pop_dict['inds'])
+
+        # Remove duplicates from the unique individual list
+        individual_list = list(set(individual_list))
+
+        # Save the model
+        models_to_return[str(model.name)] = model
+
+    # Store the  unique individuals within the ModelFile object
+    models_to_return.assign_inds(individual_list)
+
+    # Return the models
+    return models_to_return
diff -r 000000000000 -r 3830d29fca6a shapeit.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/shapeit.py	Mon Oct 15 18:15:47 2018 -0400
@@ -0,0 +1,160 @@
+import os
+import sys
+import subprocess
+import shutil
+import argparse
+import glob
+import logging
+
+sys.path.insert(0, os.path.abspath(os.path.join(os.pardir, 'jared')))
+
+from vcf_reader_func import checkFormat
+from logging_module import initLogger, logArgs
+from plink import convert_haps_to_vcf
+#from vcftools import bgzip_decompress_vcfgz
+#from bcftools import convert_to_bcf, check_for_index, create_index
+
+def check_shapeit_for_errors (shapeit_stdout, output_prefix):
+    '''
+        Checks the shapeit stdout for errors
+
+        Parameters
+        ----------
+        shapeit_stdout : str
+            shapeit stdout
+        output_prefix : str
+            Output filename prefix
+
+        Raises
+        ------
+        Exception
+            If shapeit stdout returns an error
+    '''
+
+    # Returns True if the job completed without error
+    if 'Running time:' in str(shapeit_stdout):
+        pass
+
+    # Print output if not completed and no error found. Unlikely to be used, but included.
+    else:
+        # Remove intermediate files before reporting the error
+        remove_intermediate_files(output_prefix, error_intermediates = True)
+        raise Exception(str(shapeit_stdout))
+
+def remove_intermediate_files (output_prefix, error_intermediates = False):
+    '''
+        Removes shapeit intermediate files
+
+        This function is used to remove the various intermediate files created
+        by shapeit. The exact intermediate files to be removed are defined by
+        the error-state of shapeit. The function will also return warnings if
+        the intermediate files were not found.
+
+        Parameters
+        ----------
+        output_prefix : str
+            Output filename prefix
+        error_intermediates : bool, optional
+            Defines if shapeit encountered an error
+
+    '''
+    if error_intermediates:
+
+        # Check that the log file was created, give a warning otherwise
+        if not os.path.isfile(output_prefix + '.phase.log'):
+            logging.warning('shapeit intermediate file %s.phase.log does not exist' % output_prefix)
+        else:
+            # Remove shapeit log file
+            os.remove(output_prefix + '.phase.log')
+
+    else:
+
+        # Check that the phase.ind.mm file was created, give a warning otherwise
+        if not os.path.isfile(output_prefix + '.phase.ind.mm'):
+            logging.warning('shapeit intermediate file %s.phase.ind.mm does not exist' % output_prefix)
+        else:
+            # Remove shapeit phase.ind.mm file
+            os.remove(output_prefix + '.phase.ind.mm')
+
+        # Check that the phase.snp.mm file was created, give a warning otherwise
+        if not os.path.isfile(output_prefix + '.phase.snp.mm'):
+            logging.warning('shapeit intermediate file %s.phase.snp.mm does not exist' % output_prefix)
+        else:
+            # Remove shapeit phase.snp.mm file
+            os.remove(output_prefix + '.phase.snp.mm')
+
+    # Check that the haps file was created, give a warning otherwise
+    if not os.path.isfile(output_prefix + '.haps'):
+        logging.warning('shapeit intermediate file %s.haps does not exist' % output_prefix)
+    else:
+        # Remove shapeit haps file
+        os.remove(output_prefix + '.haps')
+
+    # Check that the sample file was created, give a warning otherwise
+    if not os.path.isfile(output_prefix + '.sample'):
+        logging.warning('shapeit intermediate file %s.sample does not exist' % output_prefix)
+    else:
+        # Remove shapeit sample file
+        os.remove(output_prefix + '.sample')
+
+    logging.info('shapeit-related files removed')
+
+def standard_shapeit_call (shapeit_call_args, output_prefix):
+    '''
+        Calls shapeit using subprocess
+
+        This function is used to call shapeit and passes the resulting stdout
+        to check_shapeit_for_errors to check for errors. The function also
+        passes output_prefix to check_shapeit_for_errors to delete shapeit
+        intermediate files if shapeit results in an error.
+
+        Parameters
+        ----------
+        shapeit_call_args : list
+            Argument list for shapeit
+        output_prefix : str
+            Output filename prefix
+
+    '''
+
+    logging.info('shapeit phasing parameters assigned')
+
+    # Phasing subprocess call
+    phase_call = subprocess.Popen(['shapeit'] + shapeit_call_args, stdout = subprocess.PIPE, stderr = subprocess.PIPE)
+    phase_stdout, phase_stderr = phase_call.communicate()
+
+    # Check if code is running in python 3
+    if sys.version_info[0] == 3:
+        # Convert bytes to string
+        phase_stdout = phase_stdout.decode()
+
+    # Check shapeit call for errors
+    check_shapeit_for_errors(phase_stdout, output_prefix)
+
+    logging.info('shapeit phasing complete (HAPS format)')
+
+def call_shapeit (shapeit_call_args, output_prefix, output_format):
+    '''
+        Calls shapeit and automates file conversions
+
+        The function is used to call shapeit and also automates conversion to
+        VCF, VCF.GZ, and BCF using plink2
+
+        Parameters
+        ----------
+        shapeit_call_args : list
+            Argument list for shapeit
+        output_prefix : str
+            Output filename prefix
+        output_format : str
+            Output file format
+
+    '''
+
+    # Standard call to beagle
+    standard_shapeit_call(shapeit_call_args, output_prefix)
+
+    # Convert haps-format to vcf
+    convert_haps_to_vcf(output_prefix, output_format)
+
+    logging.info('HAPS conversion to VCF complete')
diff -r 000000000000 -r 3830d29fca6a vcf_phase.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/vcf_phase.py	Mon Oct 15 18:15:47 2018 -0400
@@ -0,0 +1,542 @@
+import os
+import sys
+import copy
+import shutil
+import argparse
+import glob
+import logging
+
+sys.path.insert(0, os.path.abspath(os.path.join(os.pardir, 'jared')))
+
+from vcf_reader_func import checkFormat
+from logging_module import initLogger, logArgs
+from model import read_model_file
+from beagle import call_beagle
+from shapeit import call_shapeit, remove_intermediate_files
+from bcftools import pipe_bcftools_to_chr, chr_subset_file, concatenate
+
+def phase_argument_parser(passed_arguments):
+    '''Phase Argument Parser - Assigns arguments for vcftools from command line.
+    Depending on the argument in question, a default value may be specified'''
+
+    def parser_confirm_file ():
+        '''Custom action to confirm file exists'''
+        class customAction(argparse.Action):
+            def __call__(self, parser, args, value, option_string=None):
+                if not os.path.isfile(value):
+                    raise IOError('%s not found' % value)
+                setattr(args, self.dest, value)
+        return customAction
+
+    def metavar_list (var_list):
+        '''Create a formmated metavar list for the help output'''
+        return '{' + ', '.join(var_list) + '}'
+
+    phase_parser = argparse.ArgumentParser()
+
+    # Input arguments
+    phase_parser.add_argument('--vcf', help = "Input VCF filename", type = str, required = True, action = parser_confirm_file())
+
+    # Model file arguments
+    phase_parser.add_argument('--model-file', help = 'Defines the model file', type = str, action = parser_confirm_file())
+    phase_parser.add_argument('--model', help = 'Defines the model to analyze', type = str)
+
+    # General arguments
+    phase_parser.add_argument('--overwrite', help = "Overwrite previous output files", action = 'store_true')
+    phase_parser.add_argument('--beagle-path', help = "Defines path to locate beagle.jar", type = str, default = 'bin/')
+
+    # Phase algorithm argument
+    phasing_list = ['beagle', 'shapeit']
+    phasing_default = 'beagle'
+    phase_parser.add_argument('--phase-algorithm', metavar = metavar_list(phasing_list), help = 'Specifies the phase algorithm to be used', type = str, choices = phasing_list, default = phasing_default)
+
+    # Common phasing arguments
+    phase_parser.add_argument('--Ne', help = 'Defines the effective population size', type = int)
+    phase_parser.add_argument('--random-seed', help="Defines the random seed value for the random number generator", type = int)
+    phase_parser.add_argument('--genetic-map', help = 'Genetic map filename', type = str, action = parser_confirm_file())
+    phase_parser.add_argument('--phase-chr', help = 'Selects a single chromosome to phase', type = str)
+    phase_parser.add_argument('--phase-from-bp', help = 'Lower bound of sites to include (Only usable with a single chromosome)', type = int)
+    phase_parser.add_argument('--phase-to-bp', help = 'Upper bound of sites to include (Only usable with a single chromosome)', type = int)
+
+    # Shapeit-specific options
+    phase_parser.add_argument('--shapeit-burn-iter', help = 'Number of the burn-in iterations (shapeit)', type = int)
+    phase_parser.add_argument('--shapeit-prune-iter', help = 'Number of pruning iterations (shapeit)', type = int)
+    phase_parser.add_argument('--shapeit-main-iter', help = 'Number of main iterations (shapeit)', type = int)
+    phase_parser.add_argument('--shapeit-states', help = 'Number of conditioning states for haplotype estimation (shapeit)', type = int)
+    phase_parser.add_argument('--shapeit-window', help = 'Model window size in Mb (shapeit)', type = float)
+
+    # Beagle-specific options
+    phase_parser.add_argument('--beagle-burn-iter', help = 'Number of the burn-in iterations (beagle)', type = int)
+    phase_parser.add_argument('--beagle-iter', help = 'Number of iterations after burn-in (beagle)', type = int)
+    phase_parser.add_argument('--beagle-states', help = 'Number of model states for genotype estimation (beagle)', type = int)
+    phase_parser.add_argument('--beagle-window', help = 'Sliding window size in cM (beagle)', type = float)
+    phase_parser.add_argument('--beagle-overlap', help = 'Overlap between neighboring sliding windows in cM (beagle)', type = float)
+    phase_parser.add_argument('--beagle-error', help = 'HMM allele mismatch probability (beagle)', type = float)
+    phase_parser.add_argument('--beagle-step', help = 'Step length in cM used for identifying short IBS segments (beagle)', type = float)
+    phase_parser.add_argument('--beagle-nsteps', help = 'Number of consecutive --beagle-steps used for identifying long IBS segments (beagle)', type = int)
+
+    # Output arguments
+    phase_parser.add_argument('--out', help = 'Defines the output filename')
+    phase_parser.add_argument('--out-prefix', help = 'Defines the output prefix (used by phasing algorithms)', default = 'out')
+
+    out_format_list = ['vcf', 'vcf.gz', 'bcf']
+    out_format_default = 'vcf.gz'
+
+    phase_parser.add_argument('--out-format', metavar = metavar_list(out_format_list), help = 'Specifies the output format.', type = str, choices = out_format_list, default = out_format_default)
+
+    if passed_arguments:
+        return phase_parser.parse_args(passed_arguments)
+    else:
+        return phase_parser.parse_args()
+
+def concatenate_logs (log_files, new_log_filename):
+    '''
+        Concatenates log files
+
+        This function is used to concatenate log files. If a single chromosome
+        is being analyzed, the function is used to concatenate the log files
+        from multiple wrappers. If multiple chromosomes are being analyzed the
+        function is used to concatenate the logs of each chromosome.
+
+        Raises
+        ------
+        IOError
+            If a log file does not exist
+    '''
+
+    # Open the new log file
+    with open(new_log_filename, 'wb') as out_file:
+        # Loop the logs to be combined
+        for log_file in log_files:
+            # Confirm the log file exists, raise an error if not
+            if not os.path.isfile(log_file):
+                raise IOError('%s not found' % log_file)
+            # Open the log file
+            with open(log_file, 'rb') as in_file:
+                # Copy the contents of the file
+                shutil.copyfileobj(in_file, out_file)
+
+    # Remove old log files
+    for log_file in log_files:
+        os.remove(log_file)
+
+def assign_vcf_extension (filename):
+    # Checks if the file is unzipped, bgzipped, or gzipped
+    vcfname_format = checkFormat(filename)
+
+    if vcfname_format == 'vcf':
+        return '.vcf'
+
+    elif vcfname_format == 'gzip' or vcfname_format == 'bgzip':
+        return '.vcf.gz'
+
+    elif vcfname_format == 'bcf':
+        raise Exception('BCF not supported in current version')
+
+    else:
+        raise Exception('Unknown file format')
+
+def run (passed_arguments = []):
+    '''
+        Phaser for VCF files.
+
+        Automates the phasing process for a specified VCF file. The function
+        allows users to select between multiple phasing algorithms: beagle
+        (default) and shapit.
+
+        Parameters
+        ----------
+        --vcf : str
+            Specifies the input VCF filename
+        --phase-algorithm : str
+            Specifies the algorithm to be used. Choices: beagle (default) and
+            shapit
+        --out : str
+            Specifies the output filename
+
+        Returns
+        -------
+        output : file
+            Phased VCF file
+
+        Raises
+        ------
+        IOError
+            Input VCF file does not exist
+        IOError
+            Output file already exists
+
+    '''
+
+    # Grab VCF arguments from command line
+    phase_args = phase_argument_parser(passed_arguments)
+
+    # Adds the arguments (i.e. parameters) to the log file
+    logArgs(phase_args, func_name = 'vcf_phase')
+
+    # Assign file extension for VCF input file
+    vcfname_ext = assign_vcf_extension(phase_args.vcf)
+
+    logging.info('Input file assigned')
+
+    # Used to confirm if the VCF input was renamed
+    vcfname_renamed = False
+
+    # Confirm input has correct file extension
+    if vcfname_ext not in phase_args.vcf:
+        vcfname_renamed = True
+        os.rename(phase_args.vcf, phase_args.vcf + vcfname_ext)
+        phase_args.vcf += vcfname_ext
+
+    # List to hold beagle/shapeit arguments
+    phase_call_args = []
+
+    # bool to check if an eclude file was created
+    exclude_file_created = False
+
+    # bool to check if an include file was created
+    include_file_created = False
+
+    # Performs actions related to --model-file and argument assignment
+    if phase_args.model_file:
+
+        # Check if a model has been specified
+        if not phase_args.model:
+            raise IOError('No selected model. Please use --model to select a model')
+
+        # Read in the model file
+        models_file = read_model_file(phase_args.model_file)
+
+        # Check that the selected model was found in the file
+        if phase_args.model not in models_file:
+            raise IOError('Selected model "%s" not found in: %s' % (phase_args.model, phase_args.model_file))
+
+        logging.info('Model assigned')
+
+        if phase_args.phase_algorithm == 'beagle':
+
+            # Get list of individuals to include
+            inds_in_model = models_file[phase_args.model].inds
+
+            # Create file with individuals to exclude (beagle has no include option)
+            models_file.create_exclude_ind_file(inds_to_include = inds_in_model, overwrite = phase_args.overwrite)
+
+            # Assign exclude file
+            phase_call_args.append('excludesamples=' + models_file.exclude_file)
+
+            # Confirm that the exclude file was created
+            exclude_file_created = True
+
+        elif phase_args.phase_algorithm == 'shapeit':
+
+            # Assign the model
+            selected_model = models_file[phase_args.model]
+
+            # Create file with individuals to exclude (beagle has no include option)
+            selected_model.create_ind_file(overwrite = phase_args.overwrite)
+
+            # Assign exclude file
+            phase_call_args.extend(['--include-ind', selected_model.ind_file])
+
+            # Confirm that the include file was created
+            include_file_created = True
+
+    # Get the list of chromosomes within the VCF
+    chrs_in_vcf = pipe_bcftools_to_chr(phase_args.vcf)
+
+    # Check if the user specified a specific chromosome
+    if phase_args.phase_chr:
+        # Check if the chromosome is within the file
+        if phase_args.phase_chr not in chrs_in_vcf:
+            raise Exception('--phase-chr %s not found in %s' % (phase_args.phase_chr, phase_args.vcf))
+
+    # Check that a chr was specified if a bp-based argument was specified
+    if (phase_args.phase_from_bp or phase_args.phase_to_bp) and not phase_args.phase_chr:
+        # Check if the file has a single chromosome
+        if len(chrs_in_vcf) == 1:
+            # Check if the beagle algorithm is being used
+            if phase_args.phase_algorithm == 'beagle':
+                # Assign the chromosome, as beagle requires it to be assigned
+                phase_args.phase_chr = chrs_in_vcf[0]
+
+                logging.info('Assigned chr %s for beagle' % chrs_in_vcf[0])
+        # Report error if multiple chromosomes are within the file
+        else:
+            raise Exception('The --phase-from-bp and --phase-to-bp arguments '
+                            'require the --phase-chr argument to function if '
+                            'multiple chrs are within %s' % phase_args.vcf)
+
+    # Assign general arguments and call beagle
+    if phase_args.phase_algorithm == 'beagle':
+
+        # Assigns the input and output arguments
+        phase_call_args.extend(['gt=' + phase_args.vcf,
+                                'out=' + phase_args.out_prefix])
+
+        # Assign the burn-in iter, if specified
+        if phase_args.beagle_burn_iter:
+            phase_call_args.append('burnin=' + phase_args.beagle_burn_iter)
+
+        # Assign the iter, if specified
+        if phase_args.beagle_iter:
+            phase_call_args.append('iterations=' + phase_args.beagle_iter)
+
+        # Assign the state count, if specified
+        if phase_args.beagle_states:
+            phase_call_args.append('phase-states=' + phase_args.beagle_states)
+
+        # Assign the window length, if specified
+        if phase_args.beagle_window:
+            phase_call_args.append('window=' + phase_args.beagle_window)
+
+        # Assign the window overlap, if specified
+        if phase_args.beagle_overlap:
+            phase_call_args.append('overlap=' + phase_args.beagle_overlap)
+
+        # Assign the HMM error, if specified
+        if phase_args.beagle_error:
+            phase_call_args.append('err=' + phase_args.beagle_error)
+
+        # Assign the step length, if specified
+        if phase_args.beagle_step:
+            phase_call_args.append('step=' + phase_args.beagle_step)
+
+        # Assign the number of steps, if specified
+        if phase_args.beagle_nsteps:
+            phase_call_args.append('nsteps=' + phase_args.beagle_nsteps)
+
+        # Assign the genetic map, if specified
+        if phase_args.genetic_map:
+            phase_call_args.append('map=' + phase_args.genetic_map)
+
+        # Assign the effective pop size, if specified
+        if phase_args.Ne:
+            phase_call_args.append('ne=' + phase_args.Ne)
+
+        # Assign the random seed, if specified
+        if phase_args.random_seed:
+            phase_call_args.append('seed=' + phase_args.random_seed)
+
+        # Assign the chromosome to phase, if specified
+        if phase_args.phase_chr:
+
+            # Store the chromosome argument
+            chr_arg = 'chrom=%s' % phase_args.phase_chr
+
+            # Check if either bp position arguments were specified
+            if phase_args.phase_from_bp or phase_args.phase_to_bp:
+
+                # List of the position arguments, in their required order
+                position_args = [':', phase_args.phase_from_bp, '-', phase_args.phase_to_bp]
+
+                # Filter the position arguments to remove empty values
+                filttered_position_args = filter(None, position_args)
+
+                # Map the arguments to str and add them to the chromosome argument
+                chr_arg += ''.join(map(str, filttered_position_args))
+
+            # Add the chromosome argument the argument list
+            phase_call_args.append(chr_arg)
+
+        # Call beagle wrapper
+        call_beagle(phase_args.beagle_path, list(map(str, phase_call_args)), phase_args.out_prefix, phase_args.out_format)
+
+        # Assign expected phased output filename
+        phased_output = '%s.%s' % (phase_args.out_prefix, phase_args.out_format)
+
+        # Rename output to phase_args.out, if specified
+        if phase_args.out:
+            shutil.move(phased_output, phase_args.out)
+            shutil.move(phase_args.out_prefix + '.log', phase_args.out + '.log')
+        # Rename log using phased_output
+        else:
+            shutil.move(phase_args.out_prefix + '.log', phased_output + '.log')
+
+        logging.info('beagle log file created')
+
+    # Assign general arguments and call shapeit
+    elif phase_args.phase_algorithm == 'shapeit':
+
+        # Assign the shapeit burn in iter, if specified
+        if phase_args.shapeit_burn_iter:
+            phase_call_args.extend(['--burn', phase_args.shapeit_burn_iter])
+
+        # Assign the shapeit prune iter, if specified
+        if phase_args.shapeit_prune_iter:
+            phase_call_args.extend(['--prune', phase_args.shapeit_prune_iter])
+
+        # Assign the shapeit main iter, if specified
+        if phase_args.shapeit_main_iter:
+            phase_call_args.extend(['--main', phase_args.shapeit_main_iter])
+
+        # Assign the number of shapeit states if specified
+        if phase_args.shapeit_states:
+            phase_call_args.extend(['--states', phase_args.shapeit_states])
+
+        # Assign the shapeit window size, if specified
+        if phase_args.shapeit_window:
+            phase_call_args.extend(['--window', phase_args.shapeit_window])
+
+        # Assign the genetic map, if specified
+        if phase_args.genetic_map:
+            phase_call_args.extend(['--input-map', phase_args.genetic_map])
+
+        # Assign the effective pop size, if specified
+        if phase_args.Ne:
+            phase_call_args.extend(['--effective-size', phase_args.Ne])
+
+        # Assign the random seed, if specified
+        if phase_args.random_seed:
+            phase_call_args.extend(['--seed', phase_args.random_seed])
+
+        # Check if only a single shapeit run is required
+        if phase_args.phase_chr or len(chrs_in_vcf) == 1:
+
+            # Holds shapeit input filename, as a temp file may be required
+            shapeit_input_vcf = None
+
+            # Check if a single chromosome is selected
+            if phase_args.phase_chr and len(chrs_in_vcf) > 1:
+                # Assign the chromosome-specific input prefix
+                chr_input_prefix = phase_args.vcf + '.' + phase_args.phase_chr
+
+                # Assign the chr-specific file as the shapeit input filename
+                shapeit_input_vcf = chr_input_prefix + '.vcf.gz'
+
+                # Create the chromosome-specific input
+                chr_subset_file(phase_args.vcf,
+                                phase_args.phase_chr,
+                                chr_input_prefix,
+                                'vcf.gz',
+                                from_bp = phase_args.phase_from_bp,
+                                to_bp = phase_args.phase_to_bp)
+
+                logging.info('Chr %s subset created' % phase_args.phase_chr)
+
+            else:
+                # Assign the shapeit input filename
+                shapeit_input_vcf = phase_args.vcf
+
+            # Assigns the input and output arguments for shapeit
+            phase_call_args.extend(['--input-vcf', shapeit_input_vcf,
+                                    '--output-max', phase_args.out_prefix,
+                                    '--output-log', phase_args.out_prefix + '.phase.log'])
+
+            # Assign the from bp position argument, if specified
+            if phase_args.phase_from_bp:
+                phase_call_args.extend(['--input-from', phase_args.phase_from_bp])
+
+            # Assign the to bp position argument, if specified
+            if phase_args.phase_to_bp:
+                phase_call_args.extend(['--input-to', phase_args.phase_to_bp])
+
+            # Call shapeit wrapper
+            call_shapeit(list(map(str, phase_call_args)), phase_args.out_prefix, phase_args.out_format)
+
+            # Assign expected phased output filename
+            phased_output = '%s.%s' % (phase_args.out_prefix, phase_args.out_format)
+
+            # Combine the log files
+            concatenate_logs([phase_args.out_prefix + '.phase.log', phase_args.out_prefix + '.log'],  phased_output + '.log')
+
+            # Rename output to phase_args.out, if specified
+            if phase_args.out:
+                shutil.move(phased_output, phase_args.out)
+                shutil.move(phased_output + '.log',  phase_args.out + '.log')
+
+            logging.info('shapeit log file created')
+
+            # Check if a chr subset file was created
+            if phase_args.phase_chr and len(chrs_in_vcf) > 1:
+                # Delete the chromosome-specific input
+                os.remove(shapeit_input_vcf)
+
+                logging.info('Chr %s subset deleted' % phase_args.phase_chr)
+
+            # Remove intermediate files created by shapeit
+            remove_intermediate_files(phase_args.out_prefix)
+
+        # Check if multiple shapeit runs are required
+        else:
+
+            # List to store the phased filenames
+            phased_filename_list = []
+
+            # List to store the phased logs
+            phased_log_list = []
+
+            logging.info('Multi-chr shapeit phasing assigned')
+
+            for chr_in_vcf in chrs_in_vcf:
+
+                logging.info('Chr %s assigned' % chr_in_vcf)
+
+                # Copy the arguments for this run
+                chr_call_args = copy.deepcopy(phase_call_args)
+
+                # Assign the chromosome-specific output prefix
+                chr_out_prefix = phase_args.out_prefix + '.' + chr_in_vcf
+
+                # Assign the expected chromosome-specific output filename
+                chr_out_filename = '%s.%s' % (chr_out_prefix, phase_args.out_format)
+
+                # Store the output filename
+                phased_filename_list.append(chr_out_filename)
+
+                # Assign the chromosome-specific input prefix
+                chr_input_prefix = phase_args.vcf + '.' + chr_in_vcf
+
+                # Assign the expected chromosome-specific input filename
+                chr_in_filename = chr_input_prefix + '.vcf.gz'
+
+                # Create the chromosome-specific input
+                chr_subset_file(phase_args.vcf,
+                                chr_in_vcf,
+                                chr_input_prefix,
+                                'vcf.gz')
+
+                # Assigns the input and output arguments for shapeit
+                chr_call_args.extend(['--input-vcf', chr_in_filename,
+                                      '--output-max', chr_out_prefix,
+                                      '--output-log', chr_out_prefix + '.phase.log'])
+
+                # Call shapeit wrapper
+                call_shapeit(list(map(str, chr_call_args)), chr_out_prefix, phase_args.out_format)
+
+                # Combine log files
+                concatenate_logs([chr_out_prefix + '.phase.log', chr_out_prefix + '.log'],  chr_out_filename + '.log')
+
+                # Store the filename of the combined logs
+                phased_log_list.append(chr_out_filename + '.log')
+
+                # Delete the chromosome-specific input
+                os.remove(chr_in_filename)
+
+                # Remove intermediate files created by shapeit
+                remove_intermediate_files(chr_out_prefix)
+
+            # Concatenate the vcf files
+            concatenate(phased_filename_list, phase_args.out_prefix, phase_args.out_format)
+
+            logging.info('Concatenated chromosomes')
+
+            # Assign expected concatenated output filename
+            phased_output = '%s.%s' % (phase_args.out_prefix, phase_args.out_format)
+
+            # Combine the log files
+            concatenate_logs(phased_log_list, phased_output + '.log')
+
+            # Rename output to phase_args.out, if specified
+            if phase_args.out:
+                shutil.move(phased_output, phase_args.out)
+                shutil.move(phased_output + '.log',  phase_args.out + '.log')
+
+            logging.info('Multi-chr shapeit log created')
+
+    # Reverts the VCF input file
+    if vcfname_renamed:
+        os.rename(phase_args.vcf, phase_args.vcf[:-len(vcfname_ext)])
+
+if __name__ == "__main__":
+    initLogger()
+    run()
diff -r 000000000000 -r 3830d29fca6a vcf_phase.xml
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/vcf_phase.xml	Mon Oct 15 18:15:47 2018 -0400
@@ -0,0 +1,166 @@
+<tool id="vcf_phase" name="Phase VCF" version="1.0.0.1">
+
+  <description>files with BEAGLE or SHAPEIT</description>
+
+  <requirements>
+      <requirement type="package" >pandas</requirement>
+      <requirement type="package" >pysam</requirement>
+      <requirement type="package" >shapeit</requirement>
+      <requirement type="package" >beagle</requirement>
+  </requirements>
+
+  <command><![CDATA[
+    #if $input.is_of_type('vcf_bgzip')
+      ln -fs $input input.vcf.gz &&
+    #end if
+    #if $input.is_of_type('vcf')
+      ln -fs $input input.vcf &&
+    #end if
+    python $__tool_directory__/vcf_phase.py
+            #if $input.is_of_type('vcf_bgzip')
+               --vcf input.vcf.gz
+            #end if
+            #if $input.is_of_type('vcf')
+               --vcf input.vcf
+            #end if
+            #if $model_file
+              --model-file $model_file
+              --model $model
+            #end if
+            --phase-algorithm $phase.phase_algorithm
+            #if $phase.beagle_burn_iter
+              --beagle-burn-iter $common.beagle_burn_iter
+            #end if
+            #if $phase.beagle_burn_iter
+              --beagle-burn-iter $phase.beagle_burn_iter
+            #end if
+            #if $phase.phase_algorithm == 'beagle'
+              --beagle-path $__tool_data_path__/shared/jars/
+              #if $phase.beagle_iter
+                --beagle-iter $phase.beagle_iter
+              #end if
+              #if $phase.beagle_states
+                --beagle-states $phase.beagle_states
+              #end if
+              #if $phase.beagle_window
+                --beagle-window $phase.beagle_window
+              #end if
+              #if $phase.beagle_overlap
+                --beagle-overlap $phase.beagle_overlap
+              #end if
+              #if $phase.beagle_error
+                --beagle-error $phase.beagle_error
+              #end if
+              #if $phase.beagle_step
+                --beagle-step $phase.beagle_step
+              #end if
+              #if $phase.beagle_nsteps
+                --beagle-nsteps $phase.beagle_nsteps
+              #end if
+            #end if
+            #if $phase.phase_algorithm == 'shapeit'
+              #if $phase.shapeit_burn_iter
+                --shapeit-burn-iter $phase.shapeit_burn_iter
+              #end if
+              #if $phase.shapeit_prune_iter
+                --shapeit-prune-iter $phase.shapeit_prune_iter
+              #end if
+              #if $phase.shapeit_main_iter
+                --shapeit-main-iter $phase.shapeit_main_iter
+              #end if
+              #if $phase.shapeit_states
+                --shapeit-states $phase.shapeit_states
+              #end if
+              #if $phase.shapeit_window
+                --shapeit-window $phase.shapeit_window
+              #end if
+            #end if
+            #if $common.genetic_map
+              --genetic-map $common.genetic_map
+            #end if
+            #if $common.ne
+              --Ne $common.ne
+            #end if
+            #if $common.random_seed
+              --random-seed $common.random_seed
+            #end if
+            #if $common.phase_chr
+              --phase-chr $common.phase_chr
+            #end if
+            #if $common.phase_from_bp
+              --phase-from-bp $common.phase_from_bp
+            #end if
+            #if $common.phase_to_bp
+              --phase-to-bp $common.phase_to_bp
+            #end if
+            --out $output
+            --out-format $out_format
+  ]]></command>
+
+  <inputs>
+
+    <param format="vcf,vcf_bgzip" name="input" type="data" label="VCF Input"/>
+
+    <param format="model" name="model_file" type="data" label="Model Input" optional="True"/>
+    <param name="model" type="select" label="Select Model" refresh_on_change="True">
+      <options>
+        <filter type="data_meta" ref="model_file" key="models"/>
+      </options>
+    </param>
+
+    <conditional name="phase">
+      <param name="phase_algorithm" type="select" label="Phase Algorithm" refresh_on_change='True'>
+        <option value="beagle" selected="True" >Beagle</option>
+        <option value="shapeit">SHAPEIT</option>
+      </param>
+      <when value="beagle">
+        <param name="beagle_burn_iter" type="integer" label="Burn-in iterations" optional="True"/>
+        <param name="beagle_iter" type="integer" label="Post burn-in iterations" optional="True"/>
+        <param name="beagle_states" type="integer" label="Model states for genotype estimation" optional="True"/>
+        <param name="beagle_window" type="float" label="Sliding window size (cM)" optional="True"/>
+        <param name="beagle_overlap" type="float" label="Overlap between neighboring sliding windows (cM)" optional="True"/>
+        <param name="beagle_error" type="float" label="HMM allele mismatch probability" optional="True"/>
+        <param name="beagle_step" type="float" label="Step length (cM)" optional="True" help="Used for identifying short IBS segments"/>
+        <param name="beagle_nsteps" type="integer" label="Number of consecutive steps" optional="True" help="Used for identifying long IBS segments"/>
+      </when>
+      <when value="shapeit">
+        <param name="shapeit_burn_iter" type="integer" label="Burn-in iterations" optional="True"/>
+        <param name="shapeit_prune_iter" type="integer" label="Pruning iterations" optional="True"/>
+        <param name="shapeit_main_iter" type="integer" label="Main iterations" optional="True"/>
+        <param name="shapeit_states" type="integer" label="Conditioning states for haplotype estimation" optional="True"/>
+        <param name="shapeit_window" type="float" label="Model window size (Mb)" optional="True"/>
+      </when>
+    </conditional>
+
+    <section name="common" title="Common Parameters" expanded="True">
+      <param name="random_seed" type="integer" label="Seed value" optional="True"/>
+      <param format="text" name="genetic_map" type="data" label="Genetic Map" optional="True"/>
+      <param name="ne" type="integer" label="Effective population size" optional="True"/>
+      <param name="phase_chr" type="text" label="Chromosome to phase" optional="True"/>
+      <param name="phase_from_bp" type="integer" label="Lower bound of sites to include" help="Only usable with a single chromosome" optional="True"/>
+      <param name="phase_to_bp" type="integer" label="Upper bound of sites to include" help="Only usable with a single chromosome" optional="True"/>
+    </section>
+
+    <param name="out_format" type="select" label="Output Format">
+      <option value="vcf">VCF File</option>
+      <option value="vcf.gz" selected="True">bgzipped-VCF File</option>
+      <option value="bcf">BCF File</option>
+    </param>
+
+  </inputs>
+  <outputs>
+
+    <data name="output" format="vcf_bgzip">
+      <change_format>
+        <when input="out_format" value="vcf" format="vcf"/>
+        <when input="out_format" value="bcf" format="bcf"/>
+      </change_format>
+    </data>
+
+  </outputs>
+
+  <help>
+      VCF Phaser Help Text
+  </help>
+
+</tool>
diff -r 000000000000 -r 3830d29fca6a vcf_reader_func.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/vcf_reader_func.py	Mon Oct 15 18:15:47 2018 -0400
@@ -0,0 +1,474 @@
+import sys
+import pysam
+import logging
+import struct
+from random import sample
+from collections import defaultdict
+import os
+import gzip
+
+def checkIfGzip(filename):
+    try:
+        gf = gzip.open(filename)
+        gl = gf.readline()
+        gf.close()
+        vcf_check = b'##fileformat=VCF'
+        if gl[0:3] == b'BCF':
+            return 'bcf'
+        elif gl[:len(vcf_check)] == vcf_check:
+            return checkHeader(filename)
+        else:
+            return 'other'
+    except:
+        return 'nozip'
+
+def checkHeader(filename):
+    f = open(filename,'rb')
+    l = f.readline()
+    f.close()
+    BGZF_HEADER=b'\x1f\x8b\x08\x04\x00\x00\x00\x00\x00\xff\x06\x00\x42\x43\x02\x00'
+    #BGZF_HEADER=b'\x1f\x8b\x08\x04\x00\x00\x00\x00\x00\xff'
+    GZF_HEADER=b'\x1f\x8b'
+    if l[:len(BGZF_HEADER)] == BGZF_HEADER:
+        return 'bgzip'
+    if l[:len(GZF_HEADER)] == GZF_HEADER:
+        return 'gzip'
+    return 'nozip'
+
+def checkFormat(vcfname):
+    """Checks header of given file for compression type
+
+
+    Given a filename, opens file and reads first line to check if
+    file has BGZF or GZIP header. May be extended to check for BCF format
+
+    Parameters
+    ----------
+    filename : str
+        Name of file to be checked
+
+    Returns
+    -------
+    extension : str {'bgzip','gzip','vcf','other'}
+        File extension as indicated by header
+
+    """
+    typ = checkIfGzip(vcfname)
+    if typ != 'nozip':
+        return typ
+    f = open(vcfname)
+    l = f.readline()
+    f.close()
+    VCF_TAG='##fileformat=VCF'
+    if l[:len(VCF_TAG)] == VCF_TAG:
+        return 'vcf'
+    return 'other'
+
+def checkIfCpG(record,fasta_ref,offset=0,add_chr=False):
+    dr = None
+    pos = record.pos
+    c = record.chrom
+    if record.alts is None:
+        return False
+    if add_chr:
+        c = 'chr'+record.chrom
+    if record.ref == 'C' and 'T' in record.alts:
+        seq = fasta_ref.fetch(c,pos-1,pos+1)
+        if seq[0].upper() != 'C':
+            logging.warning('%s %d has bad base %s' % (record.chrom,record.pos,seq[0]))
+            #raise Exception("checkIfCpG function not lining up properly")
+        if seq[1].upper() == 'G':
+            return True
+        return False
+    elif record.ref == 'G' and 'A' in record.alts:
+        seq = fasta_ref.fetch(c,pos-2,pos)
+        if seq[1].upper() != 'G':
+            logging.warning('%s %d has bad base %s' % (record.chrom,record.pos,seq[1]))
+            #raise Exception("checkIfCpg function not lining up on negative strand")
+        if seq[0].upper() == 'C':
+            return True
+        return False
+    return False
+
+def checkForDuplicates(rec_list,pass_list):
+    for i in range(len(rec_list)-1):
+        if rec_list[i].pos == rec_list[i+1].pos:
+            pass_list[i] = False
+            pass_list[i+1] = False
+
+def checkForMultiallele(rec_list,pass_list):
+    for i in range(len(rec_list)):
+        if i != len(rec_list)-1 and rec_list[i].pos == rec_list[i+1].pos:
+            pass_list[i] = False
+            pass_list[i+1] = False
+        if len(rec_list[i].alleles) > 2:
+            pass_list[i] = False
+
+def flipChrom(chrom):
+    if chrom[0:3] == 'chr':
+        return chrom[0:3]
+    return 'chr'+chrom
+
+def getAlleleCountDict(rec):
+    alleles = defaultdict(int)
+    total_sites = 0
+    missing_inds = 0
+    for j in range(len(rec.samples)):
+        samp = rec.samples[j]
+        if None in samp.alleles:
+            missing_inds += 1
+        for k in range(len(samp.alleles)):
+            b = samp.alleles[k]
+            if b is not None:
+                alleles[b] += 1
+            total_sites+=1
+    return alleles, total_sites, missing_inds
+
+def isInformative(rec, mincount=2, alleles=None):
+    count = 0
+    if alleles is None:
+        alleles, total_sites, missing_inds = getAlleleCountDict(rec)
+    if len(alleles) != 2:
+        return False
+    i1,i2 = alleles.keys()
+    return (alleles[i1] >= mincount and alleles[i2] >= mincount)
+
+def getPassSites(record_list, remove_cpg=False, remove_indels=True,
+                 remove_multiallele=True, remove_missing=0,
+                 inform_level=2, fasta_ref=None):
+    pass_list = [True for r in record_list]
+    if remove_cpg == True and fasta_ref is None:
+        raise Exception("CpG removal requires a reference")
+    if inform_level > 2 or inform_level < 0:
+        raise Exception("Inform level %d must be between 0 and 2" % inform_level)
+    if remove_multiallele:
+        checkForMultiallele(record_list,pass_list)
+    for i in range(len(record_list)):
+        rec = record_list[i]
+        logging.info(rec.pos)
+        if remove_indels and not checkRecordIsSnp(rec):
+            pass_list[i] = False
+        if remove_cpg and checkIfCpG(rec,fasta_ref):
+            pass_list[i] = False
+        alleles,total_sites,missing_inds = getAlleleCountDict(rec)
+        if remove_missing != -1 and missing_inds > remove_missing:
+            pass_list[i] = False
+        if inform_level != 0 and not isInformative(rec,mincount=inform_level,alleles=alleles):
+            pass_list[i] = False
+    #pp = zip([r.pos for r in record_list],pass_list)
+    #for ppp in pp:
+    #    logging.info(ppp)
+    return pass_list
+
+
+def filterSites(record_list, remove_cpg=False, remove_indels=True,
+                remove_multiallele=True, remove_missing=0, inform_level=2,
+                fasta_ref=None):
+    pass_list = getPassSites(record_list,remove_cpg,remove_indels,remove_multiallele,remove_missing,inform_level,fasta_ref)
+    out_list = []
+    for i in range(len(pass_list)):
+        if pass_list[i]:
+            out_list.append(record_list[i])
+    return out_list
+
+
+
+class VcfReader():
+    def __init__(self, vcfname, compress_flag=False, subsamp_num=None,
+                 subsamp_fn=None, subsamp_list=None, index=None, popmodel=None, use_allpop=False):
+
+        ext = checkFormat(vcfname)
+        if ext in ['gzip','other'] :
+            raise Exception(('Input file %s is gzip-formatted, must be either '
+                             'uncompressed or zipped with bgzip' % vcfname))
+        self.file_uncompressed = (ext == 'vcf')
+        self.reader_uncompressed = (self.file_uncompressed and not compress_flag)
+        self.popmodel = None
+        self.popkeys = None
+        if popmodel is not None and use_allpop:
+            raise Exception("Popmodel and allpop cannot both be specified")
+        if compress_flag and file_uncompressed:
+            vcfname = compressVcf(vcfname)
+        if subsamp_num is not None:
+            if subsamp_list is not None:
+                raise Exception('Multiple subsampling options called in getVcfReader')
+            subsamp_list = getSubsampleList(vcfname, subsamp_num)
+        elif subsamp_fn is not None:
+            if subsamp_list is not None:
+                raise Exception('Multiple subsampling options called in getVcfReader')
+            subsamp_file = open(subsamp_fn,'r')
+            subsamp_list = [l.strip() for l in subsamp_file.readlines()]
+            subsamp_file.close()
+
+        if index is None:
+            self.reader = pysam.VariantFile(vcfname)
+        else:
+            self.reader = pysam.VariantFile(vcfname, index_filename=index)
+        if popmodel is not None:
+            self.popmodel = popmodel
+            popsamp_list = popmodel.inds
+            self.reader.subset_samples(popsamp_list)
+            self.setPopIdx()
+        if use_allpop:
+            self.setAllPop()
+        if subsamp_list is not None:
+            logging.debug('Subsampling %d individuals from VCF file' %
+            (len(subsamp_list)))
+            self.reader.subset_samples(subsamp_list)
+        self.prev_last_rec = next(self.reader)
+        self.chr_in_chrom = (self.prev_last_rec.chrom[0:3] == 'chr')
+
+    def fetch(self, chrom=None, start=None, end=None):
+        return self.reader.fetch(chrom, start, end)
+
+    def getRecordList(self, region=None, chrom=None, start=None,
+                      end=None):
+        if self.reader_uncompressed:
+            ret, self.prev_last_rec = getRecordListUnzipped(self.reader, self.prev_last_rec, region, add_chr=self.chr_in_chrom)
+            return ret
+        else:
+            return getRecordList(self.reader, region, chrom, start, end, self.chr_in_chrom)
+
+    def setPopIdx(self):
+        self.popkeys = {}
+        sample_names = [l for l in self.reader.header.samples]
+        for p in self.popmodel.pop_list:
+            self.popkeys[p] = []
+            for ind in self.popmodel.ind_dict[p]:
+                self.popkeys[p].append(sample_names.index(ind))
+
+    def close(self):
+        self.reader.close()
+
+    def setAllPop(self):
+        self.popkeys = {'ALL':[]}
+        for i in range(len(self.reader.header.samples)):
+            self.popkeys['ALL'].append(i)
+
+def modChrom(c,vcf_chr):
+    if c is None:
+        return None
+    if vcf_chr and c[:3] != 'chr':
+        return 'chr'+c
+    if not vcf_chr and c[:3] == 'chr':
+        return c[3:]
+    return c
+
+def getRecordList(vcf_reader, region=None, chrom=None, start=None,
+                  end=None, add_chr=False):
+    """Returns list for use in subsampling from input file"""
+    if region is not None:
+        c = modChrom(region.chrom,add_chr)
+        var_sites = vcf_reader.fetch(c, region.start, region.end)
+    else:
+        c = modChrom(chrom,add_chr)
+        var_sites = vcf_reader.fetch(c, start, end)
+    lst = []
+    for rec in var_sites:
+        lst.append(rec)
+    return lst
+
+
+def getRecordListUnzipped(vcf_reader, prev_last_rec, region=None, chrom=None,
+                          start=None, end=None, add_chr=False):
+    """Method for getting record list from unzipped VCF file.
+
+    This method will sequentially look through a VCF file until it finds
+    the given `start` position on `chrom`, then add all records to a list
+    until the `end` position has been reached. Note that `prev_last_rec`
+    must be kept track of externally to ensure that if consecutive regions
+    are called, the record of the first variant outside the first region
+    is not lost.
+
+    Parameters
+    ----------
+    vcf_reader : pysam VariantFile object
+        VCF reader initialized from other function
+    region : region object
+        Region object with start and end coordinates of region of interest.
+    prev_last_rec : VariantRecord object
+        Variable with last record read from VcfReader. Stored here so that
+        if two adjacent regions are called, the overflow record from the
+        first region is still included in the next region
+
+    Returns
+    -------
+    lst : list
+        List of records in given gene region
+    prev_last_rec : VariantRecord object
+        First record after target region, for use in next call
+
+    """
+    lst = []
+    if region is None:
+        lst.append(prev_last_rec)
+        for rec in vcf_reader:
+            lst.append(rec)
+        return lst, lst[-1]
+
+    if (prev_last_rec is not None and
+        region.containsRecord(prev_last_rec) == 'in'):
+        lst.append(prev_last_rec)
+    elif (prev_last_rec is not None and
+         region.containsRecord(prev_last_rec) == 'after'):
+        return []
+    rec = next(vcf_reader,None)
+    if rec is None:
+        return lst,None
+    place = region.containsRecord(rec)
+    while rec is not None and place != 'after':
+        if place == 'in':
+            lst.append(rec)
+        rec = next(vcf_reader,None)
+        if rec is None:
+            break
+        place = region.containsRecord(rec)
+    prev_last_rec = rec
+    return lst, prev_last_rec
+
+
+def checkRecordIsSnp(rec):
+    """Checks if this record is a single nucleotide variant, returns bool."""
+    if len(rec.ref) != 1:
+        return False
+    if rec.alts is None:
+        return False
+    for allele in rec.alts:
+        if len(allele) != 1:
+            return False
+    return True
+
+
+def getSubsampleList(vcfname, ss_count):
+    """Returns a list of the first `ss_count` individuals in `vcfname`
+
+    """
+
+    vcf_o = pysam.VariantFile(vcfname)
+    rec = next(vcf_o)
+    vcf_o.close()
+    lst = []
+    for samp in rec.samples:
+        lst.append(samp)
+    return lst[:int(ss_count)]
+
+
+def compressVcf(vcfname,forceflag=False,remove=False):
+    """Runs bgzip and tabix on input VCF file.
+
+    Using the pysam library, this function runs the bgzip and tabix utilities
+    on the given input file. By default, this will not overwrite an existing
+    zipped file, but will overwrite an existing index. `remove` can be set to
+    delete the unzipped file.
+
+    Parameters
+    ----------
+    vcfname : str
+        Name of uncompressed VCF file
+    forceflag : bool (False)
+        If true, will overwrite (vcfname).gz if it exists
+    remove : bool (False)
+        If true, will delete uncompressed source file
+
+    Returns
+    -------
+    cvcfname : str
+        Filepath to compressed VCF file
+    """
+    cvcfname = vcfname+".gz"
+    pysam.tabix_compress(vcfname,cvcfname,force=forceflag)
+    pysam.tabix_index(cvcfname,preset="vcf",force=True)
+    if remove:
+        os.remove(vcfname)
+    return cvcfname
+
+def vcfRegionName(prefix, region, ext, oneidx=False,
+                  halfopen=True, sep='-'):
+    chrom = region.toStr(halfopen, oneidx, sep)
+    return prefix+'_'+chrom+'.'+ext
+
+def getRecordsInRegion(region, record_list):
+    sub_list = []
+    for i in range(len(record_list)):
+        loc = region.containsRecord(record_list[i])
+        if loc == "in":
+            sub_list.append(record_list[i])
+        elif loc == "after":
+            break
+    return sub_list
+
+
+
+
+
+#def getVcfReader(args):
+def getVcfReader(vcfname, compress_flag=False, subsamp_num=None,
+                 subsamp_fn=None, subsamp_list=None, index=None):
+    """Returns a reader for a given input VCF file.
+
+    Given a filename, filetype, compression option, and optional Subsampling
+    options, will return a pysam.VariantFile object for iteration and
+    a flag as to whether this file is compressed or uncompressed.
+
+    Parameters
+    ----------
+    vcfname : str
+        Filename for VCF file. The extension of this file will be used to
+        determine whether it is compressed or not unless `var_ext` is set.
+    var_ext : str (None)
+        Extension for VCF file if it is not included in the filename.
+    compress_flag : bool (False)
+        If filetype is uncompressed and this is set to true, will run
+        compressVcf function.
+    subsamp_num : int (None)
+        If set, will randomly select `subsamp_num` individuals (not
+        genotypes) from the input VCF file and return a reader with
+        only those data.
+    subsamp_fn : str (None)
+        If set, will return a reader with only data from the samples listed
+        in the file provided. Cannot be used with other subsampling options.
+    subsamp_list : list (None)
+        If set, will return reader with records containing only
+        individuals named in the list. Cannot be used with other subsampling
+        options.
+
+    Returns
+    -------
+    vcf_reader : pysam.VariantFile
+        A reader that can be iterated through for variant records. If
+        compressed, it will be able to use the pysam fetch method, otherwise
+        it must be read through sequentially
+    reader_uncompressed : bool
+        If True, VCF reader is uncompressed. This means the fetch method
+        cannot be used and region access must be done using the
+        "getRecordListUnzipped" method.
+
+    """
+    ext = checkFormat(vcfname)
+    if ext in ['gzip','other'] :
+        raise Exception(('Input file %s is gzip-formatted, must be either '
+                         'uncompressed or zipped with bgzip' % vcfname))
+    file_uncompressed = (ext == 'vcf')
+    reader_uncompressed = (file_uncompressed and not compress_flag)
+    if compress_flag and file_uncompressed:
+        vcfname = compressVcf(vcfname)
+    #subsamp_list = None
+    if subsamp_num is not None:
+        if subsamp_list is not None:
+            raise Exception('Multiple subsampling options called in getVcfReader')
+        subsamp_list = getSubsampleList(vcfname, subsamp_num)
+    elif subsamp_fn is not None:
+        if subsamp_list is not None:
+            raise Exception('Multiple subsampling options called in getVcfReader')
+        subsamp_file = open(subsamp_fn,'r')
+        subsamp_list = [l.strip() for l in subsamp_file.readlines()]
+        subsamp_file.close()
+    if index is None:
+        vcf_reader = pysam.VariantFile(vcfname)
+    else:
+        vcf_reader = pysam.VariantFile(vcfname, index_filename=index)
+    if subsamp_list is not None:
+        logging.debug('Subsampling %d individuals from VCF file' %
+        (len(subsamp_list)))
+        vcf_reader.subset_samples(subsamp_list)
+    return vcf_reader, reader_uncompressed