ExomeFlow

A production-quality Python pipeline for Whole Exome Sequencing analysis.

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 ANNOVAR 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, ANNOVAR, 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
• Getting Help
• License
• Citation

Features

Core workflow

• Zero-config first run: auto-detects or provisions GATK, ANNOVAR, 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 → ANNOVAR 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
• ANNOVAR 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, ANNOVAR buildver, 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, and upgrading to a version with new pipeline steps reprocesses only what's new
• 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
• Modular architecture: pipeline steps compose from a pluggable registry, so new modes gate in and out without touching unrelated code

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  ANNOVAR         Functional annotation          │
│           refGene · ClinVar · gnomAD · dbNSFP · ExAC     │
│           → multianno.vcf  +  multianno.txt              │
└──────────────────────────┬──────────────────────────────┘
                           │
                           ▼
┌─────────────────────────────────────────────────────────┐
│  Step 11  HPO + ACMG      Clinical annotation enrichment │
│           → multianno.hpo.txt                            │
└─────────────────────────────────────────────────────────┘

Benchmarks

Benchmarked on NA12878 (HG001) whole-exome sequencing data (Hg38), default germline single-sample mode. Accuracy evaluated against the GIAB NISTv4.2.1 .

Performance

Variant Quality (PASS variants)

Accuracy (vs. GIAB NISTv4.2.1, PASS-only)

Recall reflects PASS-only evaluation under conservative hard filters. PASS-only extraction is unconditional in ExomeFlow the raw pre-filter VCF is always retained alongside it. This benchmark predates a fix that removed a non-standard DP < 10 site filter — a threshold that was never part of GATK's own hard-filter recommendation and could discard real variants at modest-but-genuine depth (e.g. near capture-kit target edges). The recall above likely understates current behavior; numbers will be refreshed against a re-run once available.

Functional Annotation (NA12878)

Installation

Option 1 Python package (recommended)

pip install exomeflow

Option 2 Docker

docker pull itsrobintomar/exomeflow:2.2.2

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.2 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.2.sif docker://itsrobintomar/exomeflow:2.2.2

# Or build from the definition file
singularity build exomeflow-2.2.2.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.2.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

#!/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.2.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

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, ANNOVAR, 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

# 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 pypi.org/project/exomeflow for details.

Getting Help

For usage questions and bug reports, contact:

Robin Kumar itsrobintomar@gmail.com

Made with ❤️ || Dr Prabudh Goel lab, AIIMS New Delhi , India