How to Interpret Your Own Whole-Genome Sequencing Data

By the end of this guide you will have taken a personal whole-genome sequencing dataset from raw reads to a small, ranked set of variants you can reason about, plus the quality metrics that tell you how much of the genome you can make claims about. Concretely: an aligned CRAM with coverage statistics, a joint-annotated VCF, a set of filtered call subsets (ClinVar pathogenic, loss-of-function in constrained genes, pharmacogenomic star alleles, polygenic scores), and a written record of what was and was not callable. You need the raw data, ideally FASTQ or CRAM rather than a genotyping array or a vendor PDF. You also need a Linux machine with 16 or more cores, 64 GB of RAM and about 1 TB of free disk for a 30x genome plus references. Finally, you need comfort with the command line and a scripting language. Budget a weekend for the first pass. The interpretation itself is not a weekend project, and nothing here is a clinical result: a variant that appears to matter medically needs confirmation in an accredited laboratory and a conversation with a clinical geneticist or genetic counselor before it informs any decision.

1. Get the right files, and check what you were given

Genomic interpretation starts with knowing which artifact you hold, because each one imposes a ceiling on what questions you can ask. A consumer genotyping array measures roughly 600,000 to 900,000 pre-selected positions, which is fine for ancestry and for polygenic scores built on imputed common variants, and useless for discovering a rare variant nobody put on the chip. Exome sequencing covers about 1 to 2 percent of the genome, the protein-coding portion, at high depth. Whole-genome sequencing at 30x mean depth covers coding and noncoding sequence, gives better copy-number and structural variant sensitivity, and produces more uniform coverage than exome capture because there is no hybridization bias ¹.

Ask your provider for three things: the FASTQ or CRAM files, the reference genome build used, and the pipeline version. Build matters enormously. GRCh38 is the right target for new analysis, and much of the published annotation and clinical literature was written against GRCh37, so you will be lifting coordinates back and forth. If you receive a VCF only, open it and read the header. The ##reference, ##source, and the ##GATKCommandLine or equivalent lines tell you what happened upstream.

bcftools view -h sample.vcf.gz | grep -E '^##(reference|source|GATKCommandLine|contig=<ID=chr1,)'
bcftools stats sample.vcf.gz > sample.stats.txt

A 30x human genome should yield roughly 4.5 to 5.5 million variants against GRCh38, of which about 4 to 4.5 million are single-nucleotide variants and 600,000 to 900,000 are small insertions and deletions. If your VCF has 3 million variants, something filtered aggressively. If it has 8 million, something is wrong with quality filtering or the sample is a mixture.

2. Align and call, or verify someone else’s calls

If you have FASTQ, we would align with BWA-MEM to the GRCh38 analysis set that includes ALT contigs and the decoy sequences, because the decoys absorb reads from sequence missing in the primary assembly and reduce false positives in paralogous regions.

bwa-mem2 mem -t 32 -R '@RG\tID:L1\tSM:SAMPLE\tPL:ILLUMINA\tLB:lib1' \
  GRCh38_full_analysis_set_plus_decoy_hla.fa R1.fastq.gz R2.fastq.gz \
  | samtools sort -@ 8 -m 2G -o sample.bam -
samtools index sample.bam
gatk MarkDuplicates -I sample.bam -O sample.md.bam -M dup_metrics.txt

For germline calling, GATK HaplotypeCaller in GVCF mode remains our default. It performs local de novo assembly of haplotypes in active regions rather than evaluating each position independently, which is what makes it handle indels and clustered variation reasonably well ². DeepVariant is the alternative we would run alongside it: a convolutional network over pileup images, generally better on indels and less dependent on hand-tuned filters. Running both and taking the intersection for high-confidence work and the union for discovery is a defensible strategy for a single sample.

gatk HaplotypeCaller -R GRCh38.fa -I sample.md.bam -O sample.g.vcf.gz -ERC GVCF
gatk GenotypeGVCFs -R GRCh38.fa -V sample.g.vcf.gz -O sample.raw.vcf.gz

