--- name: bio-workflows-fastq-to-variants description: Orchestrates the end-to-end germline short-variant pipeline from FASTQ to a filtered, normalized, benchmarked VCF, chaining QC/trim, BWA-MEM2 alignment, duplicate marking, optional BQSR, calling (bcftools/GATK HaplotypeCaller/DeepVariant/DRAGEN), normalization, site+genotype filtering, annotation, and hap.py/vcfeval benchmarking. Use when deciding the pipeline-wide reference-genome commitment (GRCh38 analysis set vs T2T, ALT/decoy handling), sequencing the steps in the defensible order (normalize BEFORE annotate, filter site- then genotype-level), choosing the calling engine and single-sample vs cohort joint-calling, picking a filtering strategy by cohort size, or benchmarking stratified within GIAB confident regions. Hands off mechanism to the variant-calling and read-alignment component skills; not a re-teach of any single step. tool_type: cli primary_tool: bcftools workflow: true depends_on: - read-qc/fastp-workflow - read-alignment/bwa-alignment - alignment-files/alignment-sorting - alignment-files/duplicate-handling - variant-calling/variant-calling - variant-calling/joint-calling - variant-calling/variant-normalization - variant-calling/filtering-best-practices - variant-calling/variant-annotation - variant-calling/vcf-statistics qc_checkpoints: - after_qc: "Q30 >85%, adapter content <1%" - after_alignment: "Mapping rate >95%, properly paired >90%" - after_dedup: "Duplication rate <30% for WGS, <50% for exome" - after_calling: "Ti/Tv ratio ~2.0-2.1 for WGS, ~3.0-3.3 for exome; dbSNP overlap >95% only after annotating the ID column" --- ## 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: ` --version` then ` --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. | Choice | Commit to it when | Consequence inherited downstream | |--------|-------------------|----------------------------------| | GRCh38 analysis set + decoys (hs38DH), ALT-aware (bwa-postalt) | Human germline, research or most clinical | Decoys soak up off-target reads; ALT-aware mapping recovers reads in MHC/segdup loci that ALT-unaware mapping force-fits to the primary, inflating false positives | | GRCh38 analysis set, ALT-unaware (primary only) | Frozen clinical pipeline needing deterministic simplicity | Simpler and validated, but loses signal in ~5 Mb of ALT-bearing loci | | Masked GRCh38 (false-duplication fix) | Calling in CBS, U2AF1, KCNE1B, KCNJ18 and other affected genes | Recovers reads whose mapQ collapsed across the phantom duplicate copy | | T2T-CHM13 | Research needing segdups/centromeres/dark genes; maximum accuracy | Reveals variants in newly resolved regions and removes GRCh38 false-duplication artifacts, but no lossless liftover to GRCh37/38, so the entire annotation/interpretation stack must be revalidated (Nurk 2022 *Science* 376:44-53; Aganezov 2022 *Science* 376:eabl3533) | 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. | Situation | Lean toward | Hand off to | |-----------|-------------|-------------| | Auditable open-source, large cohort, joint calling | GATK HaplotypeCaller GVCF -> GenomicsDBImport -> GenotypeGVCFs | variant-calling/gatk-variant-calling, variant-calling/joint-calling | | Best indel/difficult-region accuracy, single or cohort | DeepVariant (+ GLnexus for cohorts) | variant-calling/deepvariant | | Quick/exploratory, non-model organism, limited compute | bcftools mpileup + call | variant-calling/variant-calling | | Maximum throughput on Illumina, hardware available | DRAGEN (or GATK `--dragen-mode` for the open equivalent) | variant-calling/variant-calling | **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 ```bash 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. ```bash 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. ```bash # 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 ```bash # 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) ```bash # 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. ```bash 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. | Cohort / data | Site-level filter | Why | |---------------|-------------------|-----| | Large WGS cohort (~30+ jointly genotyped) | VQSR (or AS_VQSR for huge cohorts) | The Gaussian-mixture model needs tens of thousands of variants and truth-resource overlap to fit; unreliable below that | | Single sample / small cohort | GATK hard filters or VETS/NVScoreVariants | VQSR is non-identifiable on few variants; a "converged" model on one exome is filtering on noise | | Exome specifically | Hard filters (do NOT use DP as a VQSR annotation) | Capture-boundary coverage cliffs break the annotation manifold | | Non-model organism | Hard filters or a bootstrapped truth set | No HapMap/Omni/Mills truth resources exist | | DeepVariant / DRAGEN output | Use the caller's own calibration; do NOT re-apply GATK hard filters | Their error modes differ; classic annotations do not describe them | 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. ## Benchmarking the pipeline **Goal:** a defensible accuracy statement, not a single number. The only rigorous way to compare a callset to truth is haplotype-aware, stratified, and confined to the truth set's confident region. Two VCFs can encode the identical haplotype with different records, so a naive `bcftools isec`/line-diff overcounts errors; use `hap.py` wrapping the `vcfeval` engine, which replays variants onto the reference and matches at the haplotype level (Krusche 2019 *Nat Biotechnol* 37:555-560; GIAB truth, Zook 2019 *Nat Biotechnol* 37:561-566). ```bash # Only meaningful when the sample IS a GIAB genome (HG001-HG007) with a truth VCF + confident BED. # -f = confident/callable region BED (TP/FP/FN counted ONLY inside it; calls outside are UNK, not FP) hap.py truth.vcf.gz query.norm.vcf.gz \ -f HG002_confident.bed \ -r reference.fa \ -o bench/hg002 \ --engine=vcfeval \ --stratification stratification.tsv # GIAB region BEDs: low-complexity, segdup, MHC, GC-extreme ``` Discipline that separates a senior benchmark from a naive one: - **Refuse the global F1.** Report SNP and INDEL separately, and show the low-complexity/segmental-duplication/MHC rows explicitly -- hiding them behind an all-regions average is the most common soft cheat. - **Do not benchmark an ML caller only on its training genome.** DeepVariant and DRAGEN-ML train on GIAB coordinates, so scoring on HG002 alone partly measures memorization; score on a held-out or semi-blinded sample (HG003/HG004). - **"We found variants GIAB missed" is almost always a category error** -- calls outside the confident region, which GIAB declined to adjudicate, not accuracy. Real hard-region claims cite CMRG or an assembly-based benchmark. - **When the sample is not a GIAB genome** (the usual case), there is no truth VCF; fall back to proxy QC -- Ti/Tv (~2.0-2.1 WGS, ~3.0-3.3 exome), dbSNP overlap, het/hom by ancestry, and trio Mendelian concordance if a family is available. These are proxies for the absence of a benchmark, not a substitute for one (variant-calling/vcf-statistics). ## QC checkpoints between steps | After | Gate | Interpretation | |-------|------|----------------| | QC/trim | Q30 >85%, adapter <1% | DNA is typically higher quality than RNA | | Alignment | Mapped >95%, properly paired >90% (`samtools flagstat`) | Low mapping rate: wrong reference or contamination | | Dedup | Duplicates <30% WGS, <50% exome | High duplication: PCR over-amplification, low input | | Calling | Ti/Tv ~2.0-2.1 WGS, ~3.0-3.3 exome; dbSNP overlap >95% | Ti/Tv sliding toward 0.5 (random) signals false-positive inflation -- filters too loose. `bcftools stats` counts known sites from the ID column, and callers leave it as `.`: run `bcftools annotate -c ID -a dbsnp.vcf.gz` first or the dbSNP-overlap line reads 0% regardless of quality | ## Common Errors | Symptom | Cause | Fix | |---------|-------|-----| | Annotation reports a known pathogenic variant as novel/absent | Annotated before normalizing; indel one base off the database coordinate | `bcftools norm -m-any -f ref` against the SAME build as the annotation DB, BEFORE annotation | | Cohort has impossible all-hom-ref genotypes at variant sites | Built the cohort by `bcftools merge` of single-sample VCFs | Emit per-sample gVCFs and joint-genotype (variant-calling/joint-calling) | | Every hom-alt site filtered out | Hand-written filter treats missing RankSum as failing | Treat missing annotation as PASS; RankSum is defined only at het sites | | VQSR "converged" on one exome but the callset is garbage | VQSR needs tens of thousands of variants across ~30+ samples | Use hard filters or VETS/NVScoreVariants for single samples/exomes | | Spurious variants in CBS/U2AF1/KCNE1B | GRCh38 false duplications collapse mapQ | Use a masked GRCh38 or T2T-CHM13; commit the reference before calling | | GATK error "sample ... has no read group" | Read groups omitted at alignment | Re-run `bwa-mem2 mem -R "@RG\t..."` (SM/ID/PL/LB) | | Different variant counts from vt vs bcftools on the same data | vt decomposes MNPs by default, bcftools does not | Standardize ONE normalization tool + flags across every cohort compared (variant-calling/variant-normalization) | ## Pipeline map (hand-offs) - read-qc/fastp-workflow -- QC/trim options and report interpretation - read-alignment/bwa-alignment -- BWA-MEM2 parameters, read groups, ALT/decoy analysis set, reference indexing - alignment-files/duplicate-handling -- the collate/fixmate/sort/markdup order and UMI cases - variant-calling/variant-calling -- engine selection (bcftools vs GATK vs DeepVariant vs DRAGEN) and ploidy - variant-calling/gatk-variant-calling -- HaplotypeCaller, GVCF, BQSR/DRAGSTR, edge cases - variant-calling/deepvariant -- CNN calling, platform models, DeepTrio, GLnexus cohorts - variant-calling/joint-calling -- per-sample gVCF -> cohort joint genotyping, the N+1 problem, scaling - variant-calling/variant-normalization -- left-align/parsimony, multiallelic split, MNP decomposition - variant-calling/filtering-best-practices -- VQSR vs hard vs ML by cohort size; site vs genotype filters - variant-calling/variant-annotation -- VEP/SnpEff on the normalized representation - variant-calling/vcf-statistics -- Ti/Tv, het/hom, contamination/relatedness identity QC The complete runnable scripts for both paths are in this skill's examples/ (`bwa_bcftools_workflow.sh`, `bwa_gatk_workflow.sh`). ## Related Skills - database-access/sra-data - Pull public FASTQ for reanalysis (ENA mirror or STRIDES cloud) - database-access/ncbi-datasets-cli - Pull reference genome assembly via Datasets v2 CLI - read-qc/fastp-workflow - Detailed QC options - sequence-io/fastq-quality - Confirm the FASTQ quality encoding (Phred+33 vs legacy Phred+64/Solexa) before trimming public or pre-2011 data - sequence-io/paired-end-fastq - Keep R1/R2 mates synchronized; independent per-mate filtering silently desyncs pairs - read-alignment/bwa-alignment - BWA-MEM2 parameters, read groups, ALT/decoy analysis set, the dedup ordering - alignment-files/duplicate-handling - Duplicate marking details - variant-calling/variant-calling - Engine selection and bcftools calling options - variant-calling/gatk-variant-calling - GATK HaplotypeCaller and DRAGEN mode - variant-calling/deepvariant - Deep-learning calling and GLnexus cohorts - variant-calling/joint-calling - Cohort joint genotyping and scaling - variant-calling/variant-normalization - Normalize before annotate/compare - variant-calling/filtering-best-practices - VQSR, hard filters, VETS - variant-calling/variant-annotation - Annotate variants with VEP - variant-calling/vcf-statistics - Ti/Tv, het/hom, and identity QC ## References - Krusche P, Trigg L, Boutros PC, et al. (GA4GH Benchmarking Team). Best practices for benchmarking germline small-variant calls in human genomes. *Nature Biotechnology* 37:555-560 (2019). DOI 10.1038/s41587-019-0054-x. Stratified haplotype-aware benchmarking (hap.py/vcfeval). - Zook JM, McDaniel J, Olson ND, et al. An open resource for accurately benchmarking small variant and reference calls. *Nature Biotechnology* 37:561-566 (2019). DOI 10.1038/s41587-019-0074-6. GIAB truth set + confident regions. - DePristo MA, Banks E, Poplin R, et al. A framework for variation discovery and genotyping using next-generation DNA sequencing data. *Nature Genetics* 43:491-498 (2011). DOI 10.1038/ng.806. GATK framework. - Van der Auwera GA, Carneiro MO, Hartl C, et al. From FastQ Data to High-Confidence Variant Calls: The Genome Analysis Toolkit Best Practices Pipeline. *Current Protocols in Bioinformatics* 43:11.10.1-11.10.33 (2013). DOI 10.1002/0471250953.bi1110s43. - Poplin R, Ruano-Rubio V, DePristo MA, et al. Scaling accurate genetic variant discovery to tens of thousands of samples. *bioRxiv* 201178 (2018). DOI 10.1101/201178. Preprint only (never journal-published); the GVCF/joint-genotyping reference. - Poplin R, Chang P-C, Alexander D, et al. A universal SNP and small-indel variant caller using deep neural networks. *Nature Biotechnology* 36:983-987 (2018). DOI 10.1038/nbt.4235. DeepVariant. - Yun T, Li H, Chang P-C, et al. Accurate, scalable cohort variant calls using DeepVariant and GLnexus. *Bioinformatics* 36:5582-5589 (2020). DOI 10.1093/bioinformatics/btaa1081. - Danecek P, Bonfield JK, Liddle J, et al. Twelve years of SAMtools and BCFtools. *GigaScience* 10:giab008 (2021). DOI 10.1093/gigascience/giab008. - Nurk S, Koren S, Rhie A, et al. The complete sequence of a human genome. *Science* 376:44-53 (2022). DOI 10.1126/science.abj6987. T2T-CHM13. - Aganezov S, Yan SM, et al. A complete reference genome improves analysis of human genetic variation. *Science* 376:eabl3533 (2022). DOI 10.1126/science.abl3533. Reference-choice variant-calling payoff.