Building a SNP Calling Pipeline for Your Own Whole Genome

Oak
wgs variant-calling gatk deepvariant bioinformatics
A glowing DNA strand passes through a row of sorting gates inside a long glass channel, separating points of light into trays.

By the end of this you will have a single-sample VCF (Variant Call Format) containing roughly four to five million single-nucleotide variants and around half a million short insertions and deletions, aligned to GRCh38, filtered, annotated, and benchmarked against a truth set so you know its precision and recall rather than guessing. You need paired-end FASTQ files from a 30x or deeper whole-genome run (typically 50 to 100 GB gzipped), a machine with at least 16 cores and 64 GB of RAM, about 1.5 TB of scratch disk, and Docker or Singularity. Everything below runs on a single workstation or one cloud instance. Cluster and cloud ports of these same stages exist and scale well, but the added complexity is mostly scheduling and data movement rather than science1.

A caution before the commands: this pipeline produces research-grade calls. If a variant in it looks medically consequential, it needs orthogonal confirmation in a clinical laboratory and interpretation with a genetics clinician. Pipelines disagree with each other at a rate that matters for single-variant decisions2.

1. Verify what you received

Start by confirming the files are intact and are what you think they are. Sequencing providers deliver either FASTQ (raw reads plus base quality strings) or an aligned CRAM. If you got CRAM and want to rerun alignment yourself, convert back to FASTQ rather than trusting someone else’s aligner settings.

md5sum -c checksums.md5
zcat sample_R1.fastq.gz | head -4
seqkit stats -a sample_R1.fastq.gz sample_R2.fastq.gz

Check three things in the seqkit output: read length (150 bp is standard on Illumina NovaSeq), total bases, and that R1 and R2 have identical read counts. Total bases divided by 3.1e9 gives your expected raw coverage. If you paid for 30x and see 78 Gbp, you have about 25x before duplicate removal, and you should say so to the provider.

fastqc -t 8 sample_R1.fastq.gz sample_R2.fastq.gz
multiqc .

Look at the per-base quality plot and the adapter content plot. Modern Illumina chemistry rarely needs trimming, and aggressive trimming costs you real alignments. We only trim when adapter content exceeds a few percent at the read ends, and then with fastp --detect_adapter_for_pe and nothing else enabled.

2. Choose and prepare the reference

Use GRCh38 with the analysis set that includes decoy contigs and excludes the false duplications on chromosome 21 and elsewhere. The specific build we use is GCA_000001405.15_GRCh38_no_alt_analysis_set.fna. The “no_alt” part matters: alternate haplotype contigs create multi-mapping reads that lower mapping quality and silently drop variants unless your aligner is ALT-aware and your caller respects it. Using the no-alt analysis set avoids that whole class of problem.

wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_GRCh38/seqs_for_alignment_pipelines.ucsc_ids/GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz
gunzip GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz
bwa-mem2 index GCA_000001405.15_GRCh38_no_alt_analysis_set.fna
samtools faidx GCA_000001405.15_GRCh38_no_alt_analysis_set.fna
gatk CreateSequenceDictionary -R GCA_000001405.15_GRCh38_no_alt_analysis_set.fna

Indexing with bwa-mem2 takes about 40 minutes and needs roughly 80 GB of RAM at peak, which is the one step where a 64 GB machine may fail. If it does, download a prebuilt index or use plain bwa for indexing. T2T-CHM13 is a better assembly, but the annotation and benchmarking ecosystem still centers on GRCh38, so we stay there and treat CHM13 as a second pass for specific hard regions.

3. Align

bwa-mem2 is the same algorithm as bwa mem with a faster implementation, roughly two to three times the throughput for identical output. Set the read group at alignment time, because GATK and most downstream tools refuse to run without it.

bwa-mem2 mem -t 30 -K 100000000 -Y \
  -R '@RG\tID:SAMPLE1_L1\tSM:SAMPLE1\tLB:lib1\tPL:ILLUMINA\tPU:HXXXX.1.NNNNNN' \
  GCA_000001405.15_GRCh38_no_alt_analysis_set.fna \
  sample_R1.fastq.gz sample_R2.fastq.gz \
| samtools sort -@ 8 -m 3G -o sample.sorted.bam -
samtools index -@ 8 sample.sorted.bam

-K 100000000 fixes the chunk size so output is deterministic regardless of thread count, which you want when comparing runs. -Y uses soft clipping for supplementary alignments, which keeps sequence available for structural variant callers later. On 30 cores this takes four to six hours for a 30x genome and writes a BAM of 60 to 80 GB.

4. Mark duplicates and check contamination

PCR and optical duplicates inflate apparent read support and bias allele fractions. Mark them rather than removing them, so the information stays recoverable.

gatk MarkDuplicatesSpark \
  -I sample.sorted.bam \
  -O sample.md.bam \
  -M sample.dup_metrics.txt \
  --conf 'spark.executor.cores=30' \
  --conf 'spark.local.dir=/scratch/tmp'

