genome-seek

#!/usr/bin/env python3
# -*- coding: UTF-8 -*-

"""
ABOUT: This is the main entry for the pipeline.
REQUIRES:
  - python>=3.6
  - snakemake   (required<8.0.0)
  - singularity (recommended==latest)
DISCLAIMER:
                    PUBLIC DOMAIN NOTICE
        NIAID Collaborative Bioinformatics Resource (NCBR)
   National Institute of Allergy and Infectious Diseases (NIAID)
This software/database is a "United  States Government Work" under
the terms of the United  States Copyright Act.  It was written as 
part of the author's official duties as a United States Government
employee and thus cannot be copyrighted. This software is freely
available to the public for use.
Although all  reasonable  efforts have been taken  to ensure  the
accuracy and reliability of the software and data, NCBR do not and
cannot warrant the performance or results that may  be obtained by 
using this software or data. NCBR and NIH disclaim all warranties,
express  or  implied,  including   warranties   of   performance, 
merchantability or fitness for any particular purpose.
Please cite the author and NIH resources like the "Biowulf Cluster" 
in any work or product based on this material.
USAGE:
  $ genome-seek <run> [OPTIONS]
EXAMPLE:
  $ genome-seek run --input *.R?.fastq.gz --output output/
"""

# Python standard library
from __future__ import print_function
import sys, os, subprocess, re, json, textwrap

# 3rd party imports from pypi
import argparse  # potential python3 3rd party package, added in python/3.5

# Local imports
from src import version 
from src.run import init, setup, bind, dryrun, runner
from src.shells import bash 
from src.utils import (
    Colors,
    err,
    exists,
    fatal,
    hashed,
    pairs,
    permissions,
    check_cache,
    require,
    tool_version
)


# Pipeline Metadata
__version__ = version
__authors__ = 'Skyler Kuhn, Justin Lack, Keyur Talsania'
__email__ = 'skyler.kuhn <AT> nih.gov, justin.lack <AT> nih.gov, keyur.talsania <AT> nih.gov'
__home__  =  os.path.dirname(os.path.abspath(__file__))
_name = os.path.basename(sys.argv[0])
_description = 'Clinical whole genome sequencing pipeline'


def unlock(sub_args):
    """Unlocks a previous runs output directory. If snakemake fails ungracefully,
    it maybe required to unlock the working directory before proceeding again.
    This is rare but it does occasionally happen. Maybe worth add a --force
    option to delete the '.snakemake/' directory in the future.
    @param sub_args <parser.parse_args() object>:
        Parsed arguments for unlock sub-command
    """
    print("Unlocking the pipeline's output directory...")
    outdir = sub_args.output

    try:
        unlock_output = subprocess.check_output([
            'snakemake', '--unlock',
            '--cores', '1',
            '--configfile=config.json'
        ], cwd = outdir,
        stderr=subprocess.STDOUT)
    except subprocess.CalledProcessError as e:
        # Unlocking process returned a non-zero exit code
        sys.exit("{}\n{}".format(e, e.output))

    print("Successfully unlocked the pipeline's working directory!")


