How to Read Your DNA Results, From Raw Files to a Single Variant

Oak
genomics vcf wgs bioinformatics variant-interpretation
A lab machine scrolls a ribbon of bead-marked bands while a robotic arm holds one lit band apart from the dimmed rest.

By the end of this piece you will be able to open your own DNA results and say what every column in a line means. You will also be able to decide whether that line is trustworthy and turn a genome-wide file into a short list of variants you can read one at a time.

You will need the raw data. That means a .vcf.gz plus its .tbi index, ideally with a .cram or .bam alongside it. A genotyping array text file is the minimum. Set aside about 50 GB of disk and work in a Unix shell. Install bcftools ≥1.18, samtools, and mosdepth. You also need either Ensembl VEP or snpEff. Everything below runs comfortably on a laptop, with two exceptions: VEP’s cache download is around 25 GB, and any step that reads the alignment file is heavier than the rest.

1. Identify which file you have

Before interpreting anything, work out what kind of file is in front of you. There are four things people call “DNA results,” and they carry wildly different amounts of information.

An array file, the kind 23andMe and AncestryDNA export, is a tab-separated text file listing roughly 600,000 to 700,000 genotyped positions:

# rsid  chromosome  position  genotype
rs4477212   1   82154   AA
rs3094315   1   752566  AG
rs4988235   2   136608646   AG

That is about 0.02% of your 3.1 billion bases, and those positions were chosen because they are common polymorphisms. Single-nucleotide polymorphisms are the dominant class of human variation, and they were mapped densely enough two decades ago to make this chip design viable, at roughly one SNP per 1.9 kb in the first genome-wide map 1. The practical consequence is simple: an array tells you the genotype at the positions on the chip and nothing whatsoever about the positions that are not.

The sequencing-based files are larger and more informative. A VCF (variant call format) from whole-genome sequencing has 4 to 5 million lines for one person, each line describing a position where the sample differs from the reference genome. A gVCF extends that idea by carrying a record for every base, including the ones that matched the reference. That is what lets you distinguish “reference” from “no data.” A CRAM or BAM is the alignment itself, about 15 to 60 GB, and it is the only file that lets you look at the actual sequencing reads.

A few commands will tell you which of these you are holding:

bcftools view -h sample.vcf.gz | grep -E '^##(source|reference|contig=<ID=chr1,)'
bcftools index -n sample.vcf.gz          # number of records
samtools quickcheck -v sample.cram && echo "cram ok"

Record two things before you do anything else: the reference build, meaning GRCh37/hg19 versus GRCh38, and whether contigs are named chr1 or 1. Almost every downstream mismatch traces back to one of those two details. Array exports are usually GRCh37, while modern whole-genome sequencing is GRCh38 or T2T-CHM13.

2. Read one VCF line, field by field

The fastest way to understand a callset is to read a single line closely. Here is a real-shaped line from a 30x germline whole-genome callset, where 30x means each base was covered by an average of thirty independent reads:

#CHROM POS       ID          REF  ALT QUAL   FILTER INFO                      FORMAT       SAMPLE
chr7   117559590 rs113993960 CTT  C   1283.6 PASS   AC=1;AF=0.5;DP=38;MQ=60   GT:AD:DP:GQ:PL 0/1:19,19:38:99:1291,0,1274

Read it left to right. chr7 117559590 is the position on the stated build. CTT → C is a two-base deletion: the caller left-aligns the event and pads it with one anchoring base, which is why the REF field starts with a base that is not part of the deletion. rs113993960 is the dbSNP identifier, which is useful as a lookup key and meaningless as evidence of anything.

The next two fields describe confidence in the call. QUAL 1283.6 is a Phred-scaled confidence that a non-reference allele exists here at all; above roughly 100 on a germline callset it stops being informative. FILTER PASS means the caller’s own hard filters, or its variant quality score recalibration step, did not flag the site. Anything other than PASS should be treated as absent unless you go and look at the reads.

The FORMAT and SAMPLE pair is where the answer to “what do the numbers mean” lives. Each key in FORMAT corresponds to a colon-separated value in SAMPLE:

  • GT 0/1 is the genotype. 0 is REF and 1 is the first ALT. So 0/0 is reference homozygous and 0/1 heterozygous, while 1/1 is homozygous alternate and ./. a no call. A | instead of / means the call is phased, so you know which parental haplotype each allele sits on.
  • AD 19,19 is the allele depth: 19 reads supporting REF and 19 supporting ALT. For a real heterozygote you want the ratio near 0.5. Under 0.25 or over 0.75 alongside a 0/1 call is a red flag for a mapping artifact or mosaicism.
  • DP 38 is the total read depth at the site. Under about 10, a heterozygote can be missed entirely by chance.
  • GQ 99 is the Phred confidence in that particular genotype call, capped at 99. Below 20 means the caller could not distinguish the reported genotype from the next most likely one.
  • PL 1291,0,1274 is the normalized likelihood of 0/0, 0/1, and 1/1. The zero marks the called genotype, and the other two numbers say how much less likely the alternatives are, in Phred units.