Single-sample variant quality score recalibration is unreliable because the model needs many samples to learn the distribution of annotations. For one genome, use GATK’s hard filters or DeepVariant’s built-in quality score, and remember that machine-learning filters trained on populations behave differently on a single exome or genome than on a cohort ³.

3. Quantify what is callable before interpreting anything

This is the step most people skip, and it is the one that determines whether an absence of findings means anything. Compute per-base coverage and then the fraction of your genes of interest that reached at least 20x, which is roughly the depth where heterozygous single-nucleotide variant sensitivity plateaus.

mosdepth --by exons_GRCh38.bed --threshold 10,20,30 -t 8 sample sample.md.bam
samtools stats sample.md.bam | grep -E '^SN' > bamstats.txt

For a 30x genome expect mean autosomal depth of 28 to 33x, duplicate rate under 10 percent, and 90 to 95 percent of the reference callable at 20x or better. The uncallable remainder is not random. Segmental duplications, the pseudoautosomal regions, GC-extreme first exons, and several clinically important genes with close paralogs are systematically hard with short reads. Those genes include PMS2, SMN1, CYP2D6 and HBA1/HBA2. If a gene you care about sits in that set, short-read whole-genome sequencing gives you a weak negative. It does not give you a negative. Long-read sequencing or a gene-specific assay is the answer there.

Also confirm sample identity and contamination. Run VerifyBamID2 and expect a contamination estimate (FREEMIX) below 0.02. Cross-check reported sex against X and Y coverage and against X heterozygosity, and check the ethnicity-matched principal components if you have them, because ancestry mismatch between your sample and the allele-frequency panels you will use downstream distorts every rarity-based filter you apply.

4. Annotate with a versioned, reproducible stack

Annotation is where most interpretation errors originate, because annotation is a moving target and the same variant can be pathogenic in one release and reclassified in the next. Pin your versions and record them. We use Ensembl VEP with a local cache, plus a small set of plugins, and we write the version string into the output filename.

vep -i sample.filtered.vcf.gz -o sample.vep.vcf.gz --vcf --compress_output bgzip \
  --cache --offline --dir_cache $VEP_CACHE --assembly GRCh38 --fasta GRCh38.fa \
  --everything --pick_allele_gene --canonical --hgvs --symbol --numbers \
  --plugin CADD,whole_genome_SNVs.tsv.gz,gnomad.genomes.indel.tsv.gz \
  --plugin SpliceAI,snv=spliceai_scores.raw.snv.hg38.vcf.gz,indel=spliceai_scores.raw.indel.hg38.vcf.gz \
  --plugin LoFtool --plugin REVEL,revel_grch38.tsv.gz \
  --custom clinvar_20260601.vcf.gz,ClinVar,vcf,exact,0,CLNSIG,CLNREVSTAT,CLNDN \
  --custom gnomad.genomes.v4.1.sites.chr%CHROM%.vcf.bgz,gnomADg,vcf,exact,0,AF,AF_popmax,nhomalt \
  --fork 16 --stats_file sample.vep.html

Two flags deserve explanation. --pick_allele_gene collapses the many transcript consequences per variant to one per gene. That is what you want for triage and it loses information you may want later, so keep the full-transcript output too. --custom ... exact means ClinVar and gnomAD annotations attach only on exact allele match. This avoids the common mistake of inheriting a pathogenic label from a different allele at the same position.

MANE Select transcripts are the right default when a gene has several. Clinical reports are increasingly written against MANE, and using a non-MANE transcript is a frequent source of mismatched HGVS nomenclature when you try to look a variant up in the literature.

5. Triage: four queries that do most of the work

You now have roughly five million annotated variants and need to get to a few dozen worth reading about. The arithmetic is unforgiving: rare disease-causing variation is rare, and in any single genome the overwhelming majority of unusual-looking variants are benign ⁴. Set your priors accordingly and treat every candidate as false until it survives filtering.

