SNP Analysis in R on Your Own Whole-Genome Data

Oak
snp r wgs bioconductor vcf population-genetics
A dusk tidal plain covered in coral pillars, each topped by a glowing colored bead, clustering into drifting lanes across wet sand.

By the end of this you will have your own genome in a GDS file that loads in under a second, a table of quality metrics you can compare against published expectations, principal components placing you against the 1000 Genomes reference panel, an annotated table of your coding variants joined to ClinVar, and at least one polygenic score expressed as a percentile rather than a meaningless raw number. You need: a WGS VCF (30x short-read is the assumption throughout, ideally with the accompanying gVCF or a callable-regions BED), R 4.3+, Bioconductor 3.18+, bcftools 1.19+, plink2, about 40 GB of free disk, and 16 GB of RAM. Most of the work is not statistics. It is making sure the coordinates, alleles, and missingness mean what you think they mean.

A note on scope before the commands. Almost every technique below was designed for cohorts. An association test on one person is undefined. What you can do with n=1 is measurement: confirm your data is what the lab says it is, place your genotypes in a reference population, annotate what is known about specific variants, and compute scores whose weights came from someone else’s cohort. Anything that looks like a clinical finding, ClinVar pathogenic calls and pharmacogene star alleles especially, goes to a clinician or a certified lab for confirmation. Research-grade WGS is not a diagnostic test.

1. Normalize the VCF before R ever sees it

Do not load a raw caller VCF into R. Fix representation first, in bcftools, where it is fast and the semantics are well defined.

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

-m -any splits multiallelic records into one row per ALT. Skip it and every downstream join on (chrom, pos, ref, alt) will silently miss the second allele at a multiallelic site. -f with the exact reference FASTA your caller used left-aligns and normalizes indels; if the FASTA is wrong, bcftools will warn about REF mismatches, and those warnings are the single most useful signal that you have a build mismatch. -d exact drops duplicate records.

Check the contig naming while you are here:

bcftools view -h sample.norm.vcf.gz | grep '^##contig' | head -3

chr1 versus 1 is the cause of roughly half of all empty-result bugs in this pipeline. Pick one convention and enforce it on every file you touch, including the ClinVar VCF and the PGS Catalog scoring file. bcftools annotate --rename-chrs is the tool.

2. Convert to GDS and stop re-parsing text

VCF is a streaming text format. Every filter you run on it re-parses it. GDS (Genomic Data Structure) is a columnar on-disk format with array indexing, and SeqArray gives you a per-variant and per-sample query interface that reads only the blocks you ask for.

library(SeqArray)
seqVCF2GDS("sample.norm.vcf.gz", "sample.gds",
           storage.option = "LZMA_RA",
           parallel = 8)

f <- seqOpen("sample.gds")
seqSummary(f)

LZMA_RA gives compression with random access blocks. A 30x single-sample WGS VCF of roughly 4.5 million variants lands around 300-500 MB as GDS. Whole-genome scans that took minutes in text take seconds.

If you plan to use SNPRelate for PCA and kinship, produce the SNP-array-style GDS too, which holds biallelic genotype dosages only:

library(SNPRelate)
seqGDS2SNP(f, "sample.snp.gds")

Keep both. SeqArray for annotation and INFO/FORMAT fields, SNPRelate for the matrix algebra.

3. Compute quality metrics and compare them to expectations

This is the step people skip, and it is the one that catches a bad sample. Do it before any interpretation.

library(SeqArray)

# variant counts by type
seqSetFilter(f, variant.sel = seqGetData(f, "$num_allele") == 2)
ref <- seqGetData(f, "$ref"); alt <- seqGetData(f, "$alt")
snv <- nchar(ref) == 1 & nchar(alt) == 1
sum(snv)                      # expect ~3.8-4.2M for 30x GRCh38 WGS

ti <- (ref == "A" & alt == "G") | (ref == "G" & alt == "A") |
      (ref == "C" & alt == "T") | (ref == "T" & alt == "C")
sum(ti[snv]) / sum(!ti[snv])  # genome-wide Ti/Tv, expect ~2.0-2.1

