How to Interpret Your Own DNA Data

By the end of this you will have a GRCh38-aligned CRAM and a quality-checked VCF with roughly 4.5–5 million variants. You will have an annotation layer from Ensembl VEP with ClinVar and gnomAD attached. You will also have a short list of variants worth a human conversation. That list runs to tens of variants and not to millions. You need your raw data, a Linux box with ~200 GB free and 16+ cores, and Docker or a conda environment. The raw data is ideally FASTQ or CRAM from a 30x whole-genome run, but this guide covers what to do if all you have is a lab VCF or an array file. Budget a weekend for the first pass and about 8–14 hours of wall clock for alignment plus calling on 16 cores. Nothing here is a diagnosis. Anything that looks clinically meaningful goes to a genetic counselor or physician for confirmation in a CLIA lab before you believe it.

1. Find out what you have

Three tiers of data get called “DNA results,” and they are not interchangeable.

An array file (23andMe, AncestryDNA) is a tab-separated text file of roughly 600,000 to 700,000 genotyped positions. That is about 0.02% of your genome, chosen because those sites are common in reference populations. It cannot see rare variants, which is where most of the high-effect biology lives. The positions are genotyped, not sequenced: the chip asks “which of these two alleles is here” and cannot report a third.

A lab VCF is a variant list relative to a reference. It tells you where you differ. It does not tell you where you were confidently checked and found identical. A position absent from a VCF might be reference, or might have had 2x coverage and been silently dropped.

FASTQ or CRAM is the real thing. Reads and base qualities, everything downstream is recomputable.

Check the reference build immediately, because build mismatch is the single most common way people generate wrong interpretations:

bcftools view -h sample.vcf.gz | grep -E '^##(reference|contig=<ID=(chr)?1,)'

If contigs are 1, 2, 3 with length 249250621, you are on GRCh37. If they are chr1, chr2 with chr1 at 248956422, you are on GRCh38. Coordinates differ by megabases in places. A GRCh37 position looked up in a GRCh38-indexed ClinVar will land on a different gene, and it will not error.

For CRAM you also need the exact reference FASTA the file was compressed against, since CRAM stores only differences from it:

samtools view -H sample.cram | grep '@SQ' | head -3
samtools quickcheck -v sample.cram

2. Understand the coordinate and strand conventions before you read anything

The 3’ and 5’ labels come from the numbered carbons of deoxyribose. The phosphodiester backbone links the 5’ carbon of one sugar to the 3’ carbon of the next, so every strand has a free 5’ phosphate at one end and a free 3’ hydroxyl at the other. DNA polymerase can only add nucleotides to a 3’ hydroxyl, so synthesis runs 5’ to 3’ and the template is therefore read 3’ to 5’. That asymmetry is why one strand is copied continuously and the other in Okazaki fragments. This chemistry is also why the first generation of sequencing methods read sequence in a defined polarity, and every convention downstream inherited it ¹.

What matters for your data: the reference genome FASTA, and therefore every coordinate in a BAM, CRAM, or VCF, is written in the 5’-to-3’ direction of the plus strand of the reference assembly. Genes on the minus strand are stored as their reverse complement. So a variant in a minus-strand gene like BRCA2… no, BRCA1 is the minus-strand one. A C>T in the VCF for a minus-strand gene appears as G>A in the coding sequence, and VEP handles that translation for you. Do not do it by hand.

This bites hardest with array data. Consumer chips have historically reported genotypes on the “dbSNP orientation,” “Illumina TOP/BOT,” or plus strand depending on the vendor and the year, and for A/T or C/G SNPs (palindromic sites) there is no way to infer the intended strand from the alleles alone. If you are converting a 23andMe file to VCF, use a tool that consults the reference base:

bcftools +fixref input.vcf.gz -Oz -o fixed.vcf.gz -- \
  -f GRCh38.fa -m flip -d
bcftools +fixref fixed.vcf.gz -- -f GRCh38.fa

The second command prints a summary. If more than a fraction of a percent remain non-reference-matching after flipping, stop and figure out the build rather than proceeding.

3. Align and call, or re-call

If you have FASTQ:

bwa-mem2 mem -t 16 -K 100000000 -Y \
  -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

Use a reference with decoy and HLA alt contigs. Reads that belong to sequence missing from the primary assembly will otherwise pile onto their nearest homolog and generate false heterozygous calls. The finished euchromatic reference covers about 99% of euchromatin at an error rate below one per 100,000 bases, but the remaining segmental duplications and satellite regions are exactly where spurious calls concentrate ².

For calling, we use DeepVariant. It is a single model, it needs no per-sample tuning, and on a 30x PCR-free genome it is more accurate on indels than the older haplotype-based callers we used before:

docker run -v "$PWD":/data google/deepvariant:1.6.1 \
  /opt/deepvariant/bin/run_deepvariant \
  --model_type=WGS \
  --ref=/data/GRCh38.fa \
  --reads=/data/sample.bam \
  --output_vcf=/data/sample.vcf.gz \
  --output_gvcf=/data/sample.g.vcf.gz \
  --num_shards=16

Ask for the gVCF. It records reference confidence in blocks, so you can later answer “was this position examined” rather than assuming absence means reference. That distinction is the difference between “you do not carry this variant” and “we did not look.”

If your lab gave you only a VCF and no reads, you can still annotate and filter. You cannot assess callability, and you should treat any negative result as unproven.

4. Quality-check before you interpret anything

Five numbers, in order. If any is off, fix it before annotating.

mosdepth -n --fast-mode -t 4 --by 1000 sample sample.bam
bcftools stats sample.vcf.gz > stats.txt
verifybamid2 --SVDPrefix resource/1000g.phase3.100k.b38.vcf.gz.dat \
  --Reference GRCh38.fa --BamFile sample.bam

Mean depth: 30x nominal should land between 28 and 35. Below ~25x, heterozygote sensitivity drops noticeably because you start getting 5-read positions where both alleles do not appear.

Coverage uniformity: look at the fraction of the genome at ≥10x and ≥20x. Above 90% at 20x is healthy for a PCR-free library. Library preparation, and fragmentation in particular, drives how evenly coverage lands across clinically relevant genes. A library with poor uniformity can leave specific exons under-covered even at a fine genome-wide mean ³. Check the genes you care about individually rather than trusting the average.

Ti/Tv ratio from bcftools stats: expect 2.0–2.1 genome-wide. Values near 1.5 mean false positives are diluting your call set.

Variant count: 4.4–5.0 million for a single genome against GRCh38, of which roughly 3.5–4 million are SNVs. Far fewer means undercalling. Far more means contamination or a build problem. For context on how deep this variation runs, the first genome-wide SNP map catalogued 1.42 million sites across a handful of individuals, and sequencing any one person now reproduces variation at roughly one site per 1,000 bases ⁴.

Contamination: verifyBamID2 FREEMIX above 0.02 means ~2% foreign DNA, which manufactures heterozygotes. Investigate rather than proceed.

Sex check: compare chrX and chrY mean depth against autosomes. A mismatch with the expected karyotype usually means a sample swap.

5. Annotate

Ensembl VEP with a local cache, offline, plus the two annotation sources that carry most of the interpretive weight:

vep -i sample.vcf.gz -o annotated.vcf.gz \
  --cache --offline --dir_cache $HOME/.vep \
  --assembly GRCh38 --fasta GRCh38.fa \
  --vcf --compress_output bgzip --fork 12 \
  --everything --pick_allele_gene \
  --custom file=clinvar_20250601.vcf.gz,short_name=ClinVar,\
format=vcf,type=exact,fields=CLNSIG%CLNREVSTAT%CLNDN \
  --custom file=gnomad.genomes.v4.1.sites.vcf.bgz,short_name=gnomADg,\
format=vcf,type=exact,fields=AF%AF_popmax \
  --plugin CADD,whole_genome_SNVs.tsv.gz \
  --plugin SpliceAI,snv=spliceai_snv.hg38.vcf.gz,indel=spliceai_indel.hg38.vcf.gz

--pick_allele_gene collapses the transcript explosion to one consequence per gene per allele. Without it a single variant in a gene with 15 transcripts yields 15 annotation records, and naive downstream counting inflates everything. If you want MANE Select transcripts specifically, add --mane_select.

gnomAD v4 aggregates roughly 807,000 individuals (about 730,000 exomes and 76,000 genomes). Population allele frequency is the strongest single filter you have. A variant claimed to cause a severe dominant disease that appears in 1% of any population is almost certainly misclassified.

6. Filter down to something a person can read

Do this in three passes, hardest filter first.

Pass one, known clinical assertions:

bcftools +split-vep annotated.vcf.gz -f '%CHROM\t%POS\t%REF\t%ALT\t%SYMBOL\t%ClinVar_CLNSIG\t%ClinVar_CLNREVSTAT\t%gnomADg_AF\t[%GT\t%DP\t%GQ]\n' \
  -i 'ClinVar_CLNSIG~"Pathogenic" && gnomADg_AF<0.005' \
  -d -A tab > clinvar_hits.tsv

Then filter that table by review status. Keep only criteria_provided,_multiple_submitters,_no_conflicts (two stars) or reviewed_by_expert_panel (three stars) or practice_guideline (four). One-star and zero-star assertions in ClinVar include a large number of entries that later reclassify, many toward benign. A typical genome produces a handful of two-star-or-better pathogenic hits, most of them recessive carrier status for conditions you will never manifest.