Taken together, the line says that one copy carries a two-base deletion at this position. That call rests on 38 reads with a clean 19/19 split and maximum genotype confidence. That is what “a DNA result” looks like at the level where the data lives.

To pull many lines into a readable table, query a region directly:

bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%QUAL\t%FILTER[\t%GT\t%AD\t%DP\t%GQ]\n' \
  -r chr7:117500000-117700000 sample.vcf.gz | column -t

3. Check the callset before trusting any single line

Individual variants are only as good as the file they came from, so it pays to look at the whole callset first. Genome-wide summary statistics catch sample swaps, contamination, and bad libraries far faster than staring at individual variants.

bcftools stats -s - sample.vcf.gz > stats.txt
grep -E '^SN|^TSTV' stats.txt

On a single-sample germline whole-genome callset, three numbers tell you most of what you need:

  • Total SNVs should be 3.5–4.8 million and indels 500k–900k. Two million SNVs means you have an exome rather than a genome, while eight million suggests contamination or a broken filter.
  • The transition-to-transversion ratio (Ti/Tv) should sit between 1.9 and 2.1 genome-wide, or around 3.0 if you subset to coding regions. A drop toward 1.5 means false positives are leaking in.
  • The ratio of heterozygous to non-reference-homozygous calls should be 1.3–2.0, varying with ancestry. Well below 1.0 suggests either loss-of-heterozygosity artifacts or a consanguinity signal worth checking against runs of homozygosity, rather than a file error.

Coverage is the other half of the picture, and you can only measure it from the alignment:

mosdepth --by exons_grch38.bed --fast-mode --threads 4 sample sample.cram
zcat sample.regions.bed.gz | awk '$4 < 20' | wc -l    # regions under 20x
cat sample.mosdepth.summary.txt

Mean depth is the number people quote and the least useful one. Coverage uniformity is what determines whether a specific gene was readable, and it depends heavily on library preparation, particularly the fragmentation method and the GC content of the target 2. A genome with a 30x mean can still have a clinically relevant exon sitting at 4x because it is GC-rich or lives in a segmental duplication. Get the per-region numbers for the genes you care about before you conclude anything from the absence of variants there.

4. Distinguish “reference” from “not sequenced”

This is the single most common misreading of DNA results, and it is worth dwelling on. A plain VCF only lists positions where a variant was called. So when a position is missing from the file, that absence has two possible meanings. Either the sample matched the reference, or there was no usable data there. A gVCF resolves the ambiguity by emitting non-variant blocks alongside the variants:

bcftools query -f '%CHROM\t%POS\t%END\t[%GT\t%DP\t%GQ]\n' \
  -r chr13:32338000-32339000 sample.g.vcf.gz
chr13  32338100  32338412  0/0  31  70
chr13  32338413  32338425  ./.  3   0

The second line is a 13-base hole at 3x depth. Any query against that interval will return “no variant found,” and that answer carries no information. Before you report a negative, confirm the region was callable at GQ ≥ 20 and depth ≥ 20.

Array data has the same problem in a harsher form, because every position not on the chip is simply unmeasured. Imputation against a reference panel can infer common variants in the gaps with high accuracy, but it cannot recover rare or novel alleles, and those are exactly the ones with large effects. Deep whole-genome sequencing of 1,070 individuals showed how much rare variation sits outside the common-SNP catalogs that arrays are built from 3.

5. Annotate

Raw coordinates are not interpretable on their own. To make sense of a variant you need its consequence on each transcript, its frequency in reference populations, and any existing clinical classification. Ensembl VEP attaches all three:

vep -i sample.vcf.gz -o sample.vep.vcf.gz \
  --cache --offline --assembly GRCh38 --vcf --compress_output bgzip \
  --everything --pick_allele_gene --mane_select \
  --custom clinvar_20250715.vcf.gz,ClinVar,vcf,exact,0,CLNSIG,CLNREVSTAT \
  --custom gnomad.genomes.v4.1.sites.vcf.bgz,gnomADg,vcf,exact,0,AF,AF_popmax \
  --fork 4

Two flags matter more than the rest. --mane_select pins consequences to the MANE Select transcript, so your HGVS notation matches what a clinical laboratory would report instead of describing some minor isoform on which the same variant looks harmless. --pick_allele_gene returns one consequence per gene per allele, which keeps the output readable. Bear in mind that population frequency is only as good as the match between the reference cohort and your own ancestry: a variant that is rare in gnomAD’s non-Finnish European subset can be common in a population that the panel undersamples 3.

6. Cut it down to a list you can read

With annotation in place, four to five million variants become about twenty lines under a few filters:

bcftools +split-vep sample.vep.vcf.gz -f '%CHROM\t%POS\t%REF\t%ALT\t%SYMBOL\t%HGVSc\t%HGVSp\t%Consequence\t%gnomADg_AF\t%ClinVar_CLNSIG[\t%GT\t%DP\t%GQ]\n' \
  -i 'FILTER="PASS" && FMT/GQ>=30 && FMT/DP>=15' -d -A tab \