gt <- seqGetData(f, "$dosage")[1, ]
het <- sum(gt == 1, na.rm = TRUE)
hom <- sum(gt == 0 | gt == 2, na.rm = TRUE)
het / hom

Genome-wide Ti/Tv below about 1.9 means false positives are leaking through your filters. In coding regions the expected ratio is higher, around 2.8-3.0, because of CpG deamination at methylated cytosines in exons. Heterozygous variant count tracks ancestry strongly: expect roughly 1.7-2.0 million het sites for a European-ancestry genome and more for African ancestry, since the reference is a poor match. A het count far below that range on a genome with normal coverage usually means an over-aggressive filter or contamination-driven allele-balance failures.

Also pull the depth distribution:

dp <- seqGetData(f, "annotation/format/DP")$data[1, ]
quantile(dp, c(0.05, 0.25, 0.5, 0.75, 0.95), na.rm = TRUE)
mean(dp < 10, na.rm = TRUE)   # fraction of called sites with thin support

Input DNA quantity and integrity drive call performance more than any downstream parameter, and degraded or low-input samples produce locus dropout that looks like homozygosity rather than missingness 1. If your median depth is under 25x or more than 5% of called sites sit below DP 10, treat everything downstream as provisional.

4. Confirm sex and check for contamination signals

seqResetFilter(f)
chrX <- seqGetData(f, "chromosome") %in% c("chrX", "X")
seqSetFilter(f, variant.sel = chrX)
pos <- seqGetData(f, "position")
# exclude pseudoautosomal regions (GRCh38)
par1 <- pos >= 10001    & pos <= 2781479
par2 <- pos >= 155701383 & pos <= 156030895
gx <- seqGetData(f, "$dosage")[1, !(par1 | par2)]
mean(gx == 1, na.rm = TRUE)   # het rate on non-PAR X

A non-PAR X heterozygosity rate near zero (under ~2%) is consistent with one X. Rates around 10-20% are consistent with two. Intermediate values, say 4-8%, are the interesting failure: that pattern shows up with sample contamination, and it is worth checking bcftools-computed allele balance at het sites genome-wide before believing anything else in the file. A clean sample has het allele fractions tightly centered on 0.5. A contaminated one has a heavy shoulder toward 0.3.

5. Ancestry PCA by projecting onto 1000 Genomes

One genome has no principal components. You need a reference panel, and you compute loadings on the panel, then project yourself.

Download the 1000 Genomes phase 3 GRCh38 calls, subset to a common, LD-pruned set, and merge:

plink2 --vcf 1kg.chr${i}.vcf.gz --set-all-var-ids '@:#:$r:$a' \
  --max-alleles 2 --snps-only --maf 0.05 --geno 0.01 \
  --make-pgen --out kg.chr${i}

Then in R:

library(SNPRelate)
kg <- snpgdsOpen("kg.merged.gds")

pruned <- snpgdsLDpruning(kg, ld.threshold = 0.2, maf = 0.05,
                          missing.rate = 0.01, slide.max.bp = 500000,
                          autosome.only = TRUE)
snpset <- unlist(pruned, use.names = FALSE)

pca <- snpgdsPCA(kg, snp.id = snpset, num.thread = 8, eigen.cnt = 20)
load <- snpgdsPCASNPLoading(pca, kg, num.thread = 8)

me <- snpgdsOpen("sample.snp.gds")
proj <- snpgdsPCASampleLoading(load, me, num.thread = 8)

Filtering thresholds of this shape, call rate above 0.99 and a minor allele frequency floor, are the standard conditioning step in SNP pipelines before any population-structure analysis, and the choice of threshold visibly moves the resulting structure 2. The LD pruning matters more than people expect. Without it, PC1 and PC2 can be dominated by a handful of long-range LD blocks (the MHC, the chromosome 8 and 17 inversions) rather than by ancestry. Exclude those regions explicitly if your PCs look strange.

