src/toil_lib/tools/variant_manipulation.py
def gatk_variant_recalibrator(job,
                              mode,
                              vcf,
                              ref_fasta, ref_fai, ref_dict,
                              annotations,
                              hapmap=None, omni=None, phase=None, dbsnp=None, mills=None,
                              max_gaussians=4,
                              unsafe_mode=False):
    """
    Runs either SNP or INDEL variant quality score recalibration using GATK VariantRecalibrator. Because the VQSR method
    models SNPs and INDELs differently, VQSR must be run separately for these variant types.

    :param JobFunctionWrappingJob job: passed automatically by Toil
    :param str mode: Determines variant recalibration mode (SNP or INDEL)
    :param str vcf: FileStoreID for input VCF file
    :param str ref_fasta: FileStoreID for reference genome fasta
    :param str ref_fai: FileStoreID for reference genome index file
    :param str ref_dict: FileStoreID for reference genome sequence dictionary file
    :param list[str] annotations: List of GATK variant annotations to filter on
    :param str hapmap: FileStoreID for HapMap resource file, required for SNP VQSR
    :param str omni: FileStoreID for Omni resource file, required for SNP VQSR
    :param str phase: FileStoreID for 1000G resource file, required for SNP VQSR
    :param str dbsnp: FilesStoreID for dbSNP resource file, required for SNP and INDEL VQSR
    :param str mills: FileStoreID for Mills resource file, required for INDEL VQSR
    :param int max_gaussians: Number of Gaussians used during training, default is 4
    :param bool unsafe_mode: If True, runs gatk UNSAFE mode: "-U ALLOW_SEQ_DICT_INCOMPATIBILITY"
    :return: FileStoreID for the variant recalibration table, tranche file, and plots file
    :rtype: tuple
    """
    mode = mode.upper()

    inputs = {'genome.fa': ref_fasta,
              'genome.fa.fai': ref_fai,
              'genome.dict': ref_dict,
              'input.vcf': vcf}

    # Refer to GATK documentation for description of recommended parameters:
    # https://software.broadinstitute.org/gatk/documentation/article?id=1259
    # https://software.broadinstitute.org/gatk/documentation/article?id=2805

    # This base command includes parameters for both INDEL and SNP VQSR.
    command = ['-T', 'VariantRecalibrator',
               '-R', 'genome.fa',
               '-input', 'input.vcf',
               '-tranche', '100.0',
               '-tranche', '99.9',
               '-tranche', '99.0',
               '-tranche', '90.0',
               '--maxGaussians', str(max_gaussians),
               '-recalFile', 'output.recal',
               '-tranchesFile', 'output.tranches',
               '-rscriptFile', 'output.plots.R']

    # Parameters and resource files for SNP VQSR.
    if mode == 'SNP':
        command.extend(
            ['-resource:hapmap,known=false,training=true,truth=true,prior=15.0', 'hapmap.vcf',
             '-resource:omni,known=false,training=true,truth=true,prior=12.0', 'omni.vcf',
             '-resource:dbsnp,known=true,training=false,truth=false,prior=2.0', 'dbsnp.vcf',
             '-resource:1000G,known=false,training=true,truth=false,prior=10.0', '1000G.vcf',
             '-mode', 'SNP'])

        inputs['hapmap.vcf'] = hapmap
        inputs['omni.vcf'] = omni
        inputs['dbsnp.vcf'] = dbsnp
        inputs['1000G.vcf'] = phase

    # Parameters and resource files for INDEL VQSR
    elif mode == 'INDEL':
        command.extend(
            ['-resource:mills,known=false,training=true,truth=true,prior=12.0', 'mills.vcf',
             '-resource:dbsnp,known=true,training=false,truth=false,prior=2.0', 'dbsnp.vcf',
             '-mode', 'INDEL'])

        inputs['mills.vcf'] = mills
        inputs['dbsnp.vcf'] = dbsnp

    else:
        raise ValueError('Variant filter modes can be SNP or INDEL, got %s' % mode)

    for annotation in annotations:
        command.extend(['-an', annotation])

    if unsafe_mode:
        command.extend(['-U', 'ALLOW_SEQ_DICT_INCOMPATIBILITY'])

    # Delay reading in files until function is configured
    work_dir = job.fileStore.getLocalTempDir()
    for name, file_store_id in inputs.iteritems():
        job.fileStore.readGlobalFile(file_store_id, os.path.join(work_dir, name))

    job.fileStore.logToMaster('Running GATK VariantRecalibrator on {mode}s using the following annotations:\n'
                              '{annotations}'.format(mode=mode, annotations='\n'.join(annotations)))

    docker_parameters = ['--rm', 'log-driver', 'none',
                         '-e', 'JAVA_OPTS=-Djava.io.tmpdir=/data/ -Xmx{}'.format(job.memory)]
    dockerCall(job=job, workDir=work_dir,
               parameters=command,
               tool='quay.io/ucsc_cgl/gatk:3.5--dba6dae49156168a909c43330350c6161dc7ecc2',
               dockerParameters=docker_parameters)

    recal_id = job.fileStore.writeGlobalFile(os.path.join(work_dir, 'output.recal'))
    tranches_id = job.fileStore.writeGlobalFile(os.path.join(work_dir, 'output.tranches'))
    plots_id = job.fileStore.writeGlobalFile(os.path.join(work_dir, 'output.plots.R'))
    return recal_id, tranches_id, plots_id
