Calling SNPs From Your Own Whole-Genome Sequencing Data

By the end of this guide you will have a single-sample VCF containing roughly four to five million variant sites relative to the human reference, of which the large majority are single-nucleotide variants, along with a measured estimate of how accurate that call set is and an annotation layer telling you which gene each variant falls in and how common it is in population databases. You need three things: your sequencing data (paired-end FASTQ, or a CRAM/BAM someone else aligned for you), a machine with at least 32 cores, 64 GB of RAM, and about 500 GB of free disk for a 30x genome, and Docker or Singularity. Everything below runs on Linux. A 30x human genome takes roughly four hours to align and two to four hours to call on that hardware, so plan for a day of wall-clock time including the mistakes.

1. Fix the vocabulary before you touch the data

A SNP is a single-nucleotide polymorphism: one position in the genome where the base differs between individuals, and where the alternate allele is common enough in the population to be called a polymorphism rather than a private change. Historically the cutoff was 1% minor allele frequency. The distinction between “mutation” and “polymorphism” has never been clean, and as personal genome sequencing became routine it became actively misleading, since “mutation” acquired a connotation of pathogenicity that frequency alone does not justify ¹. Most working pipelines sidestep the argument entirely and use “variant”: a difference between your sequence and the reference at a given position. That is what the V in VCF stands for, and the term SNV (single-nucleotide variant) is the frequency-agnostic version of SNP.

Two other distinctions matter for choosing tools. Germline variant calling assumes you are looking at a diploid individual whose variants are present in every cell at either 50% or 100% allele fraction, so the caller can use an explicit genotype model with three states. Somatic variant calling drops that assumption: a tumor subclone or a mosaic tissue can carry a variant at 3% allele fraction, which means the caller must distinguish real low-fraction signal from sequencing error, usually with a matched normal sample for comparison. The tools are different, the coverage requirements are much higher, and nothing in this guide is designed for somatic work. Finally, SNPs arise from point substitutions, which fall into the familiar categories of synonymous (no amino acid change), missense (one amino acid swapped), and nonsense (a premature stop), with insertions and deletions forming the separate indel class that most callers emit in the same VCF.

The density of SNPs across the genome is not uniform, and neither is the linkage between them, because local variation reflects the history of mutation and recombination at that locus rather than a uniform process ². This matters practically: your call set will have regions of high heterozygosity where calling is hard and regions that are essentially invariant, and quality metrics computed genome-wide hide both.

2. Inventory what you were given

Before computing anything, verify the files. Sequencing vendors ship corrupted transfers more often than you would expect, and a truncated FASTQ produces a call set that looks plausible until you benchmark it.

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

You are checking three numbers. Read count should match between R1 and R2 exactly. Read length should be what you paid for, typically 150 bp. Total bases divided by 3.1e9 gives your raw coverage, and you want that comfortably above 30 for germline calling, because post-alignment and post-duplicate-removal coverage will be 10 to 20 percent lower. If you were given a CRAM instead, run samtools quickcheck -v on it and confirm the reference it was aligned against by reading the @SQ lines in samtools view -H, since a CRAM without its exact reference is unreadable.

3. Pick a reference build and commit to it

Use the GRCh38 analysis set with decoy contigs, distributed by the Broad as Homo_sapiens_assembly38.fasta. The decoys absorb reads from sequences absent from the primary assembly that would otherwise pile up as false positives in paralogous regions, and the analysis set has the problematic PAR regions on chrY hard-masked so that pseudoautosomal reads map once rather than twice.

The alternative worth knowing about is T2T-CHM13, which closes the centromeres and acrogenic short arms. It gives better calls in those regions but breaks compatibility with essentially every annotation resource and benchmark set, so we would call against GRCh38 first and treat T2T as a second pass if you care about segmental duplications. A third option is a graph reference, which encodes known alternate alleles as paths rather than forcing every read onto a single linear sequence, and this measurably reduces reference bias in variable regions ³. Graph pipelines are more work to run and harder to compare against published call sets, so we would reach for one only when the linear pipeline has visibly failed in a region you care about.

