Whole Genome SNP Analysis: From FASTQ to an Annotated, Queryable Variant Table

By the end of this you will have a normalized, filtered, annotated VCF of roughly 3.5 to 4.5 million single-nucleotide variants from one person’s genome. You will also have a QC report that tells you whether to trust it, and a columnar table you can query with SQL or hand to an agent. You need: paired-end FASTQ or a CRAM from a 30x PCR-free Illumina run (or ONT/PacBio, noted where it differs), a Linux box with 16+ cores and 64 GB RAM, about 500 GB of free disk for intermediates, and Docker or Singularity. Budget six to twelve hours of wall time for a 30x genome on 32 cores. Everything below assumes GRCh38.

A note on scale before you start. An individual genome differs from the reference at roughly one site per thousand bases, which is where the 3 to 4 million SNV figure comes from. High-coverage sequencing of single individuals has landed in that range consistently across populations and platforms ¹². The 1000 Genomes pilot catalogued about 15 million SNPs across 179 individuals and found that each person carries roughly 250 to 300 loss-of-function variants in annotated genes ³. So “how common are SNPs” has two answers: about 1 in 1000 positions differ in you specifically, and across humanity the great majority of the genome is polymorphic somewhere.

1. Verify the data you were given before touching it

Checksums first, then read structure.

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

You want read length 150 (or 151), read count in the 700M–950M pair range for 30x, and mean quality above Q30. If seqkit reports a read length distribution rather than a fixed length, the vendor already trimmed adapters, which changes nothing downstream but tells you the FASTQ is not raw.

Then check for adapter contamination and overrepresented sequences:

fastp -i sample_R1.fastq.gz -I sample_R2.fastq.gz \
  --json fastp.json --html fastp.html --thread 8 \
  --disable_adapter_trimming --disable_quality_filtering \
  --disable_length_filtering -o /dev/null -O /dev/null

We run fastp in report-only mode. Aggressive pre-trimming before a soft-clipping aligner costs more than it gains. If duplication rate in the fastp report is above ~10% on a PCR-free library, the library was over-amplified or under-loaded, and your effective coverage is lower than the nominal number.

If you were handed a CRAM instead, skip to step 4 but confirm which reference it was compressed against, because CRAM is reference-dependent and will silently fail or corrupt bases with the wrong FASTA:

samtools view -H sample.cram | grep -E '^@SQ|^@PG|^@RG' | head -40

2. Pick the reference build and stick to it

Use the GRCh38 no-alt analysis set with decoys (GCA_000001405.15_GRCh38_no_alt_analysis_set.fna). Alt contigs help only if your aligner and caller are alt-aware end to end, and most downstream annotation resources assume primary chromosomes. Decoy sequences matter: without them, reads from unplaced repeats pile onto homologous primary loci and manufacture false heterozygous SNVs.

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
samtools faidx GCA_000001405.15_GRCh38_no_alt_analysis_set.fna
bwa-mem2 index GCA_000001405.15_GRCh38_no_alt_analysis_set.fna

T2T-CHM13v2.0 is a better reference biologically and resolves regions GRCh38 cannot, but ClinVar, gnomAD, dbSNP, and every GWAS summary statistic you will want to join against are keyed to GRCh38. Align to GRCh38 for the analysis you can act on. Run CHM13 as a second pass if you care about segmental duplications, acrocentric arms, or the extra ~200 Mb of sequence.

3. Align

bwa-mem2 mem -t 32 -K 100000000 -Y \
  -R '@RG\tID:L1\tSM:SAMPLE1\tPL:ILLUMINA\tLB:lib1\tPU:HXXXX.1' \
  GCA_000001405.15_GRCh38_no_alt_analysis_set.fna \
  sample_R1.fastq.gz sample_R2.fastq.gz \
| samtools fixmate -m -u -@ 4 - - \
| samtools sort -@ 8 -m 3G -T /scratch/sort -u -o - \
| samtools markdup -@ 8 --reference GCA_000001405.15_GRCh38_no_alt_analysis_set.fna \
  -f markdup.stats - SAMPLE1.cram
samtools index SAMPLE1.cram

-K 100000000 fixes the chunk size so output is deterministic regardless of thread count, which you want if you ever need to reproduce a result. -Y uses soft clipping for supplementary alignments, which structural variant callers need later. Write CRAM, not BAM: roughly 15–20 GB instead of 60–80 GB, with no information loss when you keep the original quality strings.

4. QC the alignment before calling anything