Pass two, predicted loss of function in genes with known haploinsufficiency: filter IMPACT=HIGH with gnomADg_AF < 0.001. High impact means stop gained, frameshift or canonical splice. You will get roughly 50–100 candidates. Most are artifacts in low-complexity regions or annotation errors on the last exon. Check each against the BAM in IGV before taking any seriously.

Pass three, splicing. SpliceAI delta scores above 0.5 on a variant that VEP called synonymous or intronic are the class of finding conventional filters miss entirely. Most of the genome is transcribed and regulatory elements are dense well outside coding exons, which is why noncoding interpretation remains the hard part ⁵. Be plain about the uncertainty. Outside of well-characterized splice and promoter positions, a noncoding variant with a high prediction score is a hypothesis. It is not a result.

Then look at every surviving call in IGV. You want read balance near 50/50 for a het and mapping quality 60. You want no strand bias, and the call should not sit in a homopolymer run. Roughly a third of what survives automated filtering will die at this step.

7. Pharmacogenes need a separate pipeline

CYP2D6, HLA-B, and the other star-allele genes involve copy number variation, hybrid genes, and long-range haplotypes that a position-by-position VCF cannot represent. Run PharmCAT for the callable genes and Aldy or Cyrius specifically for CYP2D6 from the BAM:

java -jar pharmcat-pipeline.jar sample.vcf.gz -o pgx_out/
python star_caller.py --manifest manifest.txt --genome 38 -o cyp2d6/

Read the output as a genotype, and nothing more. Diplotype-to-phenotype mapping and any consequence for medication belong to a clinician and a pharmacist working from CPIC guidance. They do not belong to you and a script.

8. Polygenic scores, with the caveat attached

plink2 --vcf sample.vcf.gz --score weights.txt 1 4 6 header \
  list-variants cols=+scoresums --out prs_out

Two things make a single-person PRS nearly uninterpretable in isolation. First, the raw sum is meaningless without a reference distribution, so you must score the same variants in a matched ancestry panel and report your percentile. Second, effect sizes are estimated mostly in European-ancestry cohorts and transfer poorly, with variance explained often dropping by half or more in other populations. Report percentiles, keep the score version pinned, and treat it as a population statistic applied to you rather than a property of you.

Common problems

Build mismatch that fails silently. Annotating a GRCh37 VCF against GRCh38 ClinVar produces plausible-looking output that is wrong. Always assert the build in code before annotating.

chr prefix inconsistency. bcftools annotate --rename-chrs fixes it. A mismatched prefix usually results in zero annotations rather than an error, which looks like “clean results.”

Treating absence as reference. A VCF only lists differences. Use the gVCF and check the reference-block GQ at any position where a negative result matters.

One-star ClinVar entries. Always read CLNREVSTAT. Single-submitter pathogenic assertions carry a materially higher reclassification rate than expert-panel ones.

Indels in homopolymers and tandem repeats. Represented inconsistently across callers. Always bcftools norm -m -any -f GRCh38.fa before comparing call sets, or the same variant will appear as two different ones.

Transcript choice. A variant can be missense on one transcript and intronic on another. Standardize on MANE Select and record the transcript ID and version next to every finding.

Paralogs. PMS2, SMN1, CYP21A2, and similar genes have near-identical pseudogenes that break short-read mapping. Standard WGS calls there are unreliable, and these genes need targeted long-read or MLPA assays.

Mitochondrial DNA. Heteroplasmy means chrM is not diploid. Use a caller that reports variant allele fraction rather than a genotype, and ignore any 0/1 call on chrM.

Acting on an unconfirmed research-grade result. Research pipelines are not CLIA. Anything that would change a medical decision needs orthogonal confirmation in an accredited lab and a conversation with a genetic counselor. Take the variant coordinates, the transcript, and the HGVS nomenclature to that appointment.

Footnotes

1. Walter Gilbert. DNA Sequencing and Gene Structure. Science, 1981. https://doi.org/10.1126/science.7313687 ↩

2. International Human Genome Sequencing Consortium. Finishing the euchromatic sequence of the human genome. Nature, 2004. https://doi.org/10.1038/nature03001 ↩

3. Vanessa Process, Madana M.R. Ambavaram, Sameer Vasantgadkar, et al. Optimization of DNA Fragmentation Techniques to Maximize Coverage Uniformity of Clinically Relevant Genes Using Whole Genome Sequencing. Diagnostics, 2025. https://doi.org/10.3390/diagnostics15182294 ↩

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

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