bwa-mem2 index Homo_sapiens_assembly38.fasta
samtools faidx Homo_sapiens_assembly38.fasta

Indexing needs about 80 GB of RAM for bwa-mem2 and takes an hour. Do it once.

4. Align, sort, and mark duplicates

bwa-mem2 mem -t 32 -K 100000000 -Y \
  -R '@RG\tID:sample1\tSM:sample1\tPL:ILLUMINA\tLB:lib1' \
  Homo_sapiens_assembly38.fasta sample_R1.fastq.gz sample_R2.fastq.gz \
| samtools fixmate -m -@ 4 - - \
| samtools sort -@ 8 -m 4G -T /tmp/sort -o sample.sorted.bam -

samtools markdup -@ 8 --reference Homo_sapiens_assembly38.fasta \
  -s -f markdup.stats.txt sample.sorted.bam sample.markdup.cram
samtools index sample.markdup.cram

Three flags earn their place. -K 100000000 fixes the batch size so that the output is deterministic regardless of thread count, which you will want the first time two runs disagree. -Y uses soft clipping for supplementary alignments, which downstream structural variant callers need and which costs nothing here. The read group line is not optional: every caller downstream reads SM to name the sample, and a missing @RG causes failures several hours later.

Read markdup.stats.txt when it finishes. A PCR-free library should show a duplicate rate under 5 percent, and a PCR library 8 to 15 percent. Above 25 percent means the library was over-amplified and your effective coverage is far below nominal. Then check coverage properly:

mosdepth -t 8 -n --fast-mode --by 1000 sample sample.markdup.cram

The sample.mosdepth.summary.txt file gives mean depth per chromosome. Chromosome X at half the autosomal depth indicates a male sample, and chrY near zero with chrX at full depth indicates a female sample. If those two disagree with what you expect, stop and work out whether samples were swapped.

5. Call germline variants

We use DeepVariant. It replaces the hand-built statistical models of earlier callers with a convolutional network over pileup images, and on short-read data it produces fewer false positives than the alternatives without per-dataset filter tuning.

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.markdup.cram \
  --output_vcf=/data/sample.dv.vcf.gz \
  --output_gvcf=/data/sample.dv.g.vcf.gz \
  --num_shards=32 \
  --intermediate_results_dir=/data/dv_tmp

Set --model_type=WES for exome or targeted panels, since the coverage distribution and error profile differ enough that the WGS model underperforms. Always emit the gVCF as well as the VCF. The gVCF records confidence at non-variant positions, which is the only way to distinguish “reference at this position” from “no data at this position,” and you will need that distinction the first time someone asks whether you carry a specific allele.

The alternative is GATK HaplotypeCaller followed by VQSR or hard filtering. It remains a reasonable choice, especially if you need to joint-call across a cohort. Expect the two callers to disagree on a meaningful fraction of sites even on the same BAM, with the discordance concentrated in low-coverage and repetitive regions rather than spread evenly ⁴. If you run both, treat the intersection as high-confidence and the symmetric difference as the set that needs manual review in IGV.

Coverage is the single largest determinant of how good this step is. Below roughly 10x, genotype calls at heterozygous sites become unreliable because one of the two alleles is frequently unobserved, and the useful approach shifts toward imputation using a reference panel and calling likelihoods rather than hard genotypes ⁵. If you are working with long reads instead, use a caller built for their error profile: PEPPER-Margin-DeepVariant phases reads with Margin and feeds haplotype-tagged pileups to DeepVariant, reaching SNV accuracy on nanopore data comparable to short reads ⁶.

6. Normalize and filter

Raw caller output is not directly comparable to anything else until it is left-aligned and split into biallelic records.

bcftools norm -m -any -f Homo_sapiens_assembly38.fasta -Ou sample.dv.vcf.gz \
| bcftools norm -d exact -Oz -o sample.norm.vcf.gz
bcftools index -t sample.norm.vcf.gz

