SNP Analysis on Your Own Whole-Genome Data

By the end of this you will have a per-sample VCF of roughly 3.5–4 million single-nucleotide variants called against GRCh38, filtered to a callset you can defend, annotated with gene consequence, gnomAD allele frequency and ClinVar assertions, and loaded into a columnar store you can query in a second or two. You need: your raw reads (FASTQ or an unaligned/aligned CRAM), about 300 GB of free disk, 32 GB of RAM, 8+ cores, and a few hours of wall clock. Everything here runs on a laptop with enough disk or on one cloud VM. Nothing here is a clinical result. If a variant looks consequential, the next step is an accredited clinical lab and a genetic counselor. A bigger script is not the answer.

A word on what “SNP analysis” means in practice. The term covers two different jobs that people conflate. The first is genotyping: you have a fixed list of known polymorphic positions and you measure which alleles you carry. That is what consumer arrays do, at 600k–900k sites. The second is discovery plus genotyping: you sequence the whole genome and determine the genotype at every position, including sites nobody has catalogued. The original human SNP map contained 1.42 million SNPs, about one per 1.9 kb of sequence, which was the state of the reference catalogue when array design started.¹ Your own 30x genome will call more variants than that entire first map. The work now is deciding which of your calls are real and which of the real ones mean anything. Finding variation is no longer the hard part.

1. Verify the raw data before you spend compute on it

Check checksums against whatever the lab gave you, then confirm the files are structurally complete. A truncated FASTQ or CRAM will happily align for two hours and then die.

md5sum -c raw.md5
samtools quickcheck -v *.cram && echo "cram ok"
seqkit stats -a -j 8 sample_R1.fastq.gz sample_R2.fastq.gz

You want read length (usually 150 bp), read count, and total yield. For 30x on a 3.1 Gb genome you need about 93 Gbp of passing bases, so roughly 310 million read pairs at 2x150. If seqkit stats reports 200 million pairs, you have ~19x, not 30x, and your heterozygous call sensitivity will be materially worse. Ask the lab for the yield per lane before you accept the delivery.

Run fastqc or fastp --json and look at two things only: per-base quality tail (Q30 fraction over the whole run should be above ~85%) and duplicate rate. Adapter content matters less than people think because modern aligners soft-clip, but a >20% duplicate rate means the library was over-amplified from too little input DNA and your effective coverage is lower than the file size suggests.

2. Align to GRCh38 with the analysis set, alt-aware

Use the full analysis set reference with decoys and HLA contigs (GCA_000001405.15_GRCh38_no_alt_analysis_set.fna if you want to skip ALT handling, or the _plus_hs38d1 version with the ALT-aware bwa postprocessing script). We use the no-ALT analysis set. ALT-aware alignment is technically better in the MHC and a handful of other regions, but it complicates every downstream tool and the gain is invisible unless you care specifically about HLA, which you should type with a dedicated tool anyway.

bwa-mem2 mem -t 16 -K 100000000 -Y \
  -R '@RG\tID:rg1\tSM:SAMPLE\tPL:ILLUMINA\tLB:lib1' \
  GRCh38_no_alt_analysis_set.fna sample_R1.fastq.gz sample_R2.fastq.gz \
| samblaster \
| samtools sort -@ 4 -m 2G -o sample.sorted.bam -
samtools index -@ 8 sample.sorted.bam
samtools view -@ 8 -T GRCh38_no_alt_analysis_set.fna -C -o sample.cram sample.sorted.bam

-K 100000000 fixes the chunk size so the output is deterministic across thread counts. -Y uses soft clipping for supplementary alignments, which structural variant callers expect. The read group is not optional: every caller downstream keys on SM.

If you were delivered a hg19/GRCh37 BAM, realign from scratch rather than lifting over. CrossMap and Picard LiftoverVcf will silently drop or flip 1–2% of variants at the coordinate level, and liftover of a VCF cannot recover variants in regions where the reference base itself changed.

