How to Build a Variant Calling Pipeline for Your Own Genome

By the end of this you will have several outputs from a 30x paired-end whole genome. The first is a CRAM aligned to the GRCh38 no-alt analysis set. Next come a gVCF and a filtered single-sample VCF of SNVs and short indels. Last are precision and recall numbers measured against a Genome in a Bottle truth set, and a VEP-annotated table you can query. Budget roughly 8-14 hours of wall clock on 32 cores and 64 GB RAM, plus about 400 GB of scratch disk for a 30x genome. The disk holds FASTQ at ~100 GB gzipped, an intermediate BAM at ~120 GB and a final CRAM at ~18 GB. You need Linux, Docker or Singularity, and the tools below.

• bwa-mem2 (or dragmap)
• samtools/bcftools ≥1.19
• mosdepth
• VerifyBamID2
• DeepVariant ≥1.6
• GATK ≥4.5
• hap.py (or rtg vcfeval)
• Ensembl VEP with a cached GRCh38 release

Write every step into a Makefile, Snakemake, or Nextflow file from the start, because you will rerun this: pipelines that are configuration-driven rather than hand-typed are the difference between a result you can reproduce next year and a directory of mystery BAMs ¹.

1. Verify what you were given

Before anything else, confirm the data is yours, complete, and not truncated.

md5sum -c raw.md5
seqkit stats -a -j 8 *.fastq.gz     # read count, length distribution, Q20/Q30

Check that read 1 and read 2 have identical read counts and that the mean length matches what you paid for (typically 150 bp). If you were delivered a CRAM instead of FASTQ, do not reheader and move on. Check the reference it was aligned against:

samtools head sample.cram | grep '^@SQ' | head -3
samtools head sample.cram | grep '^@RG'

Vendors ship alignments against at least four incompatible GRCh38 flavors (with alts, without, with or without the HLA and decoy contigs, and occasionally hs37d5). If the @SQ lines list chr1_KI270706v1_random and HLA contigs, you have an alt-containing build. If you plan to benchmark against GIAB or use a public annotation resource, realign. Converting CRAM back to FASTQ is fine as long as you restore pairing:

samtools collate -@16 -u -O sample.cram | \
  samtools fastq -@8 -1 r1.fq.gz -2 r2.fq.gz -0 /dev/null -s /dev/null -n

2. Pick the reference and pick it once

Use the GRCh38 no-alt analysis set: GCA_000001405.15_GRCh38_no_alt_analysis_set.fna. It masks the false duplications on chr21 and the GRC-flagged errors, includes the decoy (chrUn_KN707.../hs38d1) and HLA contigs, and is what DeepVariant’s released models and the GIAB benchmark BEDs assume. Alt-aware alignment with alt contigs is a real thing bwa supports, but it changes MAPQ behavior and most downstream callers and filters were not tuned for it.

bwa-mem2 index -p GRCh38 GCA_000001405.15_GRCh38_no_alt_analysis_set.fna
samtools faidx GCA_000001405.15_GRCh38_no_alt_analysis_set.fna

Indexing needs about 80 GB of RAM for bwa-mem2. If you do not have it, use bwa mem (slower by roughly 2x) or build the index on a large cloud instance once and copy it.

We would also, eventually, call against T2T-CHM13v2.0 as a second pass. It resolves regions GRCh38 simply does not represent, and you get better results in segmental duplications and the acrocentric arms. But keep GRCh38 as the primary, because nearly every annotation database is keyed to it.

3. Align, fix mates, mark duplicates

One streamed command, no intermediate BAM on disk:

RG='@RG\tID:L001\tSM:SUBJECT01\tPL:ILLUMINA\tPU:HXXXXX.1.NNNNNN\tLB:lib1'

bwa-mem2 mem -t 30 -K 100000000 -Y -R "$RG" GRCh38 r1.fq.gz r2.fq.gz \
  | samtools fixmate -m -u -@4 - - \
  | samtools sort -@6 -m 3G -T /scratch/sort -u - \
  | samtools markdup -@4 --reference GRCh38.fa -f markdup.stats -O cram -o sample.cram -
samtools index -@8 sample.cram

