GPU-accelerated WGS: Alignment + DeepVariant (Parabricks)

A self-built whole-genome-sequencing pipeline: FASTQ → aligned BAM → variant VCF, GPU-accelerated end to end with NVIDIA Clara Parabricks (fq2bam + deepvariant), running in rootless containers. Two reference genomes are supported:

• hg38 (GRCh38, in hg38/)
• T2T CHM13v2.0 (in t2t/, optional auto-download)

Two entry points, same pipeline:

Script	Target	Container runtime
run_parabricks_wgs.py	Multi-GPU host / shared cluster	rootless Podman (NVIDIA CDI, or a site GPU wrapper)
run_local.py	Single-GPU laptop/workstation	Docker (--gpus all), tuned for 12 GB VRAM / 32 GB RAM

I built this before I knew nf-core/Sarek existed — essentially a from-scratch, GPU-first take on the same idea. See the comparison at the bottom for what Sarek does that this doesn't (yet).

No genomic data in this repo — only the scripts. FASTQ/BAM/VCF and the large reference files stay local (see .gitignore).

Why this exists

Consumer genotyping services (e.g. tellmeGen) hand you a finished VCF — but it's a SNP microarray: ~700k pre-selected positions, called against an older reference build (typically GRCh37/hg19). A tiny, fixed keyhole next to whole-genome sequencing (~4–5M variants, any position). Some services let you download your raw data, so I took mine and built the "serious" version myself: a full WGS alignment + variant call against modern references — including the bleeding-edge T2T-CHM13 assembly that didn't even exist when most array chips were designed. Equal parts curiosity and "because I can."

---

1) One-time setup

On a shared multi-GPU host (rootless Podman):

# rootless podman must be set up once for your user (site-specific);
# then verify:
podman version
podman info
podman run -it --rm alpine /bin/sh

GPU passthrough is auto-detected: if SITE_GPU_WRAPPER points at a cluster GPU-allocator script, its output is used; otherwise the script uses NVIDIA CDI (--device nvidia.com/gpu=all, rootless-friendly).

On a local laptop/workstation (Docker): install Docker + the NVIDIA Container Toolkit, then docker run --rm --gpus all nvidia/cuda:12.2.0-base-ubuntu22.04 nvidia-smi should list your GPU.

You also need the Parabricks image (one-time, ~20 GB, NVIDIA NGC account):

podman pull nvcr.io/nvidia/clara/clara-parabricks:4.5.0-1   # or: docker pull ...

---

2) Layout & data

Work inside the project directory; the scripts bind it into the container as /work and resolve everything relative to it:

gen/
├── fastq/                 # *_R1/_R2.fq.gz  (your reads — NOT in repo)
├── hg38/                  # reference FASTA + indices (NOT in repo)
├── t2t/                   # reference FASTA + indices (NOT in repo)
└── output/parabricks/     # BAM / VCF results (NOT in repo)

Copy reads + references onto the target host with rsync as usual; only the scripts themselves live in git.

---

3) Run it

Cluster / multi-GPU (Podman):

python3 run_parabricks_wgs.py --ref hg38 --gpus 2 --cpus 40

Interactive when parameters are missing; anything you pass as a flag is not prompted. Prints a short resource overview and a checkpoint/resume status.

With T2T (auto-download the reference):

python3 run_parabricks_wgs.py --ref t2t --download-t2t --gpus 2 --cpus 40

Prebuilt T2T index download instead of building locally:

python3 run_parabricks_wgs.py --ref t2t --download-t2t \
  --t2t-index-mode download --t2t-index-url https://<your-index-host>/t2t \
  --gpus 2 --cpus 40

Local laptop (Docker, single GPU):

python3 run_local.py --ref hg38            # auto-detects FASTQs, sensible defaults

run_local.py is tuned for a 12 GB-VRAM / 32 GB-RAM laptop: --low-memory, a RAM-capped sort, reduced BWA queue capacity and a single DeepVariant stream.

---

4) Output

Written to output/parabricks/:

• <sample>.bam (+ .bai)
• <sample>.recal — BQSR table (cluster script)
• <sample>.metrics — duplicate metrics
• <sample>.deepvariant.vcf.gz (+ .tbi)
• optional <sample>.deepvariant.g.vcf.gz (--emit-gvcf)

Clean up large intermediates when done — a WGS BAM is ~30–60 GB.

---

5) Resume / checkpointing

run_parabricks_wgs.py scans existing outputs, validates them, removes half-written files from aborted runs, and resumes at the last completed step (reference prep → fq2bam → DeepVariant). Re-running a finished sample is a no-op unless --force is given.

---

6) Parameter reference

Required: none (interactive). Useful options:

• --r1, --r2 (when several FASTQs are present)
• --gpus, --cpus
• --ref hg38|t2t, --sex male|female (T2T: male includes the Y chromosome)
• --download-t2t, --t2t-index-mode build|download, --t2t-index-url <URL>
• --emit-gvcf
• --podman-userns keep-id (if output files end up root-owned)
• --prepare-ref-only (download/index the reference, then exit)