Store the final artifact as CRAM. A 30x genome is ~75 GB as BAM and ~17 GB as CRAM referenced to GRCh38, with no information loss if you keep all quality bins.

3. Call variants

Two reasonable choices.

DeepVariant is what we use for a single sample. It is a CNN over pileup images, it needs no base quality score recalibration, and on the GIAB HG002 benchmark it sits above 99.6% F1 for SNVs at 30x. It has one configuration knob that matters (--model_type), which removes an entire class of user error.

run_deepvariant \
  --model_type=WGS \
  --ref=GRCh38_no_alt_analysis_set.fna \
  --reads=sample.cram \
  --output_vcf=sample.dv.vcf.gz \
  --output_gvcf=sample.dv.g.vcf.gz \
  --num_shards=16 \
  --intermediate_results_dir=/scratch/dv

GATK HaplotypeCaller is the alternative, and it is the right choice if you are calling a family jointly or want to fit into an existing GATK-based workflow. Its architecture is local de novo assembly of active regions followed by a PairHMM likelihood model, with base quality recalibration upstream and variant quality recalibration downstream. That is the framework that defined how this is done.²

gatk --java-options "-Xmx24g" HaplotypeCaller \
  -R GRCh38_no_alt_analysis_set.fna -I sample.cram \
  -O sample.g.vcf.gz -ERC GVCF \
  --native-pair-hmm-threads 8
gatk GenotypeGVCFs -R GRCh38_no_alt_analysis_set.fna \
  -V sample.g.vcf.gz -O sample.gatk.vcf.gz

Emit a gVCF either way. The gVCF records reference confidence in non-variant blocks, which is the only way to distinguish “you are homozygous reference here” from “we had no coverage here.” That distinction decides whether a negative result means anything.

Expect 4.3–5.0 million variants total against GRCh38 for one human, of which ~3.5–4.1 million are SNVs and ~600–900k are indels. Roughly 10–12k of the SNVs will be missense in a protein-coding gene, and a few hundred will be nonsense, frameshift, or canonical splice-site.

4. Check the callset before you believe any single variant

This is the step people skip, and it is the step that catches sample swaps, contamination, and a mis-specified reference.

bcftools stats -F GRCh38_no_alt_analysis_set.fna sample.dv.vcf.gz > sample.stats
plot-vcfstats -p plots/ sample.stats

Four numbers:

• Ti/Tv ratio, all SNVs, genome-wide: expect 1.95–2.10. Below 1.9 means false positives are diluting your callset. Restricted to coding exons it should be near 3.0.
• Heterozygous/homozygous-alt ratio: 1.4–1.7 for most individuals, higher in African-ancestry genomes, lower with recent consanguinity. A ratio above 2.5 is the signature of contamination from a second sample.
• Total SNV count against the expected range above. Ten percent low means coverage or filtering problems. Twenty percent high means noise.
• Mean depth from samtools coverage, plus the fraction of the genome at ≥10x. You want ≥95% of the non-N reference at 10x or better.

Run a contamination estimate explicitly. VerifyBamID2 gives you an alpha, and anything above 0.02 makes low-allele-fraction calls untrustworthy:

VerifyBamID --SVDPrefix resource/1000g.phase3.100k.b38.vcf.gz.dat \
  --Reference GRCh38_no_alt_analysis_set.fna --BamFile sample.cram

If you also have an array genotype file from a consumer test, use it. Concordance between array and WGS at shared sites should exceed 99%. bcftools gtcheck will tell you in a minute whether the two files are the same person:

bcftools gtcheck -g array.vcf.gz sample.dv.vcf.gz

5. Filter to a callset you can defend

For DeepVariant, start from FILTER=PASS and then add depth and genotype quality thresholds. The defaults are calibrated, so do not invent hard filters on QD and FS the way you would for an unrecalibrated GATK callset.

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

The upper depth bound matters more than the lower one. Pileups above roughly 3x the genome mean are collapsed repeats and segmental duplications where reads from two paralogous loci stack on one reference copy, and they generate heterozygous-looking calls that are not variants at all.