def gatk_variant_recalibrator(job,
                              mode,
                              vcf,
                              ref_fasta, ref_fai, ref_dict,
                              annotations,
                              hapmap=None, omni=None, phase=None, dbsnp=None, mills=None,
                              max_gaussians=4,
                              unsafe_mode=False):
    """
    Runs either SNP or INDEL variant quality score recalibration using GATK VariantRecalibrator. Because the VQSR method
    models SNPs and INDELs differently, VQSR must be run separately for these variant types.

    :param JobFunctionWrappingJob job: passed automatically by Toil
    :param str mode: Determines variant recalibration mode (SNP or INDEL)
    :param str vcf: FileStoreID for input VCF file
    :param str ref_fasta: FileStoreID for reference genome fasta
    :param str ref_fai: FileStoreID for reference genome index file
    :param str ref_dict: FileStoreID for reference genome sequence dictionary file
    :param list[str] annotations: List of GATK variant annotations to filter on
    :param str hapmap: FileStoreID for HapMap resource file, required for SNP VQSR
    :param str omni: FileStoreID for Omni resource file, required for SNP VQSR
    :param str phase: FileStoreID for 1000G resource file, required for SNP VQSR
    :param str dbsnp: FilesStoreID for dbSNP resource file, required for SNP and INDEL VQSR
    :param str mills: FileStoreID for Mills resource file, required for INDEL VQSR
    :param int max_gaussians: Number of Gaussians used during training, default is 4
    :param bool unsafe_mode: If True, runs gatk UNSAFE mode: "-U ALLOW_SEQ_DICT_INCOMPATIBILITY"
    :return: FileStoreID for the variant recalibration table, tranche file, and plots file
    :rtype: tuple
    """
    mode = mode.upper()

    inputs = {'genome.fa': ref_fasta,
              'genome.fa.fai': ref_fai,
              'genome.dict': ref_dict,
              'input.vcf': vcf}

    # Refer to GATK documentation for description of recommended parameters:
    # https://software.broadinstitute.org/gatk/documentation/article?id=1259
    # https://software.broadinstitute.org/gatk/documentation/article?id=2805

    # This base command includes parameters for both INDEL and SNP VQSR.
    command = ['-T', 'VariantRecalibrator',
               '-R', 'genome.fa',
               '-input', 'input.vcf',
               '-tranche', '100.0',
               '-tranche', '99.9',
               '-tranche', '99.0',
               '-tranche', '90.0',
               '--maxGaussians', str(max_gaussians),
               '-recalFile', 'output.recal',
               '-tranchesFile', 'output.tranches',
               '-rscriptFile', 'output.plots.R']

    # Parameters and resource files for SNP VQSR.
    if mode == 'SNP':
        command.extend(
            ['-resource:hapmap,known=false,training=true,truth=true,prior=15.0', 'hapmap.vcf',
             '-resource:omni,known=false,training=true,truth=true,prior=12.0', 'omni.vcf',
             '-resource:dbsnp,known=true,training=false,truth=false,prior=2.0', 'dbsnp.vcf',
             '-resource:1000G,known=false,training=true,truth=false,prior=10.0', '1000G.vcf',
             '-mode', 'SNP'])

        inputs['hapmap.vcf'] = hapmap
        inputs['omni.vcf'] = omni
        inputs['dbsnp.vcf'] = dbsnp
        inputs['1000G.vcf'] = phase

    # Parameters and resource files for INDEL VQSR
    elif mode == 'INDEL':
        command.extend(
            ['-resource:mills,known=false,training=true,truth=true,prior=12.0', 'mills.vcf',
             '-resource:dbsnp,known=true,training=false,truth=false,prior=2.0', 'dbsnp.vcf',
             '-mode', 'INDEL'])

        inputs['mills.vcf'] = mills
        inputs['dbsnp.vcf'] = dbsnp

    else:
        raise ValueError('Variant filter modes can be SNP or INDEL, got %s' % mode)

    for annotation in annotations:
        command.extend(['-an', annotation])

    if unsafe_mode:
        command.extend(['-U', 'ALLOW_SEQ_DICT_INCOMPATIBILITY'])

    # Delay reading in files until function is configured
    work_dir = job.fileStore.getLocalTempDir()
    for name, file_store_id in inputs.iteritems():
        job.fileStore.readGlobalFile(file_store_id, os.path.join(work_dir, name))

    job.fileStore.logToMaster('Running GATK VariantRecalibrator on {mode}s using the following annotations:\n'
                              '{annotations}'.format(mode=mode, annotations='\n'.join(annotations)))

    docker_parameters = ['--rm', 'log-driver', 'none',
                         '-e', 'JAVA_OPTS=-Djava.io.tmpdir=/data/ -Xmx{}'.format(job.memory)]
    dockerCall(job=job, workDir=work_dir,
               parameters=command,
               tool='quay.io/ucsc_cgl/gatk:3.5--dba6dae49156168a909c43330350c6161dc7ecc2',
               dockerParameters=docker_parameters)

    recal_id = job.fileStore.writeGlobalFile(os.path.join(work_dir, 'output.recal'))
    tranches_id = job.fileStore.writeGlobalFile(os.path.join(work_dir, 'output.tranches'))
    plots_id = job.fileStore.writeGlobalFile(os.path.join(work_dir, 'output.plots.R'))
    return recal_id, tranches_id, plots_id