| awk -F'\t' '($9=="" || $9+0 < 0.001) &&
              ($8 ~ /frameshift|stop_gained|splice_acceptor|splice_donor|start_lost/ || $10 ~ /athogenic/)' \
| sort -k5,5 | column -t

Each cut has a reason behind it. Requiring FILTER=PASS, GQ≥30, and DP≥15 removes most technical noise. An allele frequency under 0.1% keeps variants rare enough to plausibly have a large effect. The consequence list keeps loss-of-function classes, where the mechanism is clear. The ClinVar branch keeps anything already classified as pathogenic regardless of predicted consequence, which catches missense and splice-region variants that the consequence filter would otherwise drop.

Once you have that short list, check review status before reading too much into any entry. A ClinVar record carrying CLNREVSTAT=no_assertion_criteria_provided is one submitter’s opinion. Two or more stars is a different quality of evidence. That means criteria_provided,_multiple_submitters,_no_conflicts or better. Filter on it.

7. What a “positive” result looks like

It is worth being precise about what a positive genotyping result is, because the word covers two separate things that routinely get collapsed together. The first is that a specific allele was observed, which is a measurement: the genotype, the depth, the allele balance. The second is that the allele carries a classification, which is an interpretation, expressed through the ACMG/AMP five-tier scheme of pathogenic, likely pathogenic, uncertain significance, likely benign, and benign. Most rare variants you find will land in “uncertain significance,” and that category means the evidence is insufficient, not that the answer is mildly bad.

There are three things a positive line does not tell you. Zygosity determines relevance for recessive conditions, so a single 0/1 in a recessive gene makes you a carrier, and carrier status for any given recessive condition is unremarkable. Penetrance varies enormously, and estimates drawn from families ascertained because they were affected will overstate risk in an unselected person. And a research-grade pipeline is not a clinical result. Clinical whole-genome reporting involves orthogonal confirmation, curated gene lists, and phenotype-driven filtering. Studies of diagnostic whole-genome sequencing in paediatrics show how much of the yield comes from that interpretation layer rather than from the sequencing itself 4.

For that reason, if a variant in your own data is classified pathogenic in a gene with medical consequences, do not act on it. Take it to a clinical geneticist or genetic counselor, who will order confirmation in an accredited laboratory. Your research pipeline’s job was to generate a hypothesis.

8. Relatedness, the other kind of DNA result

Paternity and relationship tests answer a different question from variant interpretation, and they use a different computation: they estimate a kinship coefficient from genotypes shared between two samples. Given two VCFs merged together, or array data, KING-style estimation is one command:

plink2 --vcf merged.vcf.gz --make-king-table --king-table-filter 0.0442 --out kin

The output gives a kinship coefficient φ. Expect roughly 0.5 for identical samples, 0.25 for parent-offspring and full siblings, 0.125 for half-siblings, uncle/niece, and grandparent-grandchild pairs, and 0.0625 for first cousins. The KING thresholds used in practice are above 0.354 for duplicates, 0.177–0.354 for first-degree relationships, and 0.0884–0.177 for second-degree.

Parent-offspring pairs separate from full siblings by the IBS0 column, which is the fraction of sites where one sample is 0/0 and the other 1/1. That value should be essentially zero for parent-offspring and clearly nonzero for siblings. Commercial reports convert the same evidence into a “combined paternity index” and a probability, which is this calculation with a prior attached.

Common problems

A handful of failure modes account for most wrong answers in genomic analysis, and they are all avoidable once you know to look for them.

Coordinates mismatch silently. A GRCh37 array file queried against GRCh38 annotations returns wrong genes for most positions. Convert deliberately with CrossMap or Picard LiftoverVcf, keep the reject file, and never assume.

Indel representation differs between callers. CTT→C and C→CTT at neighboring positions can describe the same event. Always run bcftools norm -f ref.fa -m -any -c w before comparing or joining callsets, otherwise you will get false negatives on intersection.

Strand ambiguity affects array data. A/T and C/G SNPs cannot be resolved by allele content alone. Vendors report on the plus strand of the reference, but merging files from two vendors without checking will scramble those sites.

Homologous regions produce phantom heterozygotes. Genes with pseudogenes or segmental duplications collect reads that belong elsewhere. The signs are an allele balance far from 0.5, an MQ below 40, and clusters of variants within a few hundred bases. Open the region in IGV before believing anything there.

Missing-because-uncovered gets read as reference. This was covered in step 4 and bears repeating, because it is the error that produces false reassurance rather than false alarm.

Multi-allelic sites break naive parsers. ALT=A,AT with GT=1/2 means two different non-reference alleles and no reference copy. Splitting with bcftools norm -m -any first avoids writing code that assumes one ALT per line.

VEP consequences can be taken from the wrong transcript. Without --mane_select or --pick_allele_gene, a variant can appear as intron_variant on one isoform and stop_gained on another, and whichever your parser grabs first becomes your answer.

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

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

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

  4. Dimitri J. Stavropoulos, Daniele Merico, Rebekah Jobling, et al. Whole-genome sequencing expands diagnostic utility and improves clinical management in paediatric medicine. npj Genomic Medicine, 2016. https://doi.org/10.1038/npjgenmed.2015.12 ↩