Then intersect with a high-confidence region mask. The GIAB v4.2.1 “not difficult” BED, or the GA4GH stratification files, remove homopolymers, tandem repeats and segmental duplications where every caller degrades:

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

Keep both files. The masked one is what you trust for individual-variant claims. The unmasked one is what you use for aggregate work like polygenic scores, where losing 8% of sites systematically is worse than admitting a small false positive rate.

6. Annotate

Ensembl VEP with a local cache is the tool. Run it offline so nothing leaves your machine.

vep --offline --cache --dir_cache $HOME/.vep --assembly GRCh38 \
  --fasta GRCh38_no_alt_analysis_set.fna \
  --input_file sample.hiconf.vcf.gz --vcf --compress_output bgzip \
  --output_file sample.vep.vcf.gz \
  --everything --pick_allele_gene --canonical --mane \
  --custom gnomad.genomes.v4.1.sites.chr%CHROM%.vcf.bgz,gnomADg,vcf,exact,0,AF,AF_nfe,AF_afr,AF_eas,AF_sas,AF_amr,nhomalt \
  --custom clinvar_20250601.vcf.gz,ClinVar,vcf,exact,0,CLNSIG,CLNREVSTAT,CLNDN \
  --plugin CADD,whole_genome_SNVs.tsv.gz \
  --plugin SpliceAI,snv=spliceai_scores.raw.snv.hg38.vcf.gz \
  --fork 8

--pick_allele_gene collapses the transcript explosion to one consequence per gene per allele, which is what you want for a first pass. Use --mane so the chosen transcript matches what a clinical lab would report. CADD and SpliceAI are prediction scores. They are not evidence. Treat a CADD PHRED of 30 as “worth reading about” and nothing more.

ClinVar review status is the field that decides how much attention a variant deserves. CLNREVSTAT=criteria_provided,_multiple_submitters,_no_conflicts (two stars) or reviewed_by_expert_panel (three stars) is real evidence. A single-submitter assertion with no assertion criteria is close to noise, and ClinVar reclassifies thousands of variants a year, mostly toward benign.

7. Load it somewhere you can query

Flatten the annotated VCF into Parquet and query with DuckDB. This is faster and less painful than writing bcftools view one-liners forever.

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

CREATE TABLE v AS SELECT * FROM read_csv_auto('sample.flat.tsv.gz',
  delim='\t', header=false, names=['chrom','pos','ref','alt','csq','gene',
  'feature','hgvsc','hgvsp','af','nhomalt','cadd','clnsig','clnrev','gt','gq','dp']);

SELECT gene, hgvsp, gt, clnsig, clnrev, af FROM v
WHERE clnsig LIKE '%athogenic%'
  AND clnrev LIKE '%multiple_submitters%'
  AND (af IS NULL OR af < 0.01)
ORDER BY af NULLS FIRST;

That query returns, for most people, somewhere between zero and a handful of rows, and most of them will be recessive carrier states with one copy. Before you read anything into a row, pull up the region in IGV and look at the reads. A real heterozygous SNV at 30x shows roughly half the reads supporting each allele (binomially, 10–20 of 30), with the alternate allele appearing on both strands and at varied positions within reads. An artifact clusters at read ends, sits on one strand, or has an allele fraction near 0.2.

8. Put allele frequency in the right population context

gnomAD’s global AF is the wrong denominator for you specifically. Use the subpopulation AF matching your genetic ancestry, and if your ancestry is not well represented in gnomAD, treat rarity as weakly informative. Deep sequencing of 1,070 Japanese individuals identified millions of variants, a large fraction of them novel to that population and absent from the then-current global catalogues.³ Earlier gene-based SNP discovery in the same population had already found 190,562 variations concentrated in and around genes, a substantial share of which were population-specific.⁴ The practical consequence: “absent from gnomAD” means “rare in the populations gnomAD sampled well,” and a variant that looks singular in one ancestry can be a common polymorphism in another.

