Bio Workflows Fastq To Variants

Version Compatibility

Reference examples tested with: BWA-MEM2 2.2.1+, GATK 4.5+, bcftools 1.19+, samtools 1.19+, fastp 0.23+, DeepVariant 1.6+, hap.py 0.3.15+, Ensembl VEP 111+

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.

Note: GenotypeGVCFs defaults (--max-alternate-alleles, --heterozygosity, --stand-call-conf), the GATK hard-filter thresholds, and DRAGEN speed/accuracy figures drift by version/vendor; confirm in-tool and against current GIAB benchmarks before quoting.

FASTQ to Variants Workflow

"Call variants from my whole-genome or exome FASTQ files" -> Chain QC/trim, alignment, duplicate marking, an engine-appropriate caller, normalization, filtering, annotation, and benchmarking into one filtered germline VCF.

• CLI: fastp -> bwa-mem2 -> samtools markdup -> (bcftools | gatk HaplotypeCaller | DeepVariant) -> bcftools norm -> filter -> VEP -> hap.py

This is a workflow skill: it owns the chaining decisions and hand-offs, not the internals of any one step. Every step below cross-references the component skill that teaches its mechanism.

The governing principle

A germline pipeline is a chain of commitments, and the two that decide whether the callset is trustworthy are made at the seams between steps, not inside them.

1. The reference is a pipeline-wide commitment made once and inherited by everything downstream. The build and analysis set chosen at alignment (step 2) fix the coordinates of every later comparison: the dbSNP/ClinVar/gnomAD records annotation matches against, the truth BED benchmarking scores within, and the cohort other samples are joint-called with must all be the SAME build. Changing it later means redoing alignment, calling, normalization, and annotation. Decide before aligning a single read.
2. The step order is not arbitrary; two orderings sink reviews. (a) Normalize BEFORE annotate. An indel that is not left-aligned to the database's canonical position is one base off, so the annotator silently misses the ClinVar/gnomAD record and reports a pathogenic variant as novel-absent -- a patient-safety failure that throws no error. (b) Never bcftools merge single-sample VCFs into a cohort. Absence of a record is then read as homozygous reference, fabricating genotypes; joint-genotype per-sample gVCFs instead so "confident hom-ref" is distinguished from "no data."
3. Filter site-level first, then genotype-level, then recompute cohort QC. Setting genotypes to no-call (./.) changes missingness, HWE, and allele frequencies, so those metrics must be computed on the genotype-filtered matrix, not before.
4. A single genome-wide F1 is nearly meaningless. Callers agree on easy SNPs (F1 > 0.999); they diverge in indels-in-repeats, segmental duplications, and MHC. Benchmark stratified, within the GIAB confident region, or do not claim accuracy (Krusche 2019 Nat Biotechnol 37:555-560).

Workflow overview

FASTQ
  | [1] QC & trim -------------------> fastp            (read-qc/fastp-workflow)
  v
  | [2] Align ----------------------> bwa-mem2          (read-alignment/bwa-alignment)
  v     ^-- reference commitment: GRCh38 analysis set / T2T, ALT/decoy handling
  | [3] Mark duplicates ------------> samtools markdup  (alignment-files/duplicate-handling)
  v
  | [4] (BQSR? optional on modern binned-quality instruments)
  v
  | [5] Call -----------------------> bcftools | GATK HaplotypeCaller | DeepVariant | DRAGEN
  v     ^-- single-sample OR per-sample gVCF -> joint-genotype (variant-calling/joint-calling)
  | [6] Normalize (BEFORE annotate) -> bcftools norm -m-any -f ref  (variant-calling/variant-normalization)
  v
  | [7] Filter: site-level THEN genotype-level         (variant-calling/filtering-best-practices)
  v
  | [8] Recompute cohort QC on the genotype-filtered matrix  (missingness, Ti/Tv, het/hom, excess-het HWE)
  v
  | [9] Annotate -------------------> VEP / SnpEff      (variant-calling/variant-annotation)
  v
  | [10] Benchmark & QC ------------> hap.py/vcfeval, bcftools stats  (variant-calling/vcf-statistics)
  v
Filtered, normalized, benchmarked VCF

Reference genome: the pipeline-wide commitment

Decision made once, before alignment; everything downstream inherits it. Mechanism of building/indexing the reference lives in read-alignment/bwa-alignment; the deeper reasoning below is what a reviewer expects justified.