Plot your projected coordinates against the panel colored by superpopulation. Two practical uses: it tells you which reference allele frequencies are appropriate when you interpret a variant, and it tells you whether a polygenic score derived in a European cohort has any calibration for you. Whole-genome data gives you far finer resolution here than a genotyping array does, and that resolution difference is the same one that makes WGS discriminate between samples that traditional genotyping calls identical 3.

6. Annotate coding variants and join to ClinVar

library(VariantAnnotation)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(BSgenome.Hsapiens.UCSC.hg38)

txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
seqlevelsStyle(txdb) <- "UCSC"

vcf <- readVcf("sample.norm.vcf.gz", "hg38")
seqlevelsStyle(vcf) <- "UCSC"

loc <- locateVariants(vcf, txdb, CodingVariants())
coding <- predictCoding(vcf, txdb, Hsapiens)
table(coding$CONSEQUENCE)

Expect roughly 20,000-25,000 coding SNVs, of which about 10,000-12,000 are nonsynonymous and 100-200 are nonsense or frameshift. Most of the loss-of-function calls in a healthy genome are annotation artifacts, variants in genes with poor transcript models, or real LoF in genes tolerant of it.

Reading the entire VCF with readVcf will use most of your RAM. Restrict first with a ScanVcfParam over a gene panel or a BED of exons, or iterate with VcfFile(yieldSize = 1e5).

For ClinVar, take the NCBI VCF and join on normalized coordinates:

clinvar <- readVcf("clinvar_GRCh38.vcf.gz", "hg38")
key <- function(v) paste(seqnames(rowRanges(v)), start(rowRanges(v)),
                         ref(v), unlist(alt(v)), sep = ":")
hits <- intersect(key(vcf), key(clinvar))

Normalize the ClinVar VCF with the same bcftools norm call and the same FASTA. Otherwise indel representation differs and you will match zero of them while matching SNVs fine, which is a confusing bug to chase.

Most non-coding variation will not appear in this table at all. Regulatory effects concentrate in distal elements whose activity varies by cell type, and the relevant annotation is epigenomic rather than transcript-based 4. Do not read an empty coding annotation as an empty genome.

7. Pull specific pharmacogene positions rather than scanning

Star-allele calling is not something to improvise in R. Use a dedicated caller (PharmCAT, Stargazer, Aldy) on the BAM, because the pharmacogenes have structural variation, gene conversion with pseudogenes, and phase-dependent haplotypes that a per-site VCF query cannot represent. What R is good for is pulling the individual defining positions so you can see the evidence.

seqResetFilter(f)
rs <- c("rs776746", "rs4244285", "rs1799853")
ids <- seqGetData(f, "annotation/id")
seqSetFilter(f, variant.sel = ids %in% rs)
data.frame(
  id  = seqGetData(f, "annotation/id"),
  pos = seqGetData(f, "position"),
  ref = seqGetData(f, "$ref"),
  alt = seqGetData(f, "$alt"),
  gt  = seqGetData(f, "$dosage")[1, ],
  dp  = seqGetData(f, "annotation/format/DP")$data[1, ]
)

rs776746 is the CYP3A5*3 splice variant, which creates a cryptic splice site in intron 3, inserts a premature stop, and is the reason most people of European ancestry express little functional CYP3A5 protein 5. Its allele frequency varies enormously across populations, which is exactly why step 5 comes before this one. If the annotation/id field is empty, your VCF was never annotated with dbSNP; run bcftools annotate -a dbsnp.vcf.gz -c ID first, or query by coordinate.

Everything in this table is measurement. Any inference about how you metabolize a drug belongs to a clinician working from a confirmed, clinical-grade genotype.

8. Compute a polygenic score as a percentile

Raw PRS values are unitless. The number only means something relative to a reference distribution computed with the same variant set and the same allele orientation.

# PGS Catalog harmonized scoring file, hmPOS_GRCh38
w <- read.delim("PGS000xxx_hmPOS_GRCh38.txt.gz", comment.char = "#")
w$key <- paste0("chr", w$hm_chr, ":", w$hm_pos, ":",
                w$other_allele, ":", w$effect_allele)