Notes on the flags that matter. -K 100000000 fixes the chunk size so output is deterministic across thread counts. -Y uses soft clipping for supplementary alignments, which structural variant callers want. -m on fixmate adds the MC/ms tags markdup needs. Set SM to a stable subject identifier, because every downstream tool keys on it. Do not use a file name.

Skip base quality score recalibration unless you are running GATK HaplotypeCaller on a non-Illumina platform or a very low-coverage sample. On modern 30x Illumina data BQSR moves single-sample F1 by a fraction of a percent and costs an hour plus a full BAM rewrite. DeepVariant does not use recalibrated qualities at all.

4. QC the alignment before you call anything

Three numbers decide whether calling is worth doing.

samtools stats -@8 --reference GRCh38.fa sample.cram > samtools.stats
mosdepth --fast-mode --by 1000 -t 4 -f GRCh38.fa sample sample.cram
VerifyBamID --SVDPrefix 1000g.phase3.100k.b38.vcf.gz.dat \
  --Reference GRCh38.fa --BamFile sample.cram --NumThread 8

What we look for: mean autosomal depth ≥28x, at least 90% of the callable genome at ≥20x, duplicate rate under 10%, insert size mode between 300 and 500 bp with a single clean peak, and FREEMIX (contamination estimate) under 0.01. A FREEMIX of 0.02-0.05 produces a characteristic signature of low-allele-fraction heterozygous calls that will survive naive filters and look like mosaicism. Above 0.05, stop and get the library rerun.

Also compute the X and Y coverage ratio from the mosdepth per-1kb output and check it against expectation. It is the cheapest sample-swap detector you have.

5. Call small variants with DeepVariant

For a single human Illumina genome, use DeepVariant. It has better indel accuracy than HaplotypeCaller out of the box, needs no recalibration step, and its confidence score is a calibrated probability rather than a heuristic annotation stack.

docker run --rm -v "$PWD":/data google/deepvariant:1.6.1 \
  /opt/deepvariant/bin/run_deepvariant \
  --model_type=WGS \
  --ref=/data/GRCh38.fa \
  --reads=/data/sample.cram \
  --output_vcf=/data/sample.dv.vcf.gz \
  --output_gvcf=/data/sample.dv.g.vcf.gz \
  --num_shards=30 \
  --intermediate_results_dir=/scratch/dv_tmp \
  --logging_dir=/data/logs

On 30 cores and 30x WGS this takes roughly 5-8 hours, with make_examples dominating. --intermediate_results_dir on fast local disk matters: the example TFRecords are hundreds of gigabytes of small writes. If you are on PacBio HiFi data use --model_type=PACBIO, and for Oxford Nanopore use the ONT_R104 model.

Emit the gVCF even for a single sample. You need it to distinguish “reference” from “no data” at any position someone later asks you about, and you need it if you ever add a family member and do joint calling through GLnexus.

Run GATK too, if only as a second opinion in regions you care about:

gatk --java-options "-Xmx16g" HaplotypeCaller \
  -R GRCh38.fa -I sample.cram -O sample.hc.g.vcf.gz \
  -ERC GVCF -G StandardAnnotation -G AS_StandardAnnotation \
  --native-pair-hmm-threads 8

Then intersect. Variants called by both callers with concordant genotypes are a much cleaner set than either alone. Variants unique to one caller are where you should be skeptical, and they cluster: homopolymers, short tandem repeats, and segmental duplications.

6. Know what your reads cannot see

Read-based calling against a linear reference has structural blind spots, and it is worth being explicit about them rather than discovering them by surprise. Assembly-based approaches call variants in regions where short-read alignment is unreliable, and direct comparisons show alignment-based and assembly-based callers agree on the bulk of SNVs while each recovering variants the other misses, with the discordance concentrated in repetitive and divergent regions ². Long reads push this further: nanopore assemblies plus long-read variant calling recover pathogenic variants in rare disease cases that short-read pipelines miss, including in segmental duplications and medically relevant genes with paralogs ³.

Practical consequence for a single short-read genome: treat calls in the regions below as unresolved until confirmed by a method designed for them.