def run(sub_args):
    """Initialize, setup, and run the pipeline.
    Calls initialize() to create output directory and copy over pipeline resources,
    setup() to create the pipeline config file, dryrun() to ensure their are no issues
    before running the pipeline, and finally run() to execute the Snakemake workflow.
    @param sub_args <parser.parse_args() object>:
        Parsed arguments for run sub-command
    """
    # Step 0. Check for required dependencies
    # The pipelines has only two requirements:
    # snakemake and singularity 
    require(['snakemake', 'singularity'], ['snakemake', 'singularity'])
    # Check the version of snakemake in the
    # user's $PATH, the pipeline supports
    # snakemake versions < 8.0.0.
    snakemake_version = tool_version('snakemake', ['snakemake', '--version'], strict = True)
    parsed_version = re.search(
        '^(?P<prefix>v)?(?P<major>0|[1-9]\d*)\.(?P<minor>0|[1-9]\d*)\.(?P<patch>0|[1-9]\d*)',
        snakemake_version.split()[-1]
    )
    # Check the parsed major version of snakemake
    if int(parsed_version.group('major')) >= 8:
        # Snakemake versions greater than or equal
        # to 8.0.0 are not supported by genome-seek.
        # This verison of snakemake introduced a set
        # of breaking changes that are not compatible
        # with the current pipeline. We will strictly 
        # enforce this until we move to profiles. 
        fatal(
            'Error: Snakemake version "{}" is not supported! Please use a version less than "8.0.0".'.format(
                snakemake_version
            )
        )

    if not sub_args.batch_id:
      # Set value of batch identifer by 
      # taking the MD5 checksum of the 
      # sorted list of input file names
      inputs = [os.path.basename(f) for f in sub_args.input] 
      batch_id = hashed(inputs)
      sub_args.batch_id = batch_id
    
    print(sub_args)

    # Step 1. Initialize working directory,
    # copy over required resources to run 
    # the pipeline
    git_repo = __home__
    input_files = init(
        repo_path = git_repo,
        output_path = sub_args.output,
        links = sub_args.input
    )

    # Step 2. Setup pipeline for execution, 
    # dynamically create config.json config
    # file from user inputs and base config
    # templates. Also override path of any 
    # reference files from default in the 
    # OpenOmics shared group area to base 
    # path provided by user via the option
    config = setup(sub_args, 
        ifiles = input_files,
        repo_path = git_repo,
        output_path = sub_args.output,
        resource_bundle = sub_args.resource_bundle
    )

    # Step 3. Resolve docker/singularity bind
    # paths from the config file.
    bindpaths = bind(
        sub_args,
        config = config
    )
    config['bindpaths'] = bindpaths

    # Step 4. Set tumor-normal pair information
    tn_pairs = {}  # No somatic variant calling 
    if sub_args.call_somatic:
        if sub_args.pairs:
            # Read in tumor, normal pairs file
            # If normal is missing for a sample,
            # then call it in tumor-only mode.
            # If a sample only has a normal, 
            # then only call germline variants.
            tn_pairs = pairs(sub_args.pairs)
        else:
            # Run all samples in tumor-only mode
            tn_pairs = {s: '' for s in config['samples']}
        
        # Check base names in the pairs 
        # file against base names from 
        # input FastQ files
        user_inputs = set(config['samples'])
        tumors = list(tn_pairs.keys())
        # Remove empty strings for the 
        # tumor-only samples 
        normals = [s for s in tn_pairs.values() if s]
        pair_inputs = set(tumors + normals)

        extra = pair_inputs - user_inputs
        if extra:
            # User provided samples in the 
            # pairs file, which were not 
            # in the --input option list
            err(
                "Error: Could not match some samples in the '--pairs' file "
                "against provided '--input' samples! Extra samples were provided.\n"
                "Set of extra samples in pairs file, {}:\n{}".format(
                    sub_args.pairs,
                    extra
                )
            )
            fatal('Please remove extra samples and try again!')

    # Add TN pair information where 
    # keys = tumor AND values = normal,
    # tumor-only samples have values set 
    # to an empty string
    config['pairs'] = tn_pairs

    # Step 5. Override default references
    if sub_args.pon: 
        # Use custom provided PON
        config['references']['PON'] = sub_args.pon
    
    if sub_args.wes_bed:
        # Use custom provided WES catpure kit
        config['references']['WES_BED'] = sub_args.wes_bed

    # Save config to output directory
    with open(os.path.join(sub_args.output, 'config.json'), 'w') as fh:
        json.dump(config, fh, indent = 4, sort_keys = True)

    # Optional Step: Dry-run pipeline
    if sub_args.dry_run:
        # Dryrun pipeline
        dryrun_output = dryrun(outdir = sub_args.output) # python3 returns byte-string representation
        print("\nDry-running {} pipeline:\n{}".format(_name, dryrun_output.decode("utf-8")))
        sys.exit(0)

    # Step 6. Orchestrate pipeline execution,
    # run pipeline in locally on a compute node
    # for debugging purposes or submit the master
    # job to the job scheduler, SLURM, and create 
    # logging file
    if not exists(os.path.join(sub_args.output, 'logfiles')):
        # Create directory for logfiles
        os.makedirs(os.path.join(sub_args.output, 'logfiles'))
    if sub_args.mode == 'local':
        log = os.path.join(sub_args.output, 'logfiles', 'snakemake.log')
    else: 
        log = os.path.join(sub_args.output, 'logfiles', 'master.log')
    logfh = open(log, 'w')
    mjob = runner(
        mode = sub_args.mode, 
        outdir = sub_args.output, 
        # additional_bind_paths = all_bind_paths,
        alt_cache = sub_args.singularity_cache,  
        threads = int(sub_args.threads),
        jobname = sub_args.job_name,
        submission_script=os.path.join(__home__, 'src', 'run.sh'),
        logger = logfh,
        additional_bind_paths = ",".join(bindpaths),
        tmp_dir = sub_args.tmp_dir, 
    )
    
    # Step 7. Wait for subprocess to complete,
    # this is blocking and not asynchronous
    if not sub_args.silent:
        print("\nRunning {} pipeline in '{}' mode...".format(_name, sub_args.mode)) 
    mjob.wait()
    logfh.close()

    # Step 7. Relay information about submission
    # of the master job or the exit code of the 
    # pipeline that ran in local mode
    if sub_args.mode == 'local':
        if int(mjob.returncode) == 0:
            print('{} pipeline has successfully completed'.format(_name))
        else:
            fatal('{} pipeline failed. Please see {} for more information.'.format(_name,
                os.path.join(sub_args.output, 'logfiles', 'snakemake.log')))
    elif sub_args.mode in ['slurm', 'uge']:
        jobid = open(os.path.join(sub_args.output, 'logfiles', 'mjobid.log')).read().strip()
        if not sub_args.silent:
            if int(mjob.returncode) == 0:
                print('Successfully submitted master job: ', end="")
            else:
                fatal('Error occurred when submitting the master job.')
        print(jobid)


