bio-variant-calling-deepvariant

What this skill does

## Version Compatibility

Reference examples tested with: DeepVariant 1.6+, GLnexus 1.4+, bcftools 1.19+

Before using code patterns, verify installed versions match. If versions differ:
- CLI: `docker run google/deepvariant:<tag> --version` to confirm build
- `bcftools --version` and `bcftools --help` to confirm flags

If code throws errors, introspect the installed container and adapt the example
to match the actual API rather than retrying.

# DeepVariant Variant Calling

## How DeepVariant Works

DeepVariant reframes variant calling as an image classification problem rather than a statistical genotyping problem. For each candidate variant site, the make_examples step encodes the local read pileup as a 100x221x6-channel image tensor. The six channels encode: read base identity, base quality, mapping quality, strand orientation, read support for the variant allele, and reference base mismatch. A convolutional neural network (CNN) then classifies each pileup image into one of three genotype classes: homozygous reference, heterozygous, or homozygous alternate.

This image-based approach is why DeepVariant outperforms statistical callers in difficult genomic contexts such as homopolymer runs, tandem repeats, and low-complexity regions -- the CNN learns visual patterns in pileup geometry that heuristic filters miss. Separate models are trained for each sequencing platform because read error profiles differ substantially (e.g., Illumina substitution errors vs. ONT indel errors in homopolymers).

DeepVariant calls germline variants only. For somatic variant calling, use DeepSomatic, a separate tool from the same team.

## Installation

```bash
docker pull google/deepvariant:1.6.1

# GPU support (NVIDIA GPU + nvidia-container-toolkit required)
docker pull google/deepvariant:1.6.1-gpu
```

Singularity alternative:

```bash
singularity pull docker://google/deepvariant:1.6.1
```

## Model Selection Guide

| Model | Data Type | Notes |
|-------|-----------|-------|
| `WGS` | Illumina short-read WGS | Default model; trained on 30-50x PCR-free data |
| `WES` | Illumina exome/targeted | Must supply `--regions` BED for efficiency; without it, wastes time scanning non-target regions |
| `PACBIO` | PacBio HiFi only | Not trained on CLR reads; CLR error profile is fundamentally different |
| `ONT_R104` | ONT R10.4+ chemistry | Accuracy lower than HiFi model; R9.4 reads perform poorly with this model |
| `HYBRID_PACBIO_ILLUMINA` | Mixed platforms | Emerging use case for samples with both HiFi and Illumina data |

Model selection has a large effect on accuracy. Using the wrong model (e.g., WGS model on HiFi data) silently degrades results because the CNN expects platform-specific error patterns in the pileup images.

## When to Use DeepVariant

**Best choice when:**
- Highest germline accuracy is the primary goal
- GPU resources are available (or CPU time is acceptable)
- Single-sample calling or small cohort with GLnexus for joint genotyping
- Data is from a supported platform (Illumina, HiFi, ONT R10.4+)

**Consider GATK HaplotypeCaller instead when:**
- Joint calling across large cohorts (GATK GenomicsDB + GenotypeGVCFs scales better)
- VQSR-based filtering is needed (DeepVariant QUAL scores are CNN confidence, not amenable to VQSR)
- Clinical pipeline requires established GATK validation and regulatory precedent

**Consider Clair3 instead when:**
- Long-read data where speed matters more than marginal accuracy gains
- ONT data specifically (Clair3 has strong ONT-specific models)
- Resource-constrained environments without GPU access

## Basic Usage

```bash
docker run -v "${PWD}:/input" -v "${PWD}/output:/output" \
    google/deepvariant:1.6.1 \
    /opt/deepvariant/bin/run_deepvariant \
    --model_type=WGS \
    --ref=/input/reference.fa \
    --reads=/input/sample.bam \
    --output_vcf=/output/sample.vcf.gz \
    --output_gvcf=/output/sample.g.vcf.gz \
    --num_shards=16
```

Always generate gVCFs (`--output_gvcf`) even for single samples -- they enable downstream joint calling without re-running DeepVariant.

## Step-by-Step Workflow

For more control over intermediate outputs, run each stage separately:

### Step 1: Make Examples

```bash
docker run -v "${PWD}:/data" google/deepvariant:1.6.1 \
    /opt/deepvariant/bin/make_examples \
    --mode calling \
    --ref /data/reference.fa \
    --reads /data/sample.bam \
    --examples /data/examples.tfrecord.gz \
    --gvcf /data/gvcf.tfrecord.gz
```