Then score yourself and every 1000 Genomes sample with the same file in plink2:

plink2 --pfile merged --score weights.txt 1 2 3 header cols=+scoresums \
  --out prs

Convert to a percentile within the reference samples closest to you on PC1-PC4. A score at the 90th percentile of a European reference and the 45th percentile of an East Asian reference is the same genotype, and both numbers are correct for their frame.

Two traps. First, palindromic SNPs (A/T and C/G) cannot be strand-resolved by allele identity, so either drop them or resolve by allele frequency, and drop them anyway if the frequency is near 0.5. Second, a variant absent from your VCF is treated as reference homozygous by --score, which is right for a well-covered site and wrong for an uncallable one. Intersect the weight file with your callable-regions BED and report how many weights survived. Losing 3% of weights is fine. Losing 30% means the score is not comparable to anything.

If you later have cohort data rather than one genome, SNPassoc is the shortest path to per-SNP association testing under codominant, dominant, recessive, overdominant, and log-additive models with the multiple-testing handling already in place 6.

Common problems

Missing means two different things. A site absent from a filtered VCF is either homozygous reference or uncalled, and the file does not distinguish them. This is the largest source of wrong answers in personal genomics. Keep the gVCF, or produce a callable-regions BED with GATK CallableLoci or mosdepth --quantize, and intersect every query against it.

Build mismatch. GRCh37 and GRCh38 coordinates differ by megabases in places. Symptoms: ClinVar joins return nothing, or return nonsense genes. Check ##contig lengths in the header. chr1 is 249,250,621 in GRCh37 and 248,956,422 in GRCh38.

Unnormalized indels. Two files can represent the same deletion three different ways. Always bcftools norm -f both sides of any join with the identical FASTA.

PLINK1 .bed drops phase and multiallelics. Fine for PCA and scoring. Not fine if you care about compound heterozygosity or haplotypes. Use pgen or stay in SeqArray.

Memory. readVcf on a whole genome will fail on 16 GB. Use ScanVcfParam(which = your_ranges) or yieldSize chunking, or do the filtering in SeqArray and only convert the surviving variants.

Ancestry-mismatched reference frequencies. gnomAD global allele frequency is not your allele frequency. Use the population-specific subset that matches your PCA projection, and note that for admixed genomes no single subset is right.

One sample, zero power. Nothing here is a hypothesis test. Treat the output as a description of one genome and a set of pointers into literature built on cohorts.

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. Steven A. Bates, Bruce Budowle, Morgan Johnson, et al. Large-scale analysis of DNA quantification metrics and SNP sequencing performance in unidentified human remains. Forensic Science International: Genetics, 2026. https://doi.org/10.1016/j.fsigen.2026.103470 ↩

  2. Bernd Gruber, Peter J. Unmack, Oliver F. Berry, et al. dartr: An r package to facilitate analysis of SNP data generated from reduced representation genome sequencing. Molecular Ecology Resources, 2018. https://doi.org/10.1111/1755-0998.12745 ↩

  3. Andreas Roetzer, Roland Diel, Thomas A. Kohl, et al. Whole Genome Sequencing versus Traditional Genotyping for Investigation of a Mycobacterium tuberculosis Outbreak: A Longitudinal Molecular Epidemiological Study. PLoS Medicine, 2013. https://doi.org/10.1371/journal.pmed.1001387 ↩

  4. Michael J. Ziller, Hongcang Gu, Fabian Müller, et al. Charting a dynamic DNA methylation landscape of the human genome. Nature, 2013. https://doi.org/10.1038/nature12433 ↩

  5. Peter Kuehl, Jiong Zhang, Yvonne Lin, et al. Sequence diversity in CYP3A promoters and characterization of the genetic basis of polymorphic CYP3A5 expression. Nature Genetics, 2001. https://doi.org/10.1038/86882 ↩

  6. Juan R. González, Lluís Armengol, Xavier Solé, et al. SNPassoc: an R package to perform whole genome association studies. Bioinformatics, 2007. https://doi.org/10.1093/bioinformatics/btm025 ↩