This is the step people skip and then spend a week chasing phantom variants.

mosdepth -n --fast-mode --by 1000 -t 4 SAMPLE1 SAMPLE1.cram
samtools stats -@ 8 --reference $REF SAMPLE1.cram > SAMPLE1.stats
verifyBamID2 --SVDPrefix resource/1000g.phase3.100k.b38.vcf.gz.dat \
  --Reference $REF --BamFile SAMPLE1.cram --NumThread 8

Thresholds we use as pass/fail:

• Mean autosomal coverage ≥ 28x, and ≥ 90% of the callable genome at ≥ 20x. Below 20x, heterozygote sensitivity drops fast because you need enough alt reads to distinguish a het from a sequencing error, and the allele-balance tail is where false negatives hide ⁴.
• FREEMIX from verifyBamID2 below 0.02. Above 0.03 you have sample contamination and your het count will be inflated by hundreds of thousands of calls.
• Duplicate rate below 10% (PCR-free), insert size mode 300–450 bp, ≥ 98% properly paired.
• Chromosome X and Y coverage ratio consistent with expected karyotype. mosdepth’s summary file gives you this directly. A mismatch here usually means a sample swap.

Long-read data changes the arithmetic. Single-molecule platforms trade per-base accuracy for read length and were the first to show that a usable individual genome could be assembled from single molecules without amplification ⁵. If you have ONT R10.4.1 or PacBio HiFi, 20–25x is enough for SNVs and buys you phasing and structural variants that short reads miss.

5. Call variants with DeepVariant

singularity run -B /data:/data docker://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/SAMPLE1.cram \
  --output_vcf=/data/SAMPLE1.dv.vcf.gz \
  --output_gvcf=/data/SAMPLE1.dv.g.vcf.gz \
  --num_shards=32 \
  --intermediate_results_dir=/scratch/dv_tmp

We recommend DeepVariant over GATK HaplotypeCaller for a single sample. HaplotypeCaller’s VQSR needs a cohort or falls back to hard filters that are tuned for someone else’s data, and the older statistical callers relied on explicit error models and base quality recalibration that a CNN learns implicitly from the pileup ⁴. DeepVariant needs no recalibration step and no truth-set-derived filter model at call time. Use --model_type=PACBIO or ONT_R104 for long reads.

Cost: about 5–8 hours on 32 cores for 30x WGS. If you want it in one hour, GPU builds exist, or use --make_examples sharded across machines.

Always produce the gVCF. It records confident reference calls with their GQ, which is the only way to distinguish “you are reference here” from “we could not call this position.” That distinction matters when you go looking for a specific variant and do not find it.

6. Normalize and filter

Raw caller output is not comparable to anything until it is left-aligned and decomposed.

bcftools norm -f $REF -m -any -Ou SAMPLE1.dv.vcf.gz \
| bcftools norm -f $REF -d exact -Ou \
| bcftools view -f PASS -Oz -o SAMPLE1.norm.vcf.gz
bcftools index -t SAMPLE1.norm.vcf.gz

-m -any splits multiallelic records into one row per ALT, which is what every annotation and join step expects. -d exact drops the duplicate rows that splitting can create. Skipping normalization is the single most common reason a variant “isn’t in gnomAD” when it plainly is.

Then mask the regions where short reads lie. Intersect with the GIAB high-confidence regions or, at minimum, exclude the ENCODE blacklist and low-mappability regions:

bcftools view -T ^GRCh38_unified_blacklist.bed -Oz -o SAMPLE1.filt.vcf.gz SAMPLE1.norm.vcf.gz

For downstream work that needs high precision, add FORMAT/GQ >= 20 && FORMAT/DP >= 10 as a filter. For discovery work, keep everything and carry GQ as a column.

7. Sanity-check the callset with ratios

bcftools stats -F $REF -s - SAMPLE1.filt.vcf.gz > SAMPLE1.vcfstats
plot-vcfstats -p plots/ SAMPLE1.vcfstats

What the numbers should look like for one 30x genome on GRCh38:

• Total SNVs: 3.5M–4.5M. European-ancestry genomes sit lower in that range, African-ancestry genomes higher, because more of their variation is absent from the reference ⁶.
• Indels: 500k–900k.
• Ti/Tv genome-wide: 2.0–2.1. If you see 1.7 or below, you have false positives, usually from contamination or an unmasked repeat region. Restricted to coding regions the expectation rises to about 3.0.
• Het/hom-alt ratio: 1.3–2.0. Much higher means contamination or duplicate-region pileups. Much lower means consanguinity, long runs of homozygosity, or a coverage problem.
• Fraction of SNVs present in dbSNP: above 95% for a well-covered genome of an ancestry represented in the databases. The novel remainder is a mix of genuinely rare variants and errors, and it was the same fraction that made early individual-genome papers hard to interpret ¹².