### Step 2: Call Variants

```bash
docker run -v "${PWD}:/data" google/deepvariant:1.6.1 \
    /opt/deepvariant/bin/call_variants \
    --outfile /data/call_variants.tfrecord.gz \
    --examples /data/examples.tfrecord.gz \
    --checkpoint /opt/models/wgs/model.ckpt
```

### Step 3: Postprocess Variants

```bash
docker run -v "${PWD}:/data" google/deepvariant:1.6.1 \
    /opt/deepvariant/bin/postprocess_variants \
    --ref /data/reference.fa \
    --infile /data/call_variants.tfrecord.gz \
    --outfile /data/output.vcf.gz \
    --gvcf_outfile /data/output.g.vcf.gz \
    --nonvariant_site_tfrecord_path /data/gvcf.tfrecord.gz
```

## GPU Acceleration

GPU acceleration primarily benefits the call_variants step (CNN inference). The make_examples and postprocess_variants steps are CPU-bound and benefit more from `--num_shards` parallelism.

```bash
docker run --gpus all -v "${PWD}:/data" \
    google/deepvariant:1.6.1-gpu \
    /opt/deepvariant/bin/run_deepvariant \
    --model_type=WGS \
    --ref=/data/reference.fa \
    --reads=/data/sample.bam \
    --output_vcf=/data/output.vcf.gz \
    --num_shards=16
```

## PacBio HiFi Calling

```bash
docker run -v "${PWD}:/data" google/deepvariant:1.6.1 \
    /opt/deepvariant/bin/run_deepvariant \
    --model_type=PACBIO \
    --ref=/data/reference.fa \
    --reads=/data/hifi_aligned.bam \
    --output_vcf=/data/hifi_variants.vcf.gz \
    --num_shards=16
```

HiFi reads achieve Q30+ per-read accuracy, giving DeepVariant cleaner pileup images.
The PACBIO model is not suitable for CLR reads (Q10-Q15 accuracy).

## ONT Calling

```bash
docker run -v "${PWD}:/data" google/deepvariant:1.6.1 \
    /opt/deepvariant/bin/run_deepvariant \
    --model_type=ONT_R104 \
    --ref=/data/reference.fa \
    --reads=/data/ont_aligned.bam \
    --output_vcf=/data/ont_variants.vcf.gz \
    --num_shards=16
```

R10.4+ chemistry substantially reduces systematic indel errors in homopolymers compared to R9.4. For R9.4 data, consider Clair3 which has dedicated R9.4 models.

## Exome/Targeted Sequencing

```bash
docker run -v "${PWD}:/data" google/deepvariant:1.6.1 \
    /opt/deepvariant/bin/run_deepvariant \
    --model_type=WES \
    --ref=/data/reference.fa \
    --reads=/data/exome.bam \
    --regions=/data/targets.bed \
    --output_vcf=/data/exome_variants.vcf.gz \
    --num_shards=8
```

The `--regions` flag is not strictly required but omitting it causes DeepVariant to scan the entire genome, wasting hours on off-target reads with minimal coverage.

## Joint Calling with GLnexus

For multi-sample cohorts, use per-sample gVCFs merged with GLnexus:

```bash
for bam in *.bam; do
    sample=$(basename $bam .bam)
    docker run -v "${PWD}:/data" google/deepvariant:1.6.1 \
        /opt/deepvariant/bin/run_deepvariant \
        --model_type=WGS \
        --ref=/data/reference.fa \
        --reads=/data/$bam \
        --output_vcf=/data/${sample}.vcf.gz \
        --output_gvcf=/data/${sample}.g.vcf.gz \
        --num_shards=16
done

docker run -v "${PWD}:/data" quay.io/mlin/glnexus:v1.4.1 \
    /usr/local/bin/glnexus_cli \
    --config DeepVariantWGS \
    /data/*.g.vcf.gz \
    | bcftools view - -Oz -o cohort.vcf.gz
```

### GLnexus Configuration

| Config | Use Case | Notes |
|--------|----------|-------|
| `DeepVariantWGS` | Illumina WGS gVCFs | Default for most WGS cohorts |
| `DeepVariantWES` | Illumina exome gVCFs | Tuned for higher-depth, narrower-region calling |
| `DeepVariant_unfiltered` | Keep all variant sites | Useful for research e