Running a Genome Analysis Pipeline on Your Own WGS Data

By the end of this you will have a CRAM aligned to GRCh38, and a single-sample gVCF and VCF with roughly 4.5–5 million small variants. You will also have a structural variant call set, an annotated TSV you can load into DuckDB or pandas, and a QC report that tells you whether any of it is trustworthy. You need: your raw reads (paired FASTQ or an unaligned BAM), about 500 GB of scratch disk, 32–64 cores, 64 GB RAM, and Docker or Singularity. Everything below runs on one machine. A 30x human genome is about 2 × 50 GB of gzipped FASTQ, roughly 70 GB as a BAM, and 15–20 GB as a CRAM. If you have 8 GB of RAM and a laptop SSD, stop here and rent a c7i.16xlarge for a day instead.

A note on what this guide does. It produces measurements: variant calls, quality metrics, annotations and allele frequencies. It does not tell you what any of it means for your health. Variant interpretation is a separate discipline with its own error modes, and clinical decisions belong with a clinician or a board-certified genetic counselor working from a CLIA-validated report ¹.

1. Verify what you were given

Before anything else, confirm the data is intact and is yours.

md5sum -c sample.md5
seqkit stats -a -j 8 sample_R1.fastq.gz sample_R2.fastq.gz

Check read length (150 bp is standard for Illumina WGS), read count, and that R1 and R2 have identical counts. At 30x coverage over a 3.1 Gb genome with 150 bp reads you expect roughly 310 million read pairs. If you got a CRAM or BAM instead of FASTQ:

samtools quickcheck -v sample.cram && echo OK
samtools view -H sample.cram | grep '^@RG'
samtools view -H sample.cram | grep '^@SQ' | head -3

The @SQ lines tell you which reference build was used. This matters more than anything else in the pipeline and is the single most common source of silent corruption downstream. chr1 with length 248956422 is GRCh38. 1 with length 249250621 is GRCh37/hg19.

2. Pick the reference, and pick it once

Use the GRCh38 full analysis set with decoy contigs and alt contigs, and use a BWA index built with alt-aware handling. The Broad’s resource bundle version is the one most tools expect:

gsutil -m cp -r gs://gcp-public-data--broad-references/hg38/v0/Homo_sapiens_assembly38.* .

That gives you the FASTA, .fai, .dict, the BWA index files, and the .alt file. Do not strip the alt contigs to “keep it simple.” Without them, reads from segmental duplications and the MHC pile onto the primary assembly and create false heterozygous calls in exactly the regions you care about.

T2T-CHM13 is a better assembly and fixes real gaps, but almost every annotation database, ClinVar coordinate, and gnomAD frequency is still keyed to GRCh38. Align to GRCh38 as primary. If you want CHM13 for the acrocentric arms and segmental duplications, run it as a second alignment rather than a replacement.

3. Align and mark duplicates

bwa-mem2 is a drop-in replacement for bwa mem with the same output and roughly 2–3x the speed. Pipe straight into sorting to avoid writing an unsorted 70 GB intermediate.

bwa-mem2 mem -t 32 -K 100000000 -Y \
  -R '@RG\tID:S1\tSM:SAMPLE1\tLB:lib1\tPL:ILLUMINA\tPU:unit1' \
  Homo_sapiens_assembly38.fasta sample_R1.fastq.gz sample_R2.fastq.gz \
| samtools sort -@ 8 -m 4G -o sample.sorted.bam -
samtools index -@ 8 sample.sorted.bam

-K 100000000 fixes the batch size so output is deterministic regardless of thread count. -Y uses soft clipping for supplementary alignments, which structural variant callers need. The read group is not optional: GATK refuses to run without SM, and downstream tools use PL to choose error models.

Then mark duplicates:

gatk --java-options "-Xmx32g" MarkDuplicates \
  -I sample.sorted.bam -O sample.md.bam -M sample.md.metrics.txt \
  --OPTICAL_DUPLICATE_PIXEL_DISTANCE 2500 \
  --CREATE_INDEX true

Expect 5–15% duplicates for a well-made PCR-free library. Above 25% means the library was over-amplified and your effective coverage is lower than the raw number suggests.

Runtime on 32 cores: about 3 hours for alignment, 1.5 hours for MarkDuplicates.

4. Decide whether you need BQSR

Base Quality Score Recalibration re-estimates per-base error rates using known variant sites as a mask. It is the most expensive step in the classic GATK pipeline relative to what it delivers, which is why there is a body of work on accelerating it in hardware ².