The first query is the ClinVar pass, restricted to assertions with real review support.

bcftools +split-vep sample.vep.vcf.gz -f '%CHROM\t%POS\t%REF\t%ALT\t%SYMBOL\t%ClinVar_CLNSIG\t%ClinVar_CLNREVSTAT\t%gnomADg_AF\t[%GT]\n' \
  -d -A tab | awk -F'\t' '$6 ~ /Pathogenic|Likely_pathogenic/ && $7 ~ /multiple_submitters|reviewed_by_expert_panel|practice_guideline/'

Expect something on the order of 10 to 60 hits in a typical genome. Most of them are low-penetrance, recessive carrier states or annotations that do not apply to your genotype. Discard anything with gnomAD allele frequency above roughly 0.5 percent for a claimed dominant, highly penetrant condition, since a truly penetrant dominant allele cannot be that common.

The second query looks for predicted loss of function in genes that tolerate it poorly. It selects stop-gained, frameshift, canonical splice-site and start-loss variants in genes with a high gnomAD loss-of-function observed/expected upper bound fraction constraint. Filter on gnomAD allele frequency below 0.1 percent and check the read pileup for every survivor in IGV, because frameshift calls in homopolymers and short tandem repeats are the single most common class of artifact you will encounter.

The third query is pharmacogenomics, where the analysis is haplotype-based rather than variant-based. Run PharmCAT or PyPGx on your VCF rather than reading individual positions, because clinically meaningful alleles are defined as star-allele haplotypes and CYP2D6 additionally requires copy-number resolution that a plain VCF does not encode.

java -jar pharmcat-pipeline.jar -vcf sample.vcf.gz -o pgx_out --reporter-save-json

Read the resulting diplotypes and the phenotype calls as measurements of metabolizer status. Any change in medication belongs to the prescribing clinician and should be based on a validated clinical pharmacogenomic test. Your own pipeline output is not a basis for it.

The fourth query is polygenic. Scores for common, complex traits come from genome-wide association summary statistics and are computed as weighted sums over hundreds to millions of variants. Use the PGS Catalog with pgsc_calc, and interpret the result only against a reference population that matches your genetic ancestry, because score distributions shift substantially across populations and a raw score without a matched reference distribution is uninterpretable. The design and interpretation of association studies, including the statistical corrections that make these scores meaningful, is well covered in the sequencing-study literature ³.

6. Read the pileup and read the literature

For every candidate that survives triage, do two things. Open the region in IGV and look at the reads. You want allele balance near 0.5 for a heterozygote, mapping quality above zero and no strand bias. You also want no soft-clip cluster suggesting a nearby structural event, and the variant should not sit in a homopolymer run. Ten minutes of pileup inspection eliminates more false candidates than any filter tuning.

Then read. Look the variant up in ClinVar with its submission-level detail, in gnomAD for population frequency and homozygote count, and in the primary literature the ClinVar submitters cite. Pay attention to whether the original evidence was a linkage study in one large family, a case-control comparison, or a functional assay. The strength of the genotype-to-phenotype claim varies enormously across those designs and many early reports have not replicated ⁵. Population frequency alone resolves a surprising number of cases: a variant reported as causal for a severe childhood condition that appears in gnomAD with 200 homozygotes is not causal for a severe childhood condition.

Also keep in mind what variation in your genome is there because of selection history rather than disease. Alleles at loci under recent selection look unusual against a global frequency panel but are common and unremarkable in the populations where they arose ⁶. These include alleles affecting lactase persistence, skin pigmentation and high-altitude physiology.

7. Version, script, and store everything

