How to Build a Local SNP Viewer for Your Own Genome
By the end of this you will have four things on your own machine: a normalized, indexed VCF you can query by rsID or coordinate in milliseconds; an annotation layer that tells you what gene and consequence a variant falls in and how common it is in population references; a way to look at the actual sequencing reads under any genotype call, which is the only way to know whether a call is real; and a DuckDB table that lets you run SQL across your whole variant set without waiting. You need: a VCF (ideally from whole-genome sequencing, gzipped and around 1–2 GB for a single sample) and, if you have it, the aligned CRAM or BAM; the matching reference FASTA; about 100 GB of free disk; bcftools and samtools 1.19 or newer; tabix (ships with htslib); IGV 2.17 or igv.js; optionally Ensembl VEP; and Python 3.11 with duckdb. Everything below runs locally. None of it requires uploading your genome anywhere.
1. Figure out what file you have and which reference build it is on
Three file types get called “SNP data” and they are not interchangeable.
An array export (23andMe, AncestryDNA) is a tab-separated text file with roughly 600,000 to 900,000 rows: rsID, chromosome, position, genotype as two letters. It reports only sites the chip was designed to type. Everything else is absent, not reference.
A VCF from sequencing contains only positions where the caller found evidence of a non-reference allele, typically 3.5–5 million variants for one human genome. A gVCF additionally contains reference blocks, so absence of a record means “not covered” rather than “reference.”
A CRAM/BAM is the alignment itself. It is the ground truth for any single lookup.
Check the build before anything else. A GRCh37 coordinate interpreted as GRCh38 is off by tens of thousands of bases in most regions, and the variant you look up will be the wrong one.
bcftools view -h sample.vcf.gz | grep -E '^##(reference|contig=<ID=(chr)?1,)'If contig 1 has length 249250621 you are on GRCh37/hg19. If it is 248956422 you are on GRCh38. Note whether contigs are named 1 or chr1; that mismatch causes more failed joins than anything else in this pipeline.
2. Normalize and index the VCF
Do not skip this. Variant callers emit multi-allelic records and right-aligned indels, and both break naive position lookups.
bcftools norm \
-m -any \
-f GRCh38_full_analysis_set_plus_decoy_hla.fa \
--check-ref w \
-Oz -o sample.norm.vcf.gz \
sample.vcf.gz
bcftools index -t sample.norm.vcf.gz-m -any splits multi-allelic sites into one record per ALT allele. -f with --check-ref w left-aligns indels against the reference and warns (rather than fails) when the REF field in your VCF disagrees with the FASTA. If you see thousands of REF mismatch warnings, your FASTA and your VCF are on different builds. Stop and fix that first.
-t builds a .tbi index. For genomes with contigs longer than 512 Mb, or if you prefer CSI, use bcftools index -c. Sanity check:
bcftools index -n sample.norm.vcf.gz # record count
bcftools stats sample.norm.vcf.gz | grep '^SN'For a single human WGS sample expect roughly 4–5 million SNVs and 500,000–900,000 indels. Substantially fewer usually means the VCF was filtered to a panel or an exome. Substantially more usually means contamination or an unfiltered caller.
3. Look up one SNP
This is the core operation a “SNP viewer” performs. With an indexed VCF it is a region query:
# by coordinate (GRCh38)
bcftools view -H -r chr9:133255935 sample.norm.vcf.gz
# formatted, one line
bcftools query -r chr9:133255935 \
-f '%CHROM\t%POS\t%ID\t%REF\t%ALT\t%QUAL\t[%GT\t%DP\t%GQ\t%AD]\n' \
sample.norm.vcf.gzRead the output in this order: GT (0/1 heterozygous, 1/1 homozygous alternate, 0/0 reference, ./. no call), then DP (total depth at the site), then AD (depth per allele), then GQ (phred-scaled confidence in the genotype). A 30x WGS sample should show DP between roughly 20 and 45 at most autosomal positions. A heterozygous call with AD of 18,2 is not a heterozygote, it is a mapping artifact or a sequencing error that squeaked past the caller.
If the VCF only contains variant sites and your query returns nothing, that means either homozygous reference or no coverage. You cannot tell which from a variants-only VCF. That is exactly what step 6 is for.
4. Attach rsIDs and population frequencies
Many pipelines leave the ID column as .. Fill it from dbSNP, then add allele frequencies from gnomAD. Both must be on the same build and use the same contig naming as your VCF.
bcftools annotate \
-a GCF_000001405.40.gz \
-c ID \
-Oz -o sample.rsid.vcf.gz \
sample.norm.vcf.gz
bcftools index -t sample.rsid.vcf.gz
bcftools annotate \
-a gnomad.genomes.v4.1.sites.chrALL.trimmed.vcf.gz \
-c INFO/AF,INFO/AF_nfe,INFO/nhomalt \
-Oz -o sample.af.vcf.gz \
sample.rsid.vcf.gz
bcftools index -t sample.af.vcf.gzThe dbSNP files distributed by NCBI use RefSeq accessions (NC_000001.11) rather than chr1. Rename first:
bcftools annotate --rename-chrs refseq_to_ucsc.txt -Oz -o dbsnp.chr.vcf.gz dbsnp.vcf.gzNow you can go from rsID to genotype without a coordinate:
bcftools view -H -i 'ID="rs1801133"' sample.af.vcf.gzThat is a linear scan of the whole file, roughly a minute. Step 7 fixes it.
Frequency matters more than most people expect when reading their own variants. The first genome-wide SNP map catalogued 1.42 million SNPs at a density of about one variant per 1.9 kb of genome1, and modern references have extended that by orders of magnitude, mostly with rare variants. Deep whole-genome sequencing of 1,070 individuals from a single population found that the large majority of discovered variants were novel relative to existing databases and present at very low frequency2. Practically: a variant you carry that is absent from gnomAD is far more likely to be a calling artifact or a population-specific common variant than a meaningful finding, and the population your reference panel was built from changes the answer.
5. Add consequence annotation
Ensembl VEP is the one we use, because its consequence terms are the Sequence Ontology standard and it handles transcript selection explicitly.
vep \
-i sample.af.vcf.gz \
-o sample.vep.vcf.gz \
--vcf --compress_output bgzip \
--cache --offline --dir_cache "$HOME/.vep" \
--assembly GRCh38 \
--species homo_sapiens \
--everything \
--pick_allele_gene \
--fasta GRCh38_full_analysis_set_plus_decoy_hla.fa \
--fork 8 \
--buffer_size 20000--pick_allele_gene gives one consequence per allele per gene, which keeps the output size sane. Drop it if you want every transcript, and expect the file to triple. --everything turns on SIFT, PolyPhen, CANONICAL, MANE, and existing variation IDs. Runtime on 4.5 million variants with 8 forks is roughly 20–40 minutes.
snpEff is a reasonable alternative and is faster, but its gene model handling is less transparent when a variant sits in overlapping transcripts.
Note what the annotation will tell you: the vast majority of your variants fall outside protein-coding exons, and coding consequence terms will apply to maybe 20,000 of them. Non-coding does not mean inert. The ENCODE pilot project, annotating 1% of the genome, found pervasive transcription and large numbers of regulatory elements outside coding sequence3. It also means a variant in an intron or intergenic region generally has no interpretation you can read off a single annotation field.
6. Verify the genotype in the reads
For any variant you care about, open the alignment. This step catches the failure modes that no annotation pipeline reports.
Terminal view, fastest:
samtools tview -d T -p chr9:133255935 sample.cram GRCh38.faOr pull the pileup numerically:
samtools mpileup -f GRCh38.fa -r chr9:133255935-133255935 -q 20 -Q 20 sample.cram-q 20 drops reads with mapping quality below 20, -Q 20 drops bases with base quality below 20. If the alternate allele evidence disappears when you add those flags, the call was built on ambiguously mapped reads.
For visual inspection, IGV. Load the CRAM, the reference, and your annotated VCF as separate tracks, then in the IGV batch language:
new
genome hg38
load sample.cram
load sample.vep.vcf.gz
goto chr9:133255835-133256035
group BASE_AT_POS chr9:133255935
snapshot rs8176719.png
exitSave that as check.bat and run igv.sh -b check.bat. What you are looking for: allele balance near 50/50 for a heterozygote, reads supporting the alternate allele in both orientations, no cluster of soft-clipped reads starting at the same base, and no pile of mismatches in the surrounding 50 bp. A heterozygous call where every alternate read is on the forward strand is strand bias, and it is usually wrong.
For homozygous reference confirmation, the question is coverage. If DP at the position is 0, you have no information about that site, whatever your viewer shows.
7. Make lookups fast with DuckDB
Region queries against a tabix index are fast. Everything else (filter by consequence, sort by frequency, join against a gene list) is not. Flatten to Parquet once.
bcftools +split-vep sample.vep.vcf.gz \
-f '%CHROM\t%POS\t%ID\t%REF\t%ALT\t%AF\t%Consequence\t%SYMBOL\t%IMPACT\t%CANONICAL\t[%GT\t%DP\t%GQ]\n' \
-d -A tab \
| gzip > sample.flat.tsv.gzimport duckdb
con = duckdb.connect("genome.duckdb")
con.execute("""
CREATE OR REPLACE TABLE variants AS
SELECT * FROM read_csv('sample.flat.tsv.gz',
delim='\t', header=false, nullstr='.',
columns={'chrom':'VARCHAR','pos':'BIGINT','rsid':'VARCHAR',
'ref':'VARCHAR','alt':'VARCHAR','af':'DOUBLE',
'consequence':'VARCHAR','symbol':'VARCHAR','impact':'VARCHAR',
'canonical':'VARCHAR','gt':'VARCHAR','dp':'INTEGER','gq':'INTEGER'});
""")
con.execute("CREATE INDEX idx_rsid ON variants(rsid)")
con.execute("CREATE INDEX idx_pos ON variants(chrom, pos)")Now the queries you want are one-liners:
-- your genotype at a known rsID
SELECT * FROM variants WHERE rsid = 'rs1801133';
-- high-impact coding variants that are rare and well-supported
SELECT symbol, rsid, consequence, gt, dp, gq, af
FROM variants
WHERE impact = 'HIGH'
AND (af IS NULL OR af < 0.001)
AND dp BETWEEN 15 AND 60
AND gq >= 40
ORDER BY symbol;
-- everything in a gene
SELECT pos, rsid, ref, alt, consequence, gt, af
FROM variants WHERE symbol = 'APOE' ORDER BY pos;On a laptop this returns in well under a second across 4.5 million rows. The dp BETWEEN 15 AND 60 filter is doing real work: very high depth at a non-repetitive locus usually means collapsed repeats, and calls there are unreliable.
8. A browser view you can point at your own files
If you want a graphical viewer without installing desktop IGV, igv.js runs against local files served over HTTP with range requests.
python -m http.server 8080 --directory ~/genome<div id="igv"></div>
<script src="https://cdn.jsdelivr.net/npm/igv@2.15.11/dist/igv.min.js"></script>
<script>
igv.createBrowser(document.getElementById("igv"), {
genome: "hg38",
locus: "chr19:44,908,684",
tracks: [
{ type: "variant", format: "vcf",
url: "http://localhost:8080/sample.vep.vcf.gz",
indexURL: "http://localhost:8080/sample.vep.vcf.gz.tbi",
name: "My variants" },
{ type: "alignment", format: "cram",
url: "http://localhost:8080/sample.cram",
indexURL: "http://localhost:8080/sample.cram.crai",
name: "Reads" }
]
});
</script>CRAM in igv.js needs the reference sequence available, so either use a hosted hg38 genome ID as above or supply fastaURL and indexURL for your own FASTA. Serve over localhost only. Do not put an open bucket of your genome on the internet.
Common problems
chr prefix mismatch. bcftools annotate silently annotates nothing when contig names differ. Always check the record count after annotating: bcftools query -f '%ID\n' out.vcf.gz | grep -c '^rs'. If that number is near zero, this is why. Fix with --rename-chrs.
Build mismatch. Array files from consumer services are often GRCh37. gnomAD v4 and current dbSNP are GRCh38. Lift over with CrossMap or bcftools +liftover and a chain file, then re-normalize. Expect to lose 0.5–2% of sites, and treat lifted coordinates in segmental duplications with suspicion.
Variants-only VCF read as a full genotype table. Absence of a record is not evidence of reference. If you need confident reference calls, work from a gVCF or check depth in the CRAM directly.
Multi-allelic records. A site with REF A and ALT G,T and genotype 1/2 will not match a naive ref='A' AND alt='G' filter. Splitting in step 2 handles this, but check that your downstream tools kept the split.
Indel representation. The same deletion can be written several ways. Left-alignment against the reference (bcftools norm -f) makes representations comparable. Without it, joins against ClinVar or gnomAD miss real matches.
Array strand ambiguity. For A/T and C/G SNPs on genotyping arrays, forward-strand and reverse-strand reports are indistinguishable without a reference panel, so a TOP/BOT convention mismatch silently flips those genotypes. Sequencing data does not have this problem.
Imputed sites treated as measured. If your array data went through an imputation server, most of the rows in the output were inferred from linkage with neighboring typed markers, not measured. Imputation quality drops sharply for rare alleles and for ancestries underrepresented in the reference panel. Check the INFO/R2 field and treat anything below about 0.8 as a guess.
Low depth or homopolymers. Calls in homopolymer runs and short tandem repeats are the least reliable part of any short-read genome. Look at the reads before you believe them.
Stale annotation. ClinVar assertions change. Record the version you annotated with and re-run periodically rather than trusting a cached call from a year ago.
One thing this pipeline does not do: tell you what a variant means for your health. Population frequency, consequence term, and read support tell you whether a call is real and how common it is. They do not tell you about penetrance, and most published variant-trait associations come from cohorts whose ancestry and environment may not match yours. If a variant looks clinically relevant, the next step is a clinician and a CLIA-certified confirmatory test, not a deeper query.
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
-
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 ↩
-
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 ↩
-
Ewan Birney, Paul Flicek, Damian Keefe, et al. Identification and analysis of functional elements in 1% of the human genome by the ENCODE pilot project. Nature, 2007. https://doi.org/10.1038/nature05874 ↩