gatk BaseRecalibrator -I sample.md.bam -R Homo_sapiens_assembly38.fasta \
  --known-sites Homo_sapiens_assembly38.dbsnp138.vcf \
  --known-sites Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
  -O recal.table
gatk ApplyBQSR -I sample.md.bam -R Homo_sapiens_assembly38.fasta \
  --bqsr-recal-file recal.table -O sample.bqsr.bam

We skip BQSR when calling with DeepVariant, which learns its own error model from the pileup image and gains nothing from recalibrated qualities. We run it when calling with GATK HaplotypeCaller, where the model assumes calibrated base qualities. Skipping it saves about 2 hours and 70 GB of disk.

5. Call small variants

Two reasonable choices for a single germline sample.

DeepVariant, which is what we use by default:

docker run -v $PWD:/data google/deepvariant:1.6.1 \
  /opt/deepvariant/bin/run_deepvariant \
  --model_type=WGS \
  --ref=/data/Homo_sapiens_assembly38.fasta \
  --reads=/data/sample.md.bam \
  --output_vcf=/data/sample.dv.vcf.gz \
  --output_gvcf=/data/sample.dv.g.vcf.gz \
  --num_shards=32

GATK HaplotypeCaller in gVCF mode, which is what you want if you will later joint-call with family members or a cohort:

gatk --java-options "-Xmx16g" HaplotypeCaller \
  -R Homo_sapiens_assembly38.fasta -I sample.bqsr.bam \
  -O sample.g.vcf.gz -ERC GVCF \
  --native-pair-hmm-threads 8 \
  -L wgs_calling_regions.hg38.interval_list

Run HaplotypeCaller in parallel across intervals (one job per chromosome, or the 50-shard scatter list from the Broad bundle) and merge with GatherVcfs. Single-threaded across the whole genome it takes well over a day.

Why DeepVariant by default: for a single sample, its QUAL scores are calibrated well enough that a simple threshold works, and it is consistently better on indels in homopolymer and tandem repeat regions. Why HaplotypeCaller anyway: the gVCF and joint genotyping machinery is the only mature path to cohort calling, and if you ever add a parent or sibling, GenotypeGVCFs across the trio recovers de novo and compound heterozygous signal that no single-sample caller can see. DeepVariant also emits gVCF, and GLnexus joint-calls it, so this is less binary than it used to be.

Runtime for DeepVariant on 32 CPU cores: 4–6 hours. On a single A100, under an hour.

6. Filter

If you used DeepVariant, filtering is one line:

bcftools view -f PASS -i 'QUAL>=20' sample.dv.vcf.gz -Oz -o sample.pass.vcf.gz

If you used HaplotypeCaller on a single sample, VQSR needs more variants than one genome provides for stable indel models. Use the modern replacement or hard filters:

gatk VariantFiltration -V sample.snps.vcf.gz -O snps.filt.vcf.gz \
  --filter-name "QD2"    --filter-expression "QD < 2.0" \
  --filter-name "QUAL30" --filter-expression "QUAL < 30.0" \
  --filter-name "SOR3"   --filter-expression "SOR > 3.0" \
  --filter-name "FS60"   --filter-expression "FS > 60.0" \
  --filter-name "MQ40"   --filter-expression "MQ < 40.0" \
  --filter-name "MQRS-12.5" --filter-expression "MQRankSum < -12.5" \
  --filter-name "RPRS-8"    --filter-expression "ReadPosRankSum < -8.0"

Indels get a different set (QD < 2.0, FS > 200.0, ReadPosRankSum < -20.0). Split SNPs and indels with SelectVariants -select-type SNP first, filter separately, then merge.

7. QC the call set before you believe any of it

This is the step people skip, and it is the one that catches sample swaps and contamination. Run all of these:

bcftools stats -F Homo_sapiens_assembly38.fasta sample.pass.vcf.gz > stats.txt
mosdepth --by 1000 -n -t 8 sample sample.md.bam
verifybamid2 --SVDPrefix 1000g.phase3.100k.b38.vcf.gz.dat \
  --Reference Homo_sapiens_assembly38.fasta --BamFile sample.md.bam

What to expect for a clean 30x PCR-free genome:

• 4.4–5.0 million small variants total, depending on ancestry (more for African-ancestry genomes, fewer for European).
• Ti/Tv ratio between 1.95 and 2.10 genome-wide. Below 1.8 means false positives are leaking in.
• Heterozygous-to-homozygous-alt ratio around 1.5–1.7.
• VerifyBamID2 FREEMIX below 0.02. Above 0.03 and you have sample contamination, which inflates heterozygous calls everywhere.
• Mean coverage within 10% of what was promised, with at least 90% of the callable genome at 20x or better.

Also check sex: compute normalized coverage on chrX and chrY from the mosdepth output and confirm it matches the expected karyotype. A mismatch is a sample swap until proven otherwise. Any variant you plan to act on should be confirmed in a CLIA-certified laboratory, because research-grade calls have not been validated for clinical use ³. Public health and clinical labs formalize this with defined acceptance thresholds for coverage, depth, and concordance against a known control, and running a Genome in a Bottle sample through your exact pipeline is the cheapest way to find out what your own error rate is ⁴.

If you have GIAB HG002, benchmark:

hap.py HG002_GRCh38_benchmark.vcf.gz sample.pass.vcf.gz \
  -f HG002_GRCh38_benchmark.bed -r Homo_sapiens_assembly38.fasta \
  -o happy_out --engine=vcfeval --threads 16

A competent short-read pipeline lands above 99.5% F1 for SNVs and above 99% for indels inside the high-confidence regions. Outside those regions, which is roughly 8–12% of the genome, you have no idea how well you are doing.

8. Structural variants and copy number

Small variants are the easy part. Short-read structural variant calling is genuinely hard, and recall for insertions above 50 bp is poor no matter which caller you use.

configManta.py --bam sample.md.bam --referenceFasta Homo_sapiens_assembly38.fasta \
  --runDir manta_out
manta_out/runWorkflow.py -j 16

smoove call --name sample --fasta Homo_sapiens_assembly38.fasta \
  --outdir smoove_out -p 8 --genotype sample.md.bam

Take the intersection of Manta and a split-read/discordant-pair caller like smoove for deletions, and treat anything called by one tool only as a hypothesis. For copy number over larger intervals, cnvpytor on the read-depth signal is fast and adequate. If structural variation is the point, pay for long reads instead. Nanopore or PacBio HiFi at 20x resolves insertions, tandem repeat expansions, and phase in a way that no amount of short-read cleverness will.

9. Annotate

Ensembl VEP with a local cache is the tool we use. Run it offline so you are not rate-limited by the REST API.

vep --offline --cache --dir_cache $HOME/.vep --assembly GRCh38 \
  --fasta Homo_sapiens_assembly38.fasta \
  -i sample.pass.vcf.gz -o sample.vep.vcf.gz --vcf --compress_output bgzip \
  --everything --pick_allele_gene --fork 16 \
  --custom clinvar.vcf.gz,ClinVar,vcf,exact,0,CLNSIG,CLNREVSTAT,CLNDN \
  --plugin CADD,whole_genome_SNVs.tsv.gz,gnomad.genomes.r4.0.indel.tsv.gz \
  --plugin SpliceAI,snv=spliceai_scores.raw.snv.hg38.vcf.gz,indel=spliceai_scores.raw.indel.hg38.vcf.gz \
  --plugin AlphaMissense,file=AlphaMissense_hg38.tsv.gz

The cache is about 25 GB, CADD whole-genome scores about 80 GB, SpliceAI about 100 GB. Budget the disk.

Then flatten to something queryable:

bcftools +split-vep sample.vep.vcf.gz -d -f \
  '%CHROM\t%POS\t%REF\t%ALT\t%QUAL\t%SYMBOL\t%Consequence\t%IMPACT\t%gnomADg_AF\t%CLNSIG\t%CADD_PHRED\t%SpliceAI_pred_DS_AG[\t%GT]\n' \
  -i 'IMPACT="HIGH" || IMPACT="MODERATE"' > annotated.tsv

Load that into DuckDB and start filtering. The usual first pass: allele frequency below 0.1% in gnomAD, predicted loss of function or ClinVar pathogenic with two or more review stars, and restricted to a gene list. The ACMG secondary findings list (around 80 genes with actionable, high-penetrance disease associations) is a reasonable starting gene set because it is the one clinical labs already agreed is worth reporting.

Expect roughly 100–200 rare predicted loss-of-function variants in a healthy genome. Nearly all are benign, in genes with no disease association, or artifacts. This is the central difficulty of personal genomics: the sequencing is solved, and interpretation is the bottleneck ⁵. Effect sizes for most common variants are small, and functional annotation of noncoding variation remains largely incomplete ⁵.