bcftools view -f PASS -i 'FORMAT/GQ>=20 && FORMAT/DP>=10' \
  -Oz -o sample.pass.vcf.gz sample.norm.vcf.gz
bcftools index -t sample.pass.vcf.gz

Then compute the sanity metrics that tell you whether the call set is sane before you look at any individual variant:

bcftools stats -F Homo_sapiens_assembly38.fasta sample.pass.vcf.gz > sample.stats
plot-vcfstats -p plots/ sample.stats

Three numbers matter. The transition/transversion ratio for genome-wide SNVs should land near 2.0 to 2.1, and for coding regions near 3.0. A Ti/Tv below 1.8 genome-wide means false positives are contaminating the set, because sequencing errors are roughly random across substitution types while real variation is strongly biased toward transitions. Heterozygous to homozygous-alternate ratio should be around 1.5 to 1.6 for a typical outbred individual, with lower values suggesting either consanguinity or allelic dropout. Total SNV count should be in the 3.5 to 4 million range against GRCh38, with substantial variation by ancestry.

7. Benchmark against a truth set

This step is the difference between a call set you trust and one you hope about. The Genome in a Bottle consortium publishes high-confidence truth VCFs and confident-region BED files for seven reference samples. You cannot benchmark your own genome directly, since no truth set exists for it, so benchmark the pipeline by running a GIAB sample’s public FASTQ through the identical commands.

docker run -v "$PWD":/data jmcdani20/hap.py:v0.3.12 /opt/hap.py/bin/hap.py \
  /data/HG002_GRCh38_1_22_v4.2.1_benchmark.vcf.gz \
  /data/hg002.pass.vcf.gz \
  -f /data/HG002_GRCh38_1_22_v4.2.1_benchmark_noinconsistent.bed \
  -r /data/Homo_sapiens_assembly38.fasta \
  --engine=vcfeval \
  --stratification /data/GRCh38-all-stratifications.tsv \
  -o /data/happy_out --threads 32

For 30x PCR-free Illumina with the pipeline above, SNV F1 inside the confident regions should exceed 0.999, and indel F1 should be in the high 0.99s. If your numbers are materially worse, the problem is in alignment or the input data, not in the caller. The stratification output is the part people skip and should not: it breaks accuracy down by region class, and you will find that performance in low-mappability regions, homopolymers longer than 6 bp, and segmental duplications is dramatically worse than the headline figure. When you later look up a specific variant, check which stratum it falls in before believing it.

8. Annotate, then be careful about what annotation means

vep --cache --offline --assembly GRCh38 --fork 16 \
  --input_file sample.pass.vcf.gz --format vcf \
  --output_file sample.vep.vcf.gz --vcf --compress_output bgzip \
  --everything --check_existing \
  --custom gnomad.genomes.v4.1.sites.vcf.gz,gnomADg,vcf,exact,0,AF,AF_popmax \
  --plugin CADD,whole_genome_SNVs.tsv.gz

--everything turns on the full consequence set including SIFT, PolyPhen, and regulatory overlap. The gnomAD custom annotation is the most useful single field: allele frequency in a large population cohort tells you immediately that a variant flagged as damaging by a prediction tool is carried by 3 percent of the population and therefore almost certainly not responsible for a rare severe phenotype.

The interpretive step is where the difficulty concentrates. A typical exome contains hundreds of rare, protein-altering variants in any individual, and the filtering strategies that narrow that to a candidate list depend heavily on assumptions about inheritance mode and phenotype specificity ⁷. Study design determines what can be concluded far more than the sequencing does, and single-sample analyses without family or cohort context are especially prone to overcalling causality ⁸. Even in carefully curated clinical settings, classifying a variant as pathogenic requires evidence beyond computational prediction: functional data, segregation, and curated database entries ⁹. If you find something in your own data that looks medically consequential, the correct next step is a clinical-grade confirmatory test through a certified laboratory and a conversation with a genetic counselor or physician. Research-grade pipelines, including this one, are not a diagnostic instrument.

9. Phase the calls