For both hg38 and T2T: empty folder → offer to download + build all indices; FASTA present but indices missing → build what's missing; everything present → continue. Reference downloads try known mirrors (NCBI/Ensembl) first, then ask for a URL.

---

7) File types — what is what?

Input

File	Where	What it is
*_R1.fq.gz / *_R2.fq.gz	fastq/	FASTQ paired-end reads — millions of short DNA reads (4 lines each: @name, sequence, +, quality). ~2 × 40–50 GB for a WGS sample. R1 = forward, R2 = reverse of the same fragment.
trimmed_r1/r2.fq.gz	fastq/	fastp-trimmed reads — adapters and low-quality ends removed.

Reference

File	What it is
*.fasta	The full human reference (~3 Gbp, chr1…chrX/Y/M), uncompressed ~3 GB.
*.fasta.fai	FASTA index (samtools faidx) — byte offsets for random access.
*.dict	Sequence dictionary (samtools/picard) — chromosome names, lengths, MD5s; required by GATK/Parabricks.
*.fasta.{bwt,pac,sa,ann,amb}	BWA / BWA-MEM2 index — Burrows-Wheeler structures that let the aligner place reads in milliseconds.

Why two references? hg38 (GRCh38) is the well-tested classic with maximum tool compatibility. T2T CHM13v2.0 is the first truly gapless human reference (2022), covering ~200 Mb that hg38 lacks (centromeres, repeats). For male samples T2T uses chm13v2.0_maskedY_rCRS (CHM13 + HG002 Y, PAR-masked + rCRS mitochondria).

Output

File	What it is
<sample>.bam (+ .bai)	BAM — every read mapped to its position in the reference; ~30–60 GB.
<sample>.recal	BQSR recalibration table — corrects systematic sequencer quality errors.
<sample>.metrics	MarkDuplicates metrics — PCR-duplicate rate (typically 5–15 % for WGS).
<sample>.deepvariant.vcf.gz (+ .tbi)	VCF — the called variants vs. the reference (SNVs + indels, genotypes, quality).
<sample>.deepvariant.g.vcf.gz	gVCF (optional) — includes non-variant positions with confidence; needed for joint cohort calling.

---

8) Programs — what does what?

Program	Step	What it does
Parabricks fq2bam	1	GPU all-in-one: BWA-MEM alignment + sort + MarkDuplicates + BQSR in one container call (replaces 4+ GATK steps). ~1–2 h for a 30× WGS on 2× A100.
Parabricks deepvariant	2	GPU variant caller using a CNN over pileup images; strong on indels and hard regions. Output: VCF (+ optional gVCF).
BWA / BWA-MEM2	ref prep	Builds the reference index (.bwt, .pac, .sa, .ann, .amb).
samtools	ref prep	faidx → .fai, dict → .dict; general BAM/FASTA tooling.
picard CreateSequenceDictionary	ref prep	Fallback for .dict when samtools dict is unavailable.
fastp	manual, pre-pipeline	Adapter/quality trimming → trimmed_r1/r2.fq.gz. Not yet integrated.
FastQC	manual, QC	Per-FASTQ quality reports.

Infrastructure

Tool	What it does
Podman	Rootless container engine for the cluster script; no root needed. GPUs via NVIDIA CDI or a site wrapper.
Docker	Used by run_local.py on the laptop (--gpus all).
Site GPU wrapper (optional)	If SITE_GPU_WRAPPER is set, the script calls it to obtain the right --device args for a shared GPU allocator. Unset → CDI passthrough.
nvidia-smi	GPU status/utilisation shown before each run.

Container image

nvcr.io/nvidia/clara/clara-parabricks:4.5.0-1 — NVIDIA's official image with Parabricks, BWA-MEM2, samtools, picard and CUDA libs (~20 GB, pulled once).

---

9) Repo layout

Only the scripts are tracked. Everything large or data-bearing is gitignored:

gen/
├── run_parabricks_wgs.py   ← cluster pipeline (Podman)
├── run_local.py            ← laptop pipeline (Docker)
├── README.md
└── .gitignore              ← excludes fastq/, hg38/, t2t/, output/

Never committed: fastq/*.fq.gz (~90 GB), hg38/ + t2t/ FASTA & indices (~8 GB each), output/parabricks/ (BAM/VCF). No genomic data leaves your host.

---

10) Roadmap / comparison with nf-core/Sarek

Covers the core pipeline (alignment + BQSR + DeepVariant). Sensible additions inspired by nf-core/sarek:

Feature	Status	Note
fastp trimming step (--trim)	missing	runs manually today; integration would be cleaner
CRAM output instead of BAM	missing	~50–80 % smaller; Parabricks fq2bam supports it
Samplesheet/CSV input	missing	currently one sample per run; Sarek does cohort batches
MultiQC reporting	missing	aggregate FastQC + samtools stats + duplicate metrics
Interval-based parallelism	n/a	Parabricks handles this internally
Annotation (VEP/SnpEff)	missing	post-DeepVariant step
Checkpoint/resume	partial	own checkpoint logic; no Nextflow needed

---

License

MIT — see LICENSE.