def cache(sub_args):
    """Caches remote resources or reference files stored on DockerHub and S3.
    Local SIFs will be created from images defined in 'config/containers/images.json'.
    @TODO: add option to cache other shared S3 resources (i.e. kraken db and fqscreen indices)
    @param sub_args <parser.parse_args() object>:
        Parsed arguments for unlock sub-command
    """
    # Check for dependencies
    require(['singularity'], ['singularity'])
    sif_cache = sub_args.sif_cache
    # Get absolute PATH to templates in exome-seek git repo
    repo_path = os.path.dirname(os.path.abspath(__file__))
    images = os.path.join(repo_path, 'config','containers.json')

    # Create image cache
    if not exists(sif_cache):
        # Pipeline output directory does not exist on filesystem
        os.makedirs(sif_cache)
    elif exists(sif_cache) and os.path.isfile(sif_cache):
        # Provided Path for pipeline output directory exists as file
        raise OSError("""\n\tFatal: Failed to create provided sif cache directory!
        User provided --sif-cache PATH already exists on the filesystem as a file.
        Please {} cache again with a different --sif-cache PATH.
        """.format(_name)
        )

    # Check if local SIFs already exist on the filesystem
    with open(images, 'r') as fh:
        data = json.load(fh)

    pull = []
    for image, uri in data['images'].items():
        sif = os.path.join(sif_cache, '{}.sif'.format(os.path.basename(uri).replace(':', '_')))
        if not exists(sif):
            # If local sif does not exist on in cache, print warning
            # and default to pulling from URI in config/containers.json
            print('Image will be pulled from "{}".'.format(uri), file=sys.stderr)
            pull.append(uri)

    if not pull:
        # Nothing to do!
        print('Singularity image cache is already up to update!')
    else:
        # There are image(s) that need to be pulled 
        if not sub_args.dry_run:
            # container cache script: src/cache.sh
            # Quote user provided values to avoid shell injections
            username = os.environ.get('USER', os.environ.get('USERNAME'))
            exitcode = bash(
                str(os.path.join(repo_path, 'src', 'cache.sh')) + 
                ' local ' +
                " -s '{}' ".format(sif_cache) +
                " -i '{}' ".format(','.join(pull)) + 
                " -t '{0}/.{1}/.singularity/' ".format(sif_cache, username)
            )
            # Check exitcode of caching script 
            if exitcode != 0:
                fatal('Fatal: Failed to pull all containers. Please try again!')
            print('Done: sucessfully pulled all software containers!')


