What a Genome Explorer Is, and How to Explore Your Own Genome

Oak
genomics wgs vcf igv variant-annotation
A tiny figure on a walkway shines a lamp on a canyon wall made of thousands of stacked glowing filaments.

A genome explorer is two things wearing one name. The first is a coordinate browser: you type chr7:117,559,590-117,559,620 and it draws the reference, the gene models, and your reads stacked underneath. The second is a query layer: you ask “which of my variants are in ClinVar with two or more review stars and a gnomAD allele frequency under 0.001?” and get a table. Consumer portals that call themselves Genome Explorer mostly do the first and hide the second behind a search box. If you can write code, you want both, running locally on your own files, because the interesting work is in the filters and a portal will not show you which filters it applied.

Below is the stack we would use on a personal 30x whole genome.

The files you should have before you start

Ask your sequencing provider for all of these. If you only get one, get the CRAM.

  • CRAM or BAM plus its index (.crai / .bai). A 30x human genome is roughly 90 Gbp of sequence. As CRAM with the reference on hand, that compresses to about 15-25 GB. The same data as BAM is 60-90 GB. CRAM is the right choice, with one catch described below.
  • VCF or gVCF plus .tbi. A single human genome called against GRCh38 typically yields 4-5 million variants, roughly 3.5-4.5M SNVs and 500-900k indels. If your VCF has 50,000 lines, you were given an exome or an array, not a genome.
  • The exact reference FASTA used for alignment, with its .fai and .dict. Not “GRCh38” as a word. The specific file.
  • QC reports: mosdepth or equivalent coverage summary, duplicate rate, contamination estimate.

Check the reference build immediately, before you build any tooling on top:

samtools view -H sample.cram | grep -m3 '^@SQ'
bcftools view -h sample.vcf.gz | grep -m1 '##reference'

GRCh37 and GRCh38 coordinates differ by hundreds of thousands to millions of bases in many regions. Annotating GRCh37 calls with a GRCh38 ClinVar file produces a table that looks fine and is wrong at nearly every row. The chr prefix is the other trap: chr1 versus 1. bcftools annotate --rename-chrs fixes it, and you should fix it once at the start rather than in five places later.

Normalize the VCF first, then annotate

Variant callers emit multiallelic sites and left-unaligned indels. Any join against ClinVar or gnomAD on position and alleles will silently miss those rows. Split and left-align once:

bcftools norm -m -any -f GRCh38.fa -c s sample.vcf.gz \
  | bcftools norm -d exact -Oz -o sample.norm.vcf.gz
bcftools index -t sample.norm.vcf.gz

Then annotate offline. We use Ensembl VEP with a local cache rather than a web API, because 4.5 million variants through a REST endpoint is a bad afternoon:

vep -i sample.norm.vcf.gz --cache --offline --assembly GRCh38 \
    --fasta GRCh38.fa --vcf --compress_output bgzip \
    --everything --fork 8 --buffer_size 20000 \
    --custom gnomad.genomes.v4.1.sites.vcf.gz,gnomADg,vcf,exact,0,AF,AF_popmax \
    -o sample.vep.vcf.gz

Add ClinVar separately, since NCBI ships it as a plain VCF and bcftools annotate is faster than a VEP custom field for this:

bcftools annotate -a clinvar.vcf.gz \
  -c INFO/CLNSIG,INFO/CLNREVSTAT,INFO/CLNDN \
  -Oz -o sample.annot.vcf.gz sample.vep.vcf.gz

Now the query layer exists. A rare-and-flagged filter looks like this:

bcftools view -i 'INFO/CLNSIG~"Pathogenic" && \
   (INFO/gnomADg_AF<0.001 || INFO/gnomADg_AF=".") && \
   INFO/CLNREVSTAT~"criteria_provided"' sample.annot.vcf.gz \
  | bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%INFO/CLNDN\t[%GT\t%DP]\n'

Expect on the order of tens of rows, most of them carrier status for recessive conditions, and expect several to be wrong. CLNREVSTAT matters more than CLNSIG: a single-submitter assertion with no assertion criteria is a claim, not a finding. Filter to two stars and above when you want signal.

Viewing reads: IGV, and when to use something else

Install IGV Desktop, not igv.org, if your CRAM is on your own disk. Load the genome from the same FASTA you checked above, then load the CRAM. IGV downsamples to 100 reads per window by default; turn that off (View > Preferences > Alignments) when you are inspecting a specific site, and turn it back on before you navigate to a 5 Mb view.

Three IGV settings that change what you see:

  • Enable “show soft-clipped bases.” Structural variant breakpoints and mismapped reads look like clean reference until you do.
  • Color alignments by insert size and pair orientation when checking a deletion or inversion call.
  • Turn on the mapping quality shading. A pile of MAPQ 0 reads in a segmental duplication is the most common reason a “variant” is not a variant.

