How to Read Your DNA Lab Results: From FASTQ to Annotated Variants

By the end of this guide you will have taken a set of DNA lab results, verified that they are internally consistent and yours, and produced an annotated variant table you can query with SQL or hand to an agent. You need the deliverables from the lab, ideally FASTQ or BAM/CRAM plus a VCF and a QC report. You also need a Linux or macOS machine with at least 8 cores and 32 GB of RAM. Plan on roughly 500 GB of free disk for a 30x whole genome and a working conda or container setup. Everything below uses open tools. For file handling and QC you need samtools, bcftools, and mosdepth. For the rest you need somalier, bwa-mem2 or dragmap, and GATK. Annotation and interval work use VEP or snpEff and bedtools. If your results came from a consumer genotyping array rather than sequencing, skip to step 7, which covers what an array file can and cannot tell you.

A note on scope before we start. Everything here is measurement and interpretation of data. Reading a variant annotation is not the same as knowing what it means for your health. Several of the categories you will encounter require a clinical geneticist or genetic counselor to interpret in the context of your history and family history. Those categories include pathogenic variants in cancer predisposition or cardiac genes, carrier status relevant to reproduction, and pharmacogenomic star alleles. Where you hit one of those, stop analyzing and make an appointment.

1. Inventory what the lab sent you

The phrase “DNA lab results” covers wildly different artifacts, and the first task is to figure out which one you have. A paternity or relationship lab sends a PDF with a table of 15 to 24 short tandem repeat (STR) markers and a combined paternity index. A clinical diagnostic lab sends a PDF report naming zero to a handful of variants with ACMG classifications, and usually nothing else unless you ask. A consumer array company sends a tab-delimited text file of about 600,000 to 1,000,000 genotypes. A sequencing provider should send you actual data files.

For whole-genome or whole-exome sequencing, the deliverable set you want starts with per-lane or merged FASTQ, meaning gzipped read files and typically 2 files per sample for paired-end. It also includes an aligned BAM or CRAM with its index and a single-sample VCF or gVCF with index. The last item is a QC report, either FastQC/MultiQC HTML or a vendor equivalent. Ask which reference genome was used. GRCh38 with ALT contigs, GRCh38 without ALTs, and T2T-CHM13 are all in circulation. Mixing coordinate systems across those is the single most common way to produce confidently wrong answers.

Check the sizes before anything else. A 30x human genome in gzipped FASTQ is roughly 80 to 110 GB total. The same data as a CRAM is roughly 15 to 25 GB, and as a BAM roughly 50 to 70 GB. A single-sample WGS VCF with all sites is tens of gigabytes as a gVCF and roughly 1 to 2 GB gzipped if it contains only variant sites. If your “whole genome” VCF is 40 MB, you have an exome or a targeted panel. You do not have a genome.

ls -lh
md5sum -c checksums.md5     # labs should provide these; if not, ask
samtools quickcheck -v *.bam *.cram && echo "containers intact"

samtools quickcheck verifies the EOF block and header, which catches truncated downloads. It does not verify content, which is what the next two steps are for.

2. Read the header and the QC numbers

The BAM/CRAM header records exactly how the file was made, and it is the most information-dense thing in the whole delivery.

samtools view -H sample.cram | grep -E '^@RG|^@PG|^@SQ' | head -30

From @RG you get the sample name and library. It also gives you the platform and sequencing center. Confirm the sample ID matches your accession. From @PG you get the aligner and its full command line, which tells you the reference path and any nonstandard flags. From @SQ you get the contig names and lengths: chr1 at length 248,956,422 means GRCh38, 1 at 249,250,621 means GRCh37/hg19, and 248,387,328 means T2T-CHM13v2.0. Write down which one you have and never guess again.

Then compute the numbers that determine what your data can support. Coverage depth and its uniformity govern which variant classes you can call reliably, and low-coverage regions are where false negatives hide ¹.

samtools stats -@ 8 sample.cram > stats.txt
grep ^SN stats.txt | cut -f2-

mosdepth --by 1000 --fast-mode -t 4 sample sample.cram
# genome-wide mean depth and coverage thresholds
cat sample.mosdepth.summary.txt | head

A usable 30x WGS should hit mean autosomal depth at or above 28x, with at least 90 percent of the callable genome at 20x or more. The duplicate rate should stay under 15 percent and the mapped-read fraction above 98 percent. For a standard PCR-free library the insert size mode should fall between 300 and 450 bp. An exome should show 90 percent or more of its target BED at 20x. If duplicate rate is above 25 percent, the library was over-amplified from too little input DNA, and your allele fractions in low-coverage regions are unreliable.

3. Verify the sample is yours and is one person

This step catches the failure modes people find most alarming: sample swaps at accessioning and cross-contamination during library prep. Both happen, and neither is visible in a variant table without a deliberate check.