Write the whole pipeline as a Nextflow or Snakemake workflow with containerized tool versions, and store the reference and annotation database checksums alongside the outputs. When you re-run in a year against a new ClinVar release and a variant’s classification has changed, you will want to know exactly what the old call was based on. Keep the CRAM (about 15 to 20 GB with reference-based compression, versus 45 to 60 GB of gzipped FASTQ) as your archival artifact, since you can always recall variants from a CRAM and cannot recover reads from a VCF. Interpretation improves faster than sequencing does, and the value of a personal genome grows as the annotation catches up to it ⁴.

Common problems

Build mismatch is the most frequent error and the hardest to notice. A GRCh37 VCF annotated against GRCh38 resources produces plausible-looking output with the wrong genes and the wrong alleles. Always confirm the contig naming convention (chr1 versus 1) and the build before annotating, and use CrossMap or Picard LiftoverVcf with an explicit reject file when you must convert, then check how many variants were rejected.

Pathogenic-by-position contamination happens when annotation joins on chromosome and position rather than the full allele. Use exact-match custom annotations and verify a handful of hits by hand against the ClinVar web record.

Absence of coverage read as absence of variant. If a gene averaged 6x, you had limited power to see a heterozygous variant in it. Report callable fraction alongside every negative finding, and treat the paralog-heavy genes named earlier as out of scope for short reads.

Structural and repeat variation missed entirely. Standard small-variant pipelines do not call copy-number changes, large deletions, or repeat expansions. Add Manta or GRIDSS for structural variants, a read-depth caller for copy number, and ExpansionHunter for known repeat loci. Accept that short-read sensitivity for these classes is meaningfully below its sensitivity for single-nucleotide variants ¹.

Overreading a single variant. Complex trait risk is distributed across many loci of small effect. The effect size attached to any one common variant in a personal genome is usually small enough that behavior, environment and family history dominate it ⁴. The reverse error also occurs in somatic contexts: tumor sequencing requires a matched normal sample and a somatic caller, and germline pipelines applied to tumor tissue produce nonsense ⁷.

Finally, a clinical finding from a research-grade pipeline is a hypothesis. Orthogonal confirmation in an accredited laboratory and interpretation with a clinical geneticist or genetic counselor come before any medical decision, and that boundary is worth holding firmly even when the evidence looks clear on your screen.

Footnotes

1. Rachel L Goldfeder, Dennis P Wall, Muin J Khoury, et al. Human Genome Sequencing at the Population Scale: A Primer on High-Throughput DNA Sequencing and Analysis. American Journal of Epidemiology, 2017. https://doi.org/10.1093/aje/kww224 ↩ ↩²

2. Aaron McKenna, Matthew Hanna, Eric Banks, et al. The Genome Analysis Toolkit: A MapReduce framework for analyzing next-generation DNA sequencing data. Genome Research, 2010. https://doi.org/10.1101/gr.107524.110 ↩

3. David B. Goldstein, Andrew Allen, Jonathan Keebler, et al. Sequencing studies in human genetics: design and interpretation. Nature Reviews Genetics, 2013. https://doi.org/10.1038/nrg3455 ↩ ↩²

4. Tuuli Lappalainen, Alexandra J. Scott, Margot Brandt, et al. Genomic Analysis in the Age of Human Genome Sequencing. Cell, 2019. https://doi.org/10.1016/j.cell.2019.02.032 ↩ ↩² ↩³

5. Hubert E Blum, Konrad Oexle. Clinical Interpretation and Implications of Whole Genome Sequencing. Journal of Next Generation Sequencing & Applications, 2014. https://doi.org/10.4172/2469-9853.1000105 ↩

6. Laura B. Scheinfeldt, Sarah A. Tishkoff. Recent human adaptation: genomic approaches, interpretation and insights. Nature Reviews Genetics, 2013. https://doi.org/10.1038/nrg3604 ↩

7. Hidewaki Nakagawa, Masashi Fujita. Whole genome sequencing analysis for cancer genomics and precision medicine. Cancer Science, 2018. https://doi.org/10.1111/cas.13505 ↩