src/toil_lib/tools/variant_annotation.py
def gatk_genotype_gvcfs(job,
                        gvcfs,
                        ref, fai, ref_dict,
                        annotations=None,
                        emit_threshold=10.0, call_threshold=30.0,
                        unsafe_mode=False):
    """
    Runs GenotypeGVCFs on one or more gVCFs generated by HaplotypeCaller.

    :param JobFunctionWrappingJob job: passed automatically by Toil
    :param dict gvcfs: Dictionary of GVCF FileStoreIDs {sample identifier: FileStoreID}
    :param str ref: FileStoreID for the reference genome fasta file
    :param str fai: FileStoreID for the reference genome index file
    :param str ref_dict: FileStoreID for the reference genome sequence dictionary
    :param list[str] annotations: Optional list of GATK variant annotations. Default: None.
    :param float emit_threshold: Minimum phred-scale confidence threshold for
                                 a variant to be emitted. GATK default: 10.0
    :param float call_threshold: Minimum phred-scale confidence threshold for
                                 a variant to be called. GATK default: 30.0
    :param bool unsafe_mode: If True, runs gatk UNSAFE mode: "-U ALLOW_SEQ_DICT_INCOMPATIBILITY"
    :return: VCF FileStoreID
    :rtype: str
    """
    inputs = {'genome.fa': ref,
              'genome.fa.fai': fai,
              'genome.dict': ref_dict}
    inputs.update(gvcfs)

    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))

    command = ['-T', 'GenotypeGVCFs',
               '-R', '/data/genome.fa',
               '--out', 'genotyped.vcf',
               '-stand_emit_conf', str(emit_threshold),
               '-stand_call_conf', str(call_threshold)]

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

    # Include all GVCFs for joint genotyping
    for uuid in gvcfs.keys():
        command.extend(['--variant', os.path.join('/data', uuid)])

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

    job.fileStore.logToMaster('Running GATK GenotypeGVCFs\n'
                              'Emit threshold: {emit_threshold}\n'
                              'Call threshold: {call_threshold}\n\n'
                              'Annotations:\n{annotations}\n\n'
                              'Samples:\n{samples}\n'.format(emit_threshold=emit_threshold,
                                                             call_threshold=call_threshold,
                                                             annotations='\n'.join(annotations) if annotations else '',
                                                             samples='\n'.join(gvcfs.keys())))

    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)

    return job.fileStore.writeGlobalFile(os.path.join(work_dir, 'genotyped.vcf'))
def gatk_genotype_gvcfs(job,
                        gvcfs,
                        ref, fai, ref_dict,
                        annotations=None,
                        emit_threshold=10.0, call_threshold=30.0,
                        unsafe_mode=False):
    """
    Runs GenotypeGVCFs on one or more gVCFs generated by HaplotypeCaller.

    :param JobFunctionWrappingJob job: passed automatically by Toil
    :param dict gvcfs: Dictionary of GVCF FileStoreIDs {sample identifier: FileStoreID}
    :param str ref: FileStoreID for the reference genome fasta file
    :param str fai: FileStoreID for the reference genome index file
    :param str ref_dict: FileStoreID for the reference genome sequence dictionary
    :param list[str] annotations: Optional list of GATK variant annotations. Default: None.
    :param float emit_threshold: Minimum phred-scale confidence threshold for
                                 a variant to be emitted. GATK default: 10.0
    :param float call_threshold: Minimum phred-scale confidence threshold for
                                 a variant to be called. GATK default: 30.0
    :param bool unsafe_mode: If True, runs gatk UNSAFE mode: "-U ALLOW_SEQ_DICT_INCOMPATIBILITY"
    :return: VCF FileStoreID
    :rtype: str
    """
    inputs = {'genome.fa': ref,
              'genome.fa.fai': fai,
              'genome.dict': ref_dict}
    inputs.update(gvcfs)

    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))

    command = ['-T', 'GenotypeGVCFs',
               '-R', '/data/genome.fa',
               '--out', 'genotyped.vcf',
               '-stand_emit_conf', str(emit_threshold),
               '-stand_call_conf', str(call_threshold)]

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

    # Include all GVCFs for joint genotyping
    for uuid in gvcfs.keys():
        command.extend(['--variant', os.path.join('/data', uuid)])

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

    job.fileStore.logToMaster('Running GATK GenotypeGVCFs\n'
                              'Emit threshold: {emit_threshold}\n'
                              'Call threshold: {call_threshold}\n\n'
                              'Annotations:\n{annotations}\n\n'
                              'Samples:\n{samples}\n'.format(emit_threshold=emit_threshold,
                                                             call_threshold=call_threshold,
                                                             annotations='\n'.join(annotations) if annotations else '',
                                                             samples='\n'.join(gvcfs.keys())))

    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)

    return job.fileStore.writeGlobalFile(os.path.join(work_dir, 'genotyped.vcf'))
Runs GenotypeGVCFs on one or more gVCFs generated by HaplotypeCaller.

    :param JobFunctionWrappingJob job: passed automatically by Toil
    :param dict gvcfs: Dictionary of GVCF FileStoreIDs {sample identifier: FileStoreID}
    :param str ref: FileStoreID for the reference genome fasta file
    :param str fai: FileStoreID for the reference genome index file
    :param str ref_dict: FileStoreID for the reference genome sequence dictionary
    :param list[str] annotations: Optional list of GATK variant annotations. Default: None.
    :param float emit_threshold: Minimum phred-scale confidence threshold for
                                 a variant to be emitted. GATK default: 10.0
    :param float call_threshold: Minimum phred-scale confidence threshold for
                                 a variant to be called. GATK default: 30.0
    :param bool unsafe_mode: If True, runs gatk UNSAFE mode: "-U ALLOW_SEQ_DICT_INCOMPATIBILITY"
    :return: VCF FileStoreID
    :rtype: str