def parsed_arguments(name, description):
    """Parses user-provided command-line arguments. Requires argparse and textwrap
    package. argparse was added to standard lib in python 3.5 and textwrap was added
    in python 3.5. To create custom help formatting for subparsers a docstring is
    used create the help message for required options. argparse does not support named
    subparser groups, which is normally what would be used to accomphish this reformatting.
    As so, the help message for require options must be suppressed. If a new required arg
    is added to a subparser, it must be added to the docstring and the usage statement
    also must be updated.
    @param name <str>:
        Name of the pipeline or command-line tool 
    @param description <str>:
        Short description of pipeline or command-line tool 
    """

    # Add styled name and description
    c = Colors
    styled_name = "{0}{1}{2}genome-{3}seek{4}".format(c.bold, c.bg_black, c.white, c.cyan, c.end)
    description = "{0}{1}{2}".format(c.bold, description, c.end)

    # Create a top-level parser
    parser = argparse.ArgumentParser(description = '{}: {}'.format(styled_name, description))

    # Adding Verison information
    parser.add_argument('--version', action = 'version', version='%(prog)s {}'.format(__version__))

    # Create sub-command parser
    subparsers = parser.add_subparsers(help='List of available sub-commands')

    # Sub-parser for the "run" sub-command
    # Grouped sub-parser arguments are currently 
    # not supported: https://bugs.python.org/issue9341
    # Here is a work around to create more useful help message for named
    # options that are required! Please note: if a required arg is added the
    # description below should be updated (i.e. update usage and add new option)
    required_run_options = textwrap.dedent("""\
        {0}: {1}

        {3}{4}Synopsis:{5}
          $ {2} run [--help] \\
                [--mode {{slurm,uge,local}}] [--job-name JOB_NAME] [--batch-id BATCH_ID] \\
                [--call-cnv] [--call-sv] [--call-hla] [--call-somatic] [--gatk-germline] \\
                [--open-cravat] [--oc-annotators OC_ANNOTATORS] [--oc-modules OC_MODULES] \\
                [--pairs PAIRS] [--pon PANEL_OF_NORMALS] [--wes-mode] [--wes-bed WES_BED] \\
                [--skip-qc] [--tmp-dir TMP_DIR] [--silent] [--sif-cache SIF_CACHE] \\
                [--singularity-cache SINGULARITY_CACHE] [--use-gpus] \\
                [--resource-bundle RESOURCE_BUNDLE] \\
                [--dry-run] [--threads THREADS] \\
                --input INPUT [INPUT ...] \\
                --output OUTPUT

        Optional arguments are shown in square brackets above.

        {3}{4}Description:{5}
          Run the whole genome clinical sequencing pipeline with your data. Please
        provide a space seperated list of FastQ (globbing is supported) and an output 
        directory to store results. Optionally call structural variation (SV), copy 
        number variation (CNV), and more! 

        {3}{4}Required arguments:{5}
          --input INPUT [INPUT ...]
                                Input FastQ file(s) to process. The pipeline does NOT
                                support single-end data. FastQ files for one or more  
                                samples can be provided. Multiple input FastQ files
                                should be seperated by a space. Globbing for multiple
                                file is also supported.
                                  Example: --input .tests/*.R?.fastq.gz
          --output OUTPUT
                                Path to an output directory. This location is where
                                the pipeline will create all of its output files, also
                                known as the pipeline's working directory. If the user
                                provided working directory has not been initialized,
                                it will be created automatically.
                                  Example: --output /data/$USER/output

        {3}{4}Analysis options:{5}
          --call-cnv            Call copy number variation (CNV).
                                  Example: --call-cnv
    
          --call-sv             Call structural variation (SV).
                                  Example: --call-sv

          --call-hla            Call HLA types.
                                  Example: --call-hla

          --call-somatic        Call somatic variants. By default, the pipeline will 
                                perform somatic variant calling for each sample in a 
                                tumor-only mode; however, if a tumor - normal pairs 
                                file is provided then the pipeline will call somatic 
                                variants using its matched normal. See the Somatic 
                                Options section below for more information about the
                                somatic variant calling pipeline and its options.
                                  Example: --call-somatic

          --gatk-germline       Call germline variants using GATK4 best pratices.
                                By default, the pipeline will call germline variants
                                using deepvariant. If this option is provided, the
                                pipeline will additionally call germline variants
                                using GATK4's set of best practices.
                                  Example: --gatk-germline
        
          --open-cravat         Run OpenCRAVAT to annotate variants. See Annotation
                                Options for more information about what modules are
                                included by default and how to include more modules. 
                                  Example: --open-cravat
          
          --skip-qc             Skips over quality control steps. When this option
                                is provided the pipeline will not run any of its QC
                                related steps. This means that only data processing 
                                steps will run like trimming, alignment, and variant 
                                calling. It is worth noting that we do not recommend 
                                skipping over QC; however, in certain scenarios it 
                                may make sense. This option is useful for testing
                                changes to data processing steps or when evaluating 
                                the overall accuracy and precision of the pipeline.
                                  Example: --skip-qc
          
          --wes-mode            Run the whole exome pipeline. By default, the whole
                                genome sequencing pipeline is run. This option allows
                                a user to process and analyze whole exome sequencing 
                                data. Please note when this mode is enabled, a sub-
                                set of the WGS rules will run. Please see the option
                                below for more information about providing a custom
                                exome targets BED file.
                                  Example: --wes-mode

          --wes-bed WES_BED
                                Path to exome targets BED file. This file can be
                                obtained from the manufacturer of the target capture
                                kit that was used. By default, a set of BED files 
                                generated from GENCODE's exon annotation for protein
                                coding gene's exon is used.
                                  Example: --wes-bed gencode_v44_protein-coding_exons.bed
        
          --batch-id BATCH_ID   
                                Unique identifer for a batch of samples. A batch 
                                identifer should be a string containing no spaces.
                                If a batch identifer is not provided, one will be
                                generated by taking the MD5 checksum of the sorted
                                list of input file names. The batch identifer is 
                                added to some output file names. This allows the 
                                pipeline to be run with extra samples without
                                overriding some of results of a previous run.
                                  Example: --batch-id WGS_2022-04-19

        {3}{4}Somatic options:{5}
          --pairs PAIRS         Tumor-normal pairs file. The tumor-normal pairs file 
                                is used to pair a tumor sample with its match normal 
                                sample. This file should only be provided when calling 
                                somatic variants. Please see the option above for more 
                                info about calling somatic mutations. By default the 
                                --call-somatic option will call somatic variants for 
                                each sample in a tumor-only mode. This tab delimited 
                                file contains two columns with the names of tumor and 
                                normal pairs, one per line. The header of the file 
                                needs to be Tumor for the tumor column and Normal for 
                                the normal column. The base name of each sample should 
                                be listed in the pairs file. The base name of a given 
                                sample can be determined by removing the extension 
                                from the sample's R1 FastQ file, eg '.R1.fastq.gz'.
                                  Contents of example pairs file:
                                    Tumor    Normal
                                    Sample4_CRL1622_S31  Sample10_ARK1_S37
                                    Sample4_CRL1622_S31  Sample11_ACI_158_S38
                                  Example: --pairs /data/$USER/pairs.tsv
          --pon PANEL_OF_NORMALS
                                Panel of normals (PON). A VCF file of containing
                                sites observed in normal tissue. Normal in this
                                context refers to samples derived from healthy
                                tissue that is NOT believed to have any somatic
                                alterations. By default, the pipeline will use a
                                PON included with its resource bundle. You can
                                provide your own PON with this option. The PON 
                                should be gzipped, AND there should be a tabix 
                                index for the PON in the same directory.
                                  Example: --pon 1000g_pon.hg38.vcf.gz

        {3}{4}Anotation options:{5}
        {6}@OpenCRAVAT{5}
          --oc-annotators OC_ANNOTATORS 
                                List of OpenCRAVAT annotators to use. Note that  
                                one or more annotators can be provided. Multiple 
                                annotators should be seperate with a space. By 
                                default the pipeline will annotate variants using 
                                modules in the resource bundle; however, a custom
                                module installation path can be defined using the 
                                --oc-modules option. The default list of OpenCRAVAT
                                annotators can be here: 'config/oc_annotators.cfg'.
                                  Example: --oc-annotators dann dann_coding
          --oc-modules OC_MODULES 
                                Set path to OpenCRAVAT's modules directory. This 
                                option overrides the default path any installed
                                modules for OpenCRAVAT. Use this option if you 
                                have your own custom installation of OpenCRAVAT
                                modules.
                                  Example: --oc-modules /data/$USER/CRAVAT/modules

        {3}{4}Orchestration options:{5}
          --mode {{slurm,uge,local}}  
                                Method of execution. Defines the mode of execution. 
                                Vaild options for this mode include: local, slurm, 
                                or uge. Additional modes of exection are coming soon, 
                                default: slurm.
                                Here is a brief description of each mode:
                                  • slurm: uses slurm execution backend. This method 
                                will submit jobs to a  cluster  using sbatch. It is 
                                recommended running the pipeline in this mode as it 
                                will be significantly faster.
                                  • uge: uses UGE execution backend. This method will
                                submit jobs to a cluster using qsub. Please set the 
                                mode to uge when running the pipeline on LOCUS. 
                                   • local: uses local method of execution. local runs 
                                will run serially on compute instance. This is useful 
                                for testing, debugging, or when a users does not have
                                access to a  high  performance  computing environment.
                                If this option is not provided, it will default to a 
                                slurm mode of execution. 
                                  Example: --mode slurm
          --job-name JOB_NAME   
                                Overrides the name of the pipeline's master job. When 
                                submitting the pipeline to a jobscheduler, this option
                                overrides the default name of the master job. This can 
                                be useful for tracking the progress or status of a run, 
                                default: pl:{2}.
                                  Example: --job-name {2}_03-14.1592
          --dry-run             
                                Does not execute anything. Only displays what steps in
                                the pipeline remain or will be run.
                                  Example: --dry-run
          --silent              
                                Silence standard output. This will reduces the amount 
                                of information displayed to standard  output  when the 
                                master job is submitted to the job scheduler. Only the 
                                job id of the master job is returned.
                                  Example: --silent
          --singularity-cache SINGULARITY_CACHE
                                Overrides the $SINGULARITY_CACHEDIR variable. Images
                                from remote registries are cached locally on the file
                                system. By default, the singularity cache is set to:
                                '/path/to/output/directory/.singularity/'. Please note
                                that this cache cannot be shared across users.
                                  Example: --singularity-cache /data/$USER
          --sif-cache SIF_CACHE
                                Path where a local cache of SIFs are stored. This cache
                                can be shared across users if permissions are properly
                                setup. If a SIF does not exist in the SIF cache, the
                                image will be pulled from Dockerhub. {2} cache
                                sub command can be used to create a local SIF cache. 
                                Please see {2} cache for more information.
                                   Example: --sif-cache /data/$USER/sifs/

          --resource-bundle RESOURCE_BUNDLE
                                Path to a resource bundle. Only provide this option
                                if you are running the pipeline outside of Biowulf. If 
                                you are running the pipeline on Biowulf, the pipeline 
                                will automatically resolve the correct path to any 
                                references files. The resource bundle contains the set
                                of required reference files for processing any data. 
                                  Example: --resource-bundle /path/to/refs/genome-seek
                                           
          --tmp-dir TMP_DIR     
                                Path on the file system for writing temporary output 
                                files. By default, the temporary directory is set to 
                                '/lscratch/$SLURM_JOBID' for backwards compatibility 
                                with the NIH's Biowulf cluster; however, if you are 
                                running the pipeline on another cluster, this option 
                                will need to be specified. Ideally, this path should 
                                point to a dedicated location on the filesystem for 
                                writing tmp files. On many systems, this location is 
                                set to somewhere in /scratch. If you need to inject a 
                                variable into this string that should NOT be expanded, 
                                please quote this options value in single quotes. 
                                  Example: --tmp-dir '/scratch/$USER/'
          --threads THREADS     
                                Max number of threads for local processes. It is
                                recommended setting this vaule to the maximum number 
                                of CPUs available on the host machine, default: 2.
                                  Example: --threads: 16
          --use-gpus
                                Utilize GPU computing for any suppported tools.
                                Several steps in the pipeline can be accelerated 
                                with GPUs. This feature is only supported on 
                                Biowulf at the NIH.
                                  Example: --use-gpus
        
        {3}{4}Misc Options:{5}
          -h, --help            Show usage information, help message, and exit.
                                  Example: --help
        """.format(styled_name, description, name, c.bold, c.url, c.end, c.italic))

    # Display example usage in epilog
    run_epilog = textwrap.dedent("""\
        {2}{3}Example:{4}
          # Step 1.) Grab an interactive node,
          # do not run on head node!
          srun -N 1 -n 1 --time=1:00:00 --mem=8gb  --cpus-per-task=2 --pty bash
          module purge
          module load singularity snakemake

          # Step 2A.) Dry-run the pipeline
          ./{0} run --input .tests/*.R?.fastq.gz \\
                         --output /data/$USER/output \\
                         --mode slurm \\
                         --dry-run

          # Step 2B.) Run the {0} pipeline
          # The slurm mode will submit jobs to 
          # the cluster. It is recommended running 
          # the pipeline in this mode.
          ./{0} run --input .tests/*.R?.fastq.gz \\
                         --output /data/$USER/output \\
                         --mode slurm

        {2}{3}Version:{4}
          {1}
        """.format(name, __version__, c.bold, c.url, c.end))

    # Supressing help message of required args to overcome no sub-parser named groups
    subparser_run = subparsers.add_parser('run',
        help = 'Run the {} pipeline with input files.'.format(name),
        usage = argparse.SUPPRESS,
        formatter_class=argparse.RawDescriptionHelpFormatter,
        description = required_run_options,
        epilog  = run_epilog,
        add_help=False
    )

    # Required Arguments
    # Input FastQ files
    subparser_run.add_argument(
        '--input',
        # Check if the file exists and if it is readable
        type = lambda file: permissions(parser, file, os.R_OK),
        required = True,
        nargs = '+',
        help = argparse.SUPPRESS
    )

    # Output Directory, i.e 
    # working directory
    subparser_run.add_argument(
        '--output',
        type = lambda option: os.path.abspath(os.path.expanduser(option)),
        required = True,
        help = argparse.SUPPRESS
    )

    # Optional Arguments
    # Add custom help message
    subparser_run.add_argument(
        '-h', '--help', 
        action='help', 
        help=argparse.SUPPRESS
    )
    
    # Analysis options
    # Call CNVs
    subparser_run.add_argument(
        '--call-cnv',
        action = 'store_true',
        required = False,
        default = False,
        help = argparse.SUPPRESS
    )

    # Call SVs
    subparser_run.add_argument(
        '--call-sv',
        action = 'store_true',
        required = False,
        default = False,
        help = argparse.SUPPRESS
    )

    # Call HLA
    subparser_run.add_argument(
        '--call-hla',
        action = 'store_true',
        required = False,
        default = False,
        help = argparse.SUPPRESS
    )

    # Call somatic variants
    subparser_run.add_argument(
        '--call-somatic',
        action = 'store_true',
        required = False,
        default = False,
        help = argparse.SUPPRESS
    )

    # Call variants using GATK4's
    # set of best practices for
    # germline variants
    subparser_run.add_argument(
        '--gatk-germline',
        action = 'store_true',
        required = False,
        default = False,
        help = argparse.SUPPRESS
    )

    # Skip QC
    subparser_run.add_argument(
        '--skip-qc',
        action = 'store_true',
        required = False,
        default = False,
        help = argparse.SUPPRESS
    )

    # Run the WES pipeline
    subparser_run.add_argument(
        '--wes-mode',
        action = 'store_true',
        required = False,
        default = False,
        help = argparse.SUPPRESS
    )

    # WES capture targets BED file
    subparser_run.add_argument(
        '--wes-bed',
        # Check if the file exists and if it is readable
        type = lambda file: permissions(parser, file, os.R_OK),
        required = False,
        default = None,
        help = argparse.SUPPRESS
    )

    # Tumor-normal pairs file
    subparser_run.add_argument(
        '--pairs',
        # Check if the file exists and if it is readable
        type = lambda file: permissions(parser, file, os.R_OK),
        required = False,
        default = None,
        help = argparse.SUPPRESS
    )

    # Custom Panel of Normals
    subparser_run.add_argument(
        '--pon',
        # Check if the file exists and if it is readable AND
        # check if its tabix index exist and it is readable
        type = lambda file: re.sub(
            '\.tbi$',
            '',
            permissions(
                parser, 
                "{}.tbi".format(permissions(parser, file, os.R_OK)), 
                os.R_OK
            )
        ),
        required = False,
        default = None,
        help = argparse.SUPPRESS
    )

    # Run OpenCRAVAT
    subparser_run.add_argument(
        '--open-cravat',
        action = 'store_true',
        required = False,
        default = False,
        help = argparse.SUPPRESS
    )

    # OpenCRAVAT modules directory
    subparser_run.add_argument(
        '--oc-modules',
        # Convert to absolute path and
        # check if user can read from 
        # directory.
        type = lambda option: permissions(
            parser, 
            os.path.abspath(os.path.expanduser(option)), 
            os.R_OK) if option else '',
        required = False,
        default = '',
        help = argparse.SUPPRESS
    )

    # OpenCRAVAT Variant Annotators
    subparser_run.add_argument(
        '--oc-annotators',
        type = str,
        nargs = '+',
        required = False,
        # Set default from config/oc_annotators.txt
        # Skips over any empty lines and any lines 
        # starting with comments, i.e. '#' 
        default = [
            l.strip() for l in open(
                os.path.join(__home__, 'config', 'oc_annotators.cfg'), 
                'r'
            ).readlines() if not l.strip().startswith('#') and l.strip()
        ],
        help = argparse.SUPPRESS
    )

    # Batch Identifer
    subparser_run.add_argument(
        '--batch-id',
        type = lambda option: option.replace(' ', '_'),
        required = False,
        default = '',
        help = argparse.SUPPRESS
    )

    # Execution Method, run locally 
    # on a compute node or submit to 
    # a supported job scheduler, etc.
    subparser_run.add_argument(
        '--mode',
        type = str,
        required = False,
        default = "slurm",
        choices = ['slurm', 'uge', 'local'],
        help = argparse.SUPPRESS
    )

    # Name of master job
    subparser_run.add_argument(
        '--job-name',
        type = str,
        required = False,
        default = 'pl:{}'.format(name),
        help = argparse.SUPPRESS
    )

    # Dry-run
    # Do not execute the workflow, 
    # prints what steps remain
    subparser_run.add_argument(
        '--dry-run',
        action = 'store_true',
        required = False,
        default = False,
        help = argparse.SUPPRESS
    )

    # Silent output mode
    subparser_run.add_argument(
        '--silent',
        action = 'store_true',
        required = False,
        default = False,
        help = argparse.SUPPRESS
    )
    
    # Singularity cache directory,
    # default uses output directory
    subparser_run.add_argument(
        '--singularity-cache',
        type = lambda option: check_cache(parser, os.path.abspath(os.path.expanduser(option))),
        required = False,
        help = argparse.SUPPRESS
    )
    
    # Local SIF cache directory, 
    # default pull from Dockerhub
    subparser_run.add_argument(
        '--sif-cache',
        type = lambda option: os.path.abspath(os.path.expanduser(option)),
        required = False,
        help = argparse.SUPPRESS
    )

    # Output Directory of downloaded resource bundle:
    #  - Biowulf: /data/OpenOmics/references/genome-seek    (default)
    #  - LOCUS:   /hpcdata/openomics/references/genome-seek
    #  - BigSky:  /data/openomics/references/genome-seek
    subparser_run.add_argument(
        '--resource-bundle',
        type = lambda option: os.path.abspath(os.path.expanduser(option)),
        required = False,
        default = None,
        help = argparse.SUPPRESS
    )

    # Base directory to write 
    # temporary/intermediate files 
    subparser_run.add_argument(
        '--tmp-dir',
        type = str,
        required = False,
        default = '/lscratch/$SLURM_JOBID/',
        help = argparse.SUPPRESS
    )

    # Number of threads for the 
    # pipeline's main proceess
    # This is only applicable for 
    # local rules or when running 
    # in local mode.
    subparser_run.add_argument(
        '--threads',
        type = int,
        required = False,
        default = 2,
        help = argparse.SUPPRESS
    )

    # Use GPU computing to speed
    # up any time consuming steps,
    # this is only supported for 
    # a few tools.
    subparser_run.add_argument(
        '--use-gpus',
        action = 'store_true',
        required = False,
        default = False,
        help = argparse.SUPPRESS
    )

    # Sub-parser for the "unlock" sub-command
    # Grouped sub-parser arguments are currently 
    # not supported: https://bugs.python.org/issue9341
    # Here is a work around to create more useful help message for named
    # options that are required! Please note: if a required arg is added the
    # description below should be updated (i.e. update usage and add new option)
    required_unlock_options = textwrap.dedent("""\
        {0}: {1}

        {3}{4}Synopsis:{5} 
          $ {2} unlock [-h] --output OUTPUT
    
        Optional arguments are shown in square brackets above.

        {3}{4}Description:{5}
          If the pipeline fails ungracefully, it maybe required to unlock
        the working directory before proceeding again. Please verify that 
        the pipeline is not running before running this command. If the 
        pipeline is still running, the workflow manager will report the 
        working directory is locked. This is normal behavior. Do NOT run 
        this command if the pipeline is still running.

        {3}{4}Required arguments:{5}
          --output OUTPUT       Path to a previous run's output directory
                                to unlock. This will remove a lock on the 
                                working directory. Please verify that the 
                                pipeline is not running before running 
                                this command.
                                  Example: --output /data/$USER/output

        {3}{4}Misc Options:{5}
          -h, --help            Show usage information, help message, 
                                and exit.
                                  Example: --help
        """.format(styled_name, description, name, c.bold, c.url, c.end))

    # Display example usage in epilog
    unlock_epilog = textwrap.dedent("""\
        {2}{3}Example:{4}
          # Unlock output directory of pipeline
          {0} unlock --output /data/$USER/output

        {2}{3}Version:{4}
          {1}
        """.format(name, __version__, c.bold, c.url, c.end))

    # Supressing help message of required args to overcome no sub-parser named groups
    subparser_unlock = subparsers.add_parser(
        'unlock',
        help = 'Unlocks a previous runs output directory.',
        usage = argparse.SUPPRESS,
        formatter_class=argparse.RawDescriptionHelpFormatter,
        description = required_unlock_options,
        epilog = unlock_epilog,
        add_help = False
    )

    # Required Arguments
    # Output Directory (analysis working directory)
    subparser_unlock.add_argument(
        '--output',
        type = str,
        required = True,
        help = argparse.SUPPRESS
    )

    # Add custom help message
    subparser_unlock.add_argument(
        '-h', '--help', 
        action='help', 
        help=argparse.SUPPRESS
    )


    # Sub-parser for the "cache" sub-command
    # Grouped sub-parser arguments are 
    # not supported: https://bugs.python.org/issue9341
    # Here is a work around to create more useful help message for named
    # options that are required! Please note: if a required arg is added the
    # description below should be updated (i.e. update usage and add new option)
    required_cache_options = textwrap.dedent("""\
        {0}: {1}

        {3}{4}Synopsis:{5} Cache software containers locally.
          $ {2} cache [-h] [--dry-run] \\
                  --sif-cache SIF_CACHE

        Create a local cache of software dependencies hosted on DockerHub.
        These containers are normally pulled onto the filesystem when the
        pipeline runs; however, due to network issues or DockerHub pull
        rate limits, it may make sense to pull the resources once so a
        shared cache can be created. It is worth noting that a singularity
        cache cannot normally be shared across users. Singularity strictly
        enforces that a cache is owned by the user. To get around this
        issue, the cache subcommand can be used to create local SIFs on
        the filesystem from images on DockerHub.

        {3}{4}Required arguments:{5}
          --sif-cache SIF_CACHE
                                Path where a local cache of SIFs will be 
                                stored. Images defined in containers.json
                                will be pulled into the local filesystem. 
                                The path provided to this option can be 
                                passed to the --sif-cache option of the
                                run sub command. Please see {2} run 
                                sub command for more information.
                                  Example: --sif-cache /data/$USER/cache
        
        {3}{4}Orchestration options:{5}
          --dry-run             Does not execute anything. Only displays 
                                what remote resources would be pulled.
                                  Example: --dry-run

        {3}{4}Misc Options:{5}
          -h, --help            Show usage information, help message, 
                                and exits.
                                  Example: --help

        """.format(styled_name, description, name, c.bold, c.url, c.end))

    # Display example usage in epilog
    cache_epilog = textwrap.dedent("""\
        {2}Example:{3}
          # Cache software containers of pipeline
          {0} cache --sif-cache /data/$USER/cache

        {2}Version:{3}
          {1}
        """.format(name, __version__, c.bold, c.end))

    # Supressing help message of required args 
    # to overcome no sub-parser named groups
    subparser_cache = subparsers.add_parser(
        'cache',
        help = 'Cache software containers locally.',
        usage = argparse.SUPPRESS,
        formatter_class=argparse.RawDescriptionHelpFormatter,
        description = required_cache_options,
        epilog = cache_epilog,
        add_help = False
    )

    # Required Arguments
    # Output Directory (analysis working directory)
    subparser_cache.add_argument(
        '--sif-cache',
        type = lambda option: os.path.abspath(os.path.expanduser(option)),
        required = True,
        help = argparse.SUPPRESS
    )

    # Optional Arguments
    # Dry-run cache command (do not pull any remote resources)
    subparser_cache.add_argument(
        '--dry-run',
        action = 'store_true',
        required = False,
        default = False,
        help=argparse.SUPPRESS
    )

    # Add custom help message
    subparser_cache.add_argument(
        '-h', '--help', 
        action='help', 
        help=argparse.SUPPRESS
    )

    # Define handlers for each sub-parser
    subparser_run.set_defaults(func = run)
    subparser_unlock.set_defaults(func = unlock)
    subparser_cache.set_defaults(func = cache)

    # Parse command-line args
    args = parser.parse_args()
    return args


def main():

    # Sanity check for usage
    if len(sys.argv) == 1:
        # Nothing was provided
        fatal('Invalid usage: {} [-h] [--version] ...'.format(_name))

    # Collect args for sub-command
    args = parsed_arguments(
        name = _name,
        description = _description
    )

    # Display version information
    err('{} ({})'.format(_name, __version__))

    # Mediator method to call sub-command's set handler function
    args.func(args)


if __name__ == '__main__':
    main()