10. Pharmacogenomics and polygenic scores

Star-allele calling for pharmacogenes needs phase and copy number, especially CYP2D6, which is subject to gene deletions, duplications, and hybrid alleles with CYP2D7. Generic VCF annotation gets this wrong.

docker run -v $PWD:/data pgkb/pharmcat \
  /pharmcat/pharmcat_pipeline /data/sample.pass.vcf.gz -o /data/pgx

PharmCAT produces diplotype calls and the corresponding CPIC guideline text. Read the report as a measurement of genotype. Any change to medication is a conversation with a prescriber, full stop.

For polygenic scores, pull scoring files from the PGS Catalog and use pgsc_calc. Interpret them with care: most scores were trained in European-ancestry cohorts and transfer poorly, and a percentile computed against the wrong reference population is meaningless.

Common problems

Reference mismatch. You align to GRCh38 and annotate with a GRCh37 ClinVar file. Coordinates shift by megabases and everything silently lands in the wrong gene. Check @SQ lengths and the VCF ##contig lines at every handoff.

Chromosome naming. chr1 versus 1 will cause VEP or bcftools to return zero annotations rather than an error. Fix with bcftools annotate --rename-chrs.

Missing read groups. GATK fails immediately, which is the good case. Some tools accept the BAM and produce quietly wrong output.

Low coverage sold as 30x. “30x” sometimes means raw yield before duplicate removal and before mapping quality filtering. Compute mean coverage on the deduplicated BAM with mosdepth and use that number.

Contamination. FREEMIX above 0.03 breaks everything downstream. Nothing in the pipeline will warn you. Run VerifyBamID2.

Believing calls in hard regions. Roughly 10% of the genome is not reliably callable with short reads: segmental duplications, the MHC, tandem repeats, centromeric and telomeric regions. Genes such as SMN1, PMS2, and CYP2D6 sit in or near these regions and need targeted methods. A “no variant found” in these genes means the method could not see. The result does not mean that nothing is there.

Treating a research pipeline as clinical. Even in clinical settings, WGS diagnostic yield depends heavily on phenotype-driven filtering by trained analysts, and the same variant list can support very different conclusions depending on the clinical question ³. The infrastructure burden of doing this properly, including storage, compute, versioned pipelines, and validated reference materials, is a real cost that clinical labs absorb and individuals usually underestimate ¹. Reproducibility comes from version-pinning every tool and database and recording the exact command lines, which is the design principle behind pipeline frameworks built for this purpose ⁶.

Cost. On a 64-core spot instance at roughly $0.70 per hour, the full path from FASTQ to annotated TSV runs about 10–14 hours of wall clock, so $10–20 of compute plus storage. The expensive part is the sequencing and the interpretation. The CPU time is not the expensive part.

Footnotes

1. Petar Brlek, Luka Bulić, Matea Bračić, et al. Implementing Whole Genome Sequencing (WGS) in Clinical Practice: Advantages, Challenges, and Future Perspectives. Cells, 2024. https://doi.org/10.3390/cells13060504 ↩ ↩²

2. Michael Lo, Zhenman Fang, Jie Wang, et al. Algorithm-Hardware Co-design for BQSR Acceleration in Genome Analysis ToolKit. 2020. https://doi.org/10.1109/fccm48280.2020.00029 ↩

3. Frederik Otzen Bagger, Line Borgwardt, Andreas Sand Jespersen, et al. Whole genome sequencing in clinical practice. BMC Medical Genomics, 2024. https://doi.org/10.1186/s12920-024-01795-w ↩ ↩²

4. Kelly F. Oakeson, Jennifer M. Wagner, Michelle Mendenhall, et al. Bioinformatic Analyses of Whole-Genome Sequence Data in a Public Health Laboratory. Emerging infectious diseases, 2017. https://doi.org/10.3201/eid2309.170416 ↩

5. Tuuli Lappalainen, Alexandra J. Scott, Margot Brandt, et al. Genomic Analysis in the Age of Human Genome Sequencing. Cell, 2019. https://doi.org/10.1016/j.cell.2019.02.032 ↩ ↩²

6. David Sims, Nicholas E. Ilott, Stephen N. Sansom, et al. CGAT: computational genomics analysis toolkit. Bioinformatics, 2014. https://doi.org/10.1093/bioinformatics/btt756 ↩