The CRAM catch: CRAM stores reference-relative differences, so samtools must find the exact reference by MD5. If you get [E::cram_get_ref] Failed to populate reference, set REF_PATH / REF_CACHE or pass -T GRCh38.fa explicitly. This is the single most common hour lost when someone opens their own genome for the first time.

For sharing a view or overlaying public tracks, the UCSC Genome Browser remains the best option, and it is free. Do not upload a 20 GB BAM. Host the CRAM/VCF and index on any HTTP server with range requests enabled and add a custom track by URL. JBrowse 2 Desktop is the better choice if you want a local browser you can script and embed.

Verify the data is yours, and that it is one person

Before drawing conclusions, run three checks:

mosdepth -n --fast-mode -t 4 sample sample.cram    # mean depth, chrX/chrY ratio
bcftools stats sample.norm.vcf.gz | grep -E '^SN|^TSTV'
verifybamid2 --SVDPrefix 1000g.phase3.100k.b38 --Reference GRCh38.fa --BamFile sample.cram

Ti/Tv for a whole genome should land near 2.0-2.1. Substantially lower means false-positive SNVs are inflating your call set, and every downstream filter will leak. Contamination (FREEMIX) above about 0.02 makes low-allele-fraction calls untrustworthy. Chromosome X and Y depth relative to autosomes is a fast identity sanity check against what you expect.

What the browser cannot tell you

A genome explorer shows you evidence, not meaning. Pathogenicity assertions in ClinVar are made about populations of patients with phenotypes, and your reads carry no phenotype. Penetrance for most flagged variants is well below 100% and for many is unknown. Pharmacogenomic star alleles (CYP2D6 especially) cannot be called reliably from short-read data because of the neighboring pseudogene, no matter what a portal reports. If a variant in your data has any clinical bearing, it needs confirmation in a CLIA/CAP laboratory and interpretation by a clinician or genetic counselor who can see your history. We are describing measurement here, not medical advice.

Whole-genome data also carries information about relatedness and ancestry that array-based tests approximate and sequencing resolves directly, including rare variants that are informative about population history 1.

Questions people also ask

How many generations until you are no longer related? Genealogically, never: pedigree ancestry keeps branching. Genetically, you stop sharing detectable DNA with most of a given ancestor’s descendants after roughly seven to ten meioses. Expected shared autosomal DNA halves each generation, but recombination is chunky, so beyond about five generations a specific ancestor may contribute zero segments to you. For converting generations into years, the best current estimate of average human generation time over the past 250,000 years is about 26.9 years, with mothers averaging shorter and fathers longer 2.

What animal is 99% human DNA? Chimpanzees and bonobos. Single-nucleotide identity in alignable regions of the genome is about 98.8%. Counting insertions and deletions drops that to roughly 95-96%. The figure you see depends entirely on which alignment and which denominator someone used.

What is the best free software for DNA analysis? For reads and variants: samtools, bcftools, and mosdepth. For visualization: IGV Desktop and JBrowse 2, with UCSC for public tracks. For annotation: Ensembl VEP with a local cache, plus the dbNSFP plugin. For calling from FASTQ: GATK or DeepVariant. For population-scale genotype work: plink2 and Hail. All open source, all with the specific flags that a portal will not expose.

How much does a genome test typically cost? Consumer genotyping arrays, which read 600k-1M pre-chosen sites rather than sequencing anything, run $50-200. A 30x whole genome with raw file delivery generally runs several hundred to a few thousand dollars depending on turnaround and whether interpretation is included. The decline that made this possible came from massively parallel sequencing displacing Sanger chemistry over the last two decades 34, with per-base cost falling by orders of magnitude as read output per instrument run grew 5. The price of generating the data is no longer the limiting factor. Storing, annotating, and reasoning over it is 6.

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. 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 ↩

  2. Richard J. Wang, Samer I. Al-Saffar, Jeffrey Rogers, et al. Human generation times across the past 250,000 years. Science Advances, 2023. https://doi.org/10.1126/sciadv.abm7047 ↩

  3. Jay Shendure, Shankar Balasubramanian, George M. Church, et al. DNA sequencing at 40: past, present and future. Nature, 2017. https://doi.org/10.1038/nature24286 ↩

  4. Elaine R. Mardis. The impact of next-generation sequencing technology on genetics. Trends in Genetics, 2008. https://doi.org/10.1016/j.tig.2007.12.007 ↩

  5. Erik Pettersson, Joakim Lundeberg, Afshin Ahmadian. Generations of sequencing technologies. Genomics, 2009. https://doi.org/10.1016/j.ygeno.2008.10.003 ↩

  6. Sang Tae Park, Jayoung Kim. Trends in Next-Generation Sequencing and a New Era for Whole Genome Sequencing. International Neurourology Journal, 2016. https://doi.org/10.5213/inj.1632742.371 ↩