ExomeFlow

Reproducible Python pipeline for whole-exome sequencing, germline variant calling, annotation, and interpretation

CI
Package
Container
Meta
Author

---

Overview

ExomeFlow takes raw paired-end FASTQ files to functionally annotated, clinically contextualized variants in a single reproducible command. It implements the GATK best-practices germline short-variant workflow — quality control, alignment, base quality recalibration, variant calling, hard filtering, and functional annotation — as a modular Python package, with cohort joint genotyping, tumor-only somatic calling, read-depth CNV detection, and HPO/ACMG clinical annotation available as opt-in extensions on top of the same pipeline.

The tool is designed to run unattended: on first invocation it detects or provisions every external dependency it needs (GATK, the annotation engine, reference genomes, annotation databases) and persists the resolved configuration, so a single exomeflow run command is sufficient for both a first-time setup and every run after it.

Table of Contents

• Overview
• Features
• Pipeline Workflow
• Benchmarks
• Installation
• System Requirements
• Python Dependencies
• Quick Start
• Commands
• Reference Files
• Input Convention
• Output Files
• Known Limitations
• License
• Getting Help

Features

Core workflow

• Zero-config first run: auto-detects or provisions GATK, the annotation engine, reference genomes, and annotation databases on first use; the resolved configuration is saved to ~/.exomeflow/config.json so every later run needs only --input-dir/--output
• GATK best-practice germline calling: fastp QC → BWA-MEM alignment → coordinate sorting → duplicate marking → BQSR → HaplotypeCaller → hard filtering → functional annotation
• Automatic sample detection: scans an input directory for paired-end FASTQ files by naming convention; no manifest required
• Per-sample output by default: any number of samples processed together still produces one separate annotated file per sample, identical to running each alone
• Functional annotation: against refGene, ClinVar, gnomAD, dbNSFP, ExAC, and avSNP150, with automatic per-build (hg38/GRCh37) database selection
• HPO and ACMG clinical enrichment: every annotated table is automatically joined with HPO gene-to-phenotype terms and ACMG/AMP pathogenicity classification (via InterVar); enrichment degrades gracefully rather than blocking a run when its dependencies aren't yet available

Cohort and advanced modes (opt-in)

• Cohort joint genotyping: (--joint-genotyping) GenomicsDBImport + GenotypeGVCFs across all samples, producing one shared cohort VCF and annotation set
• Somatic tumor-only calling: (--mode somatic) Mutect2 with contamination-aware filtering (tumor-normal pairing is not yet supported)
• Read-depth CNV calling: (--cnv) GATK CollectReadCounts/DenoiseReadCounts/ PlotDenoisedCopyRatios, no panel-of-normals required
• GRCh37 or hg38: (--genome-build) reference bundle, annotation build, and database selection all follow the requested build automatically
• Cohort-wide MultiQC rollup: aggregating fastp, alignment, and GATK metrics, generated automatically at the end of every run

Operational reliability

• Checkpointed and resumable: every completed step is recorded per sample; an interrupted run resumes exactly where it left off
• Parallel cohort processing: via --max-workers
• Structured logging: per-sample and pipeline-wide logs with INFO / WARNING / ERROR / SUCCESS levels
• Automatic dependency management: every external tool and database is auto-detected and, if missing, installed or downloaded without a separate setup step

Pipeline Workflow

Raw FASTQ
    │
    ▼
┌─────────────────────────────────────────────────────────┐
│  Step 1   fastp         Quality control & adapter trim   │
│           length ≥ 50 bp · base quality ≥ Q30            │
└──────────────────────────┬──────────────────────────────┘
                           │
                           ▼
┌─────────────────────────────────────────────────────────┐
│  Step 2   BWA MEM        Read alignment to hg38          │
│           -Y -K 100000000 · read-group tags set          │
└──────────────────────────┬──────────────────────────────┘
                           │
                           ▼
┌─────────────────────────────────────────────────────────┐
│  Step 3   GATK SortSam   Coordinate-sort BAM             │
│  Step 4   samtools       Flagstat alignment QC           │
│  Step 5   GATK MarkDuplicates   PCR duplicate removal    │
│  Step 6   GATK BuildBamIndex    BAI index                │
└──────────────────────────┬──────────────────────────────┘
                           │
                           ▼
┌─────────────────────────────────────────────────────────┐
│  Step 7   GATK BQSR      BaseRecalibrator + ApplyBQSR    │
│           Known sites: dbSNP · Mills · known indels      │
│           → recalibrated.bam  (IGV-ready)                │
└──────────────────────────┬──────────────────────────────┘
                           │
                           ▼
┌─────────────────────────────────────────────────────────┐
│  Step 8   GATK HaplotypeCaller   Variant calling         │
│           Exome intervals + padding · dbSNP annotation   │
└──────────────────────────┬──────────────────────────────┘
                           │
                    ┌──────┴──────┐
                    ▼             ▼
               SNP filters   INDEL filters
                    └──────┬──────┘
                           │  MergeVcfs
                           ▼
┌─────────────────────────────────────────────────────────┐
│  Step 9   SelectVariants  Extract PASS-only variants     │
└──────────────────────────┬──────────────────────────────┘
                           │
                           ▼
