GATK GenotypeGVCFs memory issue when including non-variant-sites and using GenomicsDB as input

[percyfal]

This is more of an issue with GATK itself but since there is a known workaround I post my attempts to resolve it, as requested by @maxulysse .

I am running Sarek on non-model organism with a smallish genome. I include non variant sites (see #1870) which requires the inclusion of the intervals file when creating the GenomicsDB (#1872). I have multiple populations and the workflow works fine for the smaller samples, but for the largest (n~100), the workflow fails due to a memory leak, despite the modest sample size. This turned out to be a known GATK issue (broadinstitute/gatk#8989). The proposed fix is

Our temporary solution until we make an updated release would be to convert imported genomicsDB instances to GVCF using

gatk SelectVariants -V gendb://instancename -O GVCF_export.g.vcf.gz -R ref.fa -L whateverintervalusedinGDBimport

Therefore, as a quick fix I have modified modules/nf-core/gatk4/genotypegvcfs/main.nf to perform two steps, namely 1. extract regions from the DB to a GVCF and 2. use GVCF as input to GenotypeGVCFs:

gatk --java-options "-Xmx${avail_mem}M -XX:-UsePerfData" \\
        SelectVariants \\
        --variant $input_command \\
        --output ${prefix}.g.vcf.gz \\
        --reference $fasta \\
        $interval_command \\
        --tmp-dir .

    gatk --java-options "-Xmx${avail_mem}M -XX:-UsePerfData" \\
        GenotypeGVCFs \\
        --variant ${prefix}.g.vcf.gz \\
        --output ${prefix}.vcf.gz \\
        --reference $fasta \\
        $interval_command \\
        $dbsnp_command \\
        --tmp-dir . \\
        $args

    rm -f ${prefix}.g.vcf.gz

I am running the workflow on a SLURM cluster. Although the fix does create the desired output file, Nextflow exits with a cryptic error along the lines of not being able to write to a trace file. I have set trace.overwrite = true (see nextflow-io/nextflow#3317), to no avail.

Since the genotyping step takes a day or two to finish on my data, I have tried to reproduce the error using one of the test data sets:

nextflow run main.nf -profile test,uppmax -params-file params.yml --project $NAISS_COMPUTE_ID --input ./tests/csv/3.0/fastq_pair.csv --output results

Omitting the uppmax profile works fine, confirming that the fix is ok, but the above command fails with the error

ERROR ~ Cannot invoke method getName() on null object

 -- Check script './workflows/sarek/../../subworkflows/local/channel_variant_calling_create_csv/main.nf' at line: 16 or see '.nextflow.log' file for more details

The full java stack trace does not really help. The line that seems to be the culprit in subworkflows/local/channel_variant_calling_create_csv/main.nf is

    main:
        // Creating csv files to restart from this step
        vcf_to_annotate.collectFile(keepHeader: true, skip: 1,sort: true, storeDir: "${outdir}/csv"){ meta, vcf ->
            patient       = meta.patient
            sample        = meta.id
            variantcaller = meta.variantcaller
            vcf = "${outdir}/variant_calling/${variantcaller}/${meta.id}/${vcf.getName()}"
            ["variantcalled.csv", "patient,sample,variantcaller,vcf\n${patient},${sample},${variantcaller},${vcf}\n"]
        }

where the vcf variable is defined.

The proper solution would be to run SelectVariants separately after the creation of GenomicsDB. I was hoping I could cut a corner, but apparently not, it seems.

[percyfal]

For reference, in my local version, I have now separated SelectVariants and GenotypeGVCFs steps. In short, I added a local module modules/local/gatk4/selectvariants (with some modifications compared to the nf-core version, notably the return of the intervals used):

Local GATK4/SelectVariants

process SELECTVARIANTS {
    tag "$meta.id"
    label 'process_single'

    conda "${moduleDir}/environment.yml"
    container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ?
        'https://community-cr-prod.seqera.io/docker/registry/v2/blobs/sha256/b2/b28daf5d9bb2f0d129dcad1b7410e0dd8a9b087aaf3ec7ced929b1f57624ad98/data':
        'community.wave.seqera.io/library/gatk4_gcnvkernel:e48d414933d188cd' }"

    input:
    tuple val(meta), path(input), path(vcf_idx), path(intervals), path(intervals_index)
    tuple val(meta2), path(fasta)
    tuple val(meta3), path(fai)
    tuple val(meta4), path(dict) // required if input is a GenomicsDB?

    output:
    tuple val(meta), path("*.g.vcf.gz")       , emit: gvcf
    tuple val(meta), path("*.g.vcf.gz.tbi")   , emit: gtbi
    tuple val(meta), path("$interval_list"), optional:true, emit: intervallist
    path "versions.yml"		                , emit: versions

    when:
    task.ext.when == null || task.ext.when

    script:
    def args = task.ext.args ?: ''
    def prefix = task.ext.prefix ?: "${meta.id}"
    def input_command = input.name.endsWith(".vcf") || input.name.endsWith(".vcf.gz") ? "$input" : "gendb://$input"
    def interval = intervals ? "--intervals ${intervals}" : ""

	interval_list = intervals ? intervals : "${prefix}.interval_list"
    def avail_mem = 3072
    if (!task.memory) {
        log.info '[GATK SelectVariants] Available memory not known - defaulting to 3GB. Specify process memory requirements to change this.'
    } else {
        avail_mem = (task.memory.mega*0.8).intValue()
    }
    """
    gatk --java-options "-Xmx${avail_mem}M -XX:-UsePerfData" \\
        SelectVariants \\
        --variant $input_command \\
        --output ${prefix}.g.vcf.gz \\
        --reference $fasta \\
        $interval \\
        --tmp-dir . \\
        $args

    cat <<-END_VERSIONS > versions.yml
    "${task.process}":
        gatk4: \$(echo \$(gatk --version 2>&1) | sed 's/^.*(GATK) v//; s/ .*\$//')
    END_VERSIONS
    """

    stub:
    def prefix = task.ext.prefix ?: "${meta.id}"
    """
    touch ${prefix}.g.vcf.gz
    touch ${prefix}.g.vcf.gz.tbi

    cat <<-END_VERSIONS > versions.yml
    "${task.process}":
        gatk4: \$(echo \$(gatk --version 2>&1) | sed 's/^.*(GATK) v//; s/ .*\$//')
    END_VERSIONS
    """
}

I modified subworkflows/local/bam_joint_calling_germline_gatk/main.nf to include the following (after GATK4_GENOMICSDB):

    // Map input for GenotypeGVCFs conditional on option --allsites;
    // currently there is a bug in GATK4 which makes it impossible to
    // use GenomicsDB as input to GenotypeGVCFs
    select_variants_in = GATK4_GENOMICSDBIMPORT.out.genomicsdb.join(
        GATK4_GENOMICSDBIMPORT.out.intervallist).map{
            meta, genomicsdb, intervallist -> [ meta, genomicsdb, [], intervallist, [] ]
    }

    // Select only the intervals that were used for GenomicsDBImport
    // and convert to gvcf for joint genotyping if --allsites is set
    SELECTVARIANTS(select_variants_in, fasta, fai, dict)

    // This step should only be needed if --allsites is set
    ch_vcf = SELECTVARIANTS.out.gvcf.collect()
    ch_tbi = SELECTVARIANTS.out.gtbi.collect()
    ch_intervallist = SELECTVARIANTS.out.intervallist.collect()
    ch_genotypegvcfs_in = ch_vcf
        .join(ch_tbi, failOnDuplicate: true, failOnMismatch: true)
        .join(ch_intervallist, failOnDuplicate: true, failOnMismatch: true)
        .map{
            meta, gvcf, gtbi, intervallist -> [ meta, gvcf, gtbi, intervallist, [] ]
        }

    // Joint genotyping performed using GenotypeGVCFs
    // Sort vcfs called by interval within each VCF
    GATK4_GENOTYPEGVCFS(ch_genotypegvcfs_in, fasta, fai, dict, dbsnp.map{ it -> [ [:], it ] }, dbsnp_tbi.map{ it -> [ [:], it ] })

Ideally, the separate SelectVariants step should only be called if a flag (e.g., --all-variants) is set.