src/toil_lib/tools/preprocessing.py
def run_gatk_preprocessing(job, bam, bai, ref, ref_dict, fai, g1k, mills, dbsnp, realign=False, unsafe=False):
    """
    GATK Preprocessing Pipeline
    0: Mark duplicates
    1: Create INDEL realignment intervals
    2: Realign INDELs
    3: Recalibrate base quality scores
    4: Apply base score recalibration

    :param JobFunctionWrappingJob job: passed automatically by Toil
    :param str bam: FileStoreID for BAM file
    :param str bai: FileStoreID for BAM index file
    :param str ref: FileStoreID for reference genome fasta file
    :param str ref_dict: FileStoreID for reference sequence dictionary file
    :param str fai: FileStoreID for reference fasta index file
    :param str g1k: FileStoreID for 1000 Genomes VCF file
    :param str mills: FileStoreID for Mills VCF file
    :param str dbsnp: FileStoreID for dbSNP VCF file
    :param bool realign: If True, then runs GATK INDEL realignment"
    :param bool unsafe: If True, runs GATK tools in UNSAFE mode: "-U ALLOW_SEQ_DICT_INCOMPATIBILITY"
    :return: FileStoreIDs for BAM and BAI files
    :rtype: tuple(str, str)
    """
    # The MarkDuplicates disk requirement depends on the input BAM and BAI files and the output
    # BAM and BAI files. The output BAM file is approximately the same size as the input BAM file.
    mdups_disk = PromisedRequirement(lambda bam_, bai_: 2 * (bam_.size + bai_.size), bam, bai)
    mdups = job.wrapJobFn(picard_mark_duplicates,
                          bam,
                          bai,
                          cores=job.cores,
                          disk=mdups_disk,
                          memory=job.memory)

    # Store input for BQSR
    bqsr_input_bam = mdups.rv(0)
    bqsr_input_bai = mdups.rv(1)

    # Get genome reference file sizes for calculating disk requirements
    genome_ref_size = ref.size + ref_dict.size + fai.size

    if realign:
        # Get INDEL resource file sizes and genome reference file sizes
        indel_ref_size = mills.size + g1k.size + genome_ref_size

        # The RealignerTargetCreator disk requirement depends on the input BAM/BAI files, the genome reference files,
        # and the output intervals file. The intervals file size is less than the reference file size, so estimate the
        # interval file size as the reference file size.
        realigner_target_disk = PromisedRequirement(lambda bam_, bai_, ref_size:
                                                    bam_.size + bai_.size + 2 * ref_size,
                                                    mdups.rv(0),
                                                    mdups.rv(1),
                                                    indel_ref_size)

        realigner_target = job.wrapJobFn(run_realigner_target_creator,
                                         mdups.rv(0),
                                         mdups.rv(1),
                                         ref, ref_dict, fai,
                                         g1k, mills,
                                         unsafe=unsafe,
                                         cores=1,  # RealignerTargetCreator is single threaded
                                         disk=realigner_target_disk,
                                         memory=job.memory)

        # The INDEL realignment disk requirement depends on the input BAM and BAI files, the intervals
        # file, the variant resource files, and the output BAM and BAI files. Here, we assume the
        # output BAM and BAI files are approximately the same size as the input BAM and BAI files.
        indel_realign_disk = PromisedRequirement(lambda bam_, bai_, intervals, ref_size:
                                                 2 * (bam_.size + bai_.size) + intervals.size + ref_size,
                                                 mdups.rv(0),
                                                 mdups.rv(1),
                                                 realigner_target.rv(),
                                                 indel_ref_size)

        indel_realign = job.wrapJobFn(run_indel_realignment,
                                      realigner_target.rv(),
                                      mdups.rv(0),
                                      mdups.rv(1),
                                      ref, ref_dict, fai,
                                      g1k, mills,
                                      unsafe=unsafe,
                                      cores=1,  # IndelRealigner is single threaded
                                      disk=indel_realign_disk,
                                      memory=job.memory)

        mdups.addChild(realigner_target)
        realigner_target.addChild(indel_realign)

        # Update input for BQSR using the realigned BAM files
        bqsr_input_bam = indel_realign.rv(0)
        bqsr_input_bai = indel_realign.rv(1)


    # Get size of BQSR databases and genome reference files
    bqsr_ref_size = dbsnp.size + mills.size + genome_ref_size

    # The BQSR disk requirement depends on the input BAM and BAI files, the reference files, and the output
    # recalibration table file. The recalibration table file size is less than the reference file sizes, so use
    # the reference file sizes to estimate the recalibration table file size.
    base_recal_disk = PromisedRequirement(lambda bam_, bai_, ref_size:
                                          bam_.size + bai_.size + 2 * ref_size,
                                          bqsr_input_bam,
                                          bqsr_input_bai,
                                          bqsr_ref_size)

    base_recal = job.wrapJobFn(run_base_recalibration,
                               bqsr_input_bam,
                               bqsr_input_bai,
                               ref, ref_dict, fai,
                               dbsnp, mills,
                               unsafe=unsafe,
                               cores=job.cores,
                               disk=base_recal_disk,
                               memory=job.memory)

    # The PrintReads disk requirement depends on the input BAM and BAI files, the recalibration table file, the
    # genome reference files, and the output BAM and BAI files. The output BAM and BAI files are approximately the
    # same size as the input BAM and BAI files.
    recalibrate_reads_disk = PromisedRequirement(lambda bam_, bai_, recal, ref_size:
                                                 2 * (bam_.size + bai_.size) + recal.size + ref_size,
                                                 bqsr_input_bam,
                                                 bqsr_input_bai,
                                                 base_recal.rv(),
                                                 genome_ref_size)

    recalibrate_reads = job.wrapJobFn(apply_bqsr_recalibration,
                                      base_recal.rv(),
                                      bqsr_input_bam,
                                      bqsr_input_bai,
                                      ref, ref_dict, fai,
                                      unsafe=unsafe,
                                      cores=job.cores,
                                      disk=recalibrate_reads_disk,
                                      memory=job.memory)

    job.addChild(mdups)
    mdups.addFollowOn(base_recal)
    base_recal.addChild(recalibrate_reads)
    return recalibrate_reads.rv(0), recalibrate_reads.rv(1)