The reason to fix this first: annotation databases, panel BEDs, benchmark truth sets, and any cohort a sample is joint-called with are all coordinate-specific. Mixing builds (e.g. normalizing to GRCh38 then annotating against a GRCh37 dbSNP) is a guaranteed silent miss. "We used GRCh38" is under-specified -- plain vs masked vs analysis-set-with-decoys materially changes results in named clinical genes.

The canonical order and why

Each step assumes the previous; the order is defensible under review (canonical preprocessing order, filtering/representation practice).

1. QC/trim -- remove adapters and low-quality tails before they corrupt alignment and duplicate detection.
2. Align -- to the committed reference, with read groups (SM/ID/PL/LB); read groups are a hard GATK requirement.
3. Mark duplicates -- PCR/optical duplicates are not independent evidence; marking (not removing) lets the caller down-weight them. Skip for amplicon/UMI data.
4. BQSR -- honestly optional on modern instruments. BQSR corrected context/cycle-dependent miscalibration on 2010-era continuous-quality Illumina. NovaSeq/NovaSeq X emit ~4 quality bins, leaving little to recalibrate; callsets are largely unchanged with vs without it. DeepVariant explicitly recommends NOT running BQSR (its CNN learned the raw-quality error model); DRAGEN handles quality internally (and uses DRAGSTR for STR/indel error modeling). Keep it for GATK-HaplotypeCaller consistency if a frozen pipeline demands it; otherwise the modern indel-accuracy lever is STR-aware error modeling (--dragen-mode), not BQSR.
5. Call -- per-sample VCF, or per-sample gVCF (-ERC GVCF) if a cohort will be joint-genotyped.
6. Normalize BEFORE annotate/compare -- bcftools norm -m-any -f ref.fa (split multiallelics, then left-align + parsimony), against the SAME reference used for annotation. Add -a (atomize) only when the downstream database is decomposed. This is the most common real ordering bug: annotate-then-normalize attaches consequences to a non-canonical representation that fails to match the database. The binding constraint is normalize-before-ANNOTATE/COMPARE, not normalize-before-filter: GATK convention runs VariantFiltration on the raw multiallelic records (its annotations are computed on that representation), then normalizes -- bwa_gatk_workflow.sh does exactly that, while bwa_bcftools_workflow.sh normalizes first. Both are correct; neither annotates before normalizing.
7. Filter site-level, then genotype-level -- site filters (VQSR/hard/ML) decide whether a site is real; genotype filters (GQ/DP/allele-balance) decide whether an individual genotype is trustworthy. SNPs and indels are filtered separately (different error processes and truth resources).
8. Recompute cohort QC on the genotype-filtered matrix -- missingness, Ti/Tv, het/hom, excess-het HWE. Doing HWE before genotype filtering lets low-GQ garbage drive spurious deviation.
9. Annotate on the normalized (and, where consequence matters, haplotype-resolved) representation.
10. Benchmark/validate after transforms (hap.py, bcftools stats).

Choosing the calling engine

Pipeline-level selection only; the mechanism and full decision table live in variant-calling/variant-calling. Hand off there to pick, then return here for chaining.

Single-sample vs cohort is a chaining decision, not a caller feature. For a cohort, emit per-sample gVCFs and joint-genotype them so a variant seen in one sample is evaluated in all (cohort rescue of low-coverage hets, squared-off genotype matrix). This is what makes the pipeline forward-compatible with new samples (the N+1 problem). Full mechanism: variant-calling/joint-calling.

Primary path: BWA-MEM2 + bcftools

Fast, dependency-light, good for exploratory work and non-model organisms; weaker on indels in homopolymers than reassembly callers.

Step 1: QC/trim with fastp

fastp -i sample_R1.fastq.gz -I sample_R2.fastq.gz \
    -o trimmed/sample_R1.fq.gz -O trimmed/sample_R2.fq.gz \
    --detect_adapter_for_pe \
    --qualified_quality_phred 20 \
    --length_required 50 \
    --html qc/sample_fastp.html

Step 2: Align with BWA-MEM2

Read groups are mandatory; add -Y (soft-clip supplementary) if structural-variant calling is downstream, and -K 100000000 for thread-count-invariant output. Reference/analysis-set choice: read-alignment/bwa-alignment.

