Bio Variant Calling Joint Calling

Version Compatibility

Reference examples tested with: GATK 4.5+, GLnexus 1.4+, bcftools 1.19+

Before using code patterns, verify installed versions match. If versions differ:

• CLI: <tool> --version then <tool> --help to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Joint Calling

"Joint genotype my cohort samples" -> Combine per-sample gVCFs into a single cohort callset with consistent genotyping across all sites, enabling cohort filtering and population-level analysis.

• CLI (GATK): gatk HaplotypeCaller -ERC GVCF -> gatk GenomicsDBImport (or CombineGVCFs) -> gatk GenotypeGVCFs
• CLI (DeepVariant cohorts): deepvariant --output_gvcf per sample -> glnexus_cli --config DeepVariantWGS

The Governing Principle: Genotype Jointly, Never Merge Callsets

Joint genotyping is not the same operation as merging single-sample VCFs, and confusing the two silently corrupts a cohort callset. Two facts drive every decision below:

• Joint genotyping rescues low-coverage het sites. Evidence is borrowed across samples: when one carrier has a confident variant, the cohort allele frequency raises the Bayesian prior at that exact site for every other sample, so a second sample with only 2-3 supporting reads (which per-sample would fall below threshold) is rescued into a confident genotype. This "borrowing information" mechanism is strongest at low coverage and is why a cohort callset is more sensitive than N independent callsets (DePristo 2011 Nat Genet 43:491; Poplin 2018 bioRxiv 201178).
• Joint genotyping produces a squared-off matrix - a genotype at every variant site for every sample. The reference-confidence <NON_REF> allele in each gVCF lets GenotypeGVCFs distinguish confident homozygous reference (0/0) from no-data/no-call (./.) at a site another sample carries.

The decision this forces: never bcftools merge single-sample callsets as if an absent record means hom-ref. A single-sample VCF omits sites where that sample looked reference, so a naive merge fills those cells with ./. (missing), NOT 0/0 - a sample genuinely hom-ref and a sample never assessed become indistinguishable, and downstream allele frequencies and association tests are wrong. ./. != 0/0 is the load-bearing distinction (see variant-calling/vcf-manipulation for merge semantics and variant-calling/vcf-basics for the genotype-field grammar). Genotype from gVCFs so every cell is filled from evidence, not assumption.

Why Joint Calling Matters

Single-sample calling discards cross-sample evidence that is critical for accurate genotyping:

• Statistical power from shared evidence - A site with 2 alt reads in one sample is borderline and would typically be missed. If 50 other samples in the cohort also show 2 alt reads at that site, the evidence is overwhelming and the variant is clearly real. Joint calling aggregates this weak-per-sample signal into strong cohort-level evidence.
• Genotype refinement via cohort priors - Individual genotype likelihoods are combined with cohort allele frequencies as a Bayesian prior. A heterozygous call at a common variant site (AF=0.3) receives more support than the same call at a site with no other carriers. This prior dramatically improves accuracy for low-coverage samples.
• Consistent site representation - All samples are genotyped at the same sites, producing homozygous-reference calls where applicable. Without joint calling, a missing genotype is ambiguous: it could mean homozygous-reference or simply insufficient coverage. This "missing = reference" assumption is a common source of false negatives in downstream analysis.
• Cohort filtering eligibility - Variant quality score recalibration (VQSR) and its successor VETS operate on the whole-cohort variant distribution and need cohort-scale variant counts: a single deep WGS genome supplies enough, but exomes are variant-poor so the ~30-sample floor applies to exomes/panels, not WGS (VQSR's Gaussian mixture needs enough variants to fit), so filtering is inherently a cohort operation, not a per-sample one (see Step 4).

The N+1 Problem and Why gVCFs Solve It

Naive joint calling re-visits every BAM whenever the cohort changes: adding one genome forces re-calling all N. The gVCF workflow decouples expensive per-sample discovery (local assembly + PairHMM likelihoods, captured once per sample in the gVCF) from cheap cohort-wide genotyping. Adding sample N+1 then requires only generating that one gVCF plus re-running the cheap consolidation and GenotypeGVCFs - the assembly work for the existing N is never repeated. The gVCF is the reusable intermediate; GenomicsDB workspaces can even be updated in place (--genomicsdb-update-workspace-path). GATK frames this as decoupling "the initial identification of potential variant sites from the genotyping step, which is the only part that really needs to be done jointly" (see variant-calling/gatk-variant-calling for per-sample gVCF generation).

Cohort Size Decision Table

GenomicsDBImport vs CombineGVCFs

Both produce a combined object that GenotypeGVCFs consumes; they differ in how they store it and how they scale.

Memory landmine specific to GenomicsDBImport: the heavy lifting runs in native C/C++ (TileDB), so cap the JVM heap (--java-options -Xmx) at ~80-90% of RAM. Over-allocating the JVM starves the native layer and causes a native out-of-memory failure that looks unrelated to heap size.

Workflow Overview

Sample BAMs
    │
    ├── HaplotypeCaller (per-sample, -ERC GVCF)
    │   └── sample1.g.vcf.gz, sample2.g.vcf.gz, ...
    │
    ├── CombineGVCFs or GenomicsDBImport
    │   └── Combine into cohort database
    │
    ├── GenotypeGVCFs
    │   └── Joint genotyping
    │
    └── VQSR or Hard Filtering
        └── Final VCF

Step 1: Per-Sample gVCF Generation

# Generate gVCF for each sample
gatk HaplotypeCaller \
    -R reference.fa \
    -I sample1.bam \
    -O sample1.g.vcf.gz \
    -ERC GVCF

# With intervals (faster)
gatk HaplotypeCaller \
    -R reference.fa \
    -I sample1.bam \
    -O sample1.g.vcf.gz \
    -ERC GVCF \
    -L intervals.bed

Batch Processing

# Process all samples
for bam in *.bam; do
    sample=$(basename $bam .bam)
    gatk HaplotypeCaller \
        -R reference.fa \
        -I $bam \
        -O ${sample}.g.vcf.gz \
        -ERC GVCF &
done
wait

Step 2a: CombineGVCFs (Small Cohorts)

For <100 samples:

gatk CombineGVCFs \
    -R reference.fa \
    -V sample1.g.vcf.gz \
    -V sample2.g.vcf.gz \
    -V sample3.g.vcf.gz \
    -O cohort.g.vcf.gz

From Sample Map

# Create sample map file
# sample1    /path/to/sample1.g.vcf.gz
# sample2    /path/to/sample2.g.vcf.gz

ls *.g.vcf.gz | while read f; do
    echo -e "$(basename $f .g.vcf.gz)\t$f"
done > sample_map.txt

# Combine with -V for each
gatk CombineGVCFs \
    -R reference.fa \
    $(cat sample_map.txt | cut -f2 | sed 's/^/-V /') \
    -O cohort.g.vcf.gz

Step 2b: GenomicsDBImport (Large Cohorts)

For >100 samples, use GenomicsDB:

# Create sample map
ls *.g.vcf.gz | while read f; do
    echo -e "$(basename $f .g.vcf.gz)\t$f"
done > sample_map.txt

# Import to GenomicsDB (per chromosome for parallelism)
gatk GenomicsDBImport \
    --sample-name-map sample_map.txt \
    --genomicsdb-workspace-path genomicsdb_chr1 \
    -L chr1 \
    --reader-threads 4

# Or all chromosomes
for chr in {1..22} X Y; do
    gatk GenomicsDBImport \
        --sample-name-map sample_map.txt \
        --genomicsdb-workspace-path genomicsdb_chr${chr} \
        -L chr${chr} &
done
wait

Update GenomicsDB with New Samples

gatk GenomicsDBImport \
    --genomicsdb-update-workspace-path genomicsdb_chr1 \
    --sample-name-map new_samples.txt \
    -L chr1

GenomicsDB Critical Caveats

GenomicsDB is powerful but has sharp edges that can cause data loss or silent failures:

• No sample replacement - Existing samples cannot be updated or overwritten. Only new samples with different names can be added. To fix a sample, the entire workspace must be recreated.
• Intervals locked at import time - The genomic intervals specified during the initial import cannot be changed on incremental updates. Adding new regions requires reimporting from scratch.
• Fragment accumulation - Each incremental batch creates a new database fragment. After thousands of incremental additions, file handle exhaustion becomes likely. Run --consolidate periodically to merge fragments.
• Corruption risk on failed adds - A failed incremental import can leave the datastore in an inconsistent state. Always backup the workspace directory before running --genomicsdb-update-workspace-path.
• Batch size for memory - Set --batch-size 50 to control memory consumption. The default is 0, which loads ALL samples in a single batch (maximum memory); a finite batch size trades a little speed for a bounded heap, so set it explicitly for large cohorts. Larger batches load more gVCFs simultaneously and can exhaust heap space.

Step 3: GenotypeGVCFs

From Combined gVCF

gatk GenotypeGVCFs \
    -R reference.fa \
    -V cohort.g.vcf.gz \
    -O cohort.vcf.gz

From GenomicsDB

gatk GenotypeGVCFs \
    -R reference.fa \
    -V gendb://genomicsdb_chr1 \
    -O chr1.vcf.gz

# All chromosomes
for chr in {1..22} X Y; do
    gatk GenotypeGVCFs \
        -R reference.fa \
        -V gendb://genomicsdb_chr${chr} \
        -O chr${chr}.vcf.gz &
done
wait

# Merge chromosomes
bcftools concat chr{1..22}.vcf.gz chrX.vcf.gz chrY.vcf.gz \
    -Oz -o cohort.vcf.gz

What GenotypeGVCFs Recomputes (and Why It Is Not a Copy)

GenotypeGVCFs re-derives genotypes jointly from the stored per-sample PL vectors (including the <NON_REF> likelihood) under a Bayesian model; it does not simply copy per-sample genotypes into a wider file. Two consequences matter when reading the output:

• GQ and PL are recomputed against the finalized cohort allele set. Once the real ALT alleles are known cohort-wide, the <NON_REF> likelihood is redistributed onto them and PLs are recomputed; GQ is then the difference of the two smallest PLs. A sample's genotype/GQ in the joint VCF can therefore differ from what its single-sample gVCF implied - this is the rescue mechanism working, not a bug.
• The allele-frequency prior comes from --heterozygosity (expected theta, ~0.001 for humans; verify in-tool) and --indel-heterozygosity, folding the cohort allele count into each sample's posterior. --stand-call-conf (~30; verify in-tool) drops sites below that QUAL.

Multiallelics and the Spanning-Deletion * Allele

Joint genotyping across a cohort surfaces two representation issues absent from single-sample calling:

• --max-alternate-alleles caps the ALT alleles genotyped per site (most-supported kept; confirm the default with gatk GenotypeGVCFs --help for the installed version). Genotyping cost scales roughly exponentially in ALT count, so GATK caps it and discourages raising it; --max-genotype-count similarly bounds genotype configurations. Highly multiallelic sites are also where GenotypeGVCFs can blow past very large RAM at scale.
• The * spanning/overlapping-deletion allele (VCF 4.3 reserved) appears at a variant position that falls inside an upstream deletion carried by some samples. It means "for a sample carrying the upstream deletion, these bases are deleted/absent" - not reference, not the local ALT. Such a sample genotypes as */A or */*, which correctly keeps deletion-carriers from being called spuriously hom-ref at the interior site. Downstream tools must special-case it: VEP/SnpEff have no ref/alt sequence to predict a consequence on, and bcftools norm decomposition often splits or filters it, so it is a frequent source of annotation surprises after joint genotyping.

With Allele-Specific Annotations

For larger cohorts where multiallelic sites are common, allele-specific annotations allow VQSR to evaluate each allele independently rather than penalizing a good allele because a co-occurring allele is poor:

gatk GenotypeGVCFs \
    -R reference.fa \
    -V gendb://genomicsdb \
    -O cohort.vcf.gz \
    -G StandardAnnotation \
    -G AS_StandardAnnotation

When allele-specific annotations are present, use -AS mode in VariantRecalibrator and ApplyVQSR for allele-level filtering.

Step 4: Filtering Is a Cohort Operation

Filtering the joint VCF is a whole-cohort step, not a per-sample one, and it must run after joint genotyping. VQSR (and its GATK successor VETS) fit a model to the cohort-wide distribution of site annotations - VQSR's Gaussian mixture needs enough variants and enough overlap with the truth resources to converge, which is why it is unreliable on a single exome or a small panel. This is the decision:

The full VariantRecalibrator/ApplyVQSR/VETS invocations, training resources, tranche levels, and hard-filter thresholds live in variant-calling/filtering-best-practices - this skill does not duplicate that mechanism. A minimal hard-filter fallback for a small cohort:

gatk VariantFiltration \
    -R reference.fa \
    -V cohort.vcf.gz \
    --filter-expression "QD < 2.0" --filter-name "QD2" \
    --filter-expression "FS > 60.0" --filter-name "FS60" \
    --filter-expression "MQ < 40.0" --filter-name "MQ40" \
    -O cohort.filtered.vcf.gz