┌─────────────────────────────────────────────────────────┐
│  Step 10  Annotation      Functional annotation          │
│           refGene · ClinVar · gnomAD · dbNSFP · ExAC     │
│           → multianno.vcf  +  multianno.txt              │
└──────────────────────────┬──────────────────────────────┘
                           │
                           ▼
┌─────────────────────────────────────────────────────────┐
│  Step 11  HPO + ACMG      Clinical annotation enrichment │
│           → multianno.hpo.txt                            │
└─────────────────────────────────────────────────────────┘

--joint-genotyping, --mode somatic, and --cnv each modify or extend this chain for the samples they apply to — see Commands.

Benchmarks

Benchmarked on NA12878 (GIAB HG001) whole-exome data (GRCh38), default germline single-sample mode, called with --intervals set to the Agilent SureSelect V8 targets. Accuracy is measured with RTG vcfeval against GIAB NISTv4.2.1, restricted to GIAB high-confidence regions ∩ capture targets.

Accuracy (vs. GIAB HG001 NISTv4.2.1)

Coverage in this public dataset is uneven: 12% of capture targets received essentially no reads, and 17% of the evaluation region is below 10x. Variants there cannot be called by any pipeline, so results are reported both over the full evaluation region and over the callable part of it (≥20x read depth), following GIAB stratification practice.

Where reads exist, ExomeFlow recovers 99.1% of true SNPs. PASS-only recall is lower because of GATK's standard hard filters; the unfiltered VCF is always kept alongside the PASS VCF. INDEL rates are based on few truth sites (130).

Variant Quality (PASS variants)

¹ Computed on the evaluation region (GIAB high-confidence ∩ capture targets), so calls in the ±100 bp padding flanks and hard-to-map regions don't skew the ratios.

Performance

BQSR is now restricted to the --intervals targets, so current runtimes are shorter than shown.

Functional Annotation (NA12878)

Installation

Option 1 — Python package (recommended)

pip install exomeflow

Option 2 — Docker

docker pull itsrobintomar/exomeflow:2.2.16

docker run --rm -it \
  -v /path/to/fastq:/data/fastq \
  -v /path/to/refs:/refs \
  -v /path/to/vcf:/vcf \
  -v /path/to/annovar:/annovar \
  -v /path/to/results:/data/results \
  itsrobintomar/exomeflow:2.2.16 run \
    --input-dir    /data/fastq \
    --output       /data/results \
    --reference    /refs/Homo_sapiens_assembly38.fasta \
    --dbsnp        /vcf/Homo_sapiens_assembly38.dbsnp138.vcf.gz \
    --mills        /vcf/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
    --known-indels /vcf/Homo_sapiens_assembly38.known_indels.vcf.gz \
    --annovar-bin  /annovar \
    --annovar-db   /annovar/humandb \
    --threads      24

ANNOVAR must be mounted — it cannot be bundled due to its licensing terms. Register and download at annovar.openbioinformatics.org.

Option 3 — Singularity (HPC clusters)

# Pull directly from Docker Hub
singularity pull exomeflow-2.2.16.sif docker://itsrobintomar/exomeflow:2.2.16

# Or build from the definition file
singularity build exomeflow-2.2.16.sif exomeflow.def

singularity exec \
  --bind /path/to/fastq:/data/fastq \
  --bind /path/to/refs:/refs \
  --bind /path/to/vcf:/vcf \
  --bind /path/to/annovar:/annovar \
  --bind /path/to/results:/data/results \
  exomeflow-2.2.16.sif exomeflow run \
    --input-dir    /data/fastq \
    --output       /data/results \
    --reference    /refs/Homo_sapiens_assembly38.fasta \
    --dbsnp        /vcf/Homo_sapiens_assembly38.dbsnp138.vcf.gz \
    --mills        /vcf/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
    --known-indels /vcf/Homo_sapiens_assembly38.known_indels.vcf.gz \
    --annovar-bin  /annovar \
    --annovar-db   /annovar/humandb \
    --threads      24

SLURM job script example

#!/bin/bash
#SBATCH --job-name=exomeflow
#SBATCH --cpus-per-task=24
#SBATCH --mem=90G
#SBATCH --time=24:00:00
#SBATCH --output=exomeflow_%j.log

singularity exec \
  --bind $FASTQ_DIR:/data/fastq \
  --bind $REFS_DIR:/refs \
  --bind $VCF_DIR:/vcf \
  --bind $ANNOVAR_DIR:/annovar \
  --bind $RESULTS_DIR:/data/results \
  exomeflow-2.2.16.sif exomeflow run \
    --input-dir    /data/fastq \
    --output       /data/results \
    --reference    /refs/Homo_sapiens_assembly38.fasta \
    --dbsnp        /vcf/Homo_sapiens_assembly38.dbsnp138.vcf.gz \
    --mills        /vcf/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
    --known-indels /vcf/Homo_sapiens_assembly38.known_indels.vcf.gz \
    --annovar-bin  /annovar \
    --annovar-db   /annovar/humandb \
    --threads      $SLURM_CPUS_PER_TASK