Expect a duplicate rate between 2 and 12 percent for PCR-free libraries at 30x. Rates above 20 percent mean low library complexity: the same molecules sequenced repeatedly, which caps your effective coverage no matter how many reads you bought. Duplicate handling becomes far more consequential when you care about low-frequency variants, where unique molecular identifiers and a barcode-aware caller are the only reliable way to separate true low-allele-fraction signal from amplification and sequencing noise3. For germline SNP calling at 30x, standard duplicate marking is sufficient.

Then check coverage and contamination:

mosdepth -t 8 -n --fast-mode --by 1000 sample sample.md.bam
verifybamid2 --SVDPrefix resource/1000g.phase3.100k.b38.vcf.gz.dat \
  --Reference GCA_000001405.15_GRCh38_no_alt_analysis_set.fna \
  --BamFile sample.md.bam --NumThread 8

Two thresholds we treat as hard gates: median autosomal coverage at or above 25x, and VerifyBamID2 FREEMIX below 0.02. FREEMIX estimates the fraction of reads from another individual. Above 2 percent, heterozygous calls start accumulating in regions where the sample is homozygous, and no downstream filter cleanly removes them.

5. Call variants

Two reasonable choices, and we recommend different ones depending on what you plan to do next.

For a single genome you will analyze on its own, use DeepVariant. It is a convolutional neural network that classifies pileup images and consistently produces higher precision and recall on Illumina WGS than the traditional statistical callers, particularly for indels.

docker run -v "$PWD":/data google/deepvariant:1.6.1 \
  /opt/deepvariant/bin/run_deepvariant \
  --model_type=WGS \
  --ref=/data/GCA_000001405.15_GRCh38_no_alt_analysis_set.fna \
  --reads=/data/sample.md.bam \
  --output_vcf=/data/sample.dv.vcf.gz \
  --output_gvcf=/data/sample.dv.g.vcf.gz \
  --num_shards=30 \
  --intermediate_results_dir=/scratch/dv_tmp

This runs eight to twelve hours on 30 CPU cores, or roughly 90 minutes with an A100. Always request the gVCF output as well as the VCF, even if you do not need it today. A gVCF records confidence at every reference position, not only at variant sites, which is what makes it possible to distinguish “reference” from “no data” later and to join your genome with others in a single genotyping step4.

For a genome you expect to joint-call with family members, use GATK HaplotypeCaller in gVCF mode:

gatk --java-options "-Xmx48g -XX:ParallelGCThreads=4" HaplotypeCaller \
  -R GCA_000001405.15_GRCh38_no_alt_analysis_set.fna \
  -I sample.md.bam \
  -O sample.g.vcf.gz \
  -ERC GVCF \
  -G StandardAnnotation -G AlleleSpecificAnnotation \
  --native-pair-hmm-threads 8

HaplotypeCaller is single-threaded in practice, so split by interval and run 20 to 30 parallel jobs, then merge. Joint calling across relatives is worth the trouble: family structure lets you identify Mendelian-inconsistent calls and estimate per-sample error rates directly from the data rather than assuming them5.

6. Filter

An unfiltered single-sample VCF has a false discovery rate high enough to ruin any downstream analysis. DeepVariant’s own quality score is well calibrated, so filtering is simple:

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

For GATK output from a single sample, Variant Quality Score Recalibration is unreliable because it needs many samples to model the annotation distributions. Use hard filters instead, applied separately to SNPs and indels:

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

Indels take looser thresholds: QD < 2.0, FS > 200.0, ReadPosRankSum < -20.0. Merge the two sets back with gatk MergeVcfs.

Beyond per-variant filters, mask regions where short-read sequencing is systematically unreliable. Certain loci in the human genome produce consistently biased base calls across platforms and pipelines, and no amount of coverage fixes them, because the error is in mapping and in the chemistry rather than in sampling6. Intersect your calls with the GIAB high-confidence regions to see which of your variants sit inside territory anyone can call well:

bcftools view -R GRCh38_notinalldifficultregions.bed.gz \
  sample.dv.pass.vcf.gz -Oz -o sample.easy.vcf.gz

Keep both files. The full callset is what you analyze; the masked one is what you trust without confirmation.

7. Benchmark the pipeline, not the sample

You cannot measure accuracy on your own genome because there is no truth set for it. Run the entire pipeline once on a Genome in a Bottle reference sample, for which a curated truth VCF exists, and carry those numbers forward as your pipeline’s expected performance.

# Download HG002 30x Illumina FASTQ, run steps 3-6 identically, then:
rtg vcfeval \
  -b HG002_GRCh38_1_22_v4.2.1_benchmark.vcf.gz \
  -c hg002.dv.pass.vcf.gz \
  -e HG002_GRCh38_1_22_v4.2.1_benchmark_noinconsistent.bed \
  -t GRCh38.sdf \
  -o vcfeval_out --vcf-score-field=QUAL