bwa-mem2 index reference.fa   # once
bwa-mem2 mem -t 8 -K 100000000 \
    -R "@RG\tID:sample\tSM:sample\tPL:ILLUMINA\tLB:lib1" \
    reference.fa trimmed/sample_R1.fq.gz trimmed/sample_R2.fq.gz \
  | samtools view -bS - > aligned/sample.bam

Step 3: Mark duplicates

Strict order (samtools convention): collate (name) -> fixmate -m -> sort (coordinate) -> markdup. Detail: alignment-files/duplicate-handling.

# collate groups mates by name; fixmate -m adds the ms/MC tags markdup needs; markdup needs coordinate order.
# Do NOT coordinate-sort before fixmate, and do NOT markdup amplicon/PCR data (use UMIs there).
samtools collate -@ 8 -O -u aligned/sample.bam \
  | samtools fixmate -m -@ 8 -u - - \
  | samtools sort -@ 8 -u - \
  | samtools markdup -@ 8 - aligned/sample.markdup.bam
samtools index aligned/sample.markdup.bam

Step 4: Call and normalize

# Single sample (mpileup passes MQ/BQ filters into the pileup)
# -a FORMAT/DP,FORMAT/AD + call -f GQ emit the per-sample DP/GQ the Step 5 genotype filter needs.
bcftools mpileup -Ou -f reference.fa -a FORMAT/DP,FORMAT/AD --max-depth 250 --min-MQ 20 --min-BQ 20 \
    aligned/sample.markdup.bam \
  | bcftools call -mv -f GQ -Oz -o variants/sample.vcf.gz

# Normalize BEFORE any annotation or cross-callset comparison, against the SAME reference
bcftools norm -m-any -f reference.fa -Oz -o variants/sample.norm.vcf.gz variants/sample.vcf.gz
bcftools index variants/sample.norm.vcf.gz

For multi-sample cohorts, bcftools can call several BAMs jointly, but the GATK/DeepVariant gVCF path is preferred at scale (variant-calling/joint-calling).

Step 5: Filter (site then genotype)

# Site-level (bcftools flags rather than removes, so failures stay auditable)
bcftools filter -Oz -s LowQual \
    -e 'QUAL<20 || INFO/DP<10 || MQ<30' \
    -o variants/sample.siteflt.vcf.gz variants/sample.norm.vcf.gz

# Genotype-level: set low-confidence genotypes to no-call (NOT 0/0)
bcftools filter -Oz -S . \
    -e 'FMT/GQ<20 | FMT/DP<8' \
    -o variants/sample.filtered.vcf.gz variants/sample.siteflt.vcf.gz
bcftools index variants/sample.filtered.vcf.gz

Alternative path: BWA-MEM2 + GATK HaplotypeCaller

Local reassembly + PairHMM; the auditable reference implementation, strong on indels. Full mechanism: variant-calling/gatk-variant-calling.

gatk CreateSequenceDictionary -R reference.fa
samtools faidx reference.fa

# DRAGEN mode: no BQSR, STR-aware indel model (DRAGSTR), improved QUAL calibration
gatk HaplotypeCaller -R reference.fa -I aligned/sample.markdup.bam \
    -O gvcf/sample.g.vcf.gz -ERC GVCF --dragen-mode

# Cohort: consolidate gVCFs then joint-genotype (NEVER bcftools merge single-sample VCFs)
gatk GenomicsDBImport --sample-name-map gvcf/map.txt \
    --genomicsdb-workspace-path genomicsdb -L intervals.bed
gatk GenotypeGVCFs -R reference.fa -V gendb://genomicsdb -O variants/cohort.vcf.gz

Filtering strategy depends on cohort size

The site-level filter is chosen by cohort size, platform, and organism; genotype-level filtering is always applied on top. Full mechanism and thresholds: variant-calling/filtering-best-practices.

SNPs and indels are filtered separately (different error processes, truth resources, abundance). Hard-filter starting points (SNPs QD<2, FS>60, MQ<40, MQRankSum<-12.5, ReadPosRankSum<-8, SOR>3; indels loosen FS>200, tighten ReadPosRankSum<-20) are lenient heuristics to tune, not universal truth. RankSum annotations are only defined at het sites -- a hand-written filter must treat a missing annotation as PASS, or every hom-alt site vanishes.