System Requirements

ExomeFlow invokes the following external tools via the command line. Each is auto-detected on first run and installed automatically where possible (see Quick Start); the table below documents the minimum versions, not a manual installation requirement.

ANNOVAR is the one dependency ExomeFlow cannot provision for you — its license requires individual registration. Everything else, including GATK itself, is detected or fetched automatically.

Python Dependencies

• typer — CLI interface
• rich — Structured, coloured terminal output
• pandas — Variant summaries and HPO/ACMG annotation joins

Installed automatically with pip install exomeflow. matplotlib (needed only for --cnv plots) is an optional extra (pip install exomeflow[viz]) that the dependency checker installs automatically the first time --cnv is used.

Quick Start

Before the first run: ExomeFlow downloads its reference resources (about 110 GB: GATK bundle ~20 GB plus annotation databases ~90 GB). Make sure the target disk has that space free, plus room for your BAMs and VCFs. See Reference Files.

pip install exomeflow

Arrange paired-end FASTQ files following the naming convention:

fastq/
├── sample1_1.fastq.gz
├── sample1_2.fastq.gz
├── sample2_1.fastq.gz
└── sample2_2.fastq.gz

Run the pipeline:

exomeflow run --input-dir fastq/ --output results/

On first run, ExomeFlow detects or provisions GATK, the annotation engine, reference data, the HPO gene-to-phenotype mapping, and InterVar, then saves the resolved configuration to ~/.exomeflow/config.json. Every later run needs nothing beyond --input-dir and --output. Add --yes to skip all interactive prompts for unattended or CI use.

To control every path explicitly instead — for example on a shared cluster where reference data already exists at known locations:

exomeflow run \
  --input-dir    fastq/ \
  --output       results/ \
  --reference    /data/references/hg38/hg38.fa \
  --dbsnp        /data/references/hg38/dbsnp.vcf.gz \
  --mills        /data/references/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
  --known-indels /data/references/hg38/Homo_sapiens_assembly38.known_indels.vcf.gz \
  --intervals    refs/Illumina_Exome_TargetedRegions_v1.2.hg38.bed \
  --annovar-bin  /opt/annovar \
  --annovar-db   /opt/annovar/humandb \
  --threads      32 \
  --max-workers  2

exomeflow setup is also available to run provisioning as its own step, independent of a pipeline run — useful for switching reference builds or refreshing databases.

Cohort, somatic, CNV, and GRCh37 modes

# Cohort joint genotyping instead of per-sample VCFs
exomeflow run --input-dir fastq/ --output results/ --joint-genotyping --intervals targets.bed

# Somatic tumor-only calling with Mutect2
exomeflow run --input-dir fastq/ --output results/ --mode somatic \
  --germline-resource af-only-gnomad.vcf.gz --panel-of-normals pon.vcf.gz

# Read-depth CNV calling alongside the standard germline workflow
exomeflow run --input-dir fastq/ --output results/ --cnv --intervals targets.bed

# GRCh37/hg19 instead of hg38
exomeflow run --input-dir fastq/ --output results/ --genome-build GRCh37

Commands

exomeflow run — Execute the WES pipeline

exomeflow run [OPTIONS]

--threads and --java-opts are sized automatically from available CPU/RAM when not set explicitly.

exomeflow setup — Run provisioning independently of a pipeline run

exomeflow setup [--refs-dir PATH] [--annovar-bin PATH] [--annovar-db PATH] [--genome-build hg38|GRCh37] [--existing-refs PATH] [--yes]

Not required before exomeflow run: first-run auto-setup covers the same ground. Useful for re-provisioning without running the pipeline itself.

Reference Files

exomeflow run/exomeflow setup download all GATK resource bundle files automatically for the requested --genome-build. Manual download, if needed:

gsutil -m cp -r gs://gcp-public-data--broad-references/hg38/v0/ /data/refs/

Input Convention

ExomeFlow detects samples automatically from paired-end FASTQ filenames:

<sample_id>_1.fastq.gz   ← Read 1
<sample_id>_2.fastq.gz   ← Read 2

sample_id can be any string sharing a common prefix between the two read files: an SRR accession, a patient ID, or anything else.

Output Files

Per-sample output (default — one complete set per sample, regardless of how many samples are in the run):

Cohort output (--joint-genotyping only — replaces the per-sample VCF/annotation files above with one shared set):

Generated at the end of every run:

Known Limitations

• Somatic mode is tumor-only. Matched tumor-normal pairing is not yet supported.
• HPO/ACMG enrichment depends on external, unversioned resources (the HPO gene-to-phenotype mapping and InterVar's own reference databases). When these aren't available, the corresponding columns are omitted rather than the run failing — check the pipeline log if HPO_terms/ACMG_classification are missing from a given sample's output.
• GRCh37/hg19 support annotates with gnomAD v2.1.1, not v4.1; gnomAD v4 was never released for that build.

License

MIT — see the LICENSE file.

Getting Help

For usage questions and bug reports, contact:

PG Lab, AIIMS New Delhi — drprabudhgoel.in/contactus

Made with ❤️ by Dr Prabudh Goel Lab, AIIMS New Delhi, India