Genotypes tell you that you carry two alternate alleles in a gene. Phasing tells you whether they sit on the same chromosome copy or opposite copies, which changes the interpretation completely.

whatshap phase -o sample.phased.vcf.gz \
  --reference Homo_sapiens_assembly38.fasta \
  --indels sample.pass.vcf.gz sample.markdup.cram
whatshap stats --gtf sample.phased.gtf sample.phased.vcf.gz

With 150 bp short reads, read-backed phasing produces short blocks, typically tens of kilobases, so most gene-level questions stay unresolved. Statistical phasing against a reference panel with SHAPEIT5 extends blocks chromosome-wide at the cost of switch errors in rare haplotypes. Long reads solve the problem properly, and this is the strongest practical argument for including them in a personal sequencing plan.

Common problems

Callers running out of memory during the make_examples stage is almost always a region of pathological coverage, usually rDNA or a mitochondrial NUMT pileup. Confirm with mosdepth and exclude the region rather than raising the memory limit.

Ti/Tv near 1.5 with a normal duplicate rate points to contamination between samples. Run VerifyBamID2 and check the estimated contamination fraction, since anything above 0.02 will inject false heterozygous calls throughout the genome.

A VCF that hap.py refuses to compare usually has multi-allelic records or unnormalized indel representations. Re-run the bcftools norm step in section 6 on both the query and the truth set before concluding the calls are wrong.

Missing calls at a position you know you carry are frequently a coverage gap rather than a reference genotype. Open the gVCF and check the reference band: if the depth is zero, the correct statement is that the position was not assayed.

Finally, disagreement between two callers on the same BAM is normal rather than a sign of failure, and the disagreement clusters in indels and low-complexity sequence ⁴. Resolve individual cases by looking at the pileup in IGV rather than by trusting whichever caller you ran last.

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. Roshan Karki, Deep Pandya, Robert C. Elston, et al. Defining “mutation” and “polymorphism” in the era of personal genomics. BMC Medical Genomics, 2015. https://doi.org/10.1186/s12920-015-0115-z ↩

2. David E. Reich, Stephen F. Schaffner, Mark J. Daly, et al. Human genome sequence variation and the influence of gene history, mutation and recombination. Nature Genetics, 2002. https://doi.org/10.1038/ng947 ↩

3. Goran Rakocevic, Vladimir Semenyuk, Wan-Ping Lee, et al. Fast and accurate genomic analyses using genome graphs. Nature Genetics, 2019. https://doi.org/10.1038/s41588-018-0316-4 ↩

4. Charles D Warden, Aaron W Adamson, Susan L Neuhausen, et al. Detailed comparison of two popular variant calling packages for exome and targeted exon studies. 2014. https://doi.org/10.7287/peerj.preprints.403v3 ↩ ↩²

5. Chris Bizon, Michael Spiegel, Scott A Chasse, et al. Variant calling in low-coverage whole genome sequencing of a Native American population sample. BMC Genomics, 2014. https://doi.org/10.1186/1471-2164-15-85 ↩

6. Kishwar Shafin, Trevor Pesout, Pi-Chuan Chang, et al. Haplotype-aware variant calling with PEPPER-Margin-DeepVariant enables high accuracy in nanopore long-reads. Nature Methods, 2021. https://doi.org/10.1038/s41592-021-01299-w ↩

7. Jacek Majewski, Jeremy Schwartzentruber, Emilie Lalonde, et al. What can exome sequencing do for you?. Journal of Medical Genetics, 2011. https://doi.org/10.1136/jmedgenet-2011-100223 ↩

8. David B. Goldstein, Andrew Allen, Jonathan Keebler, et al. Sequencing studies in human genetics: design and interpretation. Nature Reviews Genetics, 2013. https://doi.org/10.1038/nrg3455 ↩

9. Nicholas M. Murphy, Tanya S. Samarasekera, Lisa Macaskill, et al. Genome sequencing of human in vitro fertilisation embryos for pathogenic variation screening. Scientific Reports, 2020. https://doi.org/10.1038/s41598-020-60704-0 ↩