With DeepVariant on 30x Illumina you should see SNV F-measure above 0.995 and indel F-measure around 0.99 inside high-confidence regions. GATK with hard filters typically lands a few tenths of a percent lower on SNVs and noticeably lower on indels. If your numbers are far off, the problem is in your pipeline, not the caller.

Then run the cheap self-consistency checks on your own sample:

bcftools stats sample.dv.pass.vcf.gz | grep -E '^SN|^TSTV'

Expect 4.3 to 5.0 million SNVs against GRCh38, a transition/transversion ratio between 2.0 and 2.1 genome-wide, a heterozygous to non-reference homozygous ratio between 1.4 and 2.0 (it varies with ancestry), and under 1 percent of SNVs absent from dbSNP. A Ti/Tv below 1.9 is the single most sensitive indicator that false positives have leaked through your filters. Coverage depth drives most of what you can and cannot see, and the gap between the raw yield you paid for and the usable, uniquely-mapped, duplicate-free coverage you have is where most surprises live7.

8. Annotate

Convert coordinates into biology with Ensembl VEP, run locally with a cached database so nothing leaves your machine.

vep -i sample.dv.pass.vcf.gz -o sample.vep.vcf.gz \
  --vcf --compress_output bgzip --offline --cache --dir_cache $HOME/.vep \
  --assembly GRCh38 --fork 16 --everything \
  --plugin CADD,whole_genome_SNVs.tsv.gz \
  --plugin SpliceAI,snv=spliceai_scores.raw.snv.hg38.vcf.gz \
  --custom gnomad.genomes.v4.1.sites.vcf.gz,gnomADg,vcf,exact,0,AF

The gnomAD allele frequency annotation is the most useful single field. Filtering to variants with population frequency below 0.001 and a predicted high-impact consequence reduces four million calls to a few hundred worth reading. Those few hundred still include many artifacts, which is why the benchmarking step matters and why clinical interpretation belongs with a clinician.

Common problems

Read group missing or malformed. GATK fails with a message about no read groups, or, worse, MarkDuplicates succeeds but the sample name in the VCF is wrong. Set -R at alignment time and verify with samtools view -H sample.md.bam | grep '^@RG'.

Contig naming mismatch. A VCF with chr1 and a reference with 1 produces empty output rather than an error in several tools. Standardize once, at the reference download, and never mix.

Running out of scratch. MarkDuplicatesSpark and DeepVariant both write large intermediates. Point spark.local.dir and --intermediate_results_dir at real fast storage with 500 GB free, not at /tmp.

Ti/Tv near 1.5. Your filters are too loose, or contamination is present. Check FREEMIX first, then tighten QUAL.

Excess heterozygosity in specific regions. Segmental duplications and the pseudoautosomal regions collect reads from paralogous sequence. Those calls are mapping artifacts. The difficult-region BED files exist for exactly this reason.

Comparing two pipelines and finding disagreement. This is expected. Platform and pipeline choices shift concordance by several percent even on the same DNA, and the discordant calls concentrate in indels and in low-complexity sequence2. Decide which pipeline you trust by benchmarking both on a reference sample rather than by counting variants.

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

  1. Alberto Mulone, Sherine Awad, Davide Chiarugi, et al. Porting the Variant Calling Pipeline for NGS data in cloud-HPC environment. 2023 IEEE 47th Annual Computers, Software, and Applications Conference (COMPSAC), 2023. https://doi.org/10.1109/compsac57700.2023.00288 ↩

  2. Yanping Sun, Xiaochao Zhao, Xue Fan, et al. Assessing the impact of sequencing platforms and analytical pipelines on whole-exome sequencing. Frontiers in Genetics, 2024. https://doi.org/10.3389/fgene.2024.1334075 ↩ ↩2

  3. Chang Xu, Mohammad R. Nezami Ranjbar, Zhong Wu, et al. Detecting very low allele fraction variants using targeted DNA sequencing and a novel molecular barcode-aware variant caller. BMC Genomics, 2017. https://doi.org/10.1186/s12864-016-3425-4 ↩

  4. Areeya Disratthakit, Licht Toyo-oka, Penpitcha Thawong, et al. An optimized genomic VCF workflow for precise identification of Mycobacterium tuberculosis cluster from cross-platform whole genome sequencing data. Infection, Genetics and Evolution, 2020. https://doi.org/10.1016/j.meegid.2019.104152 ↩

  5. Kelley Paskov, Jae-Yoon Jung, Brianna Chrisman, et al. Estimating sequencing error rates using families. BioData Mining, 2021. https://doi.org/10.1186/s13040-021-00259-6 ↩

  6. Timothy M. Freeman, Genomics England Research Consortium, Dennis Wang, et al. Genomic loci susceptible to systematic sequencing bias in clinical whole genomes. Genome Research, 2020. https://doi.org/10.1101/gr.255349.119 ↩

  7. D. C. Koboldt, L. Ding, E. R. Mardis, et al. Challenges of sequencing human genomes. Briefings in Bioinformatics, 2010. https://doi.org/10.1093/bib/bbq016 ↩