def run_gatk_preprocessing(job, bam, bai, ref, ref_dict, fai, g1k, mills, dbsnp, realign=False, unsafe=False):
    """
    GATK Preprocessing Pipeline
    0: Mark duplicates
    1: Create INDEL realignment intervals
    2: Realign INDELs
    3: Recalibrate base quality scores
    4: Apply base score recalibration

    :param JobFunctionWrappingJob job: passed automatically by Toil
    :param str bam: FileStoreID for BAM file
    :param str bai: FileStoreID for BAM index file
    :param str ref: FileStoreID for reference genome fasta file
    :param str ref_dict: FileStoreID for reference sequence dictionary file
    :param str fai: FileStoreID for reference fasta index file
    :param str g1k: FileStoreID for 1000 Genomes VCF file
    :param str mills: FileStoreID for Mills VCF file
    :param str dbsnp: FileStoreID for dbSNP VCF file
    :param bool realign: If True, then runs GATK INDEL realignment"
    :param bool unsafe: If True, runs GATK tools in UNSAFE mode: "-U ALLOW_SEQ_DICT_INCOMPATIBILITY"
    :return: FileStoreIDs for BAM and BAI files
    :rtype: tuple(str, str)
    """
    # The MarkDuplicates disk requirement depends on the input BAM and BAI files and the output
    # BAM and BAI files. The output BAM file is approximately the same size as the input BAM file.
    mdups_disk = PromisedRequirement(lambda bam_, bai_: 2 * (bam_.size + bai_.size), bam, bai)
    mdups = job.wrapJobFn(picard_mark_duplicates,
                          bam,
                          bai,
                          cores=job.cores,
                          disk=mdups_disk,
                          memory=job.memory)

    # Store input for BQSR
    bqsr_input_bam = mdups.rv(0)
    bqsr_input_bai = mdups.rv(1)

    # Get genome reference file sizes for calculating disk requirements
    genome_ref_size = ref.size + ref_dict.size + fai.size

    if realign:
        # Get INDEL resource file sizes and genome reference file sizes
        indel_ref_size = mills.size + g1k.size + genome_ref_size

        # The RealignerTargetCreator disk requirement depends on the input BAM/BAI files, the genome reference files,
        # and the output intervals file. The intervals file size is less than the reference file size, so estimate the
        # interval file size as the reference file size.
        realigner_target_disk = PromisedRequirement(lambda bam_, bai_, ref_size:
                                                    bam_.size + bai_.size + 2 * ref_size,
                                                    mdups.rv(0),
                                                    mdups.rv(1),
                                                    indel_ref_size)

        realigner_target = job.wrapJobFn(run_realigner_target_creator,
                                         mdups.rv(0),
                                         mdups.rv(1),
                                         ref, ref_dict, fai,
                                         g1k, mills,
                                         unsafe=unsafe,
                                         cores=1,  # RealignerTargetCreator is single threaded
                                         disk=realigner_target_disk,
                                         memory=job.memory)

        # The INDEL realignment disk requirement depends on the input BAM and BAI files, the intervals
        # file, the variant resource files, and the output BAM and BAI files. Here, we assume the
        # output BAM and BAI files are approximately the same size as the input BAM and BAI files.
        indel_realign_disk = PromisedRequirement(lambda bam_, bai_, intervals, ref_size:
                                                 2 * (bam_.size + bai_.size) + intervals.size + ref_size,
                                                 mdups.rv(0),
                                                 mdups.rv(1),
                                                 realigner_target.rv(),
                                                 indel_ref_size)

        indel_realign = job.wrapJobFn(run_indel_realignment,
                                      realigner_target.rv(),
                                      mdups.rv(0),
                                      mdups.rv(1),
                                      ref, ref_dict, fai,
                                      g1k, mills,
                                      unsafe=unsafe,
                                      cores=1,  # IndelRealigner is single threaded
                                      disk=indel_realign_disk,
                                      memory=job.memory)

        mdups.addChild(realigner_target)
        realigner_target.addChild(indel_realign)

        # Update input for BQSR using the realigned BAM files
        bqsr_input_bam = indel_realign.rv(0)
        bqsr_input_bai = indel_realign.rv(1)


    # Get size of BQSR databases and genome reference files
    bqsr_ref_size = dbsnp.size + mills.size + genome_ref_size

    # The BQSR disk requirement depends on the input BAM and BAI files, the reference files, and the output
    # recalibration table file. The recalibration table file size is less than the reference file sizes, so use
    # the reference file sizes to estimate the recalibration table file size.
    base_recal_disk = PromisedRequirement(lambda bam_, bai_, ref_size:
                                          bam_.size + bai_.size + 2 * ref_size,
                                          bqsr_input_bam,
                                          bqsr_input_bai,
                                          bqsr_ref_size)

    base_recal = job.wrapJobFn(run_base_recalibration,
                               bqsr_input_bam,
                               bqsr_input_bai,
                               ref, ref_dict, fai,
                               dbsnp, mills,
                               unsafe=unsafe,
                               cores=job.cores,
                               disk=base_recal_disk,
                               memory=job.memory)

    # The PrintReads disk requirement depends on the input BAM and BAI files, the recalibration table file, the
    # genome reference files, and the output BAM and BAI files. The output BAM and BAI files are approximately the
    # same size as the input BAM and BAI files.
    recalibrate_reads_disk = PromisedRequirement(lambda bam_, bai_, recal, ref_size:
                                                 2 * (bam_.size + bai_.size) + recal.size + ref_size,
                                                 bqsr_input_bam,
                                                 bqsr_input_bai,
                                                 base_recal.rv(),
                                                 genome_ref_size)

    recalibrate_reads = job.wrapJobFn(apply_bqsr_recalibration,
                                      base_recal.rv(),
                                      bqsr_input_bam,
                                      bqsr_input_bai,
                                      ref, ref_dict, fai,
                                      unsafe=unsafe,
                                      cores=job.cores,
                                      disk=recalibrate_reads_disk,
                                      memory=job.memory)

    job.addChild(mdups)
    mdups.addFollowOn(base_recal)
    base_recal.addChild(recalibrate_reads)
    return recalibrate_reads.rv(0), recalibrate_reads.rv(1)
GATK Preprocessing Pipeline
    0: Mark duplicates
    1: Create INDEL realignment intervals
    2: Realign INDELs
    3: Recalibrate base quality scores
    4: Apply base score recalibration

    :param JobFunctionWrappingJob job: passed automatically by Toil
    :param str bam: FileStoreID for BAM file
    :param str bai: FileStoreID for BAM index file
    :param str ref: FileStoreID for reference genome fasta file
    :param str ref_dict: FileStoreID for reference sequence dictionary file
    :param str fai: FileStoreID for reference fasta index file
    :param str g1k: FileStoreID for 1000 Genomes VCF file
    :param str mills: FileStoreID for Mills VCF file
    :param str dbsnp: FileStoreID for dbSNP VCF file
    :param bool realign: If True, then runs GATK INDEL realignment"
    :param bool unsafe: If True, runs GATK tools in UNSAFE mode: "-U ALLOW_SEQ_DICT_INCOMPATIBILITY"
    :return: FileStoreIDs for BAM and BAI files
    :rtype: tuple(str, str)