If you have a GIAB sample (HG001–HG007) sequenced alongside as a control, run hap.py against the v4.2.1 truth set. Expect SNV F1 above 0.995 for a healthy 30x pipeline. Without a control, the ratios above are your only evidence.

8. Annotate

vep --offline --cache --dir_cache /data/vep --assembly GRCh38 \
  --fasta $REF --species homo_sapiens \
  --input_file SAMPLE1.filt.vcf.gz --format vcf \
  --vcf --compress_output bgzip --output_file SAMPLE1.vep.vcf.gz \
  --fork 16 --buffer_size 50000 \
  --everything --pick_allele_gene \
  --custom file=/data/gnomad.genomes.v4.1.sites.vcf.bgz,short_name=gnomADg,\
format=vcf,type=exact,fields=AF%AF_nfe%AF_afr%AF_eas%nhomalt \
  --custom file=/data/clinvar_20250608.vcf.gz,short_name=ClinVar,\
format=vcf,type=exact,fields=CLNSIG%CLNREVSTAT%CLNDN \
  --plugin dbNSFP,/data/dbNSFP5.0a_grch38.gz,REVEL_score,CADD_phred,AlphaMissense_score,\
SpliceAI_pred_DS_AG,SpliceAI_pred_DS_AL

--pick_allele_gene collapses VEP’s per-transcript output to one consequence per allele per gene, which keeps the table manageable. Use --flag_pick_allele_gene instead if you want all transcripts retained with the canonical one flagged, which is the better choice if you plan to interrogate splice regions.

Run the whole thing offline with a local cache. Remote VEP over the API will take days for 4 million variants.

9. Make it queryable

Flatten to a columnar file and stop treating a VCF as a database.

bcftools +split-vep SAMPLE1.vep.vcf.gz -d -f \
'%CHROM\t%POS\t%REF\t%ALT\t%QUAL\t[%GT\t%DP\t%GQ\t]%SYMBOL\t%Consequence\t%IMPACT\t\
%gnomADg_AF\t%gnomADg_AF_nfe\t%ClinVar_CLNSIG\t%CADD_phred\t%REVEL_score\n' \
-o SAMPLE1.tsv

import duckdb
con = duckdb.connect("SAMPLE1.duckdb")
con.execute("""
CREATE TABLE variants AS
SELECT * FROM read_csv('SAMPLE1.tsv', delim='\t', header=false,
  columns={'chrom':'VARCHAR','pos':'BIGINT','ref':'VARCHAR','alt':'VARCHAR',
           'qual':'DOUBLE','gt':'VARCHAR','dp':'INTEGER','gq':'INTEGER',
           'symbol':'VARCHAR','csq':'VARCHAR','impact':'VARCHAR',
           'af':'DOUBLE','af_nfe':'DOUBLE','clnsig':'VARCHAR',
           'cadd':'DOUBLE','revel':'DOUBLE'});
""")

The three queries worth running first:

-- Rare, high-impact, well-covered
SELECT * FROM variants
WHERE impact='HIGH' AND (af IS NULL OR af < 0.001)
  AND gq >= 30 AND dp BETWEEN 15 AND 80;

-- Anything ClinVar calls pathogenic with multi-submitter review status
SELECT * FROM variants WHERE clnsig LIKE '%athogenic%';

-- Missense with strong in-silico support
SELECT * FROM variants
WHERE csq LIKE '%missense%' AND revel > 0.75 AND (af IS NULL OR af < 0.01);

The first query typically returns 100–300 rows for one genome, most of which are annotation artifacts in the last exon or in genes with no established phenotype. The second usually returns a handful of real carrier findings for recessive conditions. Nothing in this table is a diagnosis. A variant that looks pathogenic here needs orthogonal confirmation in a CLIA-certified laboratory and interpretation by a clinical geneticist or genetic counselor before it means anything about your health.

10. Ancestry, relatedness, and polygenic scores

Convert to PLINK2 format and project onto a reference panel:

plink2 --vcf SAMPLE1.filt.vcf.gz --max-alleles 2 --snps-only just-acgt \
  --set-all-var-ids '@:#:$r:$a' --new-id-max-allele-len 60 \
  --make-pgen --out SAMPLE1
plink2 --pfile SAMPLE1 --score 1000g_pca_loadings.tsv 2 5 header-read \
  variance-standardize cols=+scoresums --score-col-nums 6-15 --out SAMPLE1.pca

Whole-genome data makes this qualitatively better than array data. Arrays genotype tag SNPs and rely on linkage disequilibrium to represent the variants they do not measure, an approach whose limits were mapped out early in the GWAS era ⁷. Sequencing observes rare variants directly, which is what makes fine-scale demographic inference from a single genome possible at all ⁶. The historical arc is worth knowing: the first genome-wide SNP maps came from reduced-representation shotgun sequencing of pooled DNA and catalogued about 1.4 million sites ⁸, and the cost trajectory that made per-person sequencing routine was being forecast two decades ago ⁹.

For polygenic scores, use plink2 --score with published weight files, and read the ancestry caveats in the source paper before believing the output. A score derived in European cohorts applied to a non-European genome is not interpretable.

Common problems

Ti/Tv of 1.6 and 6 million SNVs. Contamination or a reference mismatch. Check FREEMIX first. If contamination is clean, confirm the FASTA used for alignment matches the one used for calling, byte for byte.

Variants missing from a region you care about. Check coverage there before concluding anything. mosdepth --by regions.bed takes seconds. Genes with paralogs are systematically miscalled by short reads and need targeted long-read or MLPA work. Examples are PMS2, SMN1/2, CYP2D6, and HBA1/2.

A variant your consumer array reported that your WGS does not show. The array probe may sit on a different strand or reference allele. Normalize both to GRCh38 with bcftools norm and compare by chrom:pos:ref:alt, not rsID. rsIDs merge, split, and get retired.

gVCF is enormous. That is correct behavior. It has a record for every base. Compress with bcftools view -Ob and keep it, because it is the only record of what you confidently did not have.

VEP dies with out-of-memory. Lower --buffer_size to 5000 and --fork to 4. The dbNSFP plugin holds large tabix blocks in memory per fork.

Everything looks fine but you cannot reproduce last month’s run. Pin container digests, not tags. google/deepvariant:1.6.1 can be repointed. google/deepvariant@sha256:... cannot.

Footnotes

1. Akihiro Fujimoto, Hidewaki Nakagawa, Naoya Hosono, et al. Whole-genome sequencing and comprehensive variant analysis of a Japanese individual using massively parallel sequencing. Nature Genetics, 2010. https://doi.org/10.1038/ng.691 ↩ ↩²

2. Hui Shen, Jian Li, Jigang Zhang, et al. Comprehensive Characterization of Human Genome Variation by High Coverage Whole-Genome Sequencing of Forty Four Caucasians. PLoS ONE, 2013. https://doi.org/10.1371/journal.pone.0059494 ↩ ↩²

3. The 1000 Genomes Project Consortium. A map of human genome variation from population-scale sequencing. Nature, 2010. https://doi.org/10.1038/nature09534 ↩

4. Ruiqiang Li, Yingrui Li, Xiaodong Fang, et al. SNP detection for massively parallel whole-genome resequencing. Genome Research, 2009. https://doi.org/10.1101/gr.088013.108 ↩ ↩²

5. Dmitry Pushkarev, Norma F Neff, Stephen R Quake. Single-molecule sequencing of an individual human genome. Nature Biotechnology, 2009. https://doi.org/10.1038/nbt.1561 ↩

6. Krishna R. Veeramah, Michael F. Hammer. The impact of whole-genome sequencing on the reconstruction of human population history. Nature Reviews Genetics, 2014. https://doi.org/10.1038/nrg3625 ↩ ↩²

7. Christopher S. Carlson, Michael A. Eberle, Mark J. Rieder, et al. Additional SNPs and linkage-disequilibrium analyses are necessary for whole-genome association studies in humans. Nature Genetics, 2003. https://doi.org/10.1038/ng1128 ↩

8. David Altshuler, Victor J. Pollara, Chris R. Cowles, et al. An SNP map of the human genome generated by reduced representation shotgun sequencing. Nature, 2000. https://doi.org/10.1038/35035083 ↩

9. Simon T Bennett, Colin Barnes, Anthony Cox, et al. Toward the $1000 Human Genome. Pharmacogenomics, 2005. https://doi.org/10.1517/14622416.6.4.373 ↩