somalier extracts genotypes at a few thousand common sites and gives you relatedness, inferred sex, and an ancestry estimate in a couple of minutes.

somalier extract -d extracted/ --sites sites.hg38.vcf.gz -f GRCh38.fa sample.cram
somalier relate -o relate extracted/*.somalier
# open relate.html; also check relate.samples.tsv for X/Y depth and het ratio

Two things to read. First, inferred sex from X heterozygosity and Y coverage should match your records. A mismatch means either a swap or a sex chromosome aneuploidy, and telling those apart requires looking at per-chromosome depth and then a clinician. Second, if you have data from a relative or an earlier test, the relatedness coefficient should land near 0.5 for a parent or full sibling and near 1.0 for the same person. A coefficient near zero against your own prior sample means one of the two files is not you.

For contamination, use VerifyBamID2 or the CHARR-style estimate that GATK reports:

VerifyBamID --SVDPrefix resource/1000g.phase3.100k.b38.vcf.gz.dat \
  --Reference GRCh38.fa --BamFile sample.cram --NumThread 8

Freemix (the contamination estimate) below 0.01 is clean. Between 0.01 and 0.03 is marginal and will inflate your apparent heterozygous count. Above 0.05, the sample is a mixture and single-sample variant calls should not be trusted at all.

4. If you only got a BAM, call variants yourself

Labs sometimes deliver alignments and a thin VCF from a pipeline you cannot inspect. Calling your own gives you control over filters and a record of exactly what produced each genotype. For germline single-nucleotide variants (SNVs) and short insertions/deletions (indels), we use DeepVariant, because it consistently outperforms hand-tuned hard filters on indels and needs no training on your part.

run_deepvariant \
  --model_type=WGS \
  --ref=GRCh38.fa \
  --reads=sample.cram \
  --output_vcf=sample.dv.vcf.gz \
  --output_gvcf=sample.dv.g.vcf.gz \
  --num_shards=16

On 16 cores a 30x genome takes roughly 6 to 12 hours. Expect about 4.5 to 5.0 million variants total for a single human genome against GRCh38, of which roughly 4 to 4.5 million are SNVs and 500,000 to 800,000 are indels. Counts far outside that range are a QC signal: too few usually means low coverage or a sample restricted to a target region, too many usually means contamination or a reference mismatch.

Structural variants need separate tools, because short reads resolve them poorly and callers disagree. Run manta and smoove (Lumpy) and treat the intersection as your high-confidence set, with the union as leads only. Short-read sequencing systematically misses large events, repeat expansions, and anything in segmental duplications. This limitation is intrinsic to read length rather than a defect of any one pipeline ²³.

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

Ts/Tv (transition-to-transversion ratio) should be about 2.0 to 2.1 genome-wide and about 3.0 in exonic regions. A genome-wide ratio below 1.8 indicates false-positive inflation.

5. Annotate, and keep the provenance

Raw variant records are coordinates and alleles. Annotation attaches gene, consequence, population frequency, and clinical assertion, and it is the step where most interpretation errors originate, because annotation depends entirely on which transcript set and which database version you used.

vep --input_file sample.dv.vcf.gz --output_file sample.vep.vcf.gz \
  --vcf --compress_output bgzip --fork 8 \
  --cache --offline --dir_cache $VEP_CACHE --assembly GRCh38 \
  --everything --pick_allele_gene --transcript_version \
  --plugin CADD,whole_genome_SNVs.tsv.gz \
  --custom clinvar.vcf.gz,ClinVar,vcf,exact,0,CLNSIG,CLNREVSTAT,CLNDN \
  --custom gnomad.genomes.v4.1.sites.vcf.gz,gnomADg,vcf,exact,0,AF,AF_popmax,nhomalt

In a text file next to the output, record the VEP version and the Ensembl release. Record the ClinVar download date and the gnomAD version there as well. ClinVar assertions change: variants get reclassified from “pathogenic” to “conflicting” or to “benign” as evidence accumulates, and a report built on a 2021 snapshot can disagree with one built today on the same DNA. Note also that --pick_allele_gene collapses to one consequence per gene, which is convenient for filtering and wrong if you care about a specific clinically relevant transcript. For a named gene of interest, pull all transcripts.

6. Load it into something queryable

A VCF is a poor interface for exploration. Flatten it into a table and query it.

bcftools +split-vep sample.vep.vcf.gz -d -f \
  '%CHROM\t%POS\t%REF\t%ALT\t%QUAL\t[%GT\t%DP\t%AD]\t%SYMBOL\t%Consequence\t%IMPACT\t%gnomADg_AF\t%ClinVar_CLNSIG\t%ClinVar_CLNREVSTAT\t%CADD_PHRED\n' \
  -A tab -i 'FILTER="PASS"' > variants.tsv

Then in DuckDB:

CREATE TABLE v AS SELECT * FROM read_csv_auto('variants.tsv', header=false);

-- rare, high-impact, coding
SELECT * FROM v
WHERE IMPACT = 'HIGH'
  AND (gnomADg_AF IS NULL OR gnomADg_AF < 0.001)
  AND DP >= 15;

That query typically returns 100 to 300 rows for one genome. Nearly all of them are irrelevant: loss-of-function variants in genes with no disease association, annotation artifacts on poorly supported transcripts, and homozygous-reference calls mislabeled by transcript edge cases. The tightest useful filter is ClinVar pathogenic or likely pathogenic with two-star or better review status, restricted to genes with established gene-disease validity. That usually returns between zero and a few rows.

Two interpretive rules keep this accurate. Rare and predicted-damaging is not the same as causal. Most rare coding variants have no detectable phenotypic effect, and the statistical burden for implicating a single variant in a single person is high ⁴. And a variant discovered incidentally in someone with no relevant symptoms or family history carries far lower predictive value than the same variant found in an affected patient, because penetrance estimates from clinically ascertained families overstate risk in the general population. A clinician is the right person to make that call.

7. If your result is an array file or a paper report

A consumer array is a genotyping assay: it interrogates a fixed set of pre-chosen positions and reports the allele at each. It does not read your genome. A typical file looks like this:

rsid	chromosome	position	genotype
rs4477212	1	82154	AA
rs3094315	1	752566	AG

Three consequences follow. Positions not on the chip produce no data, so absence of a variant in an array file is not evidence of its absence in your genome. Rare-variant calls from arrays have poor positive predictive value, because the clustering algorithms are tuned on common alleles, and a single miscalled rare genotype looks identical to a real one in the output file. Any clinically meaningful finding from an array needs orthogonal confirmation by Sanger sequencing, the capillary electrophoresis method that remains the reference standard for single-locus confirmation ⁵.

Paternity and relationship reports are a different measurement entirely. They compare STR marker profiles and report a combined paternity index and a probability, which is calculated against population allele frequencies for the assumed ancestry group. A “positive” result is conventionally a probability of 99.9 percent or higher with no marker exclusions. An “exclusion” is typically two or more markers where the alleged father shares no allele with the child. Those reports do not contain health information and do not generalize.

On the privacy question people ask: results are not public record. A clinical result lives in your medical record under HIPAA, and a legal paternity result becomes discoverable only if entered into court proceedings. What matters more in practice is your own choices. Uploading genotype data to third-party matching services makes it visible to relatives and, depending on the platform’s terms, to law enforcement. Read the terms before uploading, and keep your primary copy on storage you control.

Common problems

Wrong reference build. A VCF at GRCh37 coordinates annotated against GRCh38 produces plausible-looking gene names that are wrong. Symptoms: implausibly high rates of novel variants, ClinVar hits in unexpected genes. Check @SQ lengths. If you must convert, use CrossMap with the UCSC chain file and then re-normalize with bcftools norm -f. Accept that a small percentage of sites will fail to lift over.

Unnormalized indels that silently fail to match databases. The same deletion can be written several ways. Always run bcftools norm -f ref.fa -m -both -c w before comparing anything to ClinVar or gnomAD, or you will get false negatives on real findings.

Low-coverage genes that look clean because nothing was called. Genes with high GC content or paralogs (for example, the first exon of several cardiac genes, and any gene with a close pseudogene) frequently drop below callable depth. Check explicitly: mosdepth --by genes.bed and look at the per-gene mean, rather than assuming a quiet region is a normal region.

Degraded or low-input DNA. Buccal swabs stored warm, or blood held too long before extraction, yield fragmented DNA. The signature is short insert sizes, elevated duplicate rate, and an excess of C-to-T transitions from deamination. The library reflects the sample condition, and the fix is a new collection rather than a filter.

Turnaround expectations. Sequencing chemistry for a 30x genome is roughly 24 to 48 hours on a current instrument, and alignment plus calling is another 6 to 24 hours on modest hardware. Real-world turnaround from sample receipt to delivered files is typically two to six weeks, because batching, accessioning, QC review, and reporting dominate. Clinical reports with curated interpretation run longer. The gap between instrument time and delivery time is process, not chemistry, and the underlying throughput has improved by orders of magnitude over four decades of development ³.

Expecting more than short reads can give. Short-read WGS is the wrong instrument if you are looking for a repeat expansion, a large structural rearrangement, or phasing across a long haplotype. A negative result there means the assay could not see it. Long-read sequencing or a targeted assay is the right next step, and which one depends on the specific question.

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. D. C. Koboldt, L. Ding, E. R. Mardis, et al. Challenges of sequencing human genomes. Briefings in Bioinformatics, 2010. https://doi.org/10.1093/bib/bbq016 ↩

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

5. Barry L. Karger, András Guttman. DNA sequencing by CE. ELECTROPHORESIS, 2009. https://doi.org/10.1002/elps.200900218 ↩