Runs either SNP or INDEL variant quality score recalibration using GATK VariantRecalibrator. Because the VQSR method
    models SNPs and INDELs differently, VQSR must be run separately for these variant types.

    :param JobFunctionWrappingJob job: passed automatically by Toil
    :param str mode: Determines variant recalibration mode (SNP or INDEL)
    :param str vcf: FileStoreID for input VCF file
    :param str ref_fasta: FileStoreID for reference genome fasta
    :param str ref_fai: FileStoreID for reference genome index file
    :param str ref_dict: FileStoreID for reference genome sequence dictionary file
    :param list[str] annotations: List of GATK variant annotations to filter on
    :param str hapmap: FileStoreID for HapMap resource file, required for SNP VQSR
    :param str omni: FileStoreID for Omni resource file, required for SNP VQSR
    :param str phase: FileStoreID for 1000G resource file, required for SNP VQSR
    :param str dbsnp: FilesStoreID for dbSNP resource file, required for SNP and INDEL VQSR
    :param str mills: FileStoreID for Mills resource file, required for INDEL VQSR
    :param int max_gaussians: Number of Gaussians used during training, default is 4
    :param bool unsafe_mode: If True, runs gatk UNSAFE mode: "-U ALLOW_SEQ_DICT_INCOMPATIBILITY"
    :return: FileStoreID for the variant recalibration table, tranche file, and plots file
    :rtype: tuple