Estimate your own ancestry components rather than assuming them. Project onto 1000 Genomes principal components with plink2:

plink2 --vcf sample.hiconf.vcf.gz --extract 1kg_pruned.snplist \
  --score 1kg_pcs.eigenvec.allele 2 5 header-read no-mean-imputation \
  --score-col-nums 6-15 --out sample_pcs

9. Polygenic scores, with the caveats attached

If you want a polygenic score, use the unmasked PASS callset, take weights from the PGS Catalog, and run:

plink2 --vcf sample.pass.vcf.gz --score PGS000123.txt 1 4 6 header-read \
  cols=+scoresums --out sample_pgs

Two things to keep in mind. A raw score is meaningless without the reference distribution it was calibrated in, so you need the percentile, computed against a population panel that matches your ancestry. And score transferability across ancestries is poor: predictive accuracy typically drops by half or more when a European-derived score is applied to an African-ancestry genome. A percentile from the wrong reference is a number with no interpretation.

Common problems

Reference mismatch. Your CRAM references one FASTA, your caller is passed another with the same contigs but different order or different decoy content. Symptom: samtools errors about missing reference sequences, or silently wrong calls at chrEBV and decoy contigs. Fix by using one FASTA and one .fai everywhere, and verify with samtools view -H sample.cram | grep '@SQ' | md5sum against the reference dictionary.

chr prefix chaos. GRCh38 analysis set uses chr1. Some annotation files use 1. bcftools annotate --rename-chrs fixes it. Until you fix it, every --custom annotation returns empty. You then conclude, incorrectly, that you carry nothing interesting.

Coverage dropouts in GC-extreme regions. PCR-based library preps lose first exons of GC-rich genes. Check with mosdepth --by targets.bed and look at the fraction of your genes of interest below 10x. If a gene is at 4x, your homozygous-reference call there is not evidence of anything.

Treating a single-sample callset like a joint one. Depth-based and quality-based filters behave differently when you have one sample versus fifty. The recalibration approach that underpins GATK’s filtering was built around learning error models from data, which needs enough variants to model.² With one genome, prefer DeepVariant’s calibrated FILTER field over hand-tuned hard filters.

Mitochondrial and repeat regions. chrM needs a caller that handles heteroplasmy (Mutect2 in mitochondria mode). Standard diploid callers report garbage there. Similarly, do not trust calls in the HLA region, the KIR locus, or any tandem repeat expansion without a specialized tool.

Over-reading a single variant. A missense change with a CADD of 28 in a gene with a scary name is, most of the time, a benign polymorphism nobody has catalogued. The base rate of consequential findings in an unselected genome is low: large studies place actionable secondary findings at roughly 1–3% of individuals. If something survives all the filters and the IGV look and a ClinVar two-star pathogenic assertion, the next step is orthogonal confirmation in a CLIA/CAP lab and a conversation with a clinical geneticist. Research-grade sequencing is not a diagnosis and cannot be used as one.

Footnotes

1. The International SNP Map Working Group, Cold Spring Harbor Laboratories:, Ravi Sachidanandam, et al. A map of human genome sequence variation containing 1.42 million single nucleotide polymorphisms. Nature, 2001. https://doi.org/10.1038/35057149 ↩

2. Mark A. DePristo, Eric Banks, Ryan Poplin, et al. A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nature Genetics, 2011. https://doi.org/10.1038/ng.806 ↩ ↩²

3. Masao Nagasaki, Jun Yasuda, Fumiki Katsuoka, et al. Rare variant discovery by deep whole-genome sequencing of 1,070 Japanese individuals. Nature Communications, 2015. https://doi.org/10.1038/ncomms9018 ↩

4. Hiroaki Haga, Ryo Yamada, Yoshifumi Ohnishi, et al. Gene-based SNP discovery as part of the Japanese Millennium Genome Project: identification of 190 562 genetic variations in the human genome. Journal of Human Genetics, 2002. https://doi.org/10.1007/s100380200092 ↩