• SMN1/SMN2
• CYP2D6
• HBA1/HBA2
• PMS2
• GBA
• IKBKG
• the HLA region

Exclude the low-complexity and segmental duplication BEDs from any broad screening query, then look at them separately and deliberately.

7. Benchmark against a truth set

This is the step most tutorials skip, and it is the step that tells you whether your numbers mean anything. If you have sequenced a GIAB sample (HG001-HG007) you can measure directly. If you have not, run your exact pipeline on a public HG002 FASTQ set once and carry those metrics as your pipeline’s characteristics.

docker run --rm -v "$PWD":/data pkrusche/hap.py /opt/hap.py/bin/hap.py \
  /data/HG002_GRCh38_1_22_v4.2.1_benchmark.vcf.gz \
  /data/sample.dv.vcf.gz \
  -f /data/HG002_GRCh38_1_22_v4.2.1_benchmark_noinconsistent.bed \
  -r /data/GRCh38.fa \
  -o /data/happy_out \
  --engine=vcfeval --engine-vcfeval-template=/data/GRCh38.sdf \
  --threads 16 --pass-only

Use --engine=vcfeval rather than the default xcmp: it does proper haplotype-aware comparison, so an indel written at a different position with the same resulting haplotype counts as a match instead of a false positive plus a false negative. Read happy_out.summary.csv. On 30x Illumina with DeepVariant you should see SNV F1 around 99.5% and indel F1 in the high 99s inside the v4.2.1 benchmark regions. If your SNV F1 is below 99%, something is wrong upstream: wrong reference flavor, contamination, or truncated FASTQ.

Two things to internalize about those numbers. First, they apply only inside the benchmark BED, which covers roughly 90% of the genome and deliberately excludes the hard parts. Your real genome-wide error rate is worse. Second, 99.5% SNV precision on ~4 million variant sites still leaves thousands of false positives, which is why single unreplicated calls should never drive a conclusion.

8. Normalize and filter

Normalize before you compare anything to anything. Multi-allelic records and non-left-aligned indels are the single largest source of spurious mismatches between your VCF and any annotation database.

bcftools norm -m -any -f GRCh38.fa -Ou sample.dv.vcf.gz \
  | bcftools norm -d exact -Ou \
  | bcftools view -f PASS -Oz -o sample.norm.vcf.gz
bcftools index -t sample.norm.vcf.gz

Then a working filter for a high-confidence set:

bcftools view -i 'FMT/GQ>=20 && FMT/DP>=10 && FMT/DP<=120' \
  -T ^GRCh38_difficult_regions.bed.gz \
  -Oz -o sample.hiconf.vcf.gz sample.norm.vcf.gz

The upper depth bound catches collapsed repeats and mitochondrial-like pileups. The -T ^ exclusion uses the GIAB stratification BEDs (GRCh38_alllowmapandsegdupregions.bed.gz is the one we reach for most). Keep the unfiltered file. Filtering is a view. It does not destroy anything.

9. Annotate

vep --offline --cache --dir_cache /ref/vep --assembly GRCh38 \
  --fasta GRCh38.fa --input_file sample.norm.vcf.gz --format vcf \
  --vcf --compress_output bgzip --output_file sample.vep.vcf.gz \
  --fork 12 --buffer_size 20000 \
  --everything --pick_allele_gene --hgvs --symbol --canonical \
  --custom file=/ref/gnomad.genomes.v4.1.sites.GRCh38.vcf.gz,short_name=gnomADg,format=vcf,type=exact,fields=AF%AF_grpmax%nhomalt \
  --custom file=/ref/clinvar_20250601.vcf.gz,short_name=ClinVar,format=vcf,type=exact,fields=CLNSIG%CLNREVSTAT%CLNDN \
  --plugin CADD,snv=/ref/whole_genome_SNVs.tsv.gz,indels=/ref/gnomad.genomes.r4.0.indel.tsv.gz \
  --plugin SpliceAI,snv=/ref/spliceai_scores.raw.snv.hg38.vcf.gz,indel=/ref/spliceai_scores.raw.indel.hg38.vcf.gz

Then flatten to something queryable:

bcftools +split-vep -d -f '%CHROM\t%POS\t%REF\t%ALT\t%SYMBOL\t%Consequence\t%HGVSc\t%HGVSp\t%gnomADg_AF\t%ClinVar_CLNSIG\t%ClinVar_CLNREVSTAT\t%CADD_PHRED\t%SpliceAI_pred_DS_AG\t[%GT\t%GQ\t%DP]\n' \
  -A tab sample.vep.vcf.gz | bgzip > sample.annot.tsv.gz

Load that into DuckDB or Parquet and query it. A 30x genome yields roughly 4-5 million SNVs and 800k-1M indels. After restricting to coding consequences with gnomAD allele frequency below 0.1% you are looking at a few hundred rows, which is a size a human can read.

10. Read the output with the right priors

Two filters do most of the interpretive work. First, ClinVar review status: CLNREVSTAT of criteria_provided,_single_submitter is a very different claim from reviewed_by_expert_panel, and single-submitter pathogenic assertions reclassify to benign often enough that you should treat them as hypotheses. Second, penetrance. Almost every effect size you will read about was estimated in ascertained cohorts (people who presented with disease), and the same variant measured in an unselected population cohort typically shows substantially lower penetrance ⁴. A pathogenic-labeled variant found by screening an asymptomatic person carries a different implication than the same variant found in someone with matching symptoms and family history.

Polygenic scores deserve the same skepticism. They are computed from the same VCF trivially, and their predictive value degrades sharply when the score was trained in a population with different ancestry than yours.

If anything in this output looks clinically consequential, the next step is a clinical-grade confirmation from an accredited laboratory and a conversation with a genetic counselor or physician. Research-grade calls, including everything produced by the pipeline above, are not a diagnosis, and we are not qualified to give you one.

Common problems

Wrong reference, silently. Your VCF has 1 instead of chr1, or your annotation cache is GRCh37. Symptom: VEP annotates almost nothing, or gnomAD AF is null everywhere. Check bcftools view -h contig names and lengths against the reference .fai before debugging anything else.

Non-normalized VCF compared to a normalized one. Symptom: hap.py reports thousands of indel false positives clustered in homopolymers. Run both sides through bcftools norm -f, and use --engine=vcfeval.

Duplicate rate above 20%. Usually a low-complexity library over-amplified during PCR. Effective coverage is much lower than nominal, and genotypes at heterozygous sites get noisy. Marking duplicates does not recover the lost information.

DeepVariant out of disk. make_examples writes far more intermediate data than the final VCF suggests. Point --intermediate_results_dir at a volume with 300+ GB free, and do not put it on a network filesystem.

Sex chromosome genotypes look wrong. DeepVariant and HaplotypeCaller handle ploidy on chrX and chrY differently, and calls in the pseudoautosomal regions are diploid in one and not the other. Check PAR boundaries explicitly before drawing any conclusion from a chrX call.

Mitochondrial calls treated like nuclear calls. chrM is high copy number and heteroplasmic. Standard diploid genotyping produces meaningless GT fields there. Use a dedicated path (GATK’s Mutect2 in mitochondria mode) and report heteroplasmy fractions. Do not report genotypes.

Non-reproducible results. Different thread counts, different tool minor versions, unpinned reference downloads. Pin container digests, record every tool version in the VCF header, and drive the whole thing from one config file ¹.

Footnotes

1. Manuel Rueda, Dietmar Fernandez Orth, Ivo G Gut. CBIcall: a configuration-driven framework for variant calling in large sequencing cohorts. Bioinformatics Advances, 2026. https://doi.org/10.1093/bioadv/vbag232 ↩ ↩²

2. Leihong Wu, Gökhan Yavaş, Huixiao Hong, et al. Direct comparison of performance of single nucleotide variant calling in human genome with alignment-based and assembly-based approaches. Scientific Reports, 2017. https://doi.org/10.1038/s41598-017-10826-9 ↩

3. Shloka Negi, Sarah L. Stenton, Seth Berger, et al. Advancing long-read nanopore genome assembly and accurate variant calling for rare disease detection. The American Journal of Human Genetics, 2025. https://doi.org/10.1016/j.ajhg.2